]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TAmpt/TAmpt.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
add option to fill or not histograms in eta gaps, and momentum imbalance in associate...
[u/mrichter/AliRoot.git] / TAmpt / TAmpt.cxx
... / ...
CommitLineData
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
69//______________________________________________________________________________
70TAmpt::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//______________________________________________________________________________
89TAmpt::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//______________________________________________________________________________
114TAmpt::~TAmpt()
115{
116 // Destructor
117}
118
119//______________________________________________________________________________
120TObjArray* 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//______________________________________________________________________________
177Int_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//______________________________________________________________________________
240Int_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//______________________________________________________________________________
310void TAmpt::SetEFRM(Float_t efrm)
311{
312 // Set the centre of mass (CMS) or lab-energy (LAB)
313 fEfrm=efrm;
314}
315
316//______________________________________________________________________________
317void TAmpt::SetFRAME(const char* frame)
318{
319 // Set the frame type ("CMS" or "LAB")
320 fFrame=frame;
321}
322
323//______________________________________________________________________________
324void TAmpt::SetPROJ(const char* proj)
325{
326 // Set the projectile type
327 fProj=proj;
328}
329
330//______________________________________________________________________________
331void TAmpt::SetTARG(const char* targ)
332{
333 // Set the target type
334 fTarg=targ;
335}
336
337//______________________________________________________________________________
338void TAmpt::SetIAP(Int_t iap)
339{
340 // Set the projectile atomic number
341 fIap=iap;
342}
343
344//______________________________________________________________________________
345void TAmpt::SetIZP(Int_t izp)
346{
347 // Set the projectile charge number
348 fIzp=izp;
349}
350
351//______________________________________________________________________________
352void TAmpt::SetIAT(Int_t iat)
353{
354 // Set the target atomic number
355 fIat=iat;
356}
357
358//______________________________________________________________________________
359void TAmpt::SetIZT(Int_t izt)
360{
361 // Set the target charge number
362 fIzt=izt;
363}
364
365//______________________________________________________________________________
366void TAmpt::SetBMIN(Float_t bmin)
367{
368 // Set the minimum impact parameter
369 fBmin=bmin;
370}
371
372//______________________________________________________________________________
373void TAmpt::SetBMAX(Float_t bmax)
374{
375 // Set the maximum impact parameter
376 fBmax=bmax;
377}
378
379//______________________________________________________________________________
380Float_t TAmpt::GetEFRM() const
381{
382 // Get the centre of mass (CMS) or lab-energy (LAB)
383 return fEfrm;
384}
385
386//______________________________________________________________________________
387const char* TAmpt::GetFRAME() const
388{
389 // Get the frame type ("CMS" or "LAB")
390 return fFrame.Data();
391}
392
393//______________________________________________________________________________
394const char* TAmpt::GetPROJ() const
395{
396 // Get the projectile type
397 return fProj;
398}
399
400//______________________________________________________________________________
401const char* TAmpt::GetTARG() const
402{
403 // Set the target type
404 return fTarg;
405}
406
407//______________________________________________________________________________
408Int_t TAmpt::GetIAP() const
409{
410 // Get the projectile atomic number
411 return fIap;
412}
413
414//______________________________________________________________________________
415Int_t TAmpt::GetIZP() const
416{
417 // Get the projectile charge number
418 return fIzp;
419}
420
421//______________________________________________________________________________
422Int_t TAmpt::GetIAT() const
423{
424 // Get the target atomic number
425 return fIat;
426}
427
428//______________________________________________________________________________
429Int_t TAmpt::GetIZT() const
430{
431 // Get the target charge number
432 return fIzt;
433}
434
435//______________________________________________________________________________
436Float_t TAmpt::GetBMIN() const
437{
438 // Get the minimum impact parameter
439 return fBmin;
440}
441
442//______________________________________________________________________________
443Float_t TAmpt::GetBMAX() const
444{
445 // Get the maximum impact parameter
446 return fBmax;
447}
448
449//====================== access to common HIPARNT ===============================
450
451//______________________________________________________________________________
452void 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//______________________________________________________________________________
464Float_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//______________________________________________________________________________
476void 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//______________________________________________________________________________
488Int_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//______________________________________________________________________________
500Float_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//______________________________________________________________________________
512Int_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//______________________________________________________________________________
526Int_t TAmpt::GetNATT() const
527{
528 // Get the number of particles produces
529 return HMAIN1.natt;
530}
531
532//______________________________________________________________________________
533Float_t TAmpt::GetEATT() const
534{
535 // Get total energy of particles
536 return HMAIN1.eatt;
537}
538
539//______________________________________________________________________________
540Int_t TAmpt::GetJATT() const
541{
542 // Get number of hard scatterings
543 return HMAIN1.jatt;
544}
545
546//______________________________________________________________________________
547Int_t TAmpt::GetNT() const
548{
549 // Get number of target participants
550 return HMAIN1.nt;
551}
552
553//______________________________________________________________________________
554Int_t TAmpt::GetNP() const
555{
556 // Get number of projectile participants
557 return HMAIN1.np;
558}
559
560//______________________________________________________________________________
561Int_t TAmpt::GetN0() const
562{
563 // Get number of N-N collisions
564 return HMAIN1.n0;
565}
566
567//______________________________________________________________________________
568Int_t TAmpt::GetN01() const
569{
570 // Get number of N-wounded collisions
571 return HMAIN1.n01;
572}
573
574//______________________________________________________________________________
575Int_t TAmpt::GetN10() const
576{
577 // Get number of wounded-N collisions
578 return HMAIN1.n10;
579}
580
581//______________________________________________________________________________
582Int_t TAmpt::GetN11() const
583{
584 // Get number of wounded-wounded collisions
585 return HMAIN1.n11;
586}
587
588//______________________________________________________________________________
589Float_t TAmpt::GetBB() const
590{
591 // Get impact parameter
592
593 return HPARNT.hint1[18];
594}
595
596//====================== access to common HMAIN2 ===============================
597
598//______________________________________________________________________________
599Int_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//______________________________________________________________________________
616Float_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//______________________________________________________________________________
633Float_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//______________________________________________________________________________
652Int_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//______________________________________________________________________________
664Int_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//______________________________________________________________________________
681Float_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//______________________________________________________________________________
698Float_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//______________________________________________________________________________
715Float_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//______________________________________________________________________________
732Float_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//______________________________________________________________________________
749Float_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//______________________________________________________________________________
766Int_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//______________________________________________________________________________
778Int_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//______________________________________________________________________________
795Float_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//______________________________________________________________________________
812Float_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//______________________________________________________________________________
829Float_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//______________________________________________________________________________
846Float_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//______________________________________________________________________________
863Float_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//______________________________________________________________________________
882Int_t TAmpt::GetNSG() const
883{
884 // Get value of NSG in common HIJJET2
885 return HJJET2.nsg;
886}
887
888//______________________________________________________________________________
889Int_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//______________________________________________________________________________
901Int_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//______________________________________________________________________________
918Int_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//______________________________________________________________________________
935Int_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//______________________________________________________________________________
952Float_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//______________________________________________________________________________
969Float_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//______________________________________________________________________________
986Float_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//______________________________________________________________________________
1003Float_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//______________________________________________________________________________
1020Float_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//______________________________________________________________________________
1039Int_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//______________________________________________________________________________
1056Float_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//______________________________________________________________________________
1073Int_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//______________________________________________________________________________
1090Float_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
1106void 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
1118void 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//______________________________________________________________________________
1132void 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//______________________________________________________________________________
1168void TAmpt::GenerateEvent()
1169{
1170 // Generates one event;
1171
1172 //printf("Next event ------------------------\n");
1173 Ampt(fFrame.Data(),fBmin,fBmax);
1174}
1175
1176//______________________________________________________________________________
1177void 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//______________________________________________________________________________
1195void 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
1209Float_t TAmpt::Profile(float b)
1210{
1211 // Call AMPT routine PROFILE
1212 return profile(b);
1213}
1214
1215//______________________________________________________________________________
1216void TAmpt::Rluget(Int_t lfn, Int_t move)
1217{
1218 // write seed to file
1219 rluget_ampt(lfn, move);
1220}
1221
1222//______________________________________________________________________________
1223void TAmpt::Rluset(Int_t lfn, Int_t move)
1224{
1225 // read seed from file
1226 rluset_ampt(lfn, move);
1227}
1228
1229//______________________________________________________________________________
1230void TAmpt::SetIsoft(Int_t i)
1231{
1232 // set ampt mode.
1233 ANIM.isoft = i;
1234}
1235
1236//______________________________________________________________________________
1237void TAmpt::SetNtMax(Int_t max)
1238{
1239 // set maximum number of timesteps
1240 INPUT2.ntmax = max;
1241}
1242
1243//______________________________________________________________________________
1244void TAmpt::SetIpop(Int_t pop)
1245{
1246 // set flag for popcorn mechanism(netbaryon stopping)
1247 POPCORN.ipop = pop;
1248}
1249
1250//______________________________________________________________________________
1251void TAmpt::SetXmu(Float_t m)
1252{
1253 // set parton screening mass in fm^-1.
1254 PARA2.xmu = m;
1255}
1256
1257//______________________________________________________________________________
1258void TAmpt::SetAlpha(Float_t alpha)
1259{
1260 // set parton screening mass in fm^-1.
1261 PARA2.alpha = alpha;
1262}
1263
1264//______________________________________________________________________________
1265void 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}