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