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