Update master to aliroot
[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
6750cf94 43 char *tHash;
44 tHash = ParameterHash();
45
46 PRINT_MESSAGE("Hash for these parameters is: " << tHash);
2e967919 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
ae89e57a 55 mFOHS = new Hypersurface(mFOHSlocation.Data()); /*MCH*/
6750cf94 56
57 free (tHash);
2e967919 58}
59
60double Integrator::CalcBE(double aX)
61{
8b039ab1 62 if (aX>200.0) return 0.0;
2e967919 63 return 1.0/(TMath::Exp(aX)-1);
64}
65
66double Integrator::CalcFD(double aX)
67{
8b039ab1 68 if (aX>200.0) return 0.0;
2e967919 69 return 1.0/(TMath::Exp(aX)+1);
70}
71
72double 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
9abb118d 113 double ttau0 = mFOHS->tau0/kFmToGev; // tau 0
2e967919 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
9abb118d 117 (*aTime) = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP); // t
2e967919 118
119 double tPU = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt*cosh((*aRap)-(*aAlfaP))-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
9abb118d 120 double tFC = tdHS*(ttau0+tdHS*sin(tZeta))*(
2e967919 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 );
9abb118d 125 if(tFC < 0.0) tFC = 0.0;
2e967919 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);
9abb118d 130
2e967919 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
9abb118d 141 double ttau0 = mFOHS->tau0/kFmToGev; // tau 0
2e967919 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
9abb118d 145 (*aTime) = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP); // t
2e967919 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
272double Integrator::GetMiu(double aIzo, double aBar, double aStr)
273{
274 return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar);
275}
276
277double 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
294double 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
323void
324Integrator::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);
c61a7285 340 Particle *tBuf=0;
2e967919 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
411void
412Integrator::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
420void
421Integrator::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());
ae89e57a 435 mFOHSlocation = sRPInstance->getPar("FOHSLocation");
2e967919 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
506char *
507Integrator::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
534void
535Integrator::ReadMultInteg(ParticleDB *aDB)
536{
537 // Make or read table with propabilities
538 int tK;
539 char *tHash;
ae89e57a 540 char tIntName[1000];
541 char tMultName[1000];
2e967919 542 ifstream *tFileIn = NULL;
543 ofstream *tFileOut = NULL;
544
545 tHash = ParameterHash();
546
ae89e57a 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 }
2e967919 561
2e967919 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
ae89e57a 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 }
2e967919 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