Update master to aliroot
[u/mrichter/AliRoot.git] / TTherminator / Therminator / Integrator.cxx
1 /******************************************************************************
2  *                      T H E R M I N A T O R                                 *
3  *                   THERMal heavy-IoN generATOR                              *
4  *                           version 1.0                                      *
5  *                                                                            *
6  * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
7  *                       Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl  *
8  * Authors of the code:  Adam Kisiel, kisiel@if.pw.edu.pl                     *
9  *                       Tomasz Taluc, ttaluc@if.pw.edu.pl                    *
10  * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski,            *
11  *                 Wojciech Florkowski                                        *
12  *                                                                            *
13  * For the detailed description of the program and furhter references         * 
14  * to the description of the model plesase refer to: nucl-th/0504047,         *
15  * accessibile at: http://www.arxiv.org/nucl-th/0504047                       *
16  *                                                                            *
17  * Homepage: http://hirg.if.pw.edu.pl/en/therminator/                         *
18  *                                                                            *
19  * This code can be freely used and redistributed. However if you decide to   *
20  * make modifications to the code, please contact the authors, especially     *
21  * if you plan to publish the results obtained with such modified code.       *
22  * Any publication of results obtained using this code must include the       *
23  * reference to nucl-th/0504047 and the published version of it, when         *
24  * available.                                                                 *
25  *                                                                            *
26  *****************************************************************************/
27 #include <fstream>
28 #include "Integrator.h"
29 #include "ParticleType.h"
30 #include "Particle.h"
31 #include <TMath.h>
32 #include <TFile.h>
33
34 extern ReadPar *sRPInstance;
35 extern int      sTables;
36 extern int      sModel;
37
38 Integrator::Integrator(int aNpart)
39 {
40   kFmToGev  = 0.197326960277;                           /*MCH updated: kFmToGev  = 0.197;*/
41   ReadParameters();
42
43   char *tHash;
44   tHash = ParameterHash();
45
46   PRINT_MESSAGE("Hash for these parameters is: " << tHash);
47   
48   mNPart = aNpart;
49   kTwoPi2 = TMath::Pi()*TMath::Pi()*2*2;                /*MCH*/
50   kTwoPi3 = TMath::Pi()*TMath::Pi()*TMath::Pi()*2*2*2;
51   mRandom = new TRandom2();
52   //  mRandom->SetSeed2(41321, 8457);
53   mRandom->SetSeed(41321);
54
55   mFOHS = new Hypersurface(mFOHSlocation.Data());                               /*MCH*/
56
57   free (tHash);
58 }
59
60 double Integrator::CalcBE(double aX)
61 {
62   if (aX>200.0) return 0.0;
63   return 1.0/(TMath::Exp(aX)-1);
64 }
65
66 double Integrator::CalcFD(double aX)
67 {
68   if (aX>200.0) return 0.0;
69   return 1.0/(TMath::Exp(aX)+1);
70 }
71
72 double Integrator::Calka(double aMass, double aMiu, 
73                          double *aRap, double *aPt, double *aPhiP, 
74                          double *aAlfaP, double *aRho, double *aPhiS, double *aTime, double aSpin,
75                          int aKeepPos)
76 {
77   // Generate momentum components
78   (*aRap) = mRandom->Rndm() * mRapRange - mRapRange/2.0;
79   double tZet   = mRandom->Rndm();
80   (*aPt) = tZet/(1-tZet);
81   (*aPhiP) = mRandom->Rndm() * TMath::Pi() * 2.0;
82   
83   // Generate position components
84   if (!aKeepPos) {
85     (*aAlfaP) = mRandom->Rndm()*mAlfaRange - mAlfaRange/2.0;
86     (*aRho) = mRandom->Rndm()*mRhoMax;
87     (*aPhiS) = mRandom->Rndm()*TMath::Pi()*2;
88     if (sModel == 8) {
89       (*aTime) = mTau - mBWDelay*TMath::Log(mRandom->Rndm());
90     }
91   }
92   
93   double tFpod = 0.0;
94   if ((sModel == 0) || (sModel == 3)) { // Single FreezeOut
95     double tPU = (TMath::Sqrt(aMass*aMass+(*aPt)*(*aPt))*
96                  cosh((*aAlfaP)-(*aRap))*
97                  TMath::Sqrt(1+(((*aRho)*(*aRho))/(mTau*mTau)))
98                  -(*aPt)*cos((*aPhiS)-(*aPhiP))*(*aRho)/mTau);
99     if (fabs(floor(aSpin) - aSpin)< 0.01) 
100       tFpod = 1/((1-tZet)*(1-tZet))*mTau*(*aPt)*(*aRho)*tPU*CalcBE((tPU- aMiu)/mTemp)/(kTwoPi3);
101     else
102       tFpod = 1/((1-tZet)*(1-tZet))*mTau*(*aPt)*(*aRho)*tPU*CalcFD((tPU- aMiu)/mTemp)/(kTwoPi3);
103   }
104 /*MCH begin*/
105   else if (sModel == 10) { // Lhyquid3D
106     double tZeta  = mRandom->Rndm()*0.5*TMath::Pi();            // angle in rho-t plane
107     double taHS   = mFOHS->fahs  ((*aPhiS),tZeta);              // velocity angle; 0.0=radial flow
108     double tvHS   = mFOHS->fvhs  ((*aPhiS),tZeta);              // velocity
109     double tdHS   = mFOHS->fdhs  ((*aPhiS),tZeta)/kFmToGev;     // distance in rho-t plane           [GeV^-1]
110     double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over (*aPhiS) [GeV^-1]
111     double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over tZeta    [GeV^-1]
112     double tTemp  = mFOHS->TFO;                                 // freeze-out temparature
113     double ttau0  = mFOHS->tau0/kFmToGev;                       // tau 0
114     double tMt    = TMath::Hypot(aMass, (*aPt));                // transverse mass
115     (*aRho)       = tdHS*cos(tZeta);                            // rho
116     double tdPt   = 1.0/((1-tZet)*(1-tZet));                    // dPt
117     (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);      // t
118
119     double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt*cosh((*aRap)-(*aAlfaP))-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
120     double tFC    = tdHS*(ttau0+tdHS*sin(tZeta))*(
121                       tdHS  *cos(tZeta)*( tMt*sin(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*cos(tZeta)*cos((*aPhiS)-(*aPhiP)))+
122                       tDzdHS*cos(tZeta)*(-tMt*cos(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*sin(tZeta)*cos((*aPhiS)-(*aPhiP)))+
123                       tDpdHS*(*aPt)*sin((*aPhiS)-(*aPhiP))
124                     );
125    if(tFC < 0.0) tFC = 0.0;
126     if (fabs(floor(aSpin)-aSpin) < 0.01)
127       tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcBE((tPU-aMiu)/tTemp);
128     else
129       tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
130
131   }
132   else if (sModel == 11) { // Lhyquid2D
133     (*aAlfaP)     = (*aRap);                                    // dirac delta (Y-ap)
134     double tZeta  = mRandom->Rndm()*0.5*TMath::Pi();            // angle in rho-t plane
135     double taHS   = mFOHS->fahs  ((*aPhiS),tZeta);              // velocity angle; 0.0=radial flow
136     double tvHS   = mFOHS->fvhs  ((*aPhiS),tZeta);              // velocity
137     double tdHS   = mFOHS->fdhs  ((*aPhiS),tZeta)/kFmToGev;     // distance in rho-t plane           [GeV^-1]
138     double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over (*aPhiS) [GeV^-1]
139     double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over tZeta    [GeV^-1]
140     double tTemp  = mFOHS->TFO;                                 // freeze-out temparature
141     double ttau0  = mFOHS->tau0/kFmToGev;                       // tau 0
142     double tMt    = TMath::Hypot(aMass, (*aPt));                // transverse mass
143     (*aRho)       = tdHS*cos(tZeta);                            // rho
144     double tdPt   = 1.0/((1-tZet)*(1-tZet));                    // dPt
145     (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);      // t
146
147     double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
148     double tFC    = tdHS*(
149                       tdHS  *cos(tZeta)*( (*aPt)/tMt*cos(tZeta)*cos((*aPhiS)-(*aPhiP))+sin(tZeta) )+
150                       tDzdHS*cos(tZeta)*( (*aPt)/tMt*sin(tZeta)*cos((*aPhiS)-(*aPhiP))-cos(tZeta) )+
151                       tDpdHS*             (*aPt)/tMt*           sin((*aPhiS)-(*aPhiP))
152                     );
153     if(tFC < 0.0) tFC = 0.0;
154     if (fabs(floor(aSpin)-aSpin) < 0.01)
155       tFpod =1.0/kTwoPi2 * tdPt*(*aPt) * tFC * CalcBE((tPU-aMiu)/tTemp);
156     else
157       tFpod =1.0/kTwoPi2 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
158   }
159 /*MCH end*/  
160   else if (sModel == 2) { // Blast-Wave with Vt
161     double tMt  = TMath::Hypot(aMass, (*aPt));
162     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
163     double tCHay = cosh((*aAlfaP)-(*aRap));
164     double tCSsp = cos((*aPhiS) - (*aPhiP));
165     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
166     double tGamma = 1.0/sqrt(1-mBWVt*mBWVt);
167     double tPU  = tMt*tGamma*tCHay-(*aPt)*mBWVt*tGamma*tCSsp - aMiu;
168     if (fabs(floor(aSpin) - aSpin)< 0.01) 
169       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
170     else
171       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
172   }
173   else if (sModel == 6) { // Blast-Wave with Vt
174     double tMt  = TMath::Hypot(aMass, (*aPt));
175     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
176     double tCHay = cosh((*aAlfaP)-(*aRap));
177     double tCSsp = cos((*aPhiS) - (*aPhiP));
178     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
179     double tGamma = 1.0/sqrt(1-mBWVt*mBWVt);
180     double tPU  = tMt*tGamma*tCHay-(*aPt)*mBWVt*tGamma*tCSsp - aMiu;
181     if (fabs(floor(aSpin) - aSpin)< 0.01) 
182       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
183     else
184       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
185   }
186   else if (sModel == 4) { // Blast-Wave with linear Vt profile
187     double tMt  = TMath::Hypot(aMass, (*aPt));
188     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
189     double tCHay = cosh((*aAlfaP)-(*aRap));
190     double tCSsp = cos((*aPhiS) - (*aPhiP));
191     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
192     double tVt = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
193     double tGamma = 1.0/sqrt(1-tVt*tVt);
194     double tPU  = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
195     if (fabs(floor(aSpin) - aSpin)< 0.01) 
196       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
197     else
198       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
199   }
200   else if (sModel == 7) { 
201     // Blast-Wave with linear Vt profile
202     // and delay in particle emission point - formation time
203     double tMt  = TMath::Hypot(aMass, (*aPt));
204     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
205     double tCHay = cosh((*aAlfaP)-(*aRap));
206     double tCSsp = cos((*aPhiS) - (*aPhiP));
207     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
208     double tVt = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
209     double tGamma = 1.0/sqrt(1-tVt*tVt);
210     double tPU  = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
211     if (fabs(floor(aSpin) - aSpin)< 0.01) 
212       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
213     else
214       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
215   }
216   else if (sModel == 8) { 
217     // Blast-Wave with linear Vt profile
218     // and delay in emission time - exponential decay
219     double tTau   = (*aTime);
220     double tMt    = TMath::Hypot(aMass, (*aPt));
221     double tPre   = (tTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
222     double tCHay  = cosh((*aAlfaP)-(*aRap));
223     double tCSsp  = cos((*aPhiS) - (*aPhiP));
224     double tBra   = tMt*tCHay - mBWA*(*aPt)*tCSsp;
225     double tVt    = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
226     double tGamma = 1.0/sqrt(1-tVt*tVt);
227     double tPU    = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
228     if (fabs(floor(aSpin) - aSpin)< 0.01) 
229       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
230     else
231       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
232   }
233 /* -=[ Testing with Mathematica ]=-*/
234 /*if ((sModel == 10)||(sModel == 11)) {
235   cout.precision(20);
236   cout << endl;
237   cout << "{phi  = " << (*aPhiS)  << ","  << endl;
238   cout << " zeta = " <<   tZeta   << ","  << endl;
239   cout << " z    = " <<   tZet    << ","  << endl;
240   cout << " phip = " << (*aPhiP)  << ","  << endl;
241   cout << " ap   = " << (*aAlfaP) << ","  << endl;
242   cout << " Y    = " << (*aRap)   << ","  << endl;
243   cout << " m    = " <<   aMass   << ","  << endl;
244   cout << " s    = " <<   aSpin   << ","  << endl;
245   cout << " mu   = " <<   aMiu    << "};" << endl;
246   cout << endl;
247   cout.precision(6);
248   cout << "TF0  = " <<   tTemp  << endl;
249   cout << "t    = " << (*aTime) << endl;
250   cout << "rho  = " << (*aRho)  << endl;
251   cout << "pT   = " << (*aPt)   << endl;
252   cout << "dpT  = " <<   tdPt   << endl;
253   cout << "mT   = " <<   tMt    << endl;
254   cout << endl;
255   cout << "aHS     = " << taHS     << endl;
256   cout << "vHS     = " << tvHS     << endl;
257   cout << "dHS     = " << tdHS     << endl;
258   cout << "DpdHS   = " << tDpdHS   << endl;
259   cout << "DzdHS   = " << tDzdHS   << endl;
260   cout << "PU      = " << tPU      << endl;
261   cout << "FC      = " << tFC      << endl;
262   cout << "part1   = " << (tPU - aMiu)/tTemp << endl;
263   cout << "part2BE = " << CalcBE( (tPU - aMiu)/tTemp ) << endl;
264   cout << "part2FD = " << CalcFD( (tPU - aMiu)/tTemp ) << endl;
265   cout << "Fpod    = " << tFpod    << endl;
266   getchar();
267 }*/
268
269   return tFpod;
270 }
271
272 double Integrator::GetMiu(double aIzo, double aBar, double aStr)
273 {
274   return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar);
275 }
276
277 double Integrator::CalcFunPodCalk(double aMass, double aIzo, double aBar, double aStr, double aSpin)
278 {
279   //  int proc = mNPart/100;
280   double tFpod;
281   double tMax = 0.0;
282   double tMiu = GetMiu(aIzo, aBar, aStr);
283   double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime; 
284   double tPt;
285
286   for (int tIpart=0; tIpart<mNPart; tIpart++)
287     {
288       tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
289       if (tFpod>tMax) tMax = tFpod;
290     }
291   return tMax;
292 }
293
294 double Integrator::Integrate(double aMass, double aIzo, double aBar, double aStr, double aSpin)
295 {
296   double tFpod=0.0;
297   double tFtest=0.0;
298   double tMiu = GetMiu(aIzo, aBar, aStr);
299   double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime; 
300   double tPt;
301
302   for (int tIpart=0; tIpart<mNPart; tIpart++)
303     {
304       tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
305
306       tFtest += tFpod;
307     }
308
309   double tCalka;
310   if (sModel == 5)
311     tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*(mTauf-mTau0)* tFtest / (mNPart*1.0);  
312 /*MCH begin*/
313   else if (sModel == 10)
314     tCalka = mRapRange*mAlfaRange*(1.0)*2.0*TMath::Pi()*2.0*TMath::Pi()*0.5*TMath::Pi()*tFtest / (mNPart*1.0);
315   else if (sModel == 11)
316     tCalka = mRapRange*(1.0)*2.0*TMath::Pi()*2.0*TMath::Pi()*0.5*TMath::Pi()*tFtest / (mNPart*1.0);
317 /*MCH end*/    
318   else
319     tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*mRhoMax* tFtest / (mNPart*1.0);  
320   return tCalka;
321 }
322
323 void 
324 Integrator::Generate(ParticleType *aPartType, int aPartCount, Particle ***oParticles)
325 {
326   int tIter = 0;
327   double tFpod;
328   double tFtest;
329   double tMiu = GetMiu(1.0*aPartType->GetI3(), 1.0*aPartType->GetBarionN(), 1.0*aPartType->GetStrangeness());
330   double tFMax;
331   double tSpin = aPartType->GetSpin();
332   double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime; 
333   double tPt;
334
335   PRINT_DEBUG_3("Gen for " << (aPartType->GetName()) << " i3 b s " << 1.0*aPartType->GetI3() << " " << 1.0*aPartType->GetBarionN() << " " << 1.0*aPartType->GetStrangeness() << " " << tMiu);
336   
337   tFMax = aPartType->GetFMax();
338
339   (*oParticles) = (Particle **) malloc(sizeof(Particle *) * aPartCount);
340   Particle *tBuf=0;
341   
342   while (tIter<aPartCount)
343     {
344       tFpod = Calka(aPartType->GetMass(),tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,tSpin);
345       tFtest = mRandom->Rndm()*tFMax;
346       if (tFtest<tFpod)
347         {
348           if ((sModel == 0) || (sModel == 3)) // Single freeze-out
349             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, TMath::Hypot(mTau, tRho), aPartType);
350 /*MCH begin*/
351           else if ((sModel == 10) || (sModel == 11)) {
352             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTime/cosh(tAlfaP), aPartType);
353           }
354 /*MCH end*/
355           else if ((sModel == 1) || (sModel == 2) || (sModel == 4)) { // Blast-wave
356             double tTau = mTau + mBWA * tRho;
357             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTau, aPartType);
358           }
359           else if (sModel == 7) { // Blast-wave
360             double tTau = mTau + mBWA * tRho;
361             double px = tPt*TMath::Cos(tPhiP);
362             double py = tPt*TMath::Sin(tPhiP);
363             double tMt = TMath::Hypot(aPartType->GetMass(),tPt);
364             double pz = tMt*TMath::SinH(tRap);
365   
366             double rx = tRho*TMath::Cos(tPhiS);
367             double ry = tRho*TMath::Sin(tPhiS);
368
369             double rz = tTau*TMath::SinH(tAlfaP);
370             double rt = tTau*TMath::CosH(tAlfaP);
371             double dt = -mBWDelay * TMath::Log(mRandom->Rndm());
372             rt += dt;
373             double en = sqrt(tMt*tMt + pz*pz);
374             rx += dt * px/en;
375             ry += dt * py/en;
376             rz += dt * pz/en;
377
378             tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
379           }
380           else if (sModel == 8) { // Blast-wave
381             double tTau = tTime + mBWA * tRho;
382             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTau, aPartType);
383           }
384           else if (sModel == 6) { // Blast-wave
385             double tTau = mTau + mBWA * tRho;
386             double px = tPt*TMath::Cos(tPhiP);
387             double py = tPt*TMath::Sin(tPhiP);
388             double tMt = TMath::Hypot(aPartType->GetMass(),tPt);
389             double pz = tMt*TMath::SinH(tRap);
390   
391             double rx = tRho*TMath::Cos(tPhiS);
392             double ry = tRho*TMath::Sin(tPhiS);
393
394             double rz = tTau*TMath::SinH(tAlfaP);
395             double rt = tTau*TMath::CosH(tAlfaP);
396             rt += -mBWDelay * TMath::Log(mRandom->Rndm());
397
398             tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
399           }
400           else if (sModel == 5) {
401             double tTau = sqrt(tTime*tTime - tAlfaP*tAlfaP);
402             double tAlfa = 0.5*log((tTime+tAlfaP)/(tTime-tAlfaP));
403             tBuf = new Particle(tRap, tPt, tPhiP, tAlfa, tRho, tPhiS, tTau, aPartType);
404           }
405           (*oParticles)[tIter] = tBuf;
406           tIter++;
407         }
408     }
409 }
410
411 void  
412 Integrator::Randomize()
413 {
414   TDatime tDat;
415
416   //  mRandom->SetSeed2(tDat.Get(), (tDat.Get() % 11) * 7 + (tDat.Get() / 7));
417   mRandom->SetSeed(tDat.Get());
418 }
419
420 void  
421 Integrator::ReadParameters()
422 {
423   STR tModel;
424   STR tTable;
425
426   // First read the model parameters
427   try {
428     
429     mRhoMax = atof(sRPInstance->getPar("RhoMax").Data()) / kFmToGev;
430     mTau    = atof(sRPInstance->getPar("Tau").Data()) / kFmToGev;
431     mTemp   = atof(sRPInstance->getPar("Temperature").Data());
432     mMiu_i  = atof(sRPInstance->getPar("MiuI").Data());
433     mMiu_s  = atof(sRPInstance->getPar("MiuS").Data());
434     mMiu_b  = atof(sRPInstance->getPar("MiuB").Data());
435     mFOHSlocation = sRPInstance->getPar("FOHSLocation");
436   }
437   catch (STR tError) {
438     PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
439     PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
440     exit(0);
441   }
442
443   // Read additional parameters for BW
444   if ((sModel == 1) || (sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
445     try {
446       mBWA  = atof(sRPInstance->getPar("BWA").Data());
447     }
448     catch (STR tError) {
449       // Using default value of 0
450       mBWA = 0.0;
451     }
452   }
453   
454   // Read additional parameters for MBWVt
455   if ((sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
456     try {
457       mBWVt  = atof(sRPInstance->getPar("BWVt").Data());
458     }
459     catch (STR tError) {
460       PRINT_DEBUG_1("Integrator::ReadParameters (BW Vt part) - Caught exception " << tError);
461       PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
462       exit(0);
463     }
464   }
465
466   // Read additional parameters for MBWVt
467   if ((sModel ==6) || (sModel == 7) || (sModel == 8))
468     try {
469       mBWDelay  = atof(sRPInstance->getPar("BWDelay").Data()) / kFmToGev;
470     }
471     catch (STR tError) {
472       PRINT_DEBUG_1("Integrator::ReadParameters (BW Vt Delay part) - Caught exception " << tError);
473       PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
474       exit(0);
475     }
476   
477   
478   // Then the integration range
479   try {
480     mAlfaRange = atof(sRPInstance->getPar("AlphaRange").Data());
481     mRapRange  = atof(sRPInstance->getPar("RapidityRange").Data());
482   }
483   catch (STR tError) {
484     PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
485     PRINT_MESSAGE("Did not find one of the neccessary integration ranges.");
486     exit(0);
487   }
488
489   // Read hydro-specific parameters
490   if (sModel == 5) {
491     try {
492       mTauf   = atof(sRPInstance->getPar("TauF").Data());
493       mTau0   = atof(sRPInstance->getPar("Tau0").Data());
494       mLambda = atof(sRPInstance->getPar("Lambda").Data());
495       mBN     = atof(sRPInstance->getPar("BN").Data());
496       mAlfa   = atof(sRPInstance->getPar("Alfa").Data());
497     }
498     catch (STR tError) {
499       PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
500       PRINT_MESSAGE("Did not find one of the neccessary hydro parameters.");
501       exit(0);
502     }
503   }
504 }
505
506 char *  
507 Integrator::ParameterHash()
508 {
509   char *tBuf;
510   
511   tBuf = (char *) malloc(sizeof(char *) * (15 * 3 + 3));
512   
513   if (sModel == 0) {
514     sprintf(tBuf, "%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mAlfaRange*100, mRapRange*100);
515   }
516 /*MCH begin*/
517   else if ((sModel == 10) || (sModel == 11)) {
518     sprintf(tBuf, "%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mAlfaRange*100, mRapRange*100);
519   }
520 /*MCH end*/
521   else if ((sModel == 2) || (sModel == 4))  {
522     sprintf(tBuf, "%1i%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sModel, sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mBWA*100, mBWVt*1000, mAlfaRange*100, mRapRange*100);
523   }
524   else if ((sModel == 6) || (sModel == 7) || (sModel == 8)) {
525     sprintf(tBuf, "%1i%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sModel, sTables, mRhoMax*100, mTau*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mBWA*100, mBWVt*1000, mAlfaRange*100, mRapRange*100, mBWDelay);
526   }
527   else if (sModel == 5)  {
528     sprintf(tBuf, "%1i%1i%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f%03.0f", sModel, sTables, mRhoMax*100, mTemp*1000, mMiu_i*10000, mMiu_s*100000, mMiu_b*10000, mTau0*10, mTauf*10, mAlfa*1000, mBN*100, mLambda*100, mAlfaRange*100, mRapRange*100);
529   }
530   
531   return tBuf;
532 }
533
534 void   
535 Integrator::ReadMultInteg(ParticleDB *aDB)
536 {
537   // Make or read table with propabilities
538   int tK;
539   char *tHash;
540   char  tIntName[1000];
541   char  tMultName[1000];
542   ifstream *tFileIn = NULL;
543   ofstream *tFileOut = NULL;
544   
545   tHash = ParameterHash();
546   
547   if (mFOHSlocation != "") {
548     strcpy(tIntName, mFOHSlocation.Data());
549     strcat(tIntName, "/");
550     strcat(tIntName, "fintegrandmax_");
551     strcat(tIntName, tHash);
552     strcat(tIntName, ".txt");
553     tFileIn = new ifstream(tIntName);
554   }
555   else if (!((tFileIn) && (tFileIn->is_open()))) {
556     strcpy(tIntName, "fintegrandmax_");
557     strcat(tIntName, tHash);
558     strcat(tIntName, ".txt");
559     tFileIn = new ifstream(tIntName);
560   }
561   
562   if ((tFileIn) && (tFileIn->is_open())) {
563     PRINT_MESSAGE("Reading Max Integrand values from " << tIntName);
564
565     char tPart[255];
566     float tRFunMax;
567     std::string tStrPart;
568     
569     while (!tFileIn->eof()) {
570       (*tFileIn) >> tPart >> tRFunMax;
571       
572       tStrPart = tPart;
573       aDB->GetParticleType(tStrPart)->SetFMax(tRFunMax);
574     }
575   }
576   else {
577     float tRFunMax;
578
579     PRINT_MESSAGE("Max Integrand file " << tIntName << " not found. Generating...");
580     
581     tFileOut = new ofstream(tIntName);
582
583     for(tK=0; tK<aDB->GetParticleTypeCount(); tK++)
584       {
585         tRFunMax = CalcFunPodCalk(aDB->GetParticleType(tK)->GetMass(),
586                                aDB->GetParticleType(tK)->GetI3(),
587                                aDB->GetParticleType(tK)->GetBarionN()*1.0,
588                                aDB->GetParticleType(tK)->GetStrangeness()*1.0,
589                                aDB->GetParticleType(tK)->GetSpin());
590         
591         (*tFileOut) << aDB->GetParticleType(tK)->GetName() << " "
592                 << tRFunMax <<endl;
593         aDB->GetParticleType(tK)->SetFMax(tRFunMax);
594         PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tK)->GetNumber()<<":"
595                       <<aDB->GetParticleType(tK)->GetName()
596                       <<" done");
597       }
598     tFileOut->close();
599   }
600   
601   // Calculate or read multiplicities
602   if (mFOHSlocation != "") {
603     strcpy(tMultName, mFOHSlocation.Data());
604     strcat(tMultName, "/");
605     strcat(tMultName, "fmultiplicity_");
606     strcat(tMultName, tHash);
607     strcat(tMultName, ".txt");
608     tFileIn = new ifstream(tMultName);
609   }
610   else if (!((tFileIn) && (tFileIn->is_open()))) {
611     strcpy(tMultName, "fmultiplicity_");
612     strcat(tMultName, tHash);
613     strcat(tMultName, ".txt");
614     tFileIn = new ifstream(tMultName);
615   }
616
617   tFileIn = new ifstream(tMultName);
618   if ((tFileIn) && (tFileIn->is_open())) {
619     PRINT_MESSAGE("Reading Multiplicities from " << tMultName);
620   }
621   else {
622     PRINT_MESSAGE("Multiplicities file " << tMultName << " not found. Generating...");
623     
624     tFileOut = new ofstream(tMultName);
625       
626     for (int tPart=0; tPart<aDB->GetParticleTypeCount(); tPart++) {
627       double tMult = Integrate(aDB->GetParticleType(tPart)->GetMass(),
628                                aDB->GetParticleType(tPart)->GetI3(),
629                                aDB->GetParticleType(tPart)->GetBarionN()*1.0,
630                                aDB->GetParticleType(tPart)->GetStrangeness()*1.0,
631                                aDB->GetParticleType(tPart)->GetSpin());
632       (*tFileOut) << (aDB->GetParticleType(tPart)->GetName()) << " " << tMult*TMath::Abs(aDB->GetParticleType(tPart)->GetSpin()*2+1) << endl;
633
634       PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tPart)->GetNumber()<<":"
635                     <<aDB->GetParticleType(tPart)->GetName()
636                     <<" done");
637     }
638     tFileOut->close();
639   }
640       
641   free (tHash);
642 }
643