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