Fix in the initialization (Solaris x86)
[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 // Author: Adam.Kisiel@cern.ch
5
6 #include <iostream>
7 #include <fstream>
8 #include <sstream>
9 #include <TClonesArray.h>
10 #include <TMCProcess.h>
11 #include <TDatabasePDG.h>
12 #include <TParticle.h>
13
14 #include "AliConst.h"
15 #include "AliDecayer.h"
16 #include "AliGenEventHeader.h"
17 #include "AliGenHijingEventHeader.h"
18 #include "AliGenTherminator.h"
19 #include "AliLog.h"
20 #include "AliRun.h"
21
22 ClassImp(AliGenTherminator)
23
24 using namespace std;
25
26 AliGenTherminator::AliGenTherminator():
27   AliGenMC(),
28   fNt(0),
29   fEventNumber(0),
30   fFileName(""),
31   fFreezeOutModel(""),
32   fTemperature(0.1656),
33   fMiuI(-0.0009),
34   fMiuS(0.0),
35   fMiuB(0.0008),
36   fAlfaRange(6.0),
37   fRapRange(4.0),
38   fRhoMax(8.0),
39   fTau(8.0),
40   fBWA(0.0),
41   fBWVt(1.41),
42   fBWDelay(0.0)
43 {
44   // Default constructor
45   fParticles = new TClonesArray("TParticle",20000);    
46 }
47 AliGenTherminator::AliGenTherminator(Int_t npart):
48   AliGenMC(npart),
49   fNt(0),
50   fEventNumber(0),
51   fFileName(""),
52   fFreezeOutModel(""),
53   fTemperature(0.1656),
54   fMiuI(-0.0009),
55   fMiuS(0.0),
56   fMiuB(0.0008),
57   fAlfaRange(6.0),
58   fRapRange(4.0),
59   fRhoMax(8.0),
60   fTau(8.0),
61   fBWA(0.0),
62   fBWVt(1.41),
63   fBWDelay(0.0)
64 {
65   // Constructor specifying the size of the particle table
66   fParticles = new TClonesArray("TParticle",20000);    
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   TParticle *iparticle;
104   Double_t evrot = gRandom->Rndm()*TMath::Pi();
105   
106   for (int i = 0; i < np; i++) {
107     iparticle = (TParticle *) fParticles->At(i);
108     Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
109     Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
110     
111     kf   = iparticle->GetPdgCode();
112     ks   = iparticle->GetStatusCode();
113     Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
114     Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
115     p[0] = arho*TMath::Cos(aphi + evrot);
116     p[1] = arho*TMath::Sin(aphi + evrot);
117 //     p[0] = iparticle->Px();
118 //     p[1] = iparticle->Py();
119     p[2] = iparticle->Pz();
120     mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
121     energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
122     origin[0] = origin0[0]+iparticle->Vx();
123     origin[1] = origin0[1]+iparticle->Vy();
124     origin[2] = origin0[2]+iparticle->Vz();
125               
126     imo = -1;
127     TParticle* mother = 0;
128     if (hasMother) {
129       imo = iparticle->GetFirstMother();
130       mother = (TParticle *) fParticles->At(imo);
131     } // if has mother   
132     Bool_t tFlag = (!hasDaughter);
133
134     printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
135     PushTrack(tFlag,imo,kf,
136               p[0],p[1],p[2],energy,
137               origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
138               polar[0],polar[1],polar[2],
139               hasMother ? kPDecay:kPNoProcess,nt);
140     
141     fNprimaries++;
142     KeepTrack(nt);
143   }
144
145   SetHighWaterMark(fNprimaries);
146
147   TArrayF eventVertex;
148   eventVertex.Set(3);
149   eventVertex[0] = origin0[0];
150   eventVertex[1] = origin0[1];
151   eventVertex[2] = origin0[2];
152
153 // Builds the event header, to be called after each event
154   AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
155
156   // Header
157   //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
158   // Event Vertex
159 //   header->SetPrimaryVertex(eventVertex);
160 //   header->SetNProduced(fNprimaries);
161
162   ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
163   ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
164   ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
165   ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
166   ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
167   ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
168   ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
169   ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
170   ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
171
172
173 // 4-momentum vectors of the triggered jets.
174 //
175 // Before final state gluon radiation.
176 //     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
177 //                                            fHijing->GetHINT1(22),
178 //                                            fHijing->GetHINT1(23),
179 //                                            fHijing->GetHINT1(24));
180
181 //     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
182 //                                            fHijing->GetHINT1(32),
183 //                                            fHijing->GetHINT1(33),
184 //                                            fHijing->GetHINT1(34));
185 // // After final state gluon radiation.
186 //     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
187 //                                            fHijing->GetHINT1(27),
188 //                                            fHijing->GetHINT1(28),
189 //                                            fHijing->GetHINT1(29));
190
191 //     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
192 //                                            fHijing->GetHINT1(37),
193 //                                            fHijing->GetHINT1(38),
194 //                                            fHijing->GetHINT1(39));
195 //     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
196 // Bookkeeping for kinematic bias
197 //     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
198 // Event Vertex
199   header->SetPrimaryVertex(fVertex);
200   AddHeader(header);
201   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
202
203   gAlice->SetGenEventHeader(header); 
204 }
205
206 void AliGenTherminator::Init()
207 {
208   // Initialize global variables and
209   // particle and decay tables
210   if (fFileName.Length() == 0)
211     fFileName = "event.out";
212   ReadShareParticleTable();
213
214   SetMC(new TTherminator());
215   
216   AliGenMC::Init();
217   ((TTherminator *) fMCEvGen)->Initialize();
218 }
219 void AliGenTherminator::SetFileName(const char *infilename)
220 {
221   // Set parameter filename
222   fFileName = infilename;
223 }
224 void AliGenTherminator::SetEventNumberInFile(int evnum)
225 {
226   // Set number of events to generate - default: 1
227   fEventNumber = evnum;
228 }
229
230 void AliGenTherminator::ReadShareParticleTable()
231 {
232   // Read in particle table from share
233   // and add missing particle type to TDatabasePDG
234
235   char str[50];
236   char str1[200];
237     
238   TDatabasePDG *tInstance = TDatabasePDG::Instance();
239   TParticlePDG *tParticleType;
240
241   AliWarning(Form("Reading particle types from particles.data"));
242   ifstream in("particles.data");
243   
244   int charge;
245     
246   int number=0;
247   if ((in) && (in.is_open()))
248     {
249       //START OF HEAD-LINE
250       in.ignore(200,'\n');
251       in.ignore(200,'\n');
252       in.ignore(200,'\n');
253       //END OF HEAD-LINE
254       
255       while (in>>str)
256         {
257           if (/*(*str == '#')||*/(*str<65)||(*str>122))
258             {
259               in.getline(str1,200);
260               continue;
261             }
262           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
263           
264           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
265           number++;
266           tParticleType = tInstance->GetParticle((int) mc);
267           if (!tParticleType) {
268             charge = 0;
269             if (strstr(str, "plu")) charge = 1;
270             if (strstr(str, "min")) charge = -1;
271             if (strstr(str, "plb")) charge = -1;
272             if (strstr(str, "mnb")) charge = 1;
273             if (strstr(str, "plp")) charge = 2;
274             if (strstr(str, "ppb")) charge = -2;
275             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
276             AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
277             //      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));
278
279           }
280         }
281       in.close();
282     }
283   CreateTherminatorInputFile();
284 }
285
286 void AliGenTherminator::CreateTherminatorInputFile()
287 {
288   // Create Therminator input file
289   ofstream *ostr = new ofstream("therminator.in");
290   (*ostr) << "NumberOfEvents = 1" << endl;
291   (*ostr) << "Randomize = 1" << endl;
292   (*ostr) << "TableType = SHARE" << endl;
293   (*ostr) << "InputDirSHARE = ." << endl;
294   (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
295   (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
296   (*ostr) << "BWVt = " << fBWVt << endl;
297   (*ostr) << "Tau = " << fTau << endl;
298   (*ostr) << "RhoMax = " << fRhoMax << endl;
299   (*ostr) << "Temperature = " << fTemperature << endl;
300   (*ostr) << "MiuI = " << fMiuI << endl;
301   (*ostr) << "MiuS = " << fMiuS << endl;
302   (*ostr) << "MiuB = " << fMiuB << endl;
303   (*ostr) << "AlphaRange = " << fAlfaRange << endl;
304   (*ostr) << "RapidityRange = " << fRapRange << endl;
305   (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
306 }
307
308 void AliGenTherminator::SetModel(const char *model)
309 {
310   // Set the freeze-out model to use
311   fFreezeOutModel = model;
312 }