Coverity bugs corrected.
[u/mrichter/AliRoot.git] / THijing / THijing.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 ////////////////////////////////////////////////////////////////////////////////
19 //                                                                            //
20 // THijing                                                                    //
21 //                                                                            //
22 // THijing is an interface class to F77 version of Hijing 1.36                //
23 // event generator, written by X.N. Wang and  M. Gyulassy.                    //
24 // For details see http://nta2.lbl.gov/~xnwang                                //
25 //                                                                            //
26 //          **************************************************                //
27 //          *     |             _______      /  ------/      *                //
28 //          *   ----- ------     |_____|     /_/     /       *                //
29 //          *    ||    /        |_____|      /    /          *                //
30 //          *    /|   /_/       /_______    /_  /    _       *                //
31 //          *   / |     / /     /  /  / |        -------     *                //
32 //          *     |    / /       /  /  |     /     |         *                //
33 //          *     |   / /       /  / _|    /   -------       *                //
34 //          *                                                *                //
35 //          **************************************************                //
36 //                                HIJING                                      //
37 //                 Heavy Ion Jet INteraction Generator                        //
38 //                                  by                                        //
39 //                     X. N. Wang  and  M. Gyulassy                           //
40 //                     Lawrence Berkeley Laboratory                           //
41 //****************************************************************************//
42
43
44 #include <TClonesArray.h>
45 #include <TObjArray.h>
46 #include <TParticle.h>
47 #include <TROOT.h>
48
49 #include "Hcommon.h"
50 #include "THijing.h"
51
52 #ifndef WIN32
53 # define hijset hijset_
54 # define hijing hijing_
55 # define profile profile_
56 # define rluget_hijing rluget_hijing_
57 # define rluset_hijing rluset_hijing_
58 # define type_of_call
59 #else
60 # define hijset HIJSET
61 # define hijing HIJING
62 # define profile PROFILE
63 # define rluget_hijing RLUGET_HIJING
64 # define rluset_hijing RLUSET_HIJING
65 # define type_of_call _stdcall
66 #endif
67
68 #ifndef WIN32
69 extern "C" void type_of_call hijset(Float_t & , const char *, 
70                                    const char *, const char *,
71                                    Int_t & , Int_t &, Int_t &,
72                                    Int_t &, const int, 
73                                    const int, const int);
74 extern "C" float type_of_call profile(Float_t &);
75 extern "C" void type_of_call hijing(const char *, Float_t  &,
76                                    Float_t &, const int);
77 extern "C" void type_of_call rluget_hijing(Int_t & lfn, Int_t & move);
78
79 extern "C" void type_of_call rluset_hijing(Int_t & lfn, Int_t & move);
80
81 #else
82 #endif
83
84
85
86 ClassImp(THijing)
87
88
89 THijing::THijing(): 
90     TGenerator("Hijing","Hijing"),
91     fEfrm(5500.),
92     fFrame("CMS"),
93     fProj("A"),
94     fTarg("A"),
95     fIap(208),
96     fIzp(82),
97     fIat(208),
98     fIzt(82),
99     fBmin(0.),
100     fBmax(5.)
101 {
102 // Default constructor 
103 }
104
105 //______________________________________________________________________________
106 THijing::THijing(Float_t efrm, const char *frame="CMS", 
107                  const char *proj="A", const char *targ="A", 
108                  Int_t iap=207, Int_t izp=82, Int_t iat=207, Int_t izt=82,
109                  Float_t bmin=0, Float_t bmax=20): 
110     TGenerator("Hijing","Hijing"),
111     fEfrm(efrm),
112     fFrame(frame),
113     fProj(proj),
114     fTarg(targ),
115     fIap(iap),
116     fIzp(izp),
117     fIat(iat),
118     fIzt(izt),
119     fBmin(bmin),
120     fBmax(bmax)
121 {
122 // THijing constructor: 
123 // Note that there may be only one functional THijing object
124 // at a time, so it's not use to create more than one instance of it.
125 }
126
127 //______________________________________________________________________________
128 THijing::~THijing()
129 {
130 // Destructor
131 }
132
133
134 TObjArray* THijing::ImportParticles(Option_t *option)
135 {
136 //
137 //  Default primary creation method. It reads the /HEPEVT/ common block which
138 //  has been filled by the GenerateEvent method. If the event generator does
139 //  not use the HEPEVT common block, This routine has to be overloaded by
140 //  the subclasses.
141 //  The function loops on the generated particles and store them in
142 //  the TClonesArray pointed by the argument particles.
143 //  The default action is to store only the stable particles (ISTHEP = 1)
144 //  This can be demanded explicitly by setting the option = "Final"
145 //  If the option = "All", all the particles are stored.
146 //
147     fParticles->Clear();
148     Int_t numpart = HIMAIN1.natt;
149     printf("\n THijing: HIJING stack contains %d particles.", numpart);
150     printf("\n THijing: Total energy:         %f           ", HIMAIN1.eatt);
151     printf("\n THijing: Number of hard scatterings: %d     ", HIMAIN1.jatt);
152     Int_t nump = 0;
153     if (!strcmp(option,"") || !strcmp(option,"Final")) {
154         for (Int_t i = 0; i < numpart; i++) {
155           
156             if (HIMAIN2.katt[3][i] == 1) {
157 //
158 //  Use the common block values for the TParticle constructor
159 //
160                 nump++;
161                 TParticle* p = new TParticle(
162                     HIMAIN2.katt[0][i], HIMAIN2.katt[3][i] ,
163                     -1, -1, -1, -1,
164                     HIMAIN2.patt[0][i], HIMAIN2.patt[1][i], HIMAIN2.patt[2][i], HIMAIN2.patt[3][i] ,
165                     HIMAIN2.vatt[0][i], HIMAIN2.vatt[1][i], HIMAIN2.vatt[2][i], HIMAIN2.vatt[3][i] 
166                   );
167                 p->SetUniqueID(HIMAIN2.katt[1][i]);
168                 fParticles->Add(p);
169             }
170         }
171     }
172     else if (!strcmp(option,"All")) {
173         nump = numpart; 
174         for (Int_t i = 0; i < numpart; i++) {
175             
176             Int_t iParent = HIMAIN2.katt[2][i]-1;
177             
178             if (iParent >= 0) {
179                 TParticle *mother = (TParticle*) (fParticles->UncheckedAt(iParent));       
180                 mother->SetLastDaughter(i);
181                 if (mother->GetFirstDaughter()==-1)
182                     mother->SetFirstDaughter(i);
183             }
184             
185             TParticle* p = new TParticle(
186                 HIMAIN2.katt[0][i], HIMAIN2.katt[3][i], iParent,
187                 -1, -1, -1,
188                 HIMAIN2.patt[0][i], HIMAIN2.patt[1][i], HIMAIN2.patt[2][i], HIMAIN2.patt[3][i] ,
189                 HIMAIN2.vatt[0][i], HIMAIN2.vatt[1][i], HIMAIN2.vatt[2][i], HIMAIN2.vatt[3][i]
190                 );
191             p->SetUniqueID(HIMAIN2.katt[1][i]);
192             fParticles->Add(p);
193         }
194     }
195     return fParticles;
196 }
197
198 Int_t THijing::ImportParticles(TClonesArray *particles, Option_t *option)
199 {
200 //
201 //  Default primary creation method. It reads the /HEPEVT/ common block which
202 //  has been filled by the GenerateEvent method. If the event generator does
203 //  not use the HEPEVT common block, This routine has to be overloaded by
204 //  the subclasses.
205 //  The function loops on the generated particles and store them in
206 //  the TClonesArray pointed by the argument particles.
207 //  The default action is to store only the stable particles (ISTHEP = 1)
208 //  This can be demanded explicitly by setting the option = "Final"
209 //  If the option = "All", all the particles are stored.
210 //
211   if (particles == 0) return 0;
212   TClonesArray &particlesR = *particles;
213   particlesR.Clear();
214   Int_t numpart = HIMAIN1.natt;
215   printf("\n THijing: HIJING stack contains %d particles.", numpart);
216   printf("\n THijing: Total energy:         %f           ", HIMAIN1.eatt);
217   printf("\n THijing: Number of hard scatterings: %d     ", HIMAIN1.jatt);
218   Int_t nump = 0;
219   if (!strcmp(option,"") || !strcmp(option,"Final")) {
220       for (Int_t i = 0; i < numpart; i++) {
221           
222           if (HIMAIN2.katt[3][i] == 1) {
223 //
224 //  Use the common block values for the TParticle constructor
225 //
226             nump++;
227             new(particlesR[i]) TParticle(
228                   HIMAIN2.katt[0][i] ,
229                   HIMAIN2.katt[3][i] ,
230                   -1 ,
231                   -1,
232                   -1,
233                   -1,
234                   
235                   HIMAIN2.patt[0][i] ,
236                   HIMAIN2.patt[1][i] ,
237                   HIMAIN2.patt[2][i] ,
238                   HIMAIN2.patt[3][i] ,
239
240                   HIMAIN2.vatt[0][i] ,
241                   HIMAIN2.vatt[1][i] ,
242                   HIMAIN2.vatt[2][i] ,
243                   HIMAIN2.vatt[3][i] 
244                   );
245             particlesR[i]->SetUniqueID(HIMAIN2.katt[1][i]);
246           }
247       }
248   }
249   else if (!strcmp(option,"All")) {
250       nump = numpart; 
251       for (Int_t i = 0; i < numpart; i++) {
252
253           Int_t iParent = HIMAIN2.katt[2][i]-1;
254           
255           if (iParent >= 0) {
256               TParticle *mother = (TParticle*) (particlesR.UncheckedAt(iParent));          
257               mother->SetLastDaughter(i);
258               if (mother->GetFirstDaughter()==-1)
259                   mother->SetFirstDaughter(i);
260           }
261
262           new(particlesR[i]) TParticle(
263               HIMAIN2.katt[0][i] ,
264               HIMAIN2.katt[3][i] ,
265               iParent,
266               -1,
267               -1,
268               -1,
269               
270               HIMAIN2.patt[0][i] ,
271               HIMAIN2.patt[1][i] ,
272               HIMAIN2.patt[2][i] ,
273               HIMAIN2.patt[3][i] ,
274               
275               HIMAIN2.vatt[0][i] ,
276               HIMAIN2.vatt[1][i] ,
277               HIMAIN2.vatt[2][i] ,
278               HIMAIN2.vatt[3][i]
279               );
280           particlesR[i]->SetUniqueID(HIMAIN2.katt[1][i]);
281       }
282   }
283   return nump;
284 }
285
286 //______________________________________________________________________________
287 void THijing::SetEFRM(Float_t efrm)
288 {
289 // Set the centre of mass (CMS) or lab-energy (LAB)
290    fEfrm=efrm;
291
292 //______________________________________________________________________________
293 void THijing::SetFRAME(const char* frame)
294 {
295 // Set the frame type ("CMS" or "LAB")
296    fFrame=frame;
297
298 //______________________________________________________________________________
299 void THijing::SetPROJ(const char* proj)
300 {
301 // Set the projectile type
302    fProj=proj;
303
304 //______________________________________________________________________________
305 void THijing::SetTARG(const char* targ)
306 {
307 // Set the target type
308    fTarg=targ;
309
310 //______________________________________________________________________________
311 void THijing::SetIAP(Int_t iap)
312 {
313 // Set the projectile atomic number
314    fIap=iap;
315
316 //______________________________________________________________________________
317 void THijing::SetIZP(Int_t izp)
318 {
319 // Set the projectile charge number
320    fIzp=izp;
321
322 //______________________________________________________________________________
323 void THijing::SetIAT(Int_t iat)
324 {
325 // Set the target atomic number
326    fIat=iat;
327
328 //______________________________________________________________________________
329 void THijing::SetIZT(Int_t izt)
330 {
331 // Set the target charge number
332    fIzt=izt;
333
334 //______________________________________________________________________________
335 void THijing::SetBMIN(Float_t bmin)
336 {
337 // Set the minimum impact parameter
338    fBmin=bmin;
339
340 //______________________________________________________________________________
341 void THijing::SetBMAX(Float_t bmax)
342 {
343 // Set the maximum impact parameter
344    fBmax=bmax;
345
346 //______________________________________________________________________________
347 Float_t THijing::GetEFRM() const
348 {
349 // Get the centre of mass (CMS) or lab-energy (LAB)
350    return fEfrm;
351
352 //______________________________________________________________________________
353 const char* THijing::GetFRAME() const
354 {
355 // Get the frame type ("CMS" or "LAB")
356    return fFrame.Data();
357
358 //______________________________________________________________________________
359 const char* THijing::GetPROJ() const
360 {
361 // Get the projectile type
362    return fProj;
363
364 //______________________________________________________________________________
365 const char* THijing::GetTARG() const
366 {
367 // Set the target type
368    return fTarg;
369
370 //______________________________________________________________________________
371 Int_t THijing::GetIAP() const
372 {
373 // Get the projectile atomic number
374    return fIap;
375
376 //______________________________________________________________________________
377 Int_t THijing::GetIZP() const
378 {
379 // Get the projectile charge number
380    return fIzp;
381
382 //______________________________________________________________________________
383 Int_t THijing::GetIAT() const
384 {
385 // Get the target atomic number
386    return fIat;
387
388 //______________________________________________________________________________
389 Int_t THijing::GetIZT() const
390 {
391 // Get the target charge number
392    return fIzt;
393
394 //______________________________________________________________________________
395 Float_t THijing::GetBMIN() const
396 {
397 // Get the minimum impact parameter
398    return fBmin;
399
400 //______________________________________________________________________________
401 Float_t THijing::GetBMAX() const
402 {
403 // Get the maximum impact parameter
404    return fBmax;
405
406
407 //====================== access to common HIPARNT ===============================
408
409 //______________________________________________________________________________
410 void THijing::SetHIPR1(Int_t key,Float_t value)
411 {
412 // Set the values of array HIPR1 in common HIPARNT
413    if ( key<1 || key>100 ) {
414       printf ("ERROR in THijing:SetHIPR1(key,value): \n ");
415       printf ("      key=%i is out of range [1..100]!\n",key);
416       return;
417    }
418
419    HIPARNT.hipr1[key-1]=value;
420
421 }
422
423 //______________________________________________________________________________
424 Float_t THijing::GetHIPR1(Int_t key) const
425 {
426 // Get the values of array HIPR1 in common HIPARNT
427    if ( key<1 || key>100 ) {
428       printf ("ERROR in THijing:GetHIPR1(key): \n ");
429       printf ("      key=%i is out of range [1..100]!\n",key);
430       return 0;
431    }
432
433    return HIPARNT.hipr1[key-1];
434
435 }
436
437 //______________________________________________________________________________
438 void THijing::SetIHPR2(Int_t key,Int_t value)
439 {
440 // Set the values of array HIPR2 in common HIPARNT
441    if ( key<1 || key>50 ) {
442       printf ("ERROR in THijing:SetIHPR2(key,value): \n ");
443       printf ("      key=%i is out of range [1..50]!\n",key);
444       return;
445    }
446
447    HIPARNT.ihpr2[key-1]=value;
448
449 }
450
451 //______________________________________________________________________________
452 Int_t THijing::GetIHPR2(Int_t key) const
453 {
454 // Get the values of array HIPR2 in common HIPARNT
455    if ( key<1 || key>50 ) {
456       printf ("ERROR in THijing:GetIHPR2(key): \n ");
457       printf ("      key=%i is out of range [1..50]!\n",key);
458       return 0;
459    }
460
461    return HIPARNT.ihpr2[key-1];
462
463 }
464
465
466 //______________________________________________________________________________
467 Float_t THijing::GetHINT1(Int_t key) const
468 {
469 // Get the values of array HINT1 in common HIPARNT
470    if ( key<1 || key>100 ) {
471       printf ("ERROR in THijing:GetHINT1(key): \n ");
472       printf ("      key=%i is out of range [1..100]!\n",key);
473       return 0;
474    }
475
476    return HIPARNT.hint1[key-1];
477
478 }
479
480
481 //______________________________________________________________________________
482 Int_t THijing::GetIHNT2(Int_t key) const
483 {
484 // Get the values of array HINT2 in common HIPARNT
485    if ( key<1 || key>50 ) {
486       printf ("ERROR in THijing:GetIHNT2(key): \n ");
487       printf ("      key=%i is out of range [1..50]!\n",key);
488       return 0;
489    }
490
491    return HIPARNT.ihnt2[key-1];
492
493 }
494
495
496 //====================== access to common HIMAIN1 ===============================
497
498 //______________________________________________________________________________
499 Int_t  THijing::GetNATT() const
500 {
501 // Get the number of particles produces
502    return HIMAIN1.natt;
503
504 }
505
506 //______________________________________________________________________________
507 Float_t  THijing::GetEATT() const
508 {
509 // Get total energy of particles
510
511    return HIMAIN1.eatt;
512
513 }
514
515 //______________________________________________________________________________
516 Int_t  THijing::GetJATT() const
517 {
518 // Get number of hard scatterings
519
520    return HIMAIN1.jatt;
521
522 }
523
524 //______________________________________________________________________________
525 Int_t  THijing::GetNT() const
526 {
527 // Get number of target participants
528
529    return HIMAIN1.nt;
530
531 }
532
533 //______________________________________________________________________________
534 Int_t  THijing::GetNP() const
535 {
536 // Get number of projectile participants
537    return HIMAIN1.np;
538
539 }
540
541
542 //______________________________________________________________________________
543 Int_t  THijing::GetN0() const
544 {
545 // Get number of N-N collisions
546    return HIMAIN1.n0;
547
548 }
549 //______________________________________________________________________________
550 Int_t  THijing::GetN01() const
551 {
552 // Get number of N-wounded collisions
553
554    return HIMAIN1.n01;
555
556 }
557
558 //______________________________________________________________________________
559 Int_t  THijing::GetN10() const
560 {
561 // Get number of wounded-N collisions
562
563    return HIMAIN1.n10;
564
565 }
566
567 //______________________________________________________________________________
568 Int_t  THijing::GetN11() const
569 {
570 // Get number of wounded-wounded collisions
571
572    return HIMAIN1.n11;
573
574 }
575
576 //______________________________________________________________________________
577 Float_t  THijing::GetBB() const
578 {
579 // Get impact parameter
580
581    return HIMAIN1.bb;
582
583 }
584
585 //====================== access to common HIMAIN2 ===============================
586
587 //______________________________________________________________________________
588 Int_t THijing::GetKATT(Int_t key1, Int_t key2) const
589 {
590 // Get values of array KATT in common HIMAIN2
591    if ( key1<1 || key1>200000 ) {
592       printf("ERROR in THijing::GetKATT(key1,key2):\n");
593       printf("      key1=%i is out of range [1..200000]\n",key1);
594       return 0;
595    }
596
597    if ( key2<1 || key2>4 ) {
598       printf("ERROR in THijing::GetKATT(key1,key2):\n");
599       printf("      key2=%i is out of range [1..4]\n",key2);
600       return 0;
601    }
602    
603    return   HIMAIN2.katt[key2-1][key1-1];
604 }
605
606 //______________________________________________________________________________
607 Float_t THijing::GetPATT(Int_t key1, Int_t key2) const
608 {
609 // Get values of array PATT in common HIMAIN2
610    if ( key1<1 || key1>200000 ) {
611       printf("ERROR in THijing::GetPATT(key1,key2):\n");
612       printf("      key1=%i is out of range [1..130000]\n",key1);
613       return 0;
614    }
615
616    if ( key2<1 || key2>4 ) {
617       printf("ERROR in THijing::GetPATT(key1,key2):\n");
618       printf("      key2=%i is out of range [1..4]\n",key2);
619       return 0;
620    }
621    
622    return   HIMAIN2.patt[key2-1][key1-1];
623 }
624
625 Float_t THijing::GetVATT(Int_t key1, Int_t key2) const
626 {
627 // Get values of array VATT in common HIMAIN2
628    if ( key1<1 || key1>200000 ) {
629       printf("ERROR in THijing::GetVATT(key1,key2):\n");
630       printf("      key1=%i is out of range [1..130000]\n",key1);
631       return 0;
632    }
633
634    if ( key2<1 || key2>4 ) {
635       printf("ERROR in THijing::GetVATT(key1,key2):\n");
636       printf("      key2=%i is out of range [1..4]\n",key2);
637       return 0;
638    }
639    
640    return   HIMAIN2.vatt[key2-1][key1-1];
641 }
642
643 //====================== access to common HIJJET1 ===============================
644
645 //______________________________________________________________________________
646 Int_t THijing::GetNPJ(Int_t key) const
647 {
648 // Get values of array NPJ of common HIJJET1
649    if ( key<1 || key>300 ) {
650       printf("ERROR in THijing::GetNPJ(key):\n");
651       printf("      key=%i is out of range [1..300]\n",key);
652       return 0;
653    }
654    return HIJJET1.npj[key-1];
655 }
656
657 //______________________________________________________________________________
658 Int_t THijing::GetKFPJ(Int_t key1, Int_t key2) const
659 {
660 // Get values of array KFPJ in common HIJJET1
661    if ( key1<1 || key1>300 ) {
662       printf("ERROR in THijing::GetKFPJ(key1):\n");
663       printf("      key1=%i is out of range [1..300]\n",key1);
664       return 0;
665    }
666    if ( key2<1 || key2>500 ) {
667       printf("ERROR in THijing::GetKFPJ(key1,key2):\n");
668       printf("      key2=%i is out of range [1..500]\n",key2);
669       return 0;
670    }
671    
672    return HIJJET1.kfpj[key2-1][key1-1];
673 }
674
675 //______________________________________________________________________________
676 Float_t THijing::GetPJPX(Int_t key1, Int_t key2) const
677 {
678 // Get values of array PJPX in common HIJJET1
679    if ( key1<1 || key1>300 ) {
680       printf("ERROR in THijing::GetPJPX(key1):\n");
681       printf("      key1=%i is out of range [1..300]\n",key1);
682       return 0;
683    }
684    if ( key2<1 || key2>500 ) {
685       printf("ERROR in THijing::GetPJPX(key1,key2):\n");
686       printf("      key2=%i is out of range [1..500]\n",key2);
687       return 0;
688    }
689    
690    return HIJJET1.pjpx[key2-1][key1-1];
691 }
692
693 //______________________________________________________________________________
694 Float_t THijing::GetPJPY(Int_t key1, Int_t key2) const
695 {
696 // Get values of array PJPY in common HIJJET1
697    if ( key1<1 || key1>300 ) {
698       printf("ERROR in THijing::GetPJPY(key1):\n");
699       printf("      key1=%i is out of range [1..300]\n",key1);
700       return 0;
701    }
702    if ( key2<1 || key2>500 ) {
703       printf("ERROR in THijing::GetPJPY(key1,key2):\n");
704       printf("      key2=%i is out of range [1..500]\n",key2);
705       return 0;
706    }
707    
708    return HIJJET1.pjpy[key2-1][key1-1];
709 }
710
711 //______________________________________________________________________________
712 Float_t THijing::GetPJPZ(Int_t key1, Int_t key2) const
713 {
714 // Get values of array PJPZ in common HIJJET1
715    if ( key1<1 || key1>300 ) {
716       printf("ERROR in THijing::GetPJPZ(key1):\n");
717       printf("      key1=%i is out of range [1..300]\n",key1);
718       return 0;
719    }
720    if ( key2<1 || key2>500 ) {
721       printf("ERROR in THijing::GetPJPZ(key1,key2):\n");
722       printf("      key2=%i is out of range [1..500]\n",key2);
723       return 0;
724    }
725    
726    return HIJJET1.pjpz[key2-1][key1-1];
727 }
728
729 //______________________________________________________________________________
730 Float_t THijing::GetPJPE(Int_t key1, Int_t key2) const
731 {
732 // Get values of array PJPE in common HIJJET1
733    if ( key1<1 || key1>300 ) {
734       printf("ERROR in THijing::GetPJPE(key1):\n");
735       printf("      key1=%i is out of range [1..300]\n",key1);
736       return 0;
737    }
738    if ( key2<1 || key2>500 ) {
739       printf("ERROR in THijing::GetPJPE(key1,key2):\n");
740       printf("      key2=%i is out of range [1..500]\n",key2);
741       return 0;
742    }
743    
744    return HIJJET1.pjpe[key2-1][key1-1];
745 }
746
747 //______________________________________________________________________________
748 Float_t THijing::GetPJPM(Int_t key1, Int_t key2) const
749 {
750 // Get values of array PJPM in common HIJJET1
751    if ( key1<1 || key1>300 ) {
752       printf("ERROR in THijing::GetPJPM(key1):\n");
753       printf("      key1=%i is out of range [1..300]\n",key1);
754       return 0;
755    }
756    if ( key2<1 || key2>500 ) {
757       printf("ERROR in THijing::GetPJPM(key1,key2):\n");
758       printf("      key2=%i is out of range [1..500]\n",key2);
759       return 0;
760    }
761    
762    return HIJJET1.pjpm[key2-1][key1-1];
763 }
764
765 //______________________________________________________________________________
766 Int_t THijing::GetNTJ(Int_t key) const
767 {
768 // Get values of array NTJ in common HIJJET1
769    if ( key<1 || key>300 ) {
770       printf("ERROR in THijing::GetNTJ(key):\n");
771       printf("      key=%i is out of range [1..300]\n",key);
772       return 0;
773    }
774    return HIJJET1.ntj[key-1];
775 }
776
777 //______________________________________________________________________________
778 Int_t THijing::GetKFTJ(Int_t key1, Int_t key2) const
779 {
780 // Get values of array KFTJ in common HIJJET1
781    if ( key1<1 || key1>300 ) {
782       printf("ERROR in THijing::GetKFTJ(key1):\n");
783       printf("      key1=%i is out of range [1..300]\n",key1);
784       return 0;
785    }
786    if ( key2<1 || key2>500 ) {
787       printf("ERROR in THijing::GetKFTJ(key1,key2):\n");
788       printf("      key2=%i is out of range [1..500]\n",key2);
789       return 0;
790    }
791    
792    return HIJJET1.kftj[key2-1][key1-1];
793 }
794
795 //______________________________________________________________________________
796 Float_t THijing::GetPJTX(Int_t key1, Int_t key2) const
797 {
798 // Get values of array PJTX in common HIJJET1
799    if ( key1<1 || key1>300 ) {
800       printf("ERROR in THijing::GetPJTX(key1):\n");
801       printf("      key1=%i is out of range [1..300]\n",key1);
802       return 0;
803    }
804    if ( key2<1 || key2>500 ) {
805       printf("ERROR in THijing::GetPJTX(key1,key2):\n");
806       printf("      key2=%i is out of range [1..500]\n",key2);
807       return 0;
808    }
809    
810    return HIJJET1.pjtx[key2-1][key1-1];
811 }
812
813 //______________________________________________________________________________
814 Float_t THijing::GetPJTY(Int_t key1, Int_t key2) const
815 {
816 // Get values of array PJTY in common HIJJET1
817    if ( key1<1 || key1>300 ) {
818       printf("ERROR in THijing::GetPJTY(key1):\n");
819       printf("      key1=%i is out of range [1..300]\n",key1);
820       return 0;
821    }
822    if ( key2<1 || key2>500 ) {
823       printf("ERROR in THijing::GetPJTY(key1,key2):\n");
824       printf("      key2=%i is out of range [1..500]\n",key2);
825       return 0;
826    }
827    
828    return HIJJET1.pjty[key2-1][key1-1];
829 }
830
831 //______________________________________________________________________________
832 Float_t THijing::GetPJTZ(Int_t key1, Int_t key2) const
833 {
834 // Get values of array PJTZ in common HIJJET1
835    if ( key1<1 || key1>300 ) {
836       printf("ERROR in THijing::GetPJTZ(key1):\n");
837       printf("      key1=%i is out of range [1..300]\n",key1);
838       return 0;
839    }
840    if ( key2<1 || key2>500 ) {
841       printf("ERROR in THijing::GetPJTZ(key1,key2):\n");
842       printf("      key2=%i is out of range [1..500]\n",key2);
843       return 0;
844    }
845    
846    return HIJJET1.pjtz[key2-1][key1-1];
847 }
848
849 //______________________________________________________________________________
850 Float_t THijing::GetPJTE(Int_t key1, Int_t key2) const
851 {
852 // Get values of array PJTE in common HIJJET1
853    if ( key1<1 || key1>300 ) {
854       printf("ERROR in THijing::GetPJTE(key1):\n");
855       printf("      key1=%i is out of range [1..300]\n",key1);
856       return 0;
857    }
858    if ( key2<1 || key2>500 ) {
859       printf("ERROR in THijing::GetPJTE(key1,key2):\n");
860       printf("      key2=%i is out of range [1..500]\n",key2);
861       return 0;
862    }
863    
864    return HIJJET1.pjte[key2-1][key1-1];
865 }
866
867 //______________________________________________________________________________
868 Float_t THijing::GetPJTM(Int_t key1, Int_t key2) const
869 {
870 // Get values of array PJTM in common HIJJET1
871    if ( key1<1 || key1>300 ) {
872       printf("ERROR in THijing::GetPJTM(key1):\n");
873       printf("      key1=%i is out of range [1..300]\n",key1);
874       return 0;
875    }
876    if ( key2<1 || key2>500 ) {
877       printf("ERROR in THijing::GetPJTM(key1,key2):\n");
878       printf("      key2=%i is out of range [1..500]\n",key2);
879       return 0;
880    }
881    
882    return HIJJET1.pjtm[key2-1][key1-1];
883 }
884
885 //====================== access to common HIJJET1 ===============================
886
887 //______________________________________________________________________________
888 Int_t THijing::GetNSG() const
889 {
890 // Get value of NSG in common HIJJET2
891    return HIJJET2.nsg;
892 }
893
894 //______________________________________________________________________________
895 Int_t THijing::GetNJSG(Int_t key) const
896 {
897 // Get values of array NJSG in common HIJJET2
898    if ( key<1 || key>900 ) {
899       printf ("ERROR in THijing:GetNJSG(key): \n ");
900       printf ("      key=%i is out of range [1..900]!\n",key);
901       return 0;
902    }
903
904    return HIJJET2.njsg[key-1];
905
906 }
907
908 //______________________________________________________________________________
909 Int_t THijing::GetIASG(Int_t key1, Int_t key2) const
910 {
911 // Get values of IASG in common HIJJET2
912    if ( key1<1 || key1>900 ) {
913       printf("ERROR in THijing::GetIASG(key1):\n");
914       printf("      key1=%i is out of range [1..900]\n",key1);
915       return 0;
916    }
917    if ( key2<1 || key2>3 ) {
918       printf("ERROR in THijing::GetIASG(key1,key2):\n");
919       printf("      key2=%i is out of range [1..3]\n",key2);
920       return 0;
921    }
922    
923    return HIJJET2.iasg[key2-1][key1-1];
924 }
925
926 //______________________________________________________________________________
927 Int_t THijing::GetK1SG(Int_t key1, Int_t key2) const
928 {
929 // Get values of K1SG in common HIJJET2
930    if ( key1<1 || key1>900 ) {
931       printf("ERROR in THijing::GetK1SG(key1):\n");
932       printf("      key1=%i is out of range [1..900]\n",key1);
933       return 0;
934    }
935    if ( key2<1 || key2>100 ) {
936       printf("ERROR in THijing::GetK1SG(key1,key2):\n");
937       printf("      key2=%i is out of range [1..100]\n",key2);
938       return 0;
939    }
940    
941    return HIJJET2.k1sg[key2-1][key1-1];
942 }
943
944 //______________________________________________________________________________
945 Int_t THijing::GetK2SG(Int_t key1, Int_t key2) const
946 {
947 // Get values of K2SG in common HIJJET2
948    if ( key1<1 || key1>900 ) {
949       printf("ERROR in THijing::GetK2SG(key1):\n");
950       printf("      key1=%i is out of range [1..900]\n",key1);
951       return 0;
952    }
953    if ( key2<1 || key2>100 ) {
954       printf("ERROR in THijing::GetK2SG(key1,key2):\n");
955       printf("      key2=%i is out of range [1..100]\n",key2);
956       return 0;
957    }
958    
959    return HIJJET2.k2sg[key2-1][key1-1];
960 }
961
962 //______________________________________________________________________________
963 Float_t THijing::GetPXSG(Int_t key1, Int_t key2) const
964 {
965 // Get values of PXSG in common HIJJET2
966    if ( key1<1 || key1>900 ) {
967       printf("ERROR in THijing::GetPXSG(key1):\n");
968       printf("      key1=%i is out of range [1..900]\n",key1);
969       return 0;
970    }
971    if ( key2<1 || key2>100 ) {
972       printf("ERROR in THijing::GetPXSG(key1,key2):\n");
973       printf("      key2=%i is out of range [1..100]\n",key2);
974       return 0;
975    }
976    
977    return HIJJET2.pxsg[key2-1][key1-1];
978 }
979
980 //______________________________________________________________________________
981 Float_t THijing::GetPYSG(Int_t key1, Int_t key2) const
982 {
983 // Get values of PYSG in common HIJJET2
984    if ( key1<1 || key1>900 ) {
985       printf("ERROR in THijing::GetPYSG(key1):\n");
986       printf("      key1=%i is out of range [1..900]\n",key1);
987       return 0;
988    }
989    if ( key2<1 || key2>100 ) {
990       printf("ERROR in THijing::GetPYSG(key1,key2):\n");
991       printf("      key2=%i is out of range [1..100]\n",key2);
992       return 0;
993    }
994    
995    return HIJJET2.pysg[key2-1][key1-1];
996 }
997
998 //______________________________________________________________________________
999 Float_t THijing::GetPZSG(Int_t key1, Int_t key2) const
1000 {
1001 // Get values of PZSG in common HIJJET2
1002    if ( key1<1 || key1>900 ) {
1003       printf("ERROR in THijing::GetPZSG(key1):\n");
1004       printf("      key1=%i is out of range [1..900]\n",key1);
1005       return 0;
1006    }
1007    if ( key2<1 || key2>100 ) {
1008       printf("ERROR in THijing::GetPZSG(key1,key2):\n");
1009       printf("      key2=%i is out of range [1..100]\n",key2);
1010       return 0;
1011    }
1012    
1013    return HIJJET2.pzsg[key2-1][key1-1];
1014 }
1015
1016 //______________________________________________________________________________
1017 Float_t THijing::GetPESG(Int_t key1, Int_t key2) const
1018 {
1019 // Get values of PESG in common HIJJET2
1020    if ( key1<1 || key1>900 ) {
1021       printf("ERROR in THijing::GetPESG(key1):\n");
1022       printf("      key1=%i is out of range [1..900]\n",key1);
1023       return 0;
1024    }
1025    if ( key2<1 || key2>100 ) {
1026       printf("ERROR in THijing::GetPESG(key1,key2):\n");
1027       printf("      key2=%i is out of range [1..100]\n",key2);
1028       return 0;
1029    }
1030    
1031    return HIJJET2.pesg[key2-1][key1-1];
1032 }
1033
1034 //______________________________________________________________________________
1035 Float_t THijing::GetPMSG(Int_t key1, Int_t key2) const
1036 {
1037 // Get values of PMSG in common HIJJET2
1038    if ( key1<1 || key1>900 ) {
1039       printf("ERROR in THijing::GetPMSG(key1):\n");
1040       printf("      key1=%i is out of range [1..900]\n",key1);
1041       return 0;
1042    }
1043    if ( key2<1 || key2>100 ) {
1044       printf("ERROR in THijing::GetPMSG(key1,key2):\n");
1045       printf("      key2=%i is out of range [1..100]\n",key2);
1046       return 0;
1047    }
1048    
1049    return HIJJET2.pmsg[key2-1][key1-1];
1050 }
1051
1052 //====================== access to common HISTRNG ===============================
1053
1054 //______________________________________________________________________________
1055 Int_t THijing::GetNFP(Int_t key1, Int_t key2) const
1056 {
1057 // Get values of array NFP in common HISTRNG
1058    if ( key1<1 || key1>300 ) {
1059       printf("ERROR in THijing::GetNFP(key1):\n");
1060       printf("      key1=%i is out of range [1..300]\n",key1);
1061       return 0;
1062    }
1063    if ( key2<1 || key2>15 ) {
1064       printf("ERROR in THijing::GetNFP(key1,key2):\n");
1065       printf("      key2=%i is out of range [1..15]\n",key2);
1066       return 0;
1067    }
1068    
1069    return HISTRNG.nfp[key2-1][key1-1];
1070 }
1071
1072 //______________________________________________________________________________
1073 Float_t THijing::GetPP(Int_t key1, Int_t key2) const
1074 {
1075 // Get values of array PP in common HISTRNG
1076    if ( key1<1 || key1>300 ) {
1077       printf("ERROR in THijing::GetPP(key1):\n");
1078       printf("      key1=%i is out of range [1..300]\n",key1);
1079       return 0;
1080    }
1081    if ( key2<1 || key2>15 ) {
1082       printf("ERROR in THijing::GetPP(key1,key2):\n");
1083       printf("      key2=%i is out of range [1..15]\n",key2);
1084       return 0;
1085    }
1086    
1087    return HISTRNG.pp[key2-1][key1-1];
1088 }
1089
1090 //______________________________________________________________________________
1091 Int_t THijing::GetNFT(Int_t key1, Int_t key2) const
1092 {
1093 // Get values of array NFT in common HISTRNG
1094    if ( key1<1 || key1>300 ) {
1095       printf("ERROR in THijing::GetNFT(key1):\n");
1096       printf("      key1=%i is out of range [1..300]\n",key1);
1097       return 0;
1098    }
1099    if ( key2<1 || key2>15 ) {
1100       printf("ERROR in THijing::GetNFT(key1,key2):\n");
1101       printf("      key2=%i is out of range [1..15]\n",key2);
1102       return 0;
1103    }
1104    
1105    return HISTRNG.nft[key2-1][key1-1];
1106 }
1107
1108 //______________________________________________________________________________
1109 Float_t THijing::GetPT(Int_t key1, Int_t key2) const
1110 {
1111 // Get values of array PT in common HISTRNG
1112    if ( key1<1 || key1>300 ) {
1113       printf("ERROR in THijing::GetPT(key1):\n");
1114       printf("      key1=%i is out of range [1..300]\n",key1);
1115       return 0;
1116    }
1117    if ( key2<1 || key2>15 ) {
1118       printf("ERROR in THijing::GetPT(key1,key2):\n");
1119       printf("      key2=%i is out of range [1..15]\n",key2);
1120       return 0;
1121    }
1122    
1123    return HISTRNG.pt[key2-1][key1-1];
1124 }
1125
1126 void THijing::SetPARJ(Int_t key, Float_t parm) 
1127 {
1128 // Set values of array PARJ in common HISTRNG
1129     if ( key < 1 || key > 200) {
1130         printf("ERROR in THijing::SetPARJ(key,parm):\n");
1131         printf("      key=%i is out of range [1..200]\n",key);
1132     } else {
1133       LUDAT1_HIJING.parj[key-1] = parm;
1134     }
1135 }
1136
1137
1138 void THijing::SetMSTJ(Int_t key, Int_t parm) 
1139 {
1140 // Set values of array MSTJ in common HISTRNG
1141     if ( key < 1 || key > 200) {
1142         printf("ERROR in THijing::SetMSTJ(key,parm):\n");
1143         printf("      key=%i is out of range [1..200]\n",key);
1144     } else {
1145       LUDAT1_HIJING.mstj[key-1] = parm;
1146     }
1147 }
1148
1149
1150 //====================== access to Hijing subroutines =========================
1151
1152
1153 //______________________________________________________________________________
1154 void THijing::Initialize()
1155 {
1156 //////////////////////////////////////////////////////////////////////////////////
1157 // Calls Hijset with the either default parameters or the ones set by the user  //
1158 // via SetEFRM, SetFRAME, SetPROJ, SetTARG, SetIAP, SetIZP, SetIAT, SetIZT      //
1159 //////////////////////////////////////////////////////////////////////////////////
1160
1161    if ( (!strcmp(fFrame.Data(), "CMS     "  )) &&
1162         (!strcmp(fFrame.Data(), "LAB     "  ))){
1163       printf("WARNING! In THijing:Initialize():\n");
1164       printf(" specified frame=%s is neither CMS or LAB\n",fFrame.Data());
1165       printf(" resetting to default \"CMS\" .");
1166       fFrame="CMS";
1167    }
1168
1169    if ( (!strcmp(fProj.Data(), "A       "     )) &&
1170         (!strcmp(fProj.Data(), "P       "     )) &&
1171         (!strcmp(fProj.Data(), "PBAR    "  ))){
1172       printf("WARNING! In THijing:Initialize():\n");
1173       printf(" specified projectile=%s is neither A, P or PBAR\n",fProj.Data());
1174       printf(" resetting to default \"A\" .");
1175       fProj="A";
1176    }
1177
1178    if ( (!strcmp(fTarg.Data(), "A       "     )) &&
1179         (!strcmp(fTarg.Data(), "P       "     )) &&
1180         (!strcmp(fTarg.Data(), "PBAR    "  ))){
1181       printf("WARNING! In THijing:Initialize():\n");
1182       printf(" specified target=%s is neither A, P or PBAR\n",fTarg.Data());
1183       printf(" resetting to default \"A\" .");
1184       fTarg="A";
1185    }
1186
1187    printf(" %s-%s at %f GeV \n",fProj.Data(),fTarg.Data(),fEfrm);
1188
1189    Hijset(fEfrm,fFrame.Data(),fProj.Data(),fTarg.Data(),fIap,fIzp,fIat,fIzt);
1190
1191    printf(" %s-%s at %f GeV \n",fProj.Data(),fTarg.Data(),fEfrm);
1192 }
1193
1194
1195 //______________________________________________________________________________
1196 void THijing::GenerateEvent()
1197 {
1198 // Generates one event;
1199
1200    Hijing(fFrame.Data(),fBmin,fBmax);
1201
1202 }
1203 //______________________________________________________________________________
1204 void THijing::Hijset(float efrm, const char *frame, const char *proj,  
1205                      const char *targ, int iap, int izp, int iat, int izt)
1206 {
1207 // Call HIJING routine HIJSET passing the parameters in a way accepted by
1208 // Fortran routines                                
1209
1210   int s1 = strlen(frame);
1211   int s2 = strlen(proj);
1212   int s3 = strlen(targ);
1213   printf("s1 = %d s2 = %d s3 = %d\n",s1,s2,s3);
1214 #ifndef WIN32 
1215   hijset(efrm, frame, proj, targ, iap, izp, iat, izt, s1, s2, s3);
1216 #else
1217   hijset(efrm, frame, s1, proj, s2, targ, s3, iap, izp, iat, izt);
1218 #endif
1219 }
1220 //______________________________________________________________________________
1221 void THijing::Hijing(const char *frame, float bmin, float bmax)
1222 {
1223 // Call HIJING routine HIJSET passing the parameters in a way accepted by
1224 // Fortran routines                                
1225
1226   int s1 = strlen(frame);
1227   
1228 #ifndef WIN32 
1229   hijing(frame, bmin, bmax, s1);
1230 #else
1231   hijing(frame, s1, bmin, bmax);
1232 #endif
1233 }
1234
1235
1236 Float_t  THijing::Profile(float b)
1237 {
1238 // Call HIJING routine PROFILE 
1239   return profile(b);
1240 }
1241
1242
1243 void  THijing::Rluget(Int_t lfn, Int_t move)
1244 {
1245 // write seed to file
1246   rluget_hijing(lfn, move);
1247 }
1248
1249
1250 void  THijing::Rluset(Int_t lfn, Int_t move)
1251 {
1252 // read seed from file 
1253   rluset_hijing(lfn, move);
1254 }
1255