| 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 | |
| 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 | |
| 55 | double 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 | |
| 61 | double 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 | |
| 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, |
| 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 | |
| 267 | double Integrator::GetMiu(double aIzo, double aBar, double aStr) |
| 268 | { |
| 269 | return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar); |
| 270 | } |
| 271 | |
| 272 | double 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 | |
| 289 | double 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 | |
| 318 | void |
| 319 | Integrator::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 | |
| 406 | void |
| 407 | Integrator::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 | |
| 415 | void |
| 416 | Integrator::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 | |
| 501 | char * |
| 502 | Integrator::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 | |
| 529 | void |
| 530 | Integrator::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 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame