Double check if SM is running added. Some redundant output removed from SM
[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 #include <TClass.h>
39
40 ReadPar *sRPInstance;
41 STR      sRPFileName;
42 int      sNumEvents;
43 int      sRunType; 
44 int      sTables;
45 int      sModel;
46 int      sIntegrateSample;
47
48 ClassImp(TTherminator)
49
50 TTherminator::TTherminator():
51   fCalka(0),
52   fEvent(0),
53   fPartDB(0)
54 {
55   // Default constructor
56   fPartDB = new ParticleDB();
57 }
58 TTherminator::TTherminator(const TTherminator & therm) :
59   TGenerator(therm), 
60   fCalka(0),
61   fEvent(0),
62   fPartDB(0)
63 {
64   // Copy constructor
65   //  fPartDB = new ParticleDB();
66   if (fCalka) delete fCalka;
67   fCalka = new Integrator(*therm.fCalka);
68   if (fEvent) delete fEvent;
69   fEvent = new Event(*therm.fEvent);
70   if (fPartDB) delete fPartDB;
71   fPartDB = new ParticleDB();
72 }
73 TTherminator& TTherminator::operator=(const TTherminator & therm)
74 {
75   if (this != &therm) {
76     fCalka = therm.fCalka;
77     fEvent = therm.fEvent;
78     delete fPartDB;
79     fPartDB = new ParticleDB();
80   }
81
82   return *this;
83 }
84
85 TTherminator::~TTherminator()
86 {
87   // Destructor
88   if (sRPInstance) delete sRPInstance;
89   if (fEvent) delete fEvent;
90   if (fCalka) delete fCalka;
91   if (fPartDB) delete fPartDB;
92 }
93
94 void TTherminator::ReadParameters()
95 {
96   // Read in global model parameters
97   STR tModel;
98   STR tTable;
99   
100   try {
101     sNumEvents = atoi(sRPInstance->getPar("NumberOfEvents").Data());
102     tModel = sRPInstance->getPar("FreezeOutModel");
103     if (tModel.Contains("SingleFreezeOut"))
104       sModel = 0;
105 /*MCH begin*/
106     else if (tModel.Contains("Lhyquid3D"))
107       sModel = 10;
108     else if (tModel.Contains("Lhyquid2D"))
109       sModel = 11;
110 /*MCH end*/      
111     else if (tModel.Contains("BlastWaveVTDelay"))
112       sModel = 6;
113     else if (tModel.Contains("BlastWaveVT"))
114       sModel = 2;
115     else if (tModel.Contains("BlastWaveVLinearFormation"))
116       sModel = 7;
117     else if (tModel.Contains("BlastWaveVLinearDelay"))
118       sModel = 8;
119     else if (tModel.Contains("BlastWaveVLinear"))
120       sModel = 4;
121     else if (tModel.Contains("FiniteHydro"))
122       sModel = 5;
123     else {
124       PRINT_MESSAGE("Unknown model type: " << tModel.Data());
125       PRINT_MESSAGE("Please provide the proper name");
126       exit(0);
127     }
128     if (atoi(sRPInstance->getPar("Randomize").Data()) == 1)
129       sRunType = 4;
130     else
131       sRunType = 3;
132     tTable = sRPInstance->getPar("TableType");
133     if (tTable.Contains("Mathematica"))
134       sTables = 0;
135     else if (tTable.Contains("SHARE"))
136       sTables = 1;
137     else {
138       PRINT_MESSAGE("Unknown table type: " << tTable.Data());
139       exit(0);
140     }
141     sIntegrateSample = atoi(sRPInstance->getPar("NumberOfIntegrateSamples").Data());
142     fInputDir = sRPInstance->getPar("InputDirSHARE");
143   }
144   catch (STR tError) {
145     PRINT_DEBUG_1("RunBFPW::ReadParameters - Caught exception " << tError);
146     PRINT_MESSAGE("Did not find one of the neccessary parameters in the parameters file.");
147     exit(0);
148   }
149 }
150
151 void        TTherminator::Initialize(){
152   // Initialize global variables
153   // and particle and decay tables
154   Parser       *tParser;
155   
156   sRPFileName = "therminator.in";
157   sRPInstance = new ReadPar(sRPFileName.Data());
158
159   ReadParameters();
160   
161   tParser = new Parser(fPartDB);
162
163   tParser->ReadShare();
164     
165   fCalka = new Integrator(sIntegrateSample);
166   fCalka->ReadMultInteg(fPartDB);
167
168
169   // Read in particle table from share
170   // and add missing particle type to TDatabasePDG
171
172   char str[50];
173   char str1[200];
174     
175   //  ParticleType *tPartBuf;
176
177   TDatabasePDG *tInstance = TDatabasePDG::Instance();
178   TParticlePDG *tParticleType;
179
180   //  AliWarning(Form("Reading particle types from particles.data"));
181
182   ifstream in((fInputDir+"/"+"particles.data").Data());
183   //  ifstream in("particles.data");
184   
185   int charge;
186     
187   int number=0;
188   if ((in) && (in.is_open()))
189     {
190       //START OF HEAD-LINE
191       in.ignore(200,'\n');
192       in.ignore(200,'\n');
193       in.ignore(200,'\n');
194       //END OF HEAD-LINE
195       
196       while (in>>str)
197         {
198           if (/*(*str == '#')||*/(*str<65)||(*str>122))
199             {
200               in.getline(str1,200);
201               continue;
202             }
203           double mass, gamma, spin, tI3, tI, q, s, aq, as, c, ac, mc;
204           
205           in>>mass>>gamma>>spin>>tI>>tI3>>q>>s>>aq>>as>>c>>ac>>mc;
206           number++;
207           tParticleType = tInstance->GetParticle((int) mc);
208           if (!tParticleType) {
209             charge = 0;
210             if (strstr(str, "plu")) charge = 1;
211             if (strstr(str, "min")) charge = -1;
212             if (strstr(str, "plb")) charge = -1;
213             if (strstr(str, "mnb")) charge = 1;
214             if (strstr(str, "plp")) charge = 2;
215             if (strstr(str, "ppb")) charge = -2;
216             tInstance->AddParticle(str, str, mass, gamma == 0.0 ? 1:0, gamma, charge , "meson", (int) mc);
217             //      printf("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge);
218             //AliWarning(Form("Added particle %s with PDG PID %d charge %d", str, (int) mc, charge));
219             //      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));
220
221           }
222         }
223       in.close();
224     }
225
226   delete tParser;
227 }
228
229 void        TTherminator::GenerateEvent()
230 {
231   // Run single event generation
232   if ((sRunType == 3) || (sRunType == 4))
233     {
234       if (sRunType==4)
235         fCalka->Randomize();
236
237       if (!fEvent) fEvent = new Event(fPartDB, fCalka);
238       if (sRunType == 4) fEvent->Randomize();
239       
240       fEvent->GenerateEvent();
241       fEvent->DecayParticles();
242       //      fEvent->WriteEvent(0);
243     }  
244 }
245
246 Int_t       TTherminator::ImportParticles(TClonesArray *particles, Option_t */*option*/)
247 {
248   // Import particles from a generated event into an external array
249   const double kFmToGev = 0.197327;
250
251   if (particles == 0) return 0;
252   TClonesArray &particlesR = *particles;
253   particlesR.Clear();
254   Int_t nump = 0;
255   if (!fEvent) return 0;
256   Int_t numpart = fEvent->GetParticleCount();
257   printf("\n TTherminator: Therminator stack contains %d particles.\n", numpart);
258   for (Int_t iPart=0; iPart<numpart; iPart++) {
259     Particle *tPart = fEvent->GetParticleOfCount(iPart);
260     Int_t tFather;
261     tFather = tPart->GetFather();
262
263     if (tFather> -1) {
264       TParticle *mother = (TParticle*) (particlesR.UncheckedAt(tFather+1));        
265       mother->SetLastDaughter(iPart);
266       if (mother->GetFirstDaughter()==-1)
267         mother->SetFirstDaughter(iPart);
268     }
269     nump++;
270     //    printf("Putting %d %d %lf %d\n", tPart->GetParticleType()->GetPDGCode(), iPart, tPart->px, tPart);
271     new (particlesR[iPart]) TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
272                                       tFather, -1, -1, -1,
273                                       tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
274                                       tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
275     particlesR[iPart]->SetUniqueID(iPart);
276   }
277   
278   return nump;
279 }
280
281 TObjArray*  TTherminator::ImportParticles(Option_t */*option*/)
282 {
283   // Import generated particles into an internal array
284   const double kFmToGev = 0.197327;
285
286   Int_t nump = 0;
287   fParticles->Clear();
288   if (!fEvent) return 0;
289   Int_t numpart = fEvent->GetParticleCount();
290   printf("\n TTherminator: Therminator stack contains %d particles.", numpart);
291   fEvent->InitParticleScan();
292   for (Int_t iPart=0; iPart<numpart; iPart++) {
293     Particle *tPart = fEvent->GetNextParticle();
294     Int_t tFather;
295     tFather = tPart->GetFather();
296     
297     if (tFather> -1) {
298       TParticle *mother = (TParticle*) (fParticles->UncheckedAt(tFather));         
299       mother->SetLastDaughter(iPart);
300       if (mother->GetFirstDaughter()==-1)
301         mother->SetFirstDaughter(iPart);
302     }
303     TParticle* p = new TParticle(tPart->GetParticleType()->GetPDGCode(), tPart->HadDecayed(),
304                                  tFather, -1, -1, -1,
305                                  tPart->px, tPart->py, tPart->pz, tPart->GetEnergy() ,
306                                  tPart->rx*1.e-13*kFmToGev, tPart->ry*1.e-13*kFmToGev, tPart->rz*1.e-13*kFmToGev, tPart->rt*1.e-13*kFmToGev/3e10);
307     p->SetUniqueID(iPart);
308     fParticles->Add(p);
309   }
310   nump = fEvent->GetParticleCount();
311   
312   return fParticles;
313   
314 }