1 /******************************************************************************
2 * T H E R M I N A T O R *
3 * THERMal heavy-IoN generATOR *
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 *
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 *
17 * Homepage: http://hirg.if.pw.edu.pl/en/therminator/ *
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 *
26 *****************************************************************************/
28 #include "Integrator.h"
29 #include "ParticleType.h"
34 extern ReadPar *sRPInstance;
38 Integrator::Integrator(int aNpart)
40 kFmToGev = 0.197326960277; /*MCH updated: kFmToGev = 0.197;*/
43 PRINT_MESSAGE("Hash for these parameters is: " << ParameterHash());
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);
52 mFOHS = new Hypersurface(mFOHSlocation.Data()); /*MCH*/
55 double Integrator::CalcBE(double aX)
57 if (aX>200.0) return 0.0;
58 return 1.0/(TMath::Exp(aX)-1);
61 double Integrator::CalcFD(double aX)
63 if (aX>200.0) return 0.0;
64 return 1.0/(TMath::Exp(aX)+1);
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,
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;
78 // Generate position components
80 (*aAlfaP) = mRandom->Rndm()*mAlfaRange - mAlfaRange/2.0;
81 (*aRho) = mRandom->Rndm()*mRhoMax;
82 (*aPhiS) = mRandom->Rndm()*TMath::Pi()*2;
84 (*aTime) = mTau - mBWDelay*TMath::Log(mRandom->Rndm());
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);
97 tFpod = 1/((1-tZet)*(1-tZet))*mTau*(*aPt)*(*aRho)*tPU*CalcFD((tPU- aMiu)/mTemp)/(kTwoPi3);
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
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))
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);
124 tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
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
142 double tPU = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
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))
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);
152 tFpod =1.0/kTwoPi2 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
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);
166 tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
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);
179 tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
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);
193 tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
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);
209 tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
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);
226 tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
228 /* -=[ Testing with Mathematica ]=-*/
229 /*if ((sModel == 10)||(sModel == 11)) {
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;
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;
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;
267 double Integrator::GetMiu(double aIzo, double aBar, double aStr)
269 return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar);
272 double Integrator::CalcFunPodCalk(double aMass, double aIzo, double aBar, double aStr, double aSpin)
274 // int proc = mNPart/100;
277 double tMiu = GetMiu(aIzo, aBar, aStr);
278 double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime;
281 for (int tIpart=0; tIpart<mNPart; tIpart++)
283 tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
284 if (tFpod>tMax) tMax = tFpod;
289 double Integrator::Integrate(double aMass, double aIzo, double aBar, double aStr, double aSpin)
293 double tMiu = GetMiu(aIzo, aBar, aStr);
294 double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime;
297 for (int tIpart=0; tIpart<mNPart; tIpart++)
299 tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
306 tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*(mTauf-mTau0)* tFtest / (mNPart*1.0);
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);
314 tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*mRhoMax* tFtest / (mNPart*1.0);
319 Integrator::Generate(ParticleType *aPartType, int aPartCount, Particle ***oParticles)
324 double tMiu = GetMiu(1.0*aPartType->GetI3(), 1.0*aPartType->GetBarionN(), 1.0*aPartType->GetStrangeness());
326 double tSpin = aPartType->GetSpin();
327 double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime;
330 PRINT_DEBUG_3("Gen for " << (aPartType->GetName()) << " i3 b s " << 1.0*aPartType->GetI3() << " " << 1.0*aPartType->GetBarionN() << " " << 1.0*aPartType->GetStrangeness() << " " << tMiu);
332 tFMax = aPartType->GetFMax();
334 (*oParticles) = (Particle **) malloc(sizeof(Particle *) * aPartCount);
337 while (tIter<aPartCount)
339 tFpod = Calka(aPartType->GetMass(),tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,tSpin);
340 tFtest = mRandom->Rndm()*tFMax;
343 if ((sModel == 0) || (sModel == 3)) // Single freeze-out
344 tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, TMath::Hypot(mTau, tRho), aPartType);
346 else if ((sModel == 10) || (sModel == 11)) {
347 tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTime/cosh(tAlfaP), aPartType);
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);
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);
361 double rx = tRho*TMath::Cos(tPhiS);
362 double ry = tRho*TMath::Sin(tPhiS);
364 double rz = tTau*TMath::SinH(tAlfaP);
365 double rt = tTau*TMath::CosH(tAlfaP);
366 double dt = -mBWDelay * TMath::Log(mRandom->Rndm());
368 double en = sqrt(tMt*tMt + pz*pz);
373 tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
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);
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);
386 double rx = tRho*TMath::Cos(tPhiS);
387 double ry = tRho*TMath::Sin(tPhiS);
389 double rz = tTau*TMath::SinH(tAlfaP);
390 double rt = tTau*TMath::CosH(tAlfaP);
391 rt += -mBWDelay * TMath::Log(mRandom->Rndm());
393 tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
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);
400 (*oParticles)[tIter] = tBuf;
407 Integrator::Randomize()
411 // mRandom->SetSeed2(tDat.Get(), (tDat.Get() % 11) * 7 + (tDat.Get() / 7));
412 mRandom->SetSeed(tDat.Get());
416 Integrator::ReadParameters()
421 // First read the model parameters
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");
433 PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
434 PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
438 // Read additional parameters for BW
439 if ((sModel == 1) || (sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
441 mBWA = atof(sRPInstance->getPar("BWA").Data());
444 // Using default value of 0
449 // Read additional parameters for MBWVt
450 if ((sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
452 mBWVt = atof(sRPInstance->getPar("BWVt").Data());
455 PRINT_DEBUG_1("Integrator::ReadParameters (BW Vt part) - Caught exception " << tError);
456 PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
461 // Read additional parameters for MBWVt
462 if ((sModel ==6) || (sModel == 7) || (sModel == 8))
464 mBWDelay = atof(sRPInstance->getPar("BWDelay").Data()) / kFmToGev;
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.");
473 // Then the integration range
475 mAlfaRange = atof(sRPInstance->getPar("AlphaRange").Data());
476 mRapRange = atof(sRPInstance->getPar("RapidityRange").Data());
479 PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
480 PRINT_MESSAGE("Did not find one of the neccessary integration ranges.");
484 // Read hydro-specific parameters
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());
494 PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
495 PRINT_MESSAGE("Did not find one of the neccessary hydro parameters.");
502 Integrator::ParameterHash()
506 tBuf = (char *) malloc(sizeof(char *) * (15 * 3 + 3));
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);
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);
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);
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);
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);
530 Integrator::ReadMultInteg(ParticleDB *aDB)
532 // Make or read table with propabilities
536 char tMultName[1000];
537 ifstream *tFileIn = NULL;
538 ofstream *tFileOut = NULL;
540 tHash = ParameterHash();
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);
550 else if (!((tFileIn) && (tFileIn->is_open()))) {
551 strcpy(tIntName, "fintegrandmax_");
552 strcat(tIntName, tHash);
553 strcat(tIntName, ".txt");
554 tFileIn = new ifstream(tIntName);
557 if ((tFileIn) && (tFileIn->is_open())) {
558 PRINT_MESSAGE("Reading Max Integrand values from " << tIntName);
562 std::string tStrPart;
564 while (!tFileIn->eof()) {
565 (*tFileIn) >> tPart >> tRFunMax;
568 aDB->GetParticleType(tStrPart)->SetFMax(tRFunMax);
574 PRINT_MESSAGE("Max Integrand file " << tIntName << " not found. Generating...");
576 tFileOut = new ofstream(tIntName);
578 for(tK=0; tK<aDB->GetParticleTypeCount(); tK++)
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());
586 (*tFileOut) << aDB->GetParticleType(tK)->GetName() << " "
588 aDB->GetParticleType(tK)->SetFMax(tRFunMax);
589 PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tK)->GetNumber()<<":"
590 <<aDB->GetParticleType(tK)->GetName()
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);
605 else if (!((tFileIn) && (tFileIn->is_open()))) {
606 strcpy(tMultName, "fmultiplicity_");
607 strcat(tMultName, tHash);
608 strcat(tMultName, ".txt");
609 tFileIn = new ifstream(tMultName);
612 tFileIn = new ifstream(tMultName);
613 if ((tFileIn) && (tFileIn->is_open())) {
614 PRINT_MESSAGE("Reading Multiplicities from " << tMultName);
617 PRINT_MESSAGE("Multiplicities file " << tMultName << " not found. Generating...");
619 tFileOut = new ofstream(tMultName);
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;
629 PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tPart)->GetNumber()<<":"
630 <<aDB->GetParticleType(tPart)->GetName()