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