]> git.uio.no Git - u/mrichter/AliRoot.git/blame - THydjet/THydjet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Disable retireval of DCS data points from AliShuttle for SDD
[u/mrichter/AliRoot.git] / THydjet / THydjet.cxx
CommitLineData
cb220f83 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
63extern "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
67extern "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
71extern "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
80ClassImp(THydjet)
81
82THydjet::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//______________________________________________________________________________
97THydjet::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//______________________________________________________________________________
114THydjet::~THydjet()
115{
116// Destructor
117}
118
119
120TObjArray* 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
165Int_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//______________________________________________________________________________
211void THydjet::SetEfrm(Float_t efrm)
212{
213// Set the centre of mass (CMS) or lab-energy (LAB)
214 fEfrm=efrm;
215}
216//______________________________________________________________________________
217void 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//______________________________________________________________________________
229void THydjet::SetTarg(const char* targ)
230{
231// Set the target type
232 fTarg=targ;
233}
234*/
235//______________________________________________________________________________
236void THydjet::SetAw(Float_t aw)
237{
238// Set the projectile-targed atomic number
239 fAw=aw;
240}
241//______________________________________________________________________________
242void THydjet::SetIfb(Int_t ifb)
243{
244// flag of type of centrality generation
245 fIfb=ifb;
246}
247//______________________________________________________________________________
248void THydjet::SetBmin(Float_t bmin)
249{
250// set minimum impact parameter in units of nucleus radius RA
251 fBmin=bmin;
252}
253//______________________________________________________________________________
254void THydjet::SetBmax(Float_t bmax)
255{
256// set maximum impact parameter in units of nucleus radius RA
257 fBmax=bmax;
258}
259//______________________________________________________________________________
260void THydjet::SetBfix(Float_t bfix)
261{
262// Set fixed impact parameter in units of nucleus radius RA
263 fBfix=bfix;
264}
265//______________________________________________________________________________
266void THydjet::SetNh(Int_t nh)
267{
268// Set mean soft hadron multiplicity in central Pb-Pb collisions
269 fNh=nh;
270}
271//______________________________________________________________________________
272Float_t THydjet::GetEfrm() const
273{
274// Get the centre of mass (CMS) or lab-energy (LAB)
275 return fEfrm;
276}
277//______________________________________________________________________________
278const 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//______________________________________________________________________________
290const char* THydjet::GetTarg() const
291{
292// Set the target type
293 return fTarg;
294}
295*/
296//______________________________________________________________________________
297Float_t THydjet::GetAw() const
298{
299// Get the projectile atomic number
300 return fAw;
301}
302//______________________________________________________________________________
303Int_t THydjet::GetIfb() const
304{
305// Get flag of type of centrality generation
306 return fIfb;
307}
308//______________________________________________________________________________
309Float_t THydjet::GetBmin() const
310{
311// Get minimum impact parameter in units of nucleus radius RA
312 return fBmin;
313}
314//______________________________________________________________________________
315Float_t THydjet::GetBmax() const
316{
317// Get maximum impact parameter in units of nucleus radius RA
318 return fBmax;
319}
320//______________________________________________________________________________
321Float_t THydjet::GetBfix() const
322{
323// Get fixed impact parameter in units of nucleus radius RA
324 return fBfix;
325}
326//______________________________________________________________________________
327Int_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//______________________________________________________________________________
336const void THydjet::SetYTFL(Float_t value) const
337{
338 HYFLOW.ytfl=value;
339}
340
341//______________________________________________________________________________
342Float_t THydjet::GetYTFL() const
343{
344 return HYFLOW.ytfl;
345}
346
347//______________________________________________________________________________
348const void THydjet::SetYLFL(Float_t value) const
349{
350 HYFLOW.ylfl=value;
351}
352
353//______________________________________________________________________________
354Float_t THydjet::GetYLFL() const
355{
356 return HYFLOW.ylfl;
357}
358
359//______________________________________________________________________________
360const void THydjet::SetFPART(Float_t value) const
361{
362 HYFLOW.fpart=value;
363}
364
365
366//______________________________________________________________________________
367Float_t THydjet::GetFPART() const
368{
369 return HYFLOW.fpart;
370}
371
372
373//====================== access to common HYJPAR ===============================
374
375//______________________________________________________________________________
376const void THydjet::SetNHSEL(Int_t value) const
377{
378 HYJPAR.nhsel=value;
379}
380
381//______________________________________________________________________________
382Int_t THydjet::GetNHSEL() const
383{
384 return HYJPAR.nhsel;
385}
386
387//______________________________________________________________________________
388const void THydjet::SetPTMIN(Float_t value) const
389{
390 HYJPAR.ptmin=value;
391}
392
393//______________________________________________________________________________
394Float_t THydjet::GetPTMIN() const
395{
396 return HYJPAR.ptmin;
397}
398
399//______________________________________________________________________________
400const void THydjet::SetNJET(Int_t value) const
401{
402 HYJPAR.njet=value;
403}
404
405//______________________________________________________________________________
406Int_t THydjet::GetNJET() const
407{
408 return HYJPAR.njet;
409}
410
411//====================== access to common HYFPAR ===============================
412
413//______________________________________________________________________________
414Float_t THydjet::GetBGEN() const
415{
416 return HYFPAR.bgen;
417}
418
419//______________________________________________________________________________
420Int_t THydjet::GetNBCOL() const
421{
422 return HYFPAR.nbcol;
423}
424
425//______________________________________________________________________________
426Int_t THydjet::GetNPART() const
427{
428 return HYFPAR.npart;
429}
430
431//______________________________________________________________________________
432Int_t THydjet::GetNPYT() const
433{
434 return HYFPAR.npyt;
435}
436
437//______________________________________________________________________________
438Int_t THydjet::GetNHYD() const
439{
440 return HYFPAR.nhyd;
441}
442
443
444//====================== access to common LUJETS ===============================
445
446//______________________________________________________________________________
447Int_t THydjet::GetN() const
448{
449 return LUJETS.n;
450}
451
452//______________________________________________________________________________
453Int_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//______________________________________________________________________________
472Float_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//______________________________________________________________________________
491Float_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//______________________________________________________________________________
512Int_t THydjet::GetNL() const
513{
514 return HYJETS.nl;
515}
516
517//______________________________________________________________________________
518Int_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//______________________________________________________________________________
537Float_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//______________________________________________________________________________
556Float_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//______________________________________________________________________________
578void 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//______________________________________________________________________________
589void 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//______________________________________________________________________________
600void 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//______________________________________________________________________________
611void 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//______________________________________________________________________________
625const void THydjet::SetMSEL(Int_t value) const
626{
627 PYSUBS.msel=value;
628}
629
630//______________________________________________________________________________
631void 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//______________________________________________________________________________
644void 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//______________________________________________________________________________
659void 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//______________________________________________________________________________
686void 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//______________________________________________________________________________
697void 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}