Removing extra symbols
[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 ClassImp(TUHKMgen)
34
35 TUHKMgen::TUHKMgen() : 
36   TGenerator("UHKM","UHKM"),
37   fInitialState(0x0),
38   fUseCharmParticles(kFALSE),
39   fMinWidth(0.0),
40   fMaxWidth(10.0),
41   fMinMass(0.0001),
42   fMaxMass(10.0)
43 {
44   cout << "TUHKMgen::TUHKMgen() IN" << endl;
45
46   // Set reasonable default values for LHC
47   
48   fHydjetParams.fSqrtS=5500; //LHC
49   fHydjetParams.fAw=207;//Au-Au
50   fHydjetParams.fIfb    = 1;      // flag of type of centrality generation (=0 is fixed by fBfix, not 0 
51                                    //impact parameter is generated in each event between fBfmin 
52                                    //and fBmax according with Glauber model (f-la 30)
53   fHydjetParams.fBfix=0.;
54   fHydjetParams.fBmin=0.;
55   fHydjetParams.fBmax=0.5; //0-5% centrality
56   fHydjetParams.fT = 0.170;
57   fHydjetParams.fMuB = 0.0;
58   fHydjetParams.fMuS = 0.0;
59   fHydjetParams.fMuI3 = 0.0;
60   fHydjetParams.fThFO = 0.130;
61   fHydjetParams.fMu_th_pip = 0.0;
62   fHydjetParams.fSeed=0;
63   fHydjetParams.fTau=10.;
64   fHydjetParams.fSigmaTau=3.;
65   fHydjetParams.fR=11.;
66   fHydjetParams.fYlmax=4.0;
67   fHydjetParams.fUmax=1.1;
68   fHydjetParams.fDelta=0.;
69   fHydjetParams.fEpsilon=0.;
70   fHydjetParams.fWeakDecay=0; //>=0 on ,-1 off
71   fHydjetParams.fEtaType=1;//gaus
72   fHydjetParams.fCorrS=1.;
73   fHydjetParams.fNhsel=2;
74   fHydjetParams.fIshad=1;
75   fHydjetParams.fPtmin=7.0;
76   fHydjetParams.fT0=0.8;
77   fHydjetParams.fTau0=0.1;
78   fHydjetParams.fNf=0;
79   fHydjetParams.fIenglu=0;
80   fHydjetParams.fIanglu=0;
81
82
83   // Set reasonable default values for RHIC
84 /*
85   fHydjetParams.fSqrtS=200; //RHIC          
86   fHydjetParams.fAw=197;//Au-Au        
87   fHydjetParams.fIfb    = 1;      
88   fHydjetParams.fBfix=0.;       
89   fHydjetParams.fBmin=0.;       
90   fHydjetParams.fBmax=0.5; //0-5% centrality
91   fHydjetParams.fT = 0.165;
92   fHydjetParams.fMuB = 0.0285;
93   fHydjetParams.fMuS = 0.007;
94   fHydjetParams.fMuI3 = -0.001;
95   fHydjetParams.fThFO = 0.100;
96   fHydjetParams.fMu_th_pip = 0.053;
97   fHydjetParams.fSeed=0;  
98   fHydjetParams.fTau=8.;
99   fHydjetParams.fSigmaTau=2.;
100   fHydjetParams.fR=10.;
101   fHydjetParams.fYlmax=3.3;
102   fHydjetParams.fUmax=1.1;
103   fHydjetParams.fDelta=0.;
104   fHydjetParams.fEpsilon=0.;
105   fHydjetParams.fWeakDecay=0;//>=0 on ,-1 off 
106   fHydjetParams.fEtaType=1;//gaus
107   fHydjetParams.fCorrS=1.;  
108   fHydjetParams.fNhsel=2;
109   fHydjetParams.fIshad=1;
110   fHydjetParams.fPtmin=3.4;
111   fHydjetParams.fT0=0.3;
112   fHydjetParams.fTau0=0.4;
113   fHydjetParams.fNf=2;      
114   fHydjetParams.fIenglu=0;     
115   fHydjetParams.fIanglu=0;  
116 */
117    
118   strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")));
119   strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")));
120   for(Int_t i=0; i<500; i++) {
121     fStableFlagPDG[i] = 0;
122     fStableFlagStatus[i] = kFALSE;
123   }
124   fStableFlagged = 0;
125   cout << "TUHKMgen::TUHKMgen() OUT" << endl;
126 }
127
128 //______________________________________________________________________________
129 TUHKMgen::~TUHKMgen()
130 {
131   cout << "TUHKMgen::~TUHKMgen() IN" << endl;
132   if(fInitialState)
133     delete fInitialState;
134   cout << "TUHKMgen::~TUHKMgen() OUT" << endl;
135 }
136
137 void TUHKMgen::SetAllParametersRHIC()
138 {
139   SetEcms(200.0);                // RHIC top energy
140   SetAw(197);                    // Au+Au
141   SetIfb(1);
142   SetBfix(0.);
143   SetBmin(0.0);                  // 0%
144   SetBmax(0.5);                  // 5%
145   SetChFrzTemperature(0.165);    // T_ch = 165 MeV
146   SetMuB(0.0285);                // mu_B = 28.5 MeV
147   SetMuS(0.007);                 // mu_S = 7 MeV
148   SetMuQ(-0.001);                // mu_Q = -1 MeV
149   SetThFrzTemperature(0.100);    // T_th = 100 MeV
150   SetMuPionThermal(0.053);       // mu_th_pion = 53 MeV
151   SetSeed(0);                    // use UNIX time
152   SetTauB(8.0);                  // tau = 8 fm/c
153   SetSigmaTau(2.0);              // sigma_tau = 2 fm/c
154   SetRmaxB(10.0);                // fR = 10 fm
155   SetYlMax(3.3);                 // fYmax = 3.3
156   SetEtaRMax(1.1);               // Umax = 1.1
157   SetMomAsymmPar(0.0);           // delta = 0.0
158   SetCoordAsymmPar(0.0);         // epsilon = 0.0
159   SetFlagWeakDecay(0.0);           // weak decay on (<0 off !!!)
160   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
161   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
162   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
163   SetPyquenShad(1);              // shadowing on (0 off)
164   SetPyquenPtmin(3.4);           // ptmin = 3.4 GeV/c
165   SetPyquenT0(0.3);              // T0 = 300 MeV
166   SetPyquenTau0(0.4);            // tau0 = 0.4 fm/c
167   SetPyquenNf(2);                // 2 flavours
168   SetPyquenIenglu(0);            // radiative and collisional energy loss
169   SetPyquenIanglu(0);            // small gluon angular distribution
170 }
171
172 void TUHKMgen::SetAllParametersLHC()
173 {
174   SetEcms(5500.0);               // LHC
175   SetAw(207);                    // Pb+Pb
176   SetIfb(1);
177   SetBfix(0.);                  // 0
178   SetBmin(0.0);                  // 0%
179   SetBmax(0.5);                  // 5%
180   SetChFrzTemperature(0.170);    // T_ch = 170 MeV
181   SetMuB(0.0);                   // mu_B = 0 MeV
182   SetMuS(0.0);                   // mu_S = 0 MeV
183   SetMuQ(0.0);                   // mu_Q = 0 MeV
184   SetThFrzTemperature(0.130);    // T_th = 130 MeV
185   SetMuPionThermal(0.0);         // mu_th_pion = 0 MeV
186   SetSeed(0);                    // use UNIX time
187   SetTauB(10.0);                 // tau = 10 fm/c
188   SetSigmaTau(3.0);              // sigma_tau = 3 fm/c
189   SetRmaxB(11.0);                // fR = 11 fm
190   SetYlMax(4.0);                 // fYmax = 4.0
191   SetEtaRMax(1.1);               // Umax = 1.1
192   SetMomAsymmPar(0.0);           // delta = 0.0
193   SetCoordAsymmPar(0.0);         // epsilon = 0.0
194   SetFlagWeakDecay(0);           // weak decay on (<0 off !!!)
195   SetEtaType(1);                 // gaus distributed with fYmax dispersion (0 means boost invariant)
196   SetGammaS(1.0);                // gammaS = 1.0 (no strangeness canonical suppresion)
197   SetPyquenNhsel(2);             // hydro on, jets on, jet quenching on
198   SetPyquenShad(1);              // shadowing on (0 off)
199   SetPyquenPtmin(7.0);           // ptmin = 7.0 GeV/c
200   SetPyquenT0(0.8);              // T0 = 800 MeV
201   SetPyquenTau0(0.1);            // tau0 = 0.4 fm/c
202   SetPyquenNf(0);                // 0 flavours
203   SetPyquenIenglu(0);            // radiative and collisional energy loss
204   SetPyquenIanglu(0);            // small gluon angular distribution
205 }
206
207 TObjArray* TUHKMgen::ImportParticles(Option_t *option)
208 {
209   cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) IN" << endl;
210   fParticles->Clear();
211   cout << "TObjArray* TUHKMgen::ImportParticles(Option_t): UHKM stack contains " << fNPsec << " particles." << endl;
212   Int_t nump = 0;
213   LPIT_t it,e;
214   Double_t Fm2Gev = 1.;//0.197327; //UHKM has coordinates in fm
215    
216   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
217     TVector3 pos(it->Pos().Vect());
218     TVector3 mom(it->Mom().Vect());
219     Float_t m1 = it->TableMass();
220     Int_t im1 = it->GetMother();
221     Int_t im2 = -1;
222     Int_t id1 = -1;
223     Int_t id2 = -1;
224     Int_t id3 = -1;
225
226     Int_t nd = it->GetNDaughters();
227     
228     Int_t type = it->GetType();  // 0-hydro, 1-jets
229
230     if (im1> -1) {
231       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(im1));     
232       mother->SetLastDaughter(nump);
233       if (mother->GetFirstDaughter()==-1)
234         mother->SetFirstDaughter(nump);
235     }
236
237     nump++;
238
239     TParticle* p = new TParticle(it->Encoding(), type,                          //pdg,stat
240                                  im1, im2, id1, id2,                                            //m1,m2,d1,d2
241                                  mom[0], mom[1], mom[2], TMath::Sqrt(mom.Mag2() + m1 * m1),     //px,py,pz,e
242                                  pos[0]*1.e-13*Fm2Gev, pos[1]*1.e-13*Fm2Gev, pos[2]*1.e-13*Fm2Gev, 
243                                  it->T()*1.e-13*Fm2Gev/3e10);                           //x,y,z,t
244   
245      p->SetUniqueID(nump);
246      fParticles->Add(p);
247   } //end for
248  
249   fAllocator.FreeList(fSourceList);
250   fAllocator.FreeList(fSecondariesList);
251   cout << "TObjArray* TUHKMgen::ImportParticles(Option_t) OUT" << endl;
252   return fParticles;
253 }
254
255 Int_t TUHKMgen::ImportParticles(TClonesArray *particles,Option_t* option)
256 {
257   //z-start import
258   cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) IN" << endl;
259   if(particles==0) return 0;
260   TClonesArray &particlesR=*particles;
261   particlesR.Clear();
262
263  
264   Int_t numprim,numsec;  numprim=numsec=0;
265   Int_t nump = 0;
266   LPIT_t it,e;
267   Double_t Fm2Gev = 1.; //UHKM has coordinates in fm //0.197327;
268
269   cout << "TUHKMgen::ImportParticles() option(All or Sec) = " << option << endl;
270   for(it = fSecondariesList.begin(), e = fSecondariesList.end(); it != e; ++it) {
271     //      cout << "TUHKMgen::ImportParticles() import particle pdg(" << it->Encoding() << ")" << endl;
272     TVector3 pos(it->Pos().Vect());
273     TVector3 mom(it->Mom().Vect());
274     Float_t m1 = it->TableMass();
275     Int_t im1 = it->GetMother();
276     Int_t nd = 0;
277     Int_t im2 = -1;
278     Int_t id1 = -1;
279     Int_t id2 = -1;
280     Int_t id3 = -1;
281     
282     nd = it->GetNDaughters();
283
284     Int_t type = it->GetType();  // 0-hydro, 1-jets
285     
286     if (im1> -1) {
287       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(im1+1));    
288       mother->SetLastDaughter(nump);
289       if (mother->GetFirstDaughter()==-1)
290         mother->SetFirstDaughter(nump);
291     }
292       
293     nump++;
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*Fm2Gev, pos[1]*1.e-13*Fm2Gev, pos[2]*1.e-13*Fm2Gev, 
299                                      it->T()*1.e-13*Fm2Gev/3e10);                       //x,y,z,t
300
301
302     particlesR[nump]->SetUniqueID(nump);
303     numsec++;
304   }//end for
305   
306   fAllocator.FreeList(fSourceList);
307   fAllocator.FreeList(fSecondariesList);
308   printf("Scan and add prim %d sec %d and all %d particles\n",
309          numprim,numsec,nump);
310   cout << "Int_t TUHKMgen::ImportParticles(TClonesArray*, Option_t) OUT" << endl;
311   return nump;
312 }
313
314 //______________________________________________________________________________
315 void TUHKMgen::Initialize()
316 {
317   cout << "TUHKMgen::Initialize() IN" << endl;
318   fInitialState = new InitialStateHydjet();
319   SetAllParameters();
320   fInitialState->LoadPDGInfo();
321   // set the stable flags
322   for(Int_t i=0; i<fStableFlagged; i++) 
323     fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
324
325   if(!fInitialState->MultIni())
326     Error("TUHKMgen::Initialize()", "Bad status return from MultIni(). Check it out!! \n");
327   cout << "TUHKMgen::Initialize() OUT" << endl;
328 }
329
330 void TUHKMgen::Print() 
331 {
332   cout << "TUHKMgen::Print() method not implemented yet!!" << endl;
333 }
334
335 void TUHKMgen::GenerateEvent()
336 {
337   cout << "TUHKMgen::GenerateEvent() IN" << endl;
338   // Generate the initial particles
339   fInitialState->Initialize(fSourceList, fAllocator);
340   cout << "TUHKMgen::fInitialState->Initialize" << endl;
341  
342   if(fSourceList.empty()) 
343     Error("TUHKMgen::GenerateEvent()", "Source particle list empty after fireball initialization!! \n");
344
345   // Run the decays
346  
347   cout << "after Initilize: get_time, weak_decay_limit" <<fInitialState->GetTime()<<fInitialState->GetWeakDecayLimit()<< endl; 
348   if(fInitialState->GetTime() >= 0.)
349     fInitialState->Evolve(fSourceList, fSecondariesList, fAllocator, fInitialState->GetWeakDecayLimit());
350   cout << "TUHKMgen::GenerateEvent() OUT" << endl;
351 }
352
353 void TUHKMgen::SetAllParameters() {
354   cout << "TUHKMgen::SetAllParameters() IN" << endl;
355
356   fInitialState->fParams.fSqrtS = fHydjetParams.fSqrtS;
357   fInitialState->fParams.fAw = fHydjetParams.fAw;
358   fInitialState->fParams.fIfb = fHydjetParams.fIfb ;
359   fInitialState->fParams.fBfix = fHydjetParams.fBfix;
360   fInitialState->fParams.fBmin = fHydjetParams.fBmin;
361   fInitialState->fParams.fBmax = fHydjetParams.fBmax;
362   fInitialState->fParams.fSeed = fHydjetParams.fSeed;
363   fInitialState->fParams.fT = fHydjetParams.fT;
364   fInitialState->fParams.fMuB = fHydjetParams.fMuB;
365   fInitialState->fParams.fMuS = fHydjetParams.fMuS;
366   fInitialState->fParams.fMuI3 = fHydjetParams.fMuI3;
367   fInitialState->fParams.fThFO = fHydjetParams.fThFO;
368   fInitialState->fParams.fMu_th_pip = fHydjetParams.fMu_th_pip;
369
370   fInitialState->fParams.fTau = fHydjetParams.fTau;
371   fInitialState->fParams.fSigmaTau = fHydjetParams.fSigmaTau;
372   fInitialState->fParams.fR = fHydjetParams.fR;
373   fInitialState->fParams.fYlmax = fHydjetParams.fYlmax;
374   fInitialState->fParams.fUmax = fHydjetParams.fUmax;
375   fInitialState->fParams.fDelta = fHydjetParams.fDelta;
376   fInitialState->fParams.fEpsilon = fHydjetParams.fEpsilon;
377
378  // fInitialState->fParams.fTime = fHydjetParams.fTime;
379   fInitialState->fParams.fWeakDecay = fHydjetParams.fWeakDecay;
380   fInitialState->fParams.fCorrS = fHydjetParams.fCorrS;
381   fInitialState->fParams.fEtaType = fHydjetParams.fEtaType;
382  
383   fInitialState->fParams.fPtmin = fHydjetParams.fPtmin;
384   fInitialState->fParams.fNhsel = fHydjetParams.fNhsel;
385   fInitialState->fParams.fIshad = fHydjetParams.fIshad;
386   fInitialState->fParams.fT0 = fHydjetParams.fT0;
387   fInitialState->fParams.fTau0 = fHydjetParams.fTau0;
388   fInitialState->fParams.fNf = fHydjetParams.fNf;      
389   fInitialState->fParams.fIenglu = fHydjetParams.fIenglu;     
390   fInitialState->fParams.fIanglu = fHydjetParams.fIanglu;  
391
392   fInitialState->SetPDGParticleFilename(fParticleFilename);
393   fInitialState->SetPDGDecayFilename(fDecayFilename);
394   fInitialState->SetUseCharmParticles(fUseCharmParticles);
395   fInitialState->SetWidthRange(fMinWidth, fMaxWidth);
396   fInitialState->SetMassRange(fMinMass, fMaxMass);
397   //for(Int_t i=0; i<fStableFlagged; i++) fInitialState->SetPDGParticleStable(fStableFlagPDG[i], fStableFlagStatus[i]);
398   cout << "TUHKMgen::SetAllParameters() OUT" << endl;
399 }
400