Removal of the dependencies on the DATE libs, which in fact was needed at all
[u/mrichter/AliRoot.git] / TTherminator / AliGenTherminator.cxx
... / ...
CommitLineData
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 "AliGenTherminator.h"
18#include "AliLog.h"
19#include "AliRun.h"
20
21ClassImp(AliGenTherminator)
22
23using namespace std;
24
25AliGenTherminator::AliGenTherminator():
26 AliGenMC(),
27 fNt(0),
28 fEventNumber(0),
29 fFileName(""),
30 fFreezeOutModel(""),
31 fTemperature(0.1656),
32 fMiuI(-0.0009),
33 fMiuS(0.0),
34 fMiuB(0.0008),
35 fAlfaRange(6.0),
36 fRapRange(4.0),
37 fRhoMax(8.0),
38 fTau(8.0),
39 fBWA(0.0),
40 fBWVt(1.41),
41 fBWDelay(0.0)
42{
43 // Default constructor
44 fParticles = new TClonesArray("TParticle",20000);
45}
46AliGenTherminator::AliGenTherminator(Int_t npart):
47 AliGenMC(npart),
48 fNt(0),
49 fEventNumber(0),
50 fFileName(""),
51 fFreezeOutModel(""),
52 fTemperature(0.1656),
53 fMiuI(-0.0009),
54 fMiuS(0.0),
55 fMiuB(0.0008),
56 fAlfaRange(6.0),
57 fRapRange(4.0),
58 fRhoMax(8.0),
59 fTau(8.0),
60 fBWA(0.0),
61 fBWVt(1.41),
62 fBWDelay(0.0)
63{
64 // Constructor specifying the size of the particle table
65 fParticles = new TClonesArray("TParticle",20000);
66 fNprimaries = 0;
67}
68
69AliGenTherminator::~AliGenTherminator()
70{
71 // AliGenMC::~AliGenMC();
72 // if (fTherminator) delete fTherminator;
73}
74
75void AliGenTherminator::Generate()
76{
77 // Run single event generation with the Therminator model
78 AliWarning("Generating event from AliGenTherminator");
79
80 Float_t polar[3] = {0,0,0};
81 Float_t origin[3] = {0,0,0};
82 Float_t origin0[3] = {0,0,0};
83 Float_t p[3];
84 Float_t mass, energy;
85
86 Int_t nt = 0;
87 Int_t j, kf, ks, imo;
88 kf = 0;
89
90 Vertex();
91 for (j=0; j < 3; j++) origin0[j] = fVertex[j];
92
93 // Generate one event
94
95 ((TTherminator *) fMCEvGen)->GenerateEvent();
96 AliWarning("Generated");
97 ((TTherminator *) fMCEvGen)->ImportParticles(fParticles);
98
99 Int_t np = fParticles->GetEntriesFast();
100 AliWarning(Form("Imported %d particles", np));
101
102 TParticle *iparticle;
103
104 for (int i = 0; i < np; i++) {
105 iparticle = (TParticle *) fParticles->At(i);
106 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
107 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
108
109 kf = iparticle->GetPdgCode();
110 ks = iparticle->GetStatusCode();
111 p[0] = iparticle->Px();
112 p[1] = iparticle->Py();
113 p[2] = iparticle->Pz();
114 mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
115 energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
116 origin[0] = origin0[0]+iparticle->Vx();
117 origin[1] = origin0[1]+iparticle->Vy();
118 origin[2] = origin0[2]+iparticle->Vz();
119
120 imo = -1;
121 TParticle* mother = 0;
122 if (hasMother) {
123 imo = iparticle->GetFirstMother();
124 mother = (TParticle *) fParticles->At(imo);
125 } // if has mother
126 Bool_t tFlag = (!hasDaughter);
127
128 printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
129 PushTrack(tFlag,imo,kf,
130 p[0],p[1],p[2],energy,
131 origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
132 polar[0],polar[1],polar[2],
133 hasMother ? kPDecay:kPNoProcess,nt);
134
135 fNprimaries++;
136 KeepTrack(nt);
137 }
138
139 SetHighWaterMark(fNprimaries);
140
141 TArrayF eventVertex;
142 eventVertex.Set(3);
143 eventVertex[0] = origin0[0];
144 eventVertex[1] = origin0[1];
145 eventVertex[2] = origin0[2];
146
147 // Header
148 AliGenEventHeader* header = new AliGenEventHeader("Therminator");
149 // Event Vertex
150 header->SetPrimaryVertex(eventVertex);
151 header->SetNProduced(fNprimaries);
152 gAlice->SetGenEventHeader(header);
153}
154
155void AliGenTherminator::Init()
156{
157 // Initialize global variables and
158 // particle and decay tables
159 if (fFileName.Length() == 0)
160 fFileName = "event.out";
161 ReadShareParticleTable();
162
163 SetMC(new TTherminator());
164
165 AliGenMC::Init();
166 ((TTherminator *) fMCEvGen)->Initialize();
167}
168void AliGenTherminator::SetFileName(const char *infilename)
169{
170 // Set parameter filename
171 fFileName = infilename;
172}
173void AliGenTherminator::SetEventNumberInFile(int evnum)
174{
175 // Set number of events to generate - default: 1
176 fEventNumber = evnum;
177}
178
179void AliGenTherminator::ReadShareParticleTable()
180{
181 // Read in particle table from share
182 // and add missing particle type to TDatabasePDG
183
184 char str[50];
185 char str1[200];
186
187 TDatabasePDG *tInstance = TDatabasePDG::Instance();
188 TParticlePDG *tParticleType;
189
190 AliWarning(Form("Reading particle types from particles.data"));
191 ifstream in("particles.data");
192
193 int charge;
194
195 int number=0;
196 if ((in) && (in.is_open()))
197 {
198 //START OF HEAD-LINE
199 in.ignore(200,'\n');
200 in.ignore(200,'\n');
201 in.ignore(200,'\n');
202 //END OF HEAD-LINE
203
204 while (in>>str)
205 {
206 if (/*(*str == '#')||*/(*str<65)||(*str>122))
207 {
208 in.getline(str1,200);
209 continue;
210 }
211 double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
212
213 in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
214 number++;
215 tParticleType = tInstance->GetParticle((int) mc);
216 if (!tParticleType) {
217 charge = 0;
218 if (strstr(str, "plu")) charge = 1;
219 if (strstr(str, "min")) charge = -1;
220 if (strstr(str, "plb")) charge = -1;
221 if (strstr(str, "mnb")) charge = 1;
222 if (strstr(str, "plp")) charge = 2;
223 if (strstr(str, "ppb")) charge = -2;
224 tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
225 AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
226 // 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));
227
228 }
229 }
230 in.close();
231 }
232 CreateTherminatorInputFile();
233}
234
235void AliGenTherminator::CreateTherminatorInputFile()
236{
237 // Create Therminator input file
238 ofstream *ostr = new ofstream("therminator.in");
239 (*ostr) << "NumberOfEvents = 1" << endl;
240 (*ostr) << "Randomize = 1" << endl;
241 (*ostr) << "TableType = SHARE" << endl;
242 (*ostr) << "InputDirSHARE = ." << endl;
243 (*ostr) << "EventOutputFile = " << fFileName.Data() << endl;
244 (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl;
245 (*ostr) << "BWVt = " << fBWVt << endl;
246 (*ostr) << "Tau = " << fTau << endl;
247 (*ostr) << "RhoMax = " << fRhoMax << endl;
248 (*ostr) << "Temperature = " << fTemperature << endl;
249 (*ostr) << "MiuI = " << fMiuI << endl;
250 (*ostr) << "MiuS = " << fMiuS << endl;
251 (*ostr) << "MiuB = " << fMiuB << endl;
252 (*ostr) << "AlphaRange = " << fAlfaRange << endl;
253 (*ostr) << "RapidityRange = " << fRapRange << endl;
254 (*ostr) << "NumberOfIntegrateSamples = 1000000" << endl;
255}
256
257void AliGenTherminator::SetModel(const char *model)
258{
259 // Set the freeze-out model to use
260 fFreezeOutModel = model;
261}