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