Check array ranges
[u/mrichter/AliRoot.git] / TTherminator / AliGenTherminator.cxx
1 // ALICE event generator based on the THERMINATOR model
2 // It reads the test output of the model and puts it onto
3 // the stack
4 // It has an option to use the Lhyquid3D input freeze-out
5 // hypersurface
6 // Author: Adam.Kisiel@cern.ch
7
8 #include <iostream>
9 #include <fstream>
10 #include <sstream>
11 #include <TClonesArray.h>
12 #include <TMCProcess.h>
13 #include <TDatabasePDG.h>
14 #include <TParticle.h>
15
16 #include "AliConst.h"
17 #include "AliDecayer.h"
18 #include "AliGenEventHeader.h"
19 #include "AliGenHijingEventHeader.h"
20 #include "AliGenTherminator.h"
21 #include "AliLog.h"
22 #include "AliRun.h"
23
24 ClassImp(AliGenTherminator)
25
26 using namespace std;
27
28 AliGenTherminator::AliGenTherminator():
29   AliGenMC(),
30   fNt(0),
31   fEventNumber(0),
32   fFileName(""),
33   fFreezeOutModel(""),
34   fTemperature(0.1656),
35   fMiuI(-0.0009),
36   fMiuS(0.0),
37   fMiuB(0.0008),
38   fAlfaRange(6.0),
39   fRapRange(4.0),
40   fRhoMax(8.0),
41   fTau(8.0),
42   fBWA(0.0),
43   fBWVt(1.41),
44   fBWDelay(0.0)
45 {
46   // Default constructor
47 }
48 AliGenTherminator::AliGenTherminator(Int_t npart):
49   AliGenMC(npart),
50   fNt(0),
51   fEventNumber(0),
52   fFileName(""),
53   fFreezeOutModel(""),
54   fTemperature(0.1656),
55   fMiuI(-0.0009),
56   fMiuS(0.0),
57   fMiuB(0.0008),
58   fAlfaRange(6.0),
59   fRapRange(4.0),
60   fRhoMax(8.0),
61   fTau(8.0),
62   fBWA(0.0),
63   fBWVt(1.41),
64   fBWDelay(0.0)
65 {
66   // Constructor specifying the size of the particle table
67   fNprimaries = 0;
68 }
69
70 AliGenTherminator::~AliGenTherminator()
71 {
72   //  AliGenMC::~AliGenMC();
73   //  if (fTherminator) delete fTherminator;
74 }
75
76 void AliGenTherminator::Generate()
77 {
78   // Run single event generation with the Therminator model
79   AliWarning("Generating event from AliGenTherminator");
80
81   Float_t polar[3]    =   {0,0,0};
82   Float_t origin[3]   =   {0,0,0};
83   Float_t origin0[3]  =   {0,0,0};
84   Float_t p[3];
85   Float_t mass, energy;
86
87   Int_t nt  = 0;
88   Int_t j, kf, ks, imo;
89   kf = 0;
90
91   Vertex();
92   for (j=0; j < 3; j++) origin0[j] = fVertex[j];
93
94   // Generate one event
95
96   ((TTherminator *) fMCEvGen)->GenerateEvent();
97   AliWarning("Generated");
98   ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);
99
100   Int_t np = fParticles.GetEntriesFast();
101   AliWarning(Form("Imported %d particles", np));
102
103   Int_t *idsOnStack;
104   idsOnStack = new Int_t[np];
105
106   TParticle *iparticle;
107   Double_t evrot = gRandom->Rndm()*TMath::Pi();
108   
109   for (int i = 0; i < np; i++) {
110     iparticle = (TParticle *) fParticles.At(i);
111     Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
112     Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
113     
114     if (hasDaughter) {
115       // This particle has decayed
116       // It will not be tracked
117       // Add it only once with coorduinates not
118       // smeared with primary vertex position
119
120       kf   = iparticle->GetPdgCode();
121       ks   = iparticle->GetStatusCode();
122       Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
123       Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
124       p[0] = arho*TMath::Cos(aphi + evrot);
125       p[1] = arho*TMath::Sin(aphi + evrot);
126       //     p[0] = iparticle->Px();
127       //     p[1] = iparticle->Py();
128       p[2] = iparticle->Pz();
129       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
130       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
131       
132       Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
133       Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
134       origin[0] = vrho*TMath::Cos(vphi + evrot);
135       origin[1] = vrho*TMath::Sin(vphi + evrot);
136       origin[2] = iparticle->Vz();
137       
138       imo = -1;
139       TParticle* mother = 0;
140       if (hasMother) {
141         imo = iparticle->GetFirstMother();
142         mother = (TParticle *) fParticles.At(imo);
143       } // if has mother   
144       Bool_t tFlag = (!hasDaughter);
145       
146       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
147       PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
148                 p[0],p[1],p[2],energy,
149                 origin[0],origin[1],origin[2],iparticle->T(),
150                 polar[0],polar[1],polar[2],
151                 hasMother ? kPDecay:kPNoProcess,nt);
152       idsOnStack[i] = nt;
153       fNprimaries++;
154       KeepTrack(nt);
155     }
156     else {
157       // This is a final state particle
158       // It will be tracked
159       // Add it TWICE to the stack !!!
160       // First time with event-wide coordicates (for femtoscopy) - 
161       //   this one will not be tracked
162       // Second time with event-wide ccordiantes and vertex smearing
163       //   this one will be tracked
164       
165       kf   = iparticle->GetPdgCode();
166       ks   = iparticle->GetStatusCode();
167       Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
168       Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
169       p[0] = arho*TMath::Cos(aphi + evrot);
170       p[1] = arho*TMath::Sin(aphi + evrot);
171       //     p[0] = iparticle->Px();
172       //     p[1] = iparticle->Py();
173       p[2] = iparticle->Pz();
174       mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
175       energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
176       
177       Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
178       Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
179       origin[0] = vrho*TMath::Cos(vphi + evrot);
180       origin[1] = vrho*TMath::Sin(vphi + evrot);
181       origin[2] = iparticle->Vz();
182       
183       imo = -1;
184       TParticle* mother = 0;
185       if (hasMother) {
186         imo = iparticle->GetFirstMother();
187         mother = (TParticle *) fParticles.At(imo);
188       } // if has mother   
189       Bool_t tFlag = (hasDaughter);
190       
191       printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo);
192       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
193       PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
194                 p[0],p[1],p[2],energy,
195                 origin[0],origin[1],origin[2],iparticle->T(),
196                 polar[0],polar[1],polar[2],
197                 hasMother ? kPDecay:kPNoProcess,nt);
198       idsOnStack[i] = nt;
199       fNprimaries++;
200       KeepTrack(nt);
201
202       origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
203       origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
204       origin[2] = origin0[2]+iparticle->Vz();
205
206       imo = nt;
207       //      mother = (TParticle *) fParticles.At(nt);
208       tFlag = (!hasDaughter);
209       
210       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
211       PushTrack(tFlag,imo,kf,
212                 p[0],p[1],p[2],energy,
213                 origin[0],origin[1],origin[2],iparticle->T(),
214                 polar[0],polar[1],polar[2],
215                 hasMother ? kPDecay:kPNoProcess,nt);
216       fNprimaries++;
217       KeepTrack(nt);
218     }
219   }
220
221
222
223   SetHighWaterMark(fNprimaries);
224
225   TArrayF eventVertex;
226   eventVertex.Set(3);
227   eventVertex[0] = origin0[0];
228   eventVertex[1] = origin0[1];
229   eventVertex[2] = origin0[2];
230
231 // Builds the event header, to be called after each event
232   AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
233
234   // Header
235   //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
236   // Event Vertex
237 //   header->SetPrimaryVertex(eventVertex);
238 //   header->SetNProduced(fNprimaries);
239
240   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
241   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
242   ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
243   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
244   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
245   ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
246   ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
247   ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
248   ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
249
250
251 // 4-momentum vectors of the triggered jets.
252 //
253 // Before final state gluon radiation.
254 //     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
255 //                                            fHijing->GetHINT1(22),
256 //                                            fHijing->GetHINT1(23),
257 //                                            fHijing->GetHINT1(24));
258
259 //     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
260 //                                            fHijing->GetHINT1(32),
261 //                                            fHijing->GetHINT1(33),
262 //                                            fHijing->GetHINT1(34));
263 // // After final state gluon radiation.
264 //     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
265 //                                            fHijing->GetHINT1(27),
266 //                                            fHijing->GetHINT1(28),
267 //                                            fHijing->GetHINT1(29));
268
269 //     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
270 //                                            fHijing->GetHINT1(37),
271 //                                            fHijing->GetHINT1(38),
272 //                                            fHijing->GetHINT1(39));
273 //     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
274 // Bookkeeping for kinematic bias
275 //     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
276 // Event Vertex
277   header->SetPrimaryVertex(fVertex);
278   AddHeader(header);
279   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
280
281   delete idsOnStack;
282
283   //  gAlice->SetGenEventHeader(header); 
284 }
285
286 void AliGenTherminator::Init()
287 {
288   // Initialize global variables and
289   // particle and decay tables
290   if (fFileName.Length() == 0)
291     fFileName = "event.out";
292   ReadShareParticleTable();
293
294   SetMC(new TTherminator());
295   
296   AliGenMC::Init();
297   ((TTherminator *) fMCEvGen)->Initialize();
298 }
299 void AliGenTherminator::SetFileName(const char *infilename)
300 {
301   // Set parameter filename
302   fFileName = infilename;
303 }
304 void AliGenTherminator::SetEventNumberInFile(int evnum)
305 {
306   // Set number of events to generate - default: 1
307   fEventNumber = evnum;
308 }
309
310 void AliGenTherminator::ReadShareParticleTable()
311 {
312   // Read in particle table from share
313   // and add missing particle type to TDatabasePDG
314
315   char str[50];
316   char str1[200];
317     
318   TDatabasePDG *tInstance = TDatabasePDG::Instance();
319   TParticlePDG *tParticleType;
320
321   AliWarning(Form("Reading particle types from particles.data"));
322
323   TString aroot = gSystem->Getenv("ALICE_ROOT");
324   ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data());
325   //  ifstream in("particles.data");
326   
327   int charge;
328     
329   int number=0;
330   if ((in) && (in.is_open()))
331     {
332       //START OF HEAD-LINE
333       in.ignore(200,'\n');
334       in.ignore(200,'\n');
335       in.ignore(200,'\n');
336       //END OF HEAD-LINE
337       
338       while (in>>str)
339         {
340           if (/*(*str == '#')||*/(*str<65)||(*str>122))
341             {
342               in.getline(str1,200);
343               continue;
344             }
345           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
346           
347           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
348           number++;
349           tParticleType = tInstance->GetParticle((int) mc);
350           if (!tParticleType) {
351             charge = 0;
352             if (strstr(str, "plu")) charge = 1;
353             if (strstr(str, "min")) charge = -1;
354             if (strstr(str, "plb")) charge = -1;
355             if (strstr(str, "mnb")) charge = 1;
356             if (strstr(str, "plp")) charge = 2;
357             if (strstr(str, "ppb")) charge = -2;
358             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
359             AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
360             //      AliWarning(Form("Quantum numbers q s c aq as ac tI3 %lf %lf %lf %lf %lf %lf %lf", q, s, c, aq, as, ac, tI3));
361
362           }
363         }
364       in.close();
365     }
366   CreateTherminatorInputFile();
367 }
368
369 void AliGenTherminator::CreateTherminatorInputFile()
370 {
371   // Create Therminator input file
372   const char *aroot = gSystem->Getenv("ALICE_ROOT");
373   ofstream *ostr = new ofstream("therminator.in");
374   (*ostr) << "NumberOfEvents = 1" << endl;
375   (*ostr) << "Randomize = 1" << endl;
376   (*ostr) << "TableType = SHARE" << endl;
377   (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl;
378   (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
379   (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl;
380   (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
381   (*ostr) << "BWVt = " << fBWVt << endl;
382   (*ostr) << "Tau = " << fTau << endl;
383   (*ostr) << "RhoMax = " << fRhoMax << endl;
384   (*ostr) << "Temperature = " << fTemperature << endl;
385   (*ostr) << "MiuI = " << fMiuI << endl;
386   (*ostr) << "MiuS = " << fMiuS << endl;
387   (*ostr) << "MiuB = " << fMiuB << endl;
388   (*ostr) << "AlphaRange = " << fAlfaRange << endl;
389   (*ostr) << "RapidityRange = " << fRapRange << endl;
390   (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
391 }
392
393 void AliGenTherminator::SetModel(const char *model)
394 {
395   // Set the freeze-out model to use
396   fFreezeOutModel = model;
397 }
398
399 void AliGenTherminator::SetLhyquidSet(const char *set)
400 {
401   // Select one of pregenerated Lhyquid hypersurfaces
402   const char *aroot = gSystem->Getenv("ALICE_ROOT");
403   if (strstr(set, "LHC500C0005")) {
404     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
405     AliWarning(Form("  Pb-Pb collisions, centrality 0-5%"));
406     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
407     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
408     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
409     fFOHSlocation = aroot;
410     fFOHSlocation += "/TTherminator/data/LHC500C0005";
411   }
412   else if (strstr(set, "LHC500C2030")) {
413     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
414     AliWarning(Form("  Pb-Pb collisions, centrality 20-30%"));
415     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
416     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
417     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
418     fFOHSlocation = aroot;
419     fFOHSlocation += "/TTherminator/data/LHC500C2030";
420   }
421   else {
422     AliWarning(Form("Did not find Lhyquid set %s", set));
423     AliWarning(Form("Reverting to default: current directory"));
424     fFOHSlocation = "";
425   }
426 }
427
428 void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)
429 {
430   // Select Your own Lhyquid hypersurface
431   fFOHSlocation = inputdir;
432 }