TTherminator/Therminator/Integrator.cxx

   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   char *tHash;
  44   tHash = ParameterHash();
  45
  46   PRINT_MESSAGE("Hash for these parameters is: " << tHash);
  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
  55   mFOHS = new Hypersurface(mFOHSlocation.Data());                               /*MCH*/
  56
  57   free (tHash);
  58 }
  59
  60 double Integrator::CalcBE(double aX)
  61 {
  62   if (aX>200.0) return 0.0;
  63   return 1.0/(TMath::Exp(aX)-1);
  64 }
  65
  66 double Integrator::CalcFD(double aX)
  67 {
  68   if (aX>200.0) return 0.0;
  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
 113     double ttau0  = mFOHS->tau0/kFmToGev;                       // tau 0
 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
 117     (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);      // t
 118
 119     double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt*cosh((*aRap)-(*aAlfaP))-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
 120     double tFC    = tdHS*(ttau0+tdHS*sin(tZeta))*(
 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                     );
 125    if(tFC < 0.0) tFC = 0.0;
 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);
 130
 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
 141     double ttau0  = mFOHS->tau0/kFmToGev;                       // tau 0
 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
 145     (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);      // t
 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);
 340   Particle *tBuf=0;
 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());
 435     mFOHSlocation = sRPInstance->getPar("FOHSLocation");
 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;
 540   char  tIntName[1000];
 541   char  tMultName[1000];
 542   ifstream *tFileIn = NULL;
 543   ofstream *tFileOut = NULL;
 544
 545   tHash = ParameterHash();
 546
 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   }
 561
 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
 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   }
 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