Adding a reminder for coders
[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   strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 255);
125   strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 255);
126   for(Int_t i=0; i<500; i++) {
127     fStableFlagPDG[i] = 0;
128     fStableFlagStatus[i] = kFALSE;
129   }
130   fStableFlagged = 0;
131   //  cout << "TUHKMgen::TUHKMgen() OUT" << endl;
132 }
133
134 //______________________________________________________________________________
135 TUHKMgen::~TUHKMgen()
136 {
137   // destructor, deletes the InitialStateHydjet object
138   if(fInitialState)
139     delete fInitialState;
140 }
141
142 void TUHKMgen::SetAllParametersRHIC()
143 {
144   // Set reasonable input parameters for 0-5% central Au+Au collisions
145   // at 200 GeV at RHIC
146   SetEcms(200.0);                // RHIC top energy
147   SetAw(197);                    // Au+Au
148   SetIfb(1);
149   SetBfix(0.);
150   SetBmin(0.0);                  // 0%
151   SetBmax(0.5);                  // 5%
152   SetChFrzTemperature(0.165);    // T_ch = 165 MeV
153   SetMuB(0.0285);                // mu_B = 28.5 MeV
154   SetMuS(0.007);                 // mu_S = 7 MeV
155   SetMuQ(-0.001);                // mu_Q = -1 MeV
156   SetThFrzTemperature(0.100);    // T_th = 100 MeV
157   SetMuPionThermal(0.053);       // mu_th_pion = 53 MeV
158   SetSeed(0);                    // use UNIX time
159   SetTauB(8.0);                  // tau = 8 fm/c
160   SetSigmaTau(2.0);              // sigma_tau = 2 fm/c
161   SetRmaxB(10.0);                // fR = 10 fm
162   SetYlMax(3.3);                 // fYmax = 3.3
163   SetEtaRMax(1.1);               // Umax = 1.1
164   SetMomAsymmPar(0.0);           // delta = 0.0
165   SetCoordAsymmPar(0.0);         // epsilon = 0.0
166   SetFlagWeakDecay(0.0);           // weak decay on (<0 off !!!)
167   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
168   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
169   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
170   SetPyquenShad(1);              // shadowing on (0 off)
171   SetPyquenPtmin(3.4);           // ptmin = 3.4 GeV/c
172   SetPyquenT0(0.3);              // T0 = 300 MeV
173   SetPyquenTau0(0.4);            // tau0 = 0.4 fm/c
174   SetPyquenNf(2);                // 2 flavours
175   SetPyquenIenglu(0);            // radiative and collisional energy loss
176   SetPyquenIanglu(0);            // small gluon angular distribution
177 }
178
179 void TUHKMgen::SetAllParametersLHC()
180 {
181   // Set reasonable input parameters for 0-5% Pb+Pb collisions at 5.5 TeV at LHC
182   SetEcms(5500.0);               // LHC
183   SetAw(207);                    // Pb+Pb
184   SetIfb(1);
185   SetBfix(0.);                  // 0
186   SetBmin(0.0);                  // 0%
187   SetBmax(0.5);                  // 5%
188   SetChFrzTemperature(0.170);    // T_ch = 170 MeV
189   SetMuB(0.0);                   // mu_B = 0 MeV
190   SetMuS(0.0);                   // mu_S = 0 MeV
191   SetMuQ(0.0);                   // mu_Q = 0 MeV
192   SetThFrzTemperature(0.130);    // T_th = 130 MeV
193   SetMuPionThermal(0.0);         // mu_th_pion = 0 MeV
194   SetSeed(0);                    // use UNIX time
195   SetTauB(10.0);                 // tau = 10 fm/c
196   SetSigmaTau(3.0);              // sigma_tau = 3 fm/c
197   SetRmaxB(11.0);                // fR = 11 fm
198   SetYlMax(4.0);                 // fYmax = 4.0
199   SetEtaRMax(1.1);               // Umax = 1.1
200   SetMomAsymmPar(0.0);           // delta = 0.0
201   SetCoordAsymmPar(0.0);         // epsilon = 0.0
202   SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
203   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
204   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
205   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
206   SetPyquenShad(1);              // shadowing on (0 off)
207   SetPyquenPtmin(7.0);           // ptmin = 7.0 GeV/c
208   SetPyquenT0(0.8);              // T0 = 800 MeV
209   SetPyquenTau0(0.1);            // tau0 = 0.4 fm/c
210   SetPyquenNf(0);                // 0 flavours
211   SetPyquenIenglu(0);            // radiative and collisional energy loss
212   SetPyquenIanglu(0);            // small gluon angular distribution
213 }
214
215 TObjArray* TUHKMgen::ImportParticles(const Option_t *)
216 {
217   // Function overloading the TGenerator::ImportParticles() member function.
218   // The particles from the local particle list (fSecondariesList) are
219   // forwarded to the TGenerator::fParticles 
220   fParticles->Clear();
221   Int_t nump = 0;
222   LPIT_t it,e;
223    
224   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
225     TVector3 pos(it->Pos().Vect());
226     TVector3 mom(it->Mom().Vect());
227     Float_t m1 = it->TableMass();
228     Int_t im1 = it->GetMother();
229     Int_t im2 = -1;
230     Int_t id1 = -1;
231     Int_t id2 = -1;
232
233     Int_t type = it->GetType();  // 0-hydro, 1-jets
234
235     if (im1> -1) {
236       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1+1));           
237       mother->SetLastDaughter(nump);
238       if (mother->GetFirstDaughter()==-1)
239         mother->SetFirstDaughter(nump+1);
240     }
241
242     nump++;
243
244     TParticle* p = new TParticle(it->Encoding(), type,                          //pdg,stat
245                                  im1, im2, id1, id2,                                            //m1,m2,d1,d2
246                                  mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1),     //px,py,pz,e
247                                  pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13, 
248                                  it->T()*1.e-13/3e10);                          //x,y,z,t
249   
250      p->SetUniqueID(nump);
251      fParticles->Add(p);
252   } //end for
253  
254   fAllocator.FreeList(fSecondariesList);
255
256   return fParticles;
257 }
258
259 Int_t TUHKMgen::ImportParticles(TClonesArray *particles, const Option_t* option)
260 {
261   // Function overloading the TGenerator::ImportParticles() member function.
262   // The particles from the local particle list (fSecondariesList) are
263   // forwarded to the TGenerator::fParticles
264   option = option;   // just to avoid the warning
265   
266
267   if(particles==0) return 0;
268   TClonesArray &particlesR=*particles;
269   particlesR.Clear();
270  
271   Int_t numprim,numsec;  numprim=numsec=0;
272   Int_t nump = 0;
273   LPIT_t it,e;
274   
275   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; it++) {
276     TVector3 pos(it->Pos().Vect());  
277     TVector3 mom(it->Mom().Vect());  
278     Float_t m1 = it->TableMass();  
279     Int_t im1 = it->GetMother();  
280     Int_t im2 = -1;
281     Int_t id1 = -1;
282     Int_t id2 = -1;
283     
284     Int_t type = it->GetType();  
285       
286     if (im1> -1) {
287       // particle not a primary -> set the daughter indexes for the mother particle"<< endl;
288       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1));
289       mother->SetLastDaughter(nump);
290       if(mother->GetFirstDaughter()==-1)
291         mother->SetFirstDaughter(nump);
292     } else 
293       ++numprim;
294     
295     new (particlesR[nump]) TParticle(it->Encoding(), type,                                              //pdg,stat
296                                      im1, im2, id1, id2,                                                //m1,m2,d1,d2
297                                      mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1), //px,py,pz,e
298                                      pos[0]*1.e-13, pos[1]*1.e-13, pos[2]*1.e-13, 
299                                      it->T()*1.e-13/3e10);                      //x,y,z,t
300     
301     particlesR[nump]->SetUniqueID(nump);
302     nump++;
303     numsec++;
304   }//end for
305   
306   fSecondariesList.clear();
307   printf("Scan and add prim %d sec %d and all %d particles\n",
308          numprim,numsec,nump);
309   return nump;
310 }
311
312 //______________________________________________________________________________
313 void TUHKMgen::Initialize()
314 {
315   // Function overloading the TGenerator::Initialize() member function.
316   // The Monte-Carlo model is initialized (input parameters are transmitted, 
317   // particle list and decay channels are loaded, average multiplicities are calculated, etc.)
318  
319   fInitialState = new InitialStateHydjet();  
320   SetAllParameters();  
321   fInitialState->LoadPDGInfo();  
322   // set the stable flags
323   for(Int_t i=0; i<fStableFlagged; i++) 
324     fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]); 
325
326   if(!fInitialState->MultIni()) 
327     Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n");  //
328 }
329
330 void TUHKMgen::Print(const Option_t*) const
331 {
332   cout << "**********************************************************************************" << endl;
333   cout << "* UHKM Generator interface to ROOT::TGenerator                                   *" << endl;
334   cout << "*  Documentation:                                                                *" << endl;
335   cout << "* I.P.Lokhtin, L.V.Malinina, S.V.Petrushanko, A.M.Snigirev, I.Arsene, K.Tywoniuk *" << endl;
336   cout << "*    Comput.Phys.Commun.180:779-799, 2009                                        *" << endl;
337   cout << "**********************************************************************************" << endl;
338 }
339
340 void TUHKMgen::GenerateEvent()
341 {
342   // Member function overloading the TGenerator::GenerateEvent()
343   // The HYDJET++ model is run and the particle lists (fSourceList and fSecondariesList) are filled
344   
345   fInitialState->Initialize(fSecondariesList, fAllocator);  
346  
347   if(fSecondariesList.empty())Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");  //
348
349   // Run the decays
350
351   if(fInitialState->RunDecays())  
352     fInitialState->Evolve(fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());  
353   
354 }
355
356 void TUHKMgen::SetAllParameters() {
357   // forward all model input parameters to the InitialStateHydjet object
358   // which will handle all the Monte-Carlo simulation using HYDJET++ model
359
360   fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
361   fInitialState->fParams.fAw = fHydjetParams.fAw;
362   fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
363   fInitialState->fParams.fBfix = fHydjetParams.fBfix;
364   fInitialState->fParams.fBmin = fHydjetParams.fBmin;
365   fInitialState->fParams.fBmax = fHydjetParams.fBmax;
366   fInitialState->fParams.fSeed = fHydjetParams.fSeed;
367   fInitialState->fParams.fT = fHydjetParams.fT;
368   fInitialState->fParams.fMuB = fHydjetParams.fMuB;
369   fInitialState->fParams.fMuS = fHydjetParams.fMuS;
370   fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
371   fInitialState->fParams.fThFO = fHydjetParams.fThFO;
372   fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
373
374   fInitialState->fParams.fTau = fHydjetParams.fTau;
375   fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
376   fInitialState->fParams.fR = fHydjetParams.fR;
377   fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
378   fInitialState->fParams.fUmax = fHydjetParams.fUmax;
379   fInitialState->fParams.fDelta = fHydjetParams.fDelta;
380   fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
381
382   fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
383   fInitialState->fParams.fDecay = 1;                   // run decays
384   fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
385   fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
386  
387   fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
388   fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
389   fInitialState->fParams.fIshad = fHydjetParams.fIshad;
390   fInitialState->fParams.fT0 = fHydjetParams.fT0;
391   fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
392   fInitialState->fParams.fNf = fHydjetParams.fNf;      
393   fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;     
394   fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;  
395
396   fInitialState->SetPDGParticleFilename(fParticleFilename);
397   fInitialState->SetPDGDecayFilename(fDecayFilename);
398   //  fInitialState->SetUseCharmParticles(fUseCharmParticles);
399   //  fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
400   //  fInitialState->SetMassRange(fMinMass, fMaxMass);
401   //  for(Int_t i=0; i<fStableFlagged; i++) 
402   //    fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
403   //  cout << "TUHKMgen::SetAllParameters() OUT" << endl;
404 }