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