1 ////////////////////////////////////////////////////////////////////////////////
5 // THydjet is an interface class to fortran version of Hydjet event generator //
7 // ------------------------------------------------------------- //
8 // HYDJET, fast MC code to simulate flow effects, jet production //
9 // and jet quenching in heavy ion AA collisions at the LHC //
10 // ------------------------------------------------------------- //
11 // This code is merging HYDRO (flow effects), PYTHIA6.4 (hard jet //
12 // production) and PYQUEN (jet quenching) //
13 // -------------------------------------------------------------- //
15 // Igor Lokhtin, SINP MSU, Moscow, RU //
16 // e-mail: Igor.Lokhtin@cern.ch //
18 // Reference for HYDJET: //
19 // I.P. Lokhtin, A.M. Snigirev, //
20 // Eur. Phys. J. C 46 (2006) 211. //
22 // References for HYDRO: //
23 // N.A.Kruglov, I.P.Lokhtin, L.I.Sarycheva, A.M.Snigirev, //
24 // Z. Phys. C 76 (1997) 99; //
25 // I.P.Lokhtin, L.I.Sarycheva, A.M.Snigirev, //
26 // Phys. Lett. B 537 (2002) 261; //
27 // I.P.Lokhtin, A.M.Snigirev, Preprint SINP MSU 2004-14/753,hep-ph/0312204.//
29 // References for PYQUEN: //
30 // I.P.Lokhtin, A.M.Snigirev, Eur.Phys.J. C16 (2000) 527; //
31 // I.P.Lokhtin, A.M.Snigirev, Preprint SINP MSU 2004-13/752, hep-ph/0406038.//
33 // References for PYTHIA: //
34 // T.Sjostrand et al., Comput.Phys.Commun. 135 (2001) 238; //
35 // T.Sjostrand, S. Mrena and P. Skands, hep-ph/0603175. //
37 // Reference for JETSET event format: //
38 // T.Sjostrand, Comput.Phys.Commun. 82 (1994) 74. //
40 // -------------------------------------------------------------- //
42 // http://cern.ch/lokhtin/hydro //
43 // -------------------------------------------------------------- //
45 //**************************************************************************** //
48 #include "TObjArray.h"
49 #include "HydCommon.h"
50 #include "TParticle.h"
54 # define pyinit pyinit_
58 # define pyinit PYINIT
60 # define type_of_call _stdcall
63 extern "C" void type_of_call hydro(float* A, int* ifb, float* bmin,
64 float* bmax, float* bfix, int* nh);
65 //extern "C" void type_of_call luedit(Int_t& medit);
67 extern "C" void type_of_call pyinit( const char *frame, const char *beam, const char *target,
68 double *win, Long_t l_frame, Long_t l_beam,
71 extern "C" void type_of_call pyinit( const char *frame, Long_t l_frame,
72 const char *beam, Long_t l_beam,
73 const char *target, Long_t l_target,
78 #include <TClonesArray.h>
83 TGenerator("Hydjet","Hydjet"),
93 // Default constructor
96 //______________________________________________________________________________
97 THydjet::THydjet(Float_t efrm, const char *frame="CMS",
98 Float_t aw=207., Int_t ifb=0, Float_t bmin=0, Float_t bmax=1, Float_t bfix=0,
100 TGenerator("Hydjet","Hydjet"),
110 // THydjet constructor:
113 //______________________________________________________________________________
120 TObjArray* THydjet::ImportParticles(Option_t *option)
123 // Default primary creation method. It reads the /LUJETS common block which
124 // has been filled by the GenerateEvent method.
125 // The function loops on the generated particles and store them in
126 // the TClonesArray pointed by the argument particles.
127 // The default action is to store only the stable particles (LUJETS.k[0][i] == 1)
128 // This can be demanded explicitly by setting the option = "Final"
129 // If the option = "All", all the particles are stored.
132 Int_t numpart = LUJETS.n;
133 printf("\n THydjet: Hydjet stack contains %d particles.", numpart);
135 if(!strcmp(option,"") || !strcmp(option,"Final")) {
136 for(Int_t i = 0; i < numpart; i++) {
137 if(LUJETS.k[0][i] == 1) {
138 //Use the common block values for the TParticle constructor
140 TParticle* p = new TParticle(
141 LUJETS.k[1][i], LUJETS.k[0][i] ,
142 LUJETS.k[2][i], -1, LUJETS.k[3][i], LUJETS.k[4][i],
143 LUJETS.p[0][i], LUJETS.p[1][i], LUJETS.p[2][i], LUJETS.p[3][i] ,
144 LUJETS.v[0][i], LUJETS.v[1][i], LUJETS.v[2][i], LUJETS.v[3][i]
150 else if(!strcmp(option,"All")) {
152 for(Int_t i = 0; i < numpart; i++){
153 TParticle* p = new TParticle(
154 LUJETS.k[1][i], LUJETS.k[0][i] ,
155 LUJETS.k[2][i], -1, LUJETS.k[3][i], LUJETS.k[4][i],
156 LUJETS.p[0][i], LUJETS.p[1][i], LUJETS.p[2][i], LUJETS.p[3][i] ,
157 LUJETS.v[0][i], LUJETS.v[1][i], LUJETS.v[2][i], LUJETS.v[3][i]
165 Int_t THydjet::ImportParticles(TClonesArray *particles, Option_t *option)
168 // Default primary creation method. It reads the /LUJETS common block which
169 // has been filled by the GenerateEvent method.
170 // The function loops on the generated particles and store them in
171 // the TClonesArray pointed by the argument particles.
172 // The default action is to store only the stable particles (LUJETS.k[0][i] == 1)
173 // This can be demanded explicitly by setting the option = "Final"
174 // If the option = "All", all the particles are stored.
176 if (particles == 0) return 0;
177 TClonesArray &particlesR = *particles;
179 Int_t numpart = LUJETS.n;
180 printf("\n THydjet: Hydjet stack contains %d particles.", numpart);
182 if(!strcmp(option,"") || !strcmp(option,"Final")) {
183 for(Int_t i = 0; i < numpart; i++) {
184 if(LUJETS.k[0][i] == 1) {
185 //Use the common block values for the TParticle constructor
187 new(particlesR[i]) TParticle(
188 LUJETS.k[1][i], LUJETS.k[0][i] ,
189 LUJETS.k[2][i], -1, LUJETS.k[3][i], LUJETS.k[4][i],
190 LUJETS.p[0][i], LUJETS.p[1][i], LUJETS.p[2][i], LUJETS.p[3][i] ,
191 LUJETS.v[0][i], LUJETS.v[1][i], LUJETS.v[2][i], LUJETS.v[3][i]
196 else if(!strcmp(option,"All")){
198 for(Int_t i = 0; i < numpart; i++){
199 new(particlesR[i]) TParticle(
200 LUJETS.k[1][i], LUJETS.k[0][i] ,
201 LUJETS.k[2][i], -1, LUJETS.k[3][i], LUJETS.k[4][i],
202 LUJETS.p[0][i], LUJETS.p[1][i], LUJETS.p[2][i], LUJETS.p[3][i] ,
203 LUJETS.v[0][i], LUJETS.v[1][i], LUJETS.v[2][i], LUJETS.v[3][i]
210 //______________________________________________________________________________
211 void THydjet::SetEfrm(Float_t efrm)
213 // Set the centre of mass (CMS) or lab-energy (LAB)
216 //______________________________________________________________________________
217 void THydjet::SetFrame(const char* frame)
219 // Set the frame type ("CMS" or "LAB")
222 //______________________________________________________________________________
223 /*void THydjet::SetProj(const char* proj)
225 // Set the projectile type
228 //______________________________________________________________________________
229 void THydjet::SetTarg(const char* targ)
231 // Set the target type
235 //______________________________________________________________________________
236 void THydjet::SetAw(Float_t aw)
238 // Set the projectile-targed atomic number
241 //______________________________________________________________________________
242 void THydjet::SetIfb(Int_t ifb)
244 // flag of type of centrality generation
247 //______________________________________________________________________________
248 void THydjet::SetBmin(Float_t bmin)
250 // set minimum impact parameter in units of nucleus radius RA
253 //______________________________________________________________________________
254 void THydjet::SetBmax(Float_t bmax)
256 // set maximum impact parameter in units of nucleus radius RA
259 //______________________________________________________________________________
260 void THydjet::SetBfix(Float_t bfix)
262 // Set fixed impact parameter in units of nucleus radius RA
265 //______________________________________________________________________________
266 void THydjet::SetNh(Int_t nh)
268 // Set mean soft hadron multiplicity in central Pb-Pb collisions
271 //______________________________________________________________________________
272 Float_t THydjet::GetEfrm() const
274 // Get the centre of mass (CMS) or lab-energy (LAB)
277 //______________________________________________________________________________
278 const char* THydjet::GetFrame() const
280 // Get the frame type ("CMS" or "LAB")
281 return fFrame.Data();
283 //______________________________________________________________________________
284 /*const char* THydjet::GetProj() const
286 // Get the projectile type
289 //______________________________________________________________________________
290 const char* THydjet::GetTarg() const
292 // Set the target type
296 //______________________________________________________________________________
297 Float_t THydjet::GetAw() const
299 // Get the projectile atomic number
302 //______________________________________________________________________________
303 Int_t THydjet::GetIfb() const
305 // Get flag of type of centrality generation
308 //______________________________________________________________________________
309 Float_t THydjet::GetBmin() const
311 // Get minimum impact parameter in units of nucleus radius RA
314 //______________________________________________________________________________
315 Float_t THydjet::GetBmax() const
317 // Get maximum impact parameter in units of nucleus radius RA
320 //______________________________________________________________________________
321 Float_t THydjet::GetBfix() const
323 // Get fixed impact parameter in units of nucleus radius RA
326 //______________________________________________________________________________
327 Int_t THydjet::GetNh() const
329 // Get mean soft hadron multiplicity in central Pb-Pb collisions
333 //====================== access to common HYFLOW ===============================
335 //______________________________________________________________________________
336 const void THydjet::SetYTFL(Float_t value) const
341 //______________________________________________________________________________
342 Float_t THydjet::GetYTFL() const
347 //______________________________________________________________________________
348 const void THydjet::SetYLFL(Float_t value) const
353 //______________________________________________________________________________
354 Float_t THydjet::GetYLFL() const
359 //______________________________________________________________________________
360 const void THydjet::SetFPART(Float_t value) const
366 //______________________________________________________________________________
367 Float_t THydjet::GetFPART() const
373 //====================== access to common HYJPAR ===============================
375 //______________________________________________________________________________
376 const void THydjet::SetNHSEL(Int_t value) const
381 //______________________________________________________________________________
382 Int_t THydjet::GetNHSEL() const
387 //______________________________________________________________________________
388 const void THydjet::SetPTMIN(Float_t value) const
393 //______________________________________________________________________________
394 Float_t THydjet::GetPTMIN() const
399 //______________________________________________________________________________
400 const void THydjet::SetNJET(Int_t value) const
405 //______________________________________________________________________________
406 Int_t THydjet::GetNJET() const
411 //====================== access to common HYFPAR ===============================
413 //______________________________________________________________________________
414 Float_t THydjet::GetBGEN() const
419 //______________________________________________________________________________
420 Int_t THydjet::GetNBCOL() const
425 //______________________________________________________________________________
426 Int_t THydjet::GetNPART() const
431 //______________________________________________________________________________
432 Int_t THydjet::GetNPYT() const
437 //______________________________________________________________________________
438 Int_t THydjet::GetNHYD() const
444 //====================== access to common LUJETS ===============================
446 //______________________________________________________________________________
447 Int_t THydjet::GetN() const
452 //______________________________________________________________________________
453 Int_t THydjet::GetK(Int_t key1, Int_t key2) const
455 // Get Particle codes information
456 if ( key1<1 || key1>150000 ) {
457 printf("ERROR in THydjet::GetK(key1,key2):\n");
458 printf(" key1=%i is out of range [1..150000]\n",key1);
462 if ( key2<1 || key2>5 ) {
463 printf("ERROR in THydjet::GetK(key1,key2):\n");
464 printf(" key2=%i is out of range [1..5]\n",key2);
468 return LUJETS.k[key2-1][key1-1];
471 //______________________________________________________________________________
472 Float_t THydjet::GetP(Int_t key1, Int_t key2) const
474 // Get Particle four momentum and mass
475 if ( key1<1 || key1>150000 ) {
476 printf("ERROR in THydjet::GetP(key1,key2):\n");
477 printf(" key1=%i is out of range [1..150000]\n",key1);
481 if ( key2<1 || key2>5 ) {
482 printf("ERROR in THydjet::GetP(key1,key2):\n");
483 printf(" key2=%i is out of range [1..5]\n",key2);
487 return LUJETS.p[key2-1][key1-1];
490 //______________________________________________________________________________
491 Float_t THydjet::GetV(Int_t key1, Int_t key2) const
493 // Get particle vertex, production time and lifetime
494 if ( key1<1 || key1>150000 ) {
495 printf("ERROR in THydjet::GetV(key1,key2):\n");
496 printf(" key1=%i is out of range [1..150000]\n",key1);
500 if ( key2<1 || key2>5 ) {
501 printf("ERROR in THydjet::GetV(key1,key2):\n");
502 printf(" key2=%i is out of range [1..5]\n",key2);
506 return LUJETS.v[key2-1][key1-1];
509 //====================== access to common HYJETS ===============================
511 //______________________________________________________________________________
512 Int_t THydjet::GetNL() const
517 //______________________________________________________________________________
518 Int_t THydjet::GetKL(Int_t key1, Int_t key2) const
520 // Get Particle codes information
521 if ( key1<1 || key1>150000 ) {
522 printf("ERROR in THydjet::GetKL(key1,key2):\n");
523 printf(" key1=%i is out of range [1..150000]\n",key1);
527 if ( key2<1 || key2>5 ) {
528 printf("ERROR in THydjet::GetKL(key1,key2):\n");
529 printf(" key2=%i is out of range [1..5]\n",key2);
533 return HYJETS.kl[key2-1][key1-1];
536 //______________________________________________________________________________
537 Float_t THydjet::GetPL(Int_t key1, Int_t key2) const
539 // Get Particle four momentum and mass
540 if ( key1<1 || key1>150000 ) {
541 printf("ERROR in THydjet::GetPL(key1,key2):\n");
542 printf(" key1=%i is out of range [1..150000]\n",key1);
546 if ( key2<1 || key2>5 ) {
547 printf("ERROR in THydjet::GetPL(key1,key2):\n");
548 printf(" key2=%i is out of range [1..5]\n",key2);
552 return HYJETS.pl[key2-1][key1-1];
555 //______________________________________________________________________________
556 Float_t THydjet::GetVL(Int_t key1, Int_t key2) const
558 // Get particle vertex, production time and lifetime
559 if ( key1<1 || key1>150000 ) {
560 printf("ERROR in THydjet::GetVL(key1,key2):\n");
561 printf(" key1=%i is out of range [1..150000]\n",key1);
565 if ( key2<1 || key2>5 ) {
566 printf("ERROR in THydjet::GetVL(key1,key2):\n");
567 printf(" key2=%i is out of range [1..5]\n",key2);
571 return HYJETS.vl[key2-1][key1-1];
575 //====================== access to common PYDAT1 ===============================
577 //______________________________________________________________________________
578 void THydjet::SetMSTU(Int_t key, Int_t value)
581 if ( key<1 || key>200 ) {
582 printf("ERROR in THydjet::SetMSTU(key,value):\n");
583 printf(" key=%i is out of range [1..200]\n",key);
585 PYDAT1.mstu[key-1] = value;
588 //______________________________________________________________________________
589 void THydjet::SetPARU(Int_t key, Double_t value)
592 if ( key<1 || key>200 ) {
593 printf("ERROR in THydjet::SetPARU(key,value):\n");
594 printf(" key=%i is out of range [1..200]\n",key);
596 PYDAT1.paru[key-1] = value;
599 //______________________________________________________________________________
600 void THydjet::SetMSTJ(Int_t key, Int_t value)
603 if ( key<1 || key>200 ) {
604 printf("ERROR in THydjet::SetMSTJ(key,value):\n");
605 printf(" key=%i is out of range [1..200]\n",key);
607 PYDAT1.mstj[key-1] = value;
610 //______________________________________________________________________________
611 void THydjet::SetPARJ(Int_t key, Double_t value)
614 if ( key<1 || key>200 ) {
615 printf("ERROR in THydjet::SetPARJ(key,value):\n");
616 printf(" key=%i is out of range [1..200]\n",key);
618 PYDAT1.parj[key-1] = value;
622 //====================== access to common PYSUBS ===============================
624 //______________________________________________________________________________
625 const void THydjet::SetMSEL(Int_t value) const
630 //______________________________________________________________________________
631 void THydjet::SetCKIN(Int_t key, Double_t value)
634 if ( key<1 || key>200 ) {
635 printf("ERROR in THydjet::SetCKIN(key,value):\n");
636 printf(" key=%i is out of range [1..200]\n",key);
638 PYSUBS.ckin[key-1] = value;
641 //====================== access to common PYPARS ===============================
643 //______________________________________________________________________________
644 void THydjet::SetMSTP(Int_t key, Int_t value)
647 if ( key<1 || key>200 ) {
648 printf("ERROR in THydjet::SetMSTP(key,value):\n");
649 printf(" key=%i is out of range [1..200]\n",key);
651 PYPARS.mstp[key-1] = value;
655 //====================== access to Hijing subroutines =========================
658 //______________________________________________________________________________
659 void THydjet::Initialize()
662 // Initialize PYTHIA for hard parton-parton scattering
663 if ( (!strcmp(fFrame.Data(), "CMS " )) &&
664 (!strcmp(fFrame.Data(), "LAB " ))){
665 printf("WARNING! In THydjet:Initialize():\n");
666 printf(" specified frame=%s is neither CMS or LAB\n",fFrame.Data());
667 printf(" resetting to default \"CMS\" .");
670 Int_t nhselflag = GetNHSEL();
672 Double_t lwin = fEfrm;
673 Long_t s1 = strlen(fFrame);
674 Long_t s2 = strlen("p");
675 Long_t s3 = strlen("p");
677 pyinit(fFrame,"p","p",&lwin,s1,s2,s3);
679 pyinit(fFrame, s1, "p" , s2, "p", s3, &lwin);
685 //______________________________________________________________________________
686 void THydjet::GenerateEvent()
688 // Generates one event;
693 hydro(&xAw,&fIfb,&xbmin,&xbmax,&xbfix,&fNh);
696 //______________________________________________________________________________
697 void THydjet::Hydro()
699 // Generates one event;
704 hydro(&xAw,&fIfb,&xbmin,&xbmax,&xbfix,&fNh);