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