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