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