Fix in the initialization (Solaris x86)
[u/mrichter/AliRoot.git] / TTherminator / TTherminator.cxx
1 ///////////////////////////////////////////////////////////////////////////////
2 //                                                                           //
3 // TTherminator: global interface class to the Therminator model.            //
4 // Initialized global variables and runs the event generation                //
5 //                                                                           //
6 ///////////////////////////////////////////////////////////////////////////////
7 /******************************************************************************
8  *                      T H E R M I N A T O R                                 *
9  *                   THERMal heavy-IoN generATOR                              *
10  *                           version 1.0                                      *
11  *                                                                            *
12  * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
13  *                       Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl  *
14  * Authors of the code:  Adam Kisiel, kisiel@if.pw.edu.pl                     *
15  *                       Tomasz Taluc, ttaluc@if.pw.edu.pl                    *
16  * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski,            *
17  *                 Wojciech Florkowski                                        *
18  *                                                                            *
19  * For the detailed description of the program and furhter references         * 
20  * to the description of the model plesase refer to: nucl-th/0504047,         *
21  * accessibile at: http://www.arxiv.org/nucl-th/0504047                       *
22  *                                                                            *
23  * Homepage: http://hirg.if.pw.edu.pl/en/therminator/                         *
24  *                                                                            *
25  * This code can be freely used and redistributed. However if you decide to   *
26  * make modifications to the code, please contact the authors, especially     *
27  * if you plan to publish the results obtained with such modified code.       *
28  * Any publication of results obtained using this code must include the       *
29  * reference to nucl-th/0504047 and the published version of it, when         *
30  * available.                                                                 *
31  *                                                                            *
32  *****************************************************************************/
33 #include <TTherminator.h>
34 #include <fstream>
35 #include <TDatabasePDG.h>
36 #include <TParticle.h>
37 #include <TClonesArray.h>
38
39 ReadPar *sRPInstance;
40 STR      sRPFileName;
41 int      sNumEvents;
42 int      sRunType; 
43 int      sTables;
44 int      sModel;
45 int      sIntegrateSample;
46
47 ClassImp(TTherminator)
48
49 TTherminator::TTherminator():
50   fCalka(0),
51   fEvent(0),
52   fPartDB(0)
53 {
54   // Default constructor
55   fPartDB = new ParticleDB();
56 }
57 TTherminator::TTherminator(const TTherminator & therm) :
58   fCalka(0),
59   fEvent(0),
60   fPartDB(0)
61 {
62   // Copy constructor
63   fPartDB = new ParticleDB();
64 }
65 TTherminator::~TTherminator()
66 {
67   // Destructor
68   if (sRPInstance) delete sRPInstance;
69   if (fEvent) delete fEvent;
70   if (fCalka) delete fCalka;
71   if (fPartDB) delete fPartDB;
72 }
73
74 void TTherminator::ReadParameters()
75 {
76   // Read in global model parameters
77   STR tModel;
78   STR tTable;
79   
80   try {
81     sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
82     tModel = sRPInstance->getPar("FreezeOutModel");
83     if (tModel.Contains("SingleFreezeOut"))
84       sModel = 0;
85 /*MCH begin*/
86     else if (tModel.Contains("Lhyquid3D"))
87       sModel = 10;
88     else if (tModel.Contains("Lhyquid2D"))
89       sModel = 11;
90 /*MCH end*/      
91     else if (tModel.Contains("BlastWaveVTDelay"))
92       sModel = 6;
93     else if (tModel.Contains("BlastWaveVT"))
94       sModel = 2;
95     else if (tModel.Contains("BlastWaveVLinearFormation"))
96       sModel = 7;
97     else if (tModel.Contains("BlastWaveVLinearDelay"))
98       sModel = 8;
99     else if (tModel.Contains("BlastWaveVLinear"))
100       sModel = 4;
101     else if (tModel.Contains("FiniteHydro"))
102       sModel = 5;
103     else {
104       PRINT_MESSAGE("Unknown model type: " << tModel.Data());
105       PRINT_MESSAGE("Please provide the proper name");
106       exit(0);
107     }
108     if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
109       sRunType = 4;
110     else
111       sRunType = 3;
112     tTable = sRPInstance->getPar("TableType");
113     if (tTable.Contains("Mathematica"))
114       sTables = 0;
115     else if (tTable.Contains("SHARE"))
116       sTables = 1;
117     else {
118       PRINT_MESSAGE("Unknown table type: " << tTable.Data());
119       exit(0);
120     }
121     sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
122   }
123   catch (STR tError) {
124     PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
125     PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
126     exit(0);
127   }
128 }
129
130 void        TTherminator::Initialize(){
131   // Initialize global variables
132   // and particle and decay tables
133   Parser       *tParser;
134   
135   sRPFileName = "therminator.in";
136   sRPInstance = new ReadPar(sRPFileName.Data());
137
138   ReadParameters();
139   
140   tParser = new Parser(fPartDB);
141
142   tParser->ReadShare();
143     
144   fCalka = new Integrator(sIntegrateSample);
145   fCalka->ReadMultInteg(fPartDB);
146
147
148   // Read in particle table from share
149   // and add missing particle type to TDatabasePDG
150
151   char str[50];
152   char str1[200];
153     
154   //  ParticleType *tPartBuf;
155
156   TDatabasePDG *tInstance = TDatabasePDG::Instance();
157   TParticlePDG *tParticleType;
158
159   //  AliWarning(Form("Reading particle types from particles.data"));
160   ifstream in("particles.data");
161   
162   int charge;
163     
164   int number=0;
165   if ((in) && (in.is_open()))
166     {
167       //START OF HEAD-LINE
168       in.ignore(200,'\n');
169       in.ignore(200,'\n');
170       in.ignore(200,'\n');
171       //END OF HEAD-LINE
172       
173       while (in>>str)
174         {
175           if (/*(*str == '#')||*/(*str<65)||(*str>122))
176             {
177               in.getline(str1,200);
178               continue;
179             }
180           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
181           
182           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
183           number++;
184           tParticleType = tInstance->GetParticle((int) mc);
185           if (!tParticleType) {
186             charge = 0;
187             if (strstr(str, "plu")) charge = 1;
188             if (strstr(str, "min")) charge = -1;
189             if (strstr(str, "plb")) charge = -1;
190             if (strstr(str, "mnb")) charge = 1;
191             if (strstr(str, "plp")) charge = 2;
192             if (strstr(str, "ppb")) charge = -2;
193             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
194             //      printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
195             //AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
196             //      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));
197
198           }
199         }
200       in.close();
201     }
202 }
203
204 void        TTherminator::GenerateEvent()
205 {
206   // Run single event generation
207   if ((sRunType == 3) || (sRunType == 4))
208     {
209       if (sRunType==4)
210         fCalka->Randomize();
211
212       if (!fEvent) fEvent = new Event(fPartDB, fCalka);
213       if (sRunType == 4) fEvent->Randomize();
214       
215       fEvent->GenerateEvent();
216       fEvent->DecayParticles();
217       //      fEvent->WriteEvent(0);
218     }  
219 }
220
221 Int_t       TTherminator::ImportParticles(TClonesArray *particles, Option_t *option)
222 {
223   // Import particles from a generated event into an external array
224   if (particles == 0) return 0;
225   TClonesArray &particlesR = *particles;
226   particlesR.Clear();
227   Int_t nump = 0;
228   if (!fEvent) return 0;
229   Int_t numpart = fEvent->GetParticleCount();
230   printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
231   for (Int_t iPart=0; iPart<numpart; iPart++) {
232     Particle *tPart = fEvent->GetParticleOfCount(iPart);
233     Int_t tFather;
234     tFather = tPart->GetFather();
235
236     if (tFather> -1) {
237       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));        
238       mother->SetLastDaughter(iPart);
239       if (mother->GetFirstDaughter()==-1)
240         mother->SetFirstDaughter(iPart);
241     }
242     nump++;
243     //    printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
244     new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
245                                       tFather, -1, -1, -1,
246                                       tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
247                                       tPart->rx*1.e-13, tPart->ry*1.e-13, tPart->rz*1.e-13, tPart->rt*1.e-13/300000000 
248                                       );
249     particlesR[iPart]->SetUniqueID(iPart);
250   }
251   
252   return nump;
253 }
254
255 TObjArray*  TTherminator::ImportParticles(Option_t *option)
256 {
257   // Import generated particles into an internal array
258   Int_t nump = 0;
259   fParticles->Clear();
260   if (!fEvent) return 0;
261   Int_t numpart = fEvent->GetParticleCount();
262   printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
263   fEvent->InitParticleScan();
264   for (Int_t iPart=0; iPart<numpart; iPart++) {
265     Particle *tPart = fEvent->GetNextParticle();
266     Int_t tFather;
267     tFather = tPart->GetFather();
268     
269     if (tFather> -1) {
270       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));         
271       mother->SetLastDaughter(iPart);
272       if (mother->GetFirstDaughter()==-1)
273         mother->SetFirstDaughter(iPart);
274     }
275     TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
276                                  tFather, -1, -1, -1,
277                                  tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
278                                  tPart->rx*1.e-13, tPart->ry*1.e-13, tPart->rz*1.e-13, tPart->rt*1.e-13/300000000 
279                                  );
280     p->SetUniqueID(iPart);
281     fParticles->Add(p);
282   }
283   nump = fEvent->GetParticleCount();
284   
285   return fParticles;
286   
287 }