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