Minor fix for compatibility with MacOS
[u/mrichter/AliRoot.git] / TTherminator / AliGenTherminator.cxx
CommitLineData
2e967919 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"
9abb118d 17#include "AliGenHijingEventHeader.h"
2e967919 18#include "AliGenTherminator.h"
19#include "AliLog.h"
20#include "AliRun.h"
21
22ClassImp(AliGenTherminator)
23
24using namespace std;
25
26AliGenTherminator::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}
47AliGenTherminator::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
70AliGenTherminator::~AliGenTherminator()
71{
72 // AliGenMC::~AliGenMC();
73 // if (fTherminator) delete fTherminator;
74}
75
76void 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;
9abb118d 104 Double_t evrot = gRandom->Rndm()*TMath::Pi();
2e967919 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();
9abb118d 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();
2e967919 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
9abb118d 153// Builds the event header, to be called after each event
154 AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
155
2e967919 156 // Header
9abb118d 157 // AliGenEventHeader* header = new AliGenEventHeader("Therminator");
2e967919 158 // Event Vertex
9abb118d 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
2e967919 203 gAlice->SetGenEventHeader(header);
204}
205
206void 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}
219void AliGenTherminator::SetFileName(const char *infilename)
220{
221 // Set parameter filename
222 fFileName = infilename;
223}
224void AliGenTherminator::SetEventNumberInFile(int evnum)
225{
226 // Set number of events to generate - default: 1
227 fEventNumber = evnum;
228}
229
230void 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
286void 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
308void AliGenTherminator::SetModel(const char *model)
309{
310 // Set the freeze-out model to use
311 fFreezeOutModel = model;
312}