2726f7c745553c47a88975e7193e931ec0252de0
[u/mrichter/AliRoot.git] / TUHKMgen / TUHKMgen.cxx
1 /////////////////////////////////////////////////////////////////////////////
2 // TUHKM is an interface to 
3 // UHKM 3.0                                                                //
4 // ( only the HYDJET++ part is implemented)                                // 
5 // temporary link:                                                         //
6 // http://lav01.sinp.msu.ru/~igor/hydjet++/hydjet++.txt                    //  
7 // The main UHKM options are accessable  through this interface.           //
8 // Uses the TUHKMgen implementation of TGenerator.                         //
9 // Author of the first implementation: Sergey Zaporozhets                  //
10 // (zaporozh@sunhe.jinr.ru)                                                // 
11 // Futhers modifications were made by                                      //
12 // Ionut Cristian Arsene (i.c.arsene@fys.uio.no)                           //
13 // & Malinina Liudmila(malinina@lav01.sinp.msu.ru) using as an example     //
14 //  AliGenTherminator.cxx created by Adam Kisiel                           //   
15 //                                                                         //
16 ////////////////////////////////////////////////////////////////////////////
17
18 #ifndef TUHKMGEN_H
19 #include "TUHKMgen.h"
20 #endif
21 #include "TObjArray.h"
22 #include "TParticle.h"
23 #include "Particle.h"
24 #include "TROOT.h"
25 #include "TError.h"
26 #include "TClonesArray.h"
27 #include "TSystem.h"
28
29 #include <string>
30 #include <iostream>
31 using namespace std;
32
33 //class TObjArray;
34
35 ClassImp(TUHKMgen)
36
37 TUHKMgen::TUHKMgen() : 
38   TGenerator("UHKM","UHKM"),
39   fInitialState(0x0),
40   fAllocator(),
41   fSecondariesList(),
42   fNPprim(0),
43   fNPsec(0),
44   fHydjetParams(),
45   fStableFlagged(0)
46 {
47   // default constructor setting reasonable defaults for initial parameters (central Pb+Pb at 5.5 TeV)
48
49   //  ParticleAllocator fAllocator;
50   //  List_t fSecondariesList;
51
52   // Set reasonable default values for LHC
53   
54   fHydjetParams.fSqrtS=5500; //LHC
55   fHydjetParams.fAw=207;//Au-Au
56   fHydjetParams.fIfb    = 1;      // flag of type of centrality generation (=0 is fixed by fBfix, not 0 
57                                    //impact parameter is generated in each event between fBfmin 
58                                    //and fBmax according with Glauber model (f-la 30)
59   fHydjetParams.fBfix=0.;
60   fHydjetParams.fBmin=0.;
61   fHydjetParams.fBmax=0.5; //0-5% centrality
62   fHydjetParams.fT = 0.170;
63   fHydjetParams.fMuB = 0.0;
64   fHydjetParams.fMuS = 0.0;
65   fHydjetParams.fMuI3 = 0.0;
66   fHydjetParams.fThFO = 0.130;
67   fHydjetParams.fMu_th_pip = 0.0;
68   fHydjetParams.fSeed=0;
69   fHydjetParams.fTau=10.;
70   fHydjetParams.fSigmaTau=3.;
71   fHydjetParams.fR=11.;
72   fHydjetParams.fYlmax=4.0;
73   fHydjetParams.fUmax=1.1;
74   fHydjetParams.fDelta=0.;
75   fHydjetParams.fEpsilon=0.;
76   fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
77   fHydjetParams.fEtaType=1;//gaus
78   fHydjetParams.fCorrS=1.;
79   fHydjetParams.fNhsel=2;
80   fHydjetParams.fIshad=1;
81   fHydjetParams.fPtmin=7.0;
82   fHydjetParams.fT0=0.8;
83   fHydjetParams.fTau0=0.1;
84   fHydjetParams.fNf=0;
85   fHydjetParams.fIenglu=0;
86   fHydjetParams.fIanglu=0;
87
88
89   // Set reasonable default values for RHIC
90 /*
91   fHydjetParams.fSqrtS=200; //RHIC          
92   fHydjetParams.fAw=197;//Au-Au        
93   fHydjetParams.fIfb    = 1;      
94   fHydjetParams.fBfix=0.;       
95   fHydjetParams.fBmin=0.;       
96   fHydjetParams.fBmax=0.5; //0-5% centrality
97   fHydjetParams.fT = 0.165;
98   fHydjetParams.fMuB = 0.0285;
99   fHydjetParams.fMuS = 0.007;
100   fHydjetParams.fMuI3 = -0.001;
101   fHydjetParams.fThFO = 0.100;
102   fHydjetParams.fMu_th_pip = 0.053;
103   fHydjetParams.fSeed=0;  
104   fHydjetParams.fTau=8.;
105   fHydjetParams.fSigmaTau=2.;
106   fHydjetParams.fR=10.;
107   fHydjetParams.fYlmax=3.3;
108   fHydjetParams.fUmax=1.1;
109   fHydjetParams.fDelta=0.;
110   fHydjetParams.fEpsilon=0.;
111   fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off 
112   fHydjetParams.fEtaType=1;//gaus
113   fHydjetParams.fCorrS=1.;  
114   fHydjetParams.fNhsel=2;
115   fHydjetParams.fIshad=1;
116   fHydjetParams.fPtmin=3.4;
117   fHydjetParams.fT0=0.3;
118   fHydjetParams.fTau0=0.4;
119   fHydjetParams.fNf=2;      
120   fHydjetParams.fIenglu=0;     
121   fHydjetParams.fIanglu=0;  
122 */
123
124   if(strlen(Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")))<255)   
125     strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 256);
126   if(strlen(Form("%s/TUHKMgen/UHKM/tabledecay.txt",gSystem->Getenv("ALICE_ROOT")))<255)
127     strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 256);
128   for(Int_t i=0; i<500; i++) {
129     fStableFlagPDG[i] = 0;
130     fStableFlagStatus[i] = kFALSE;
131   }
132   fStableFlagged = 0;
133   //  cout << "TUHKMgen::TUHKMgen() OUT" << endl;
134 }
135
136 //______________________________________________________________________________
137 TUHKMgen::~TUHKMgen()
138 {
139   // destructor, deletes the InitialStateHydjet object
140   if(fInitialState)
141     delete fInitialState;
142 }
143
144 void TUHKMgen::SetAllParametersRHIC()
145 {
146   // Set reasonable input parameters for 0-5% central Au+Au collisions
147   // at 200 GeV at RHIC
148   SetEcms(200.0);                // RHIC top energy
149   SetAw(197);                    // Au+Au
150   SetIfb(1);
151   SetBfix(0.);
152   SetBmin(0.0);                  // 0%
153   SetBmax(0.5);                  // 5%
154   SetChFrzTemperature(0.165);    // T_ch = 165 MeV
155   SetMuB(0.0285);                // mu_B = 28.5 MeV
156   SetMuS(0.007);                 // mu_S = 7 MeV
157   SetMuQ(-0.001);                // mu_Q = -1 MeV
158   SetThFrzTemperature(0.100);    // T_th = 100 MeV
159   SetMuPionThermal(0.053);       // mu_th_pion = 53 MeV
160   SetSeed(0);                    // use UNIX time
161   SetTauB(8.0);                  // tau = 8 fm/c
162   SetSigmaTau(2.0);              // sigma_tau = 2 fm/c
163   SetRmaxB(10.0);                // fR = 10 fm
164   SetYlMax(3.3);                 // fYmax = 3.3
165   SetEtaRMax(1.1);               // Umax = 1.1
166   SetMomAsymmPar(0.0);           // delta = 0.0
167   SetCoordAsymmPar(0.0);         // epsilon = 0.0
168   SetFlagWeakDecay(0.0);           // weak decay on (<0 off !!!)
169   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
170   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
171   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
172   SetPyquenShad(1);              // shadowing on (0 off)
173   SetPyquenPtmin(3.4);           // ptmin = 3.4 GeV/c
174   SetPyquenT0(0.3);              // T0 = 300 MeV
175   SetPyquenTau0(0.4);            // tau0 = 0.4 fm/c
176   SetPyquenNf(2);                // 2 flavours
177   SetPyquenIenglu(0);            // radiative and collisional energy loss
178   SetPyquenIanglu(0);            // small gluon angular distribution
179 }
180
181 void TUHKMgen::SetAllParametersLHC()
182 {
183   // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
184   SetEcms(5500.0);               // LHC
185   SetAw(207);                    // Pb+Pb
186   SetIfb(1);
187   SetBfix(0.);                  // 0
188   SetBmin(0.0);                  // 0%
189   SetBmax(0.5);                  // 5%
190   SetChFrzTemperature(0.170);    // T_ch = 170 MeV
191   SetMuB(0.0);                   // mu_B = 0 MeV
192   SetMuS(0.0);                   // mu_S = 0 MeV
193   SetMuQ(0.0);                   // mu_Q = 0 MeV
194   SetThFrzTemperature(0.130);    // T_th = 130 MeV
195   SetMuPionThermal(0.0);         // mu_th_pion = 0 MeV
196   SetSeed(0);                    // use UNIX time
197   SetTauB(10.0);                 // tau = 10 fm/c
198   SetSigmaTau(3.0);              // sigma_tau = 3 fm/c
199   SetRmaxB(11.0);                // fR = 11 fm
200   SetYlMax(4.0);                 // fYmax = 4.0
201   SetEtaRMax(1.1);               // Umax = 1.1
202   SetMomAsymmPar(0.0);           // delta = 0.0
203   SetCoordAsymmPar(0.0);         // epsilon = 0.0
204   SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
205   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
206   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
207   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
208   SetPyquenShad(1);              // shadowing on (0 off)
209   SetPyquenPtmin(7.0);           // ptmin = 7.0 GeV/c
210   SetPyquenT0(0.8);              // T0 = 800 MeV
211   SetPyquenTau0(0.1);            // tau0 = 0.4 fm/c
212   SetPyquenNf(0);                // 0 flavours
213   SetPyquenIenglu(0);            // radiative and collisional energy loss
214   SetPyquenIanglu(0);            // small gluon angular distribution
215 }
216
217 TObjArray* TUHKMgen::ImportParticles(const Option_t *)
218 {
219   // Function overloading the TGenerator::ImportParticles() member function.
220   // The particles from the local particle list (fSecondariesList) are
221   // forwarded to the TGenerator::fParticles 
222   fParticles->Clear();
223   Int_t nump = 0;
224   LPIT_t it,e;
225    
226   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
227     TVector3 pos(it->Pos().Vect());
228     TVector3 mom(it->Mom().Vect());
229     Float_t m1 = it->TableMass();
230     Int_t im1 = it->GetMother();
231     Int_t im2 = -1;
232     Int_t id1 = -1;
233     Int_t id2 = -1;
234
235     Int_t type = it->GetType();  // 0-hydro, 1-jets
236
237     if (im1> -1) {
238       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1+1));           
239       mother->SetLastDaughter(nump);
240       if (mother->GetFirstDaughter()==-1)
241         mother->SetFirstDaughter(nump+1);
242     }
243
244     nump++;
245
246     TParticle* p = new TParticle(it->Encoding(), type,                          //pdg,stat
247                                  im1, im2, id1, id2,                                            //m1,m2,d1,d2
248                                  mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1),     //px,py,pz,e
249                                  pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13, 
250                                  it->T()*1.e-13/3e10);                          //x,y,z,t
251   
252      p->SetUniqueID(nump);
253      fParticles->Add(p);
254   } //end for
255  
256   fAllocator.FreeList(fSecondariesList);
257
258   return fParticles;
259 }
260
261 Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
262 {
263   // Function overloading the TGenerator::ImportParticles() member function.
264   // The particles from the local particle list (fSecondariesList) are
265   // forwarded to the TGenerator::fParticles
266   option = option;   // just to avoid the warning
267   
268
269   if(particles==0) return 0;
270   TClonesArray &particlesR=*particles;
271   particlesR.Clear();
272  
273   Int_t numprim,numsec;  numprim=numsec=0;
274   Int_t nump = 0;
275   LPIT_t it,e;
276   
277   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
278     TVector3 pos(it->Pos().Vect());  
279     TVector3 mom(it->Mom().Vect());  
280     Float_t m1 = it->TableMass();  
281     Int_t im1 = it->GetMother();  
282     Int_t im2 = -1;
283     Int_t id1 = -1;
284     Int_t id2 = -1;
285     
286     Int_t type = it->GetType();  
287       
288     if (im1> -1) {
289       // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
290       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
291       mother->SetLastDaughter(nump);
292       if(mother->GetFirstDaughter()==-1)
293         mother->SetFirstDaughter(nump);
294     } else 
295       ++numprim;
296     
297     new (particlesR[nump]) TParticle(it->Encoding(), type,                                              //pdg,stat
298                                      im1, im2, id1, id2,                                                //m1,m2,d1,d2
299                                      mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
300                                      pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13, 
301                                      it->T()*1.e-13/3e10);                      //x,y,z,t
302     
303     particlesR[nump]->SetUniqueID(nump);
304     nump++;
305     numsec++;
306   }//end for
307   
308   fSecondariesList.clear();
309   printf("Scan and add prim %d sec %d and all %d particles\n",
310          numprim,numsec,nump);
311   return nump;
312 }
313
314 //______________________________________________________________________________
315 void TUHKMgen::Initialize()
316 {
317   // Function overloading the TGenerator::Initialize() member function.
318   // The Monte-Carlo model is initialized (input parameters are transmitted, 
319   // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
320  
321   fInitialState = new InitialStateHydjet();  
322   SetAllParameters();  
323   fInitialState->LoadPDGInfo();  
324   // set the stable flags
325   for(Int_t i=0; i<fStableFlagged; i++) 
326     fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]); 
327
328   if(!fInitialState->MultIni()) 
329     Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n");  //
330 }
331
332 void TUHKMgen::Print(const Option_t*) const
333 {
334   cout << "**********************************************************************************" << endl;
335   cout << "* UHKM Generator interface to ROOT::TGenerator                                   *" << endl;
336   cout << "*  Documentation:                                                                *" << endl;
337   cout << "* I.P.Lokhtin, L.V.Malinina, S.V.Petrushanko, A.M.Snigirev, I.Arsene, K.Tywoniuk *" << endl;
338   cout << "*    Comput.Phys.Commun.180:779-799, 2009                                        *" << endl;
339   cout << "**********************************************************************************" << endl;
340 }
341
342 void TUHKMgen::GenerateEvent()
343 {
344   // Member function overloading the TGenerator::GenerateEvent()
345   // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
346   
347   fInitialState->Initialize(fSecondariesList, fAllocator);  
348  
349   if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");  //
350
351   // Run the decays
352
353   if(fInitialState->RunDecays())  
354     fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());  
355   
356 }
357
358 void TUHKMgen::SetAllParameters() {
359   // forward all model input parameters to the InitialStateHydjet object
360   // which will handle all the Monte-Carlo simulation using HYDJET++ model
361
362   fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
363   fInitialState->fParams.fAw = fHydjetParams.fAw;
364   fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
365   fInitialState->fParams.fBfix = fHydjetParams.fBfix;
366   fInitialState->fParams.fBmin = fHydjetParams.fBmin;
367   fInitialState->fParams.fBmax = fHydjetParams.fBmax;
368   fInitialState->fParams.fSeed = fHydjetParams.fSeed;
369   fInitialState->fParams.fT = fHydjetParams.fT;
370   fInitialState->fParams.fMuB = fHydjetParams.fMuB;
371   fInitialState->fParams.fMuS = fHydjetParams.fMuS;
372   fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
373   fInitialState->fParams.fThFO = fHydjetParams.fThFO;
374   fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
375
376   fInitialState->fParams.fTau = fHydjetParams.fTau;
377   fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
378   fInitialState->fParams.fR = fHydjetParams.fR;
379   fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
380   fInitialState->fParams.fUmax = fHydjetParams.fUmax;
381   fInitialState->fParams.fDelta = fHydjetParams.fDelta;
382   fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
383
384   fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
385   fInitialState->fParams.fDecay = 1;                   // run decays
386   fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
387   fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
388  
389   fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
390   fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
391   fInitialState->fParams.fIshad = fHydjetParams.fIshad;
392   fInitialState->fParams.fT0 = fHydjetParams.fT0;
393   fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
394   fInitialState->fParams.fNf = fHydjetParams.fNf;      
395   fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;     
396   fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;  
397
398   fInitialState->SetPDGParticleFilename(fParticleFilename);
399   fInitialState->SetPDGDecayFilename(fDecayFilename);
400   //  fInitialState->SetUseCharmParticles(fUseCharmParticles);
401   //  fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
402   //  fInitialState->SetMassRange(fMinMass, fMaxMass);
403   //  for(Int_t i=0; i<fStableFlagged; i++) 
404   //    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
405   //  cout << "TUHKMgen::SetAllParameters() OUT" << endl;
406 }