Debug msg
[u/mrichter/AliRoot.git] / THydjet / THydjet.cxx
1 ////////////////////////////////////////////////////////////////////////////////
2 //                                                                            //
3 // THydjet                                                                    //
4 //                                                                            //
5 // THydjet is an interface class to fortran version of Hydjet event generator //
6 //                                                                            //
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 //      --------------------------------------------------------------        //
14 //                                                                            //
15 //      Igor Lokhtin, SINP MSU, Moscow, RU                                    //
16 //        e-mail: Igor.Lokhtin@cern.ch                                        // 
17 //                                                                            //
18 //      Reference for HYDJET:                                                 // 
19 //      I.P. Lokhtin, A.M. Snigirev,                                          // 
20 //      Eur. Phys. J. C 46 (2006) 211.                                        //
21 //                                                                            //
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.//
28 //                                                                            //
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.//
32 //                                                                             //
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.                   //
36 //                                                                             //
37 //      Reference for JETSET event format:                                     //
38 //      T.Sjostrand, Comput.Phys.Commun. 82 (1994) 74.                         //
39 //                                                                             // 
40 //      --------------------------------------------------------------         //
41 //      Web-page:                                                              //
42 //      http://cern.ch/lokhtin/hydro                                           //
43 //      --------------------------------------------------------------         //
44 //                                                                             //
45 //**************************************************************************** //
46
47 #include "THydjet.h"
48 #include "TObjArray.h"
49 #include "HydCommon.h"
50 #include "TParticle.h"
51 #include "TROOT.h"
52
53 #ifndef WIN32
54 # define pyinit pyinit_
55 # define hydro hydro_
56 # define type_of_call
57 #else
58 # define pyinit PYINIT
59 # define hydro HYDRO
60 # define type_of_call _stdcall
61 #endif
62
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);
66 #ifndef WIN32
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,
69                                      Long_t l_target);
70 #else
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,
74                                     double *win
75                                     );
76 #endif
77
78 #include <TClonesArray.h>
79
80 ClassImp(THydjet)
81
82 THydjet::THydjet() :
83   TGenerator("Hydjet","Hydjet"),
84   fEfrm(5500),
85   fFrame("CMS"),
86   fAw(207),
87   fIfb(0),
88   fBmin(0.),
89   fBmax(1.),
90   fBfix(0.),
91   fNh(20000)
92 {
93 // Default constructor
94 }
95
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,
99                  Int_t nh=20000) :
100   TGenerator("Hydjet","Hydjet"),
101   fEfrm(efrm),
102   fFrame(frame),
103   fAw(aw),
104   fIfb(ifb),
105   fBmin(bmin),
106   fBmax(bmax),
107   fBfix(bfix),
108   fNh(nh)
109 {
110 // THydjet constructor:
111 }
112
113 //______________________________________________________________________________
114 THydjet::~THydjet()
115 {
116 // Destructor
117 }
118
119
120 TObjArray* THydjet::ImportParticles(Option_t *option)
121 {
122 //
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.
130 //
131     fParticles->Clear();
132     Int_t numpart = LUJETS.n;
133     printf("\n THydjet: Hydjet stack contains %d particles.", numpart);
134     Int_t nump = 0;
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
139               nump++;
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]
145               );
146               fParticles->Add(p);
147                  }
148              }
149     }
150     else if(!strcmp(option,"All")) {
151                nump = numpart;
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]
158               );
159               fParticles->Add(p);
160           }
161     }
162     return fParticles;
163 }
164
165 Int_t THydjet::ImportParticles(TClonesArray *particles, Option_t *option)
166 {
167 //
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.
175 //
176   if (particles == 0) return 0;
177   TClonesArray &particlesR = *particles;
178   particlesR.Clear();
179   Int_t numpart = LUJETS.n;
180   printf("\n THydjet: Hydjet stack contains %d particles.", numpart);
181   Int_t nump = 0;
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
186               nump++;
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]
192               );
193                  }
194              }
195     }
196     else if(!strcmp(option,"All")){
197                nump = numpart;
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]
204               );
205           }
206     }
207     return nump;
208 }
209
210 //______________________________________________________________________________
211 void THydjet::SetEfrm(Float_t efrm)
212 {
213 // Set the centre of mass (CMS) or lab-energy (LAB)
214    fEfrm=efrm;
215 }
216 //______________________________________________________________________________
217 void THydjet::SetFrame(const char* frame)
218 {
219 // Set the frame type ("CMS" or "LAB")
220    fFrame=frame;
221 }
222 //______________________________________________________________________________
223 /*void THydjet::SetProj(const char* proj)
224 {
225 // Set the projectile type
226    fProj=proj;
227 }
228 //______________________________________________________________________________
229 void THydjet::SetTarg(const char* targ)
230 {
231 // Set the target type
232    fTarg=targ;
233 }
234 */
235 //______________________________________________________________________________
236 void THydjet::SetAw(Float_t aw)
237 {
238 // Set the projectile-targed atomic number
239    fAw=aw;
240 }
241 //______________________________________________________________________________
242 void THydjet::SetIfb(Int_t ifb)
243 {
244 // flag of type of centrality generation
245    fIfb=ifb;
246
247 //______________________________________________________________________________
248 void THydjet::SetBmin(Float_t bmin)
249 {
250 // set minimum impact parameter in units of nucleus radius RA
251    fBmin=bmin;
252 }
253 //______________________________________________________________________________
254 void THydjet::SetBmax(Float_t bmax)
255 {
256 // set maximum impact parameter in units of nucleus radius RA
257    fBmax=bmax;
258
259 //______________________________________________________________________________
260 void THydjet::SetBfix(Float_t bfix)
261 {
262 // Set fixed impact parameter in units of nucleus radius RA
263    fBfix=bfix;
264
265 //______________________________________________________________________________
266 void THydjet::SetNh(Int_t nh)
267 {
268 // Set mean soft hadron multiplicity in central Pb-Pb collisions
269    fNh=nh;
270
271 //______________________________________________________________________________
272 Float_t THydjet::GetEfrm() const
273 {
274 // Get the centre of mass (CMS) or lab-energy (LAB)
275    return fEfrm;
276 }
277 //______________________________________________________________________________
278 const char* THydjet::GetFrame() const
279 {
280 // Get the frame type ("CMS" or "LAB")
281    return fFrame.Data();
282 }
283 //______________________________________________________________________________
284 /*const char* THydjet::GetProj() const
285 {
286 // Get the projectile type
287    return fProj;
288 }
289 //______________________________________________________________________________
290 const char* THydjet::GetTarg() const
291 {
292 // Set the target type
293    return fTarg;
294 }
295 */
296 //______________________________________________________________________________
297 Float_t THydjet::GetAw() const
298 {
299 // Get the projectile atomic number
300    return fAw;
301 }
302 //______________________________________________________________________________
303 Int_t THydjet::GetIfb() const
304 {
305 // Get flag of type of centrality generation
306    return fIfb;
307 }
308 //______________________________________________________________________________
309 Float_t THydjet::GetBmin() const
310 {
311 // Get minimum impact parameter in units of nucleus radius RA
312    return fBmin;
313
314 //______________________________________________________________________________
315 Float_t THydjet::GetBmax() const
316 {
317 // Get maximum impact parameter in units of nucleus radius RA
318    return fBmax;
319 }
320 //______________________________________________________________________________
321 Float_t THydjet::GetBfix() const
322 {
323 // Get fixed impact parameter in units of nucleus radius RA
324    return fBfix;
325 }
326 //______________________________________________________________________________
327 Int_t THydjet::GetNh() const
328 {
329 // Get mean soft hadron multiplicity in central Pb-Pb collisions
330    return fNh;
331
332
333 //====================== access to common HYFLOW ===============================
334
335 //______________________________________________________________________________
336 const void THydjet::SetYTFL(Float_t value) const
337 {
338    HYFLOW.ytfl=value;
339 }
340
341 //______________________________________________________________________________
342 Float_t THydjet::GetYTFL() const
343 {
344    return HYFLOW.ytfl;
345 }
346
347 //______________________________________________________________________________
348 const void THydjet::SetYLFL(Float_t value) const
349 {
350    HYFLOW.ylfl=value;
351 }
352
353 //______________________________________________________________________________
354 Float_t THydjet::GetYLFL() const
355 {
356    return HYFLOW.ylfl;
357 }
358
359 //______________________________________________________________________________
360 const void THydjet::SetFPART(Float_t value) const
361 {
362    HYFLOW.fpart=value;
363 }
364
365
366 //______________________________________________________________________________
367 Float_t THydjet::GetFPART() const
368 {
369    return HYFLOW.fpart;
370 }
371
372
373 //====================== access to common HYJPAR ===============================
374
375 //______________________________________________________________________________
376 const void THydjet::SetNHSEL(Int_t value) const
377 {
378    HYJPAR.nhsel=value;
379 }
380
381 //______________________________________________________________________________
382 Int_t THydjet::GetNHSEL() const
383 {
384    return HYJPAR.nhsel;
385 }
386
387 //______________________________________________________________________________
388 const void THydjet::SetPTMIN(Float_t value) const
389 {
390    HYJPAR.ptmin=value;
391 }
392
393 //______________________________________________________________________________
394 Float_t THydjet::GetPTMIN() const
395 {
396    return HYJPAR.ptmin;
397 }
398
399 //______________________________________________________________________________
400 const void THydjet::SetNJET(Int_t value) const
401 {
402   HYJPAR.njet=value;
403 }
404
405 //______________________________________________________________________________
406 Int_t  THydjet::GetNJET() const
407 {
408    return HYJPAR.njet;
409 }
410
411 //====================== access to common HYFPAR ===============================
412
413 //______________________________________________________________________________
414 Float_t THydjet::GetBGEN() const
415 {
416    return HYFPAR.bgen;
417 }
418
419 //______________________________________________________________________________
420 Int_t THydjet::GetNBCOL() const
421 {
422    return HYFPAR.nbcol;
423 }
424
425 //______________________________________________________________________________
426 Int_t THydjet::GetNPART() const
427 {
428    return HYFPAR.npart;
429 }
430
431 //______________________________________________________________________________
432 Int_t THydjet::GetNPYT() const
433 {
434    return HYFPAR.npyt;
435 }
436
437 //______________________________________________________________________________
438 Int_t THydjet::GetNHYD() const
439 {
440    return HYFPAR.nhyd;
441 }
442
443
444 //====================== access to common LUJETS ===============================
445
446 //______________________________________________________________________________
447 Int_t THydjet::GetN() const
448 {
449    return LUJETS.n;
450 }
451
452 //______________________________________________________________________________
453 Int_t THydjet::GetK(Int_t key1, Int_t key2) const
454 {
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);
459       return 0;
460    }
461
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);
465       return 0;
466    }
467
468    return   LUJETS.k[key2-1][key1-1];
469 }
470
471 //______________________________________________________________________________
472 Float_t THydjet::GetP(Int_t key1, Int_t key2) const
473 {
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);
478       return 0;
479    }
480
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);
484       return 0;
485    }
486
487    return   LUJETS.p[key2-1][key1-1];
488 }
489
490 //______________________________________________________________________________
491 Float_t THydjet::GetV(Int_t key1, Int_t key2) const
492 {
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);
497       return 0;
498    }
499
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);
503       return 0;
504    }
505
506    return   LUJETS.v[key2-1][key1-1];
507 }
508
509 //====================== access to common HYJETS ===============================
510
511 //______________________________________________________________________________
512 Int_t THydjet::GetNL() const
513 {
514    return HYJETS.nl;
515 }
516
517 //______________________________________________________________________________
518 Int_t THydjet::GetKL(Int_t key1, Int_t key2) const
519 {
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);
524       return 0;
525    }
526
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);
530       return 0;
531    }
532
533    return   HYJETS.kl[key2-1][key1-1];
534 }
535
536 //______________________________________________________________________________
537 Float_t THydjet::GetPL(Int_t key1, Int_t key2) const
538 {
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);
543       return 0;
544    }
545
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);
549       return 0;
550    }
551
552    return   HYJETS.pl[key2-1][key1-1];
553 }
554
555 //______________________________________________________________________________
556 Float_t THydjet::GetVL(Int_t key1, Int_t key2) const
557 {
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);
562       return 0;
563    }
564
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);
568       return 0;
569    }
570
571    return   HYJETS.vl[key2-1][key1-1];
572 }
573
574
575 //====================== access to common PYDAT1 ===============================
576
577 //______________________________________________________________________________
578 void THydjet::SetMSTU(Int_t key, Int_t value)
579 {
580    //Set MSTU in Pythia 
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);
584    }
585    PYDAT1.mstu[key-1] = value;
586 }
587
588 //______________________________________________________________________________
589 void THydjet::SetPARU(Int_t key, Double_t value)
590 {
591    //Set PARU in Pythia 
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);
595    }
596    PYDAT1.paru[key-1] = value;
597 }
598
599 //______________________________________________________________________________
600 void THydjet::SetMSTJ(Int_t key, Int_t value)
601 {
602    //Set MSTJ in Pythia 
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);
606    }
607    PYDAT1.mstj[key-1] = value;
608 }
609
610 //______________________________________________________________________________
611 void THydjet::SetPARJ(Int_t key, Double_t value)
612 {
613    //Set PARJ in Pythia 
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);
617    }
618    PYDAT1.parj[key-1] = value;
619 }
620
621
622 //====================== access to common PYSUBS ===============================
623
624 //______________________________________________________________________________
625 const void THydjet::SetMSEL(Int_t value) const
626 {
627   PYSUBS.msel=value;
628 }
629
630 //______________________________________________________________________________
631 void THydjet::SetCKIN(Int_t key, Double_t value)
632 {
633    //Set CKIN in Pythia 
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);
637    }
638    PYSUBS.ckin[key-1] = value;
639 }
640
641 //====================== access to common PYPARS ===============================
642
643 //______________________________________________________________________________
644 void THydjet::SetMSTP(Int_t key, Int_t value)
645 {
646    //Set MSTP in Pythia 
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);
650    }
651    PYPARS.mstp[key-1] = value;
652 }
653
654
655 //====================== access to Hijing subroutines =========================
656
657
658 //______________________________________________________________________________
659 void THydjet::Initialize()
660 {
661
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\" .");
668       fFrame="CMS";
669    }
670    Int_t nhselflag = GetNHSEL();
671    if(nhselflag != 0) {
672       Double_t lwin = fEfrm;
673       Long_t  s1    = strlen(fFrame);
674       Long_t  s2    = strlen("p");
675       Long_t  s3    = strlen("p");
676 #ifndef WIN32
677      pyinit(fFrame,"p","p",&lwin,s1,s2,s3);
678 #else
679      pyinit(fFrame, s1, "p" , s2, "p", s3, &lwin);
680 #endif
681    }
682 }
683
684
685 //______________________________________________________________________________
686 void THydjet::GenerateEvent()
687 {
688 // Generates one event;
689   float xbmin = fBmin;
690   float xbmax = fBmax;
691   float xbfix = fBfix;
692   float xAw   = fAw;
693   hydro(&xAw,&fIfb,&xbmin,&xbmax,&xbfix,&fNh);
694
695 }
696 //______________________________________________________________________________
697 void THydjet::Hydro()
698 {
699   // Generates one event;
700   float xbmin = fBmin;
701   float xbmax = fBmax;
702   float xbfix = fBfix;
703   float xAw   = fAw;
704   hydro(&xAw,&fIfb,&xbmin,&xbmax,&xbfix,&fNh);
705 }