| 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 | |
| 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 | |
| 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 | |
| 60 | double 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 | |
| 66 | double 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 | |
| 72 | double 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 | |
| 272 | double Integrator::GetMiu(double aIzo, double aBar, double aStr) |
| 273 | { |
| 274 | return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar); |
| 275 | } |
| 276 | |
| 277 | double 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 | |
| 294 | double 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 | |
| 323 | void |
| 324 | Integrator::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 | |
| 411 | void |
| 412 | Integrator::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 | |
| 420 | void |
| 421 | Integrator::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 | |
| 506 | char * |
| 507 | Integrator::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 | |
| 534 | void |
| 535 | Integrator::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 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame