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