Fixing error on gcc 4.5.1
[u/mrichter/AliRoot.git] / TTherminator / Therminator / Integrator.cxx
CommitLineData
2e967919 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
34extern ReadPar *sRPInstance;
35extern int sTables;
36extern int sModel;
37
38Integrator::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
ae89e57a 52 mFOHS = new Hypersurface(mFOHSlocation.Data()); /*MCH*/
2e967919 53}
54
55double Integrator::CalcBE(double aX)
56{
8b039ab1 57 if (aX>200.0) return 0.0;
2e967919 58 return 1.0/(TMath::Exp(aX)-1);
59}
60
61double Integrator::CalcFD(double aX)
62{
8b039ab1 63 if (aX>200.0) return 0.0;
2e967919 64 return 1.0/(TMath::Exp(aX)+1);
65}
66
67double 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
9abb118d 108 double ttau0 = mFOHS->tau0/kFmToGev; // tau 0
2e967919 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
9abb118d 112 (*aTime) = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP); // t
2e967919 113
114 double tPU = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt*cosh((*aRap)-(*aAlfaP))-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
9abb118d 115 double tFC = tdHS*(ttau0+tdHS*sin(tZeta))*(
2e967919 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 );
9abb118d 120 if(tFC < 0.0) tFC = 0.0;
2e967919 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);
9abb118d 125
2e967919 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
9abb118d 136 double ttau0 = mFOHS->tau0/kFmToGev; // tau 0
2e967919 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
9abb118d 140 (*aTime) = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP); // t
2e967919 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
267double Integrator::GetMiu(double aIzo, double aBar, double aStr)
268{
269 return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar);
270}
271
272double 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
289double 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
318void
319Integrator::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);
c61a7285 335 Particle *tBuf=0;
2e967919 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
406void
407Integrator::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
415void
416Integrator::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());
ae89e57a 430 mFOHSlocation = sRPInstance->getPar("FOHSLocation");
2e967919 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
501char *
502Integrator::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
529void
530Integrator::ReadMultInteg(ParticleDB *aDB)
531{
532 // Make or read table with propabilities
533 int tK;
534 char *tHash;
ae89e57a 535 char tIntName[1000];
536 char tMultName[1000];
2e967919 537 ifstream *tFileIn = NULL;
538 ofstream *tFileOut = NULL;
539
540 tHash = ParameterHash();
541
ae89e57a 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 }
2e967919 556
2e967919 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
ae89e57a 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 }
2e967919 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