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   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
  52   mFOHS = new Hypersurface();                           /*MCH*/
  53 }
  54
  55 double Integrator::CalcBE(double aX)
  56 {
  57   return 1.0/(TMath::Exp(aX)-1);
  58 }
  59
  60 double Integrator::CalcFD(double aX)
  61 {
  62   return 1.0/(TMath::Exp(aX)+1);
  63 }
  64
  65 double Integrator::Calka(double aMass, double aMiu,
  66                          double *aRap, double *aPt, double *aPhiP,
  67                          double *aAlfaP, double *aRho, double *aPhiS, double *aTime, double aSpin,
  68                          int aKeepPos)
  69 {
  70   // Generate momentum components
  71   (*aRap) = mRandom->Rndm() * mRapRange - mRapRange/2.0;
  72   double tZet   = mRandom->Rndm();
  73   (*aPt) = tZet/(1-tZet);
  74   (*aPhiP) = mRandom->Rndm() * TMath::Pi() * 2.0;
  75
  76   // Generate position components
  77   if (!aKeepPos) {
  78     (*aAlfaP) = mRandom->Rndm()*mAlfaRange - mAlfaRange/2.0;
  79     (*aRho) = mRandom->Rndm()*mRhoMax;
  80     (*aPhiS) = mRandom->Rndm()*TMath::Pi()*2;
  81     if (sModel == 8) {
  82       (*aTime) = mTau - mBWDelay*TMath::Log(mRandom->Rndm());
  83     }
  84   }
  85
  86   double tFpod = 0.0;
  87   if ((sModel == 0) || (sModel == 3)) { // Single FreezeOut
  88     double tPU = (TMath::Sqrt(aMass*aMass+(*aPt)*(*aPt))*
  89                  cosh((*aAlfaP)-(*aRap))*
  90                  TMath::Sqrt(1+(((*aRho)*(*aRho))/(mTau*mTau)))
  91                  -(*aPt)*cos((*aPhiS)-(*aPhiP))*(*aRho)/mTau);
  92     if (fabs(floor(aSpin) - aSpin)< 0.01)
  93       tFpod = 1/((1-tZet)*(1-tZet))*mTau*(*aPt)*(*aRho)*tPU*CalcBE((tPU- aMiu)/mTemp)/(kTwoPi3);
  94     else
  95       tFpod = 1/((1-tZet)*(1-tZet))*mTau*(*aPt)*(*aRho)*tPU*CalcFD((tPU- aMiu)/mTemp)/(kTwoPi3);
  96   }
  97 /*MCH begin*/
  98   else if (sModel == 10) { // Lhyquid3D
  99     double tZeta  = mRandom->Rndm()*0.5*TMath::Pi();            // angle in rho-t plane
 100     double taHS   = mFOHS->fahs  ((*aPhiS),tZeta);              // velocity angle; 0.0=radial flow
 101     double tvHS   = mFOHS->fvhs  ((*aPhiS),tZeta);              // velocity
 102     double tdHS   = mFOHS->fdhs  ((*aPhiS),tZeta)/kFmToGev;     // distance in rho-t plane           [GeV^-1]
 103     double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over (*aPhiS) [GeV^-1]
 104     double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over tZeta    [GeV^-1]
 105     double tTemp  = mFOHS->TFO;                                 // freeze-out temparature
 106     double ttau0  = mFOHS->tau0/kFmToGev;                       // tau 0
 107     double tMt    = TMath::Hypot(aMass, (*aPt));                // transverse mass
 108     (*aRho)       = tdHS*cos(tZeta);                            // rho
 109     double tdPt   = 1.0/((1-tZet)*(1-tZet));                    // dPt
 110     (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);      // t
 111
 112     double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt*cosh((*aRap)-(*aAlfaP))-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
 113     double tFC    = tdHS*(ttau0+tdHS*sin(tZeta))*(
 114                       tdHS  *cos(tZeta)*( tMt*sin(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*cos(tZeta)*cos((*aPhiS)-(*aPhiP)))+
 115                       tDzdHS*cos(tZeta)*(-tMt*cos(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*sin(tZeta)*cos((*aPhiS)-(*aPhiP)))+
 116                       tDpdHS*(*aPt)*sin((*aPhiS)-(*aPhiP))
 117                     );
 118    if(tFC < 0.0) tFC = 0.0;
 119     if (fabs(floor(aSpin)-aSpin) < 0.01)
 120       tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcBE((tPU-aMiu)/tTemp);
 121     else
 122       tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
 123
 124   }
 125   else if (sModel == 11) { // Lhyquid2D
 126     (*aAlfaP)     = (*aRap);                                    // dirac delta (Y-ap)
 127     double tZeta  = mRandom->Rndm()*0.5*TMath::Pi();            // angle in rho-t plane
 128     double taHS   = mFOHS->fahs  ((*aPhiS),tZeta);              // velocity angle; 0.0=radial flow
 129     double tvHS   = mFOHS->fvhs  ((*aPhiS),tZeta);              // velocity
 130     double tdHS   = mFOHS->fdhs  ((*aPhiS),tZeta)/kFmToGev;     // distance in rho-t plane           [GeV^-1]
 131     double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over (*aPhiS) [GeV^-1]
 132     double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;     // distance derivative over tZeta    [GeV^-1]
 133     double tTemp  = mFOHS->TFO;                                 // freeze-out temparature
 134     double ttau0  = mFOHS->tau0/kFmToGev;                       // tau 0
 135     double tMt    = TMath::Hypot(aMass, (*aPt));                // transverse mass
 136     (*aRho)       = tdHS*cos(tZeta);                            // rho
 137     double tdPt   = 1.0/((1-tZet)*(1-tZet));                    // dPt
 138     (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);      // t
 139
 140     double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
 141     double tFC    = tdHS*(
 142                       tdHS  *cos(tZeta)*( (*aPt)/tMt*cos(tZeta)*cos((*aPhiS)-(*aPhiP))+sin(tZeta) )+
 143                       tDzdHS*cos(tZeta)*( (*aPt)/tMt*sin(tZeta)*cos((*aPhiS)-(*aPhiP))-cos(tZeta) )+
 144                       tDpdHS*             (*aPt)/tMt*           sin((*aPhiS)-(*aPhiP))
 145                     );
 146     if(tFC < 0.0) tFC = 0.0;
 147     if (fabs(floor(aSpin)-aSpin) < 0.01)
 148       tFpod =1.0/kTwoPi2 * tdPt*(*aPt) * tFC * CalcBE((tPU-aMiu)/tTemp);
 149     else
 150       tFpod =1.0/kTwoPi2 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
 151   }
 152 /*MCH end*/
 153   else if (sModel == 2) { // Blast-Wave with Vt
 154     double tMt  = TMath::Hypot(aMass, (*aPt));
 155     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
 156     double tCHay = cosh((*aAlfaP)-(*aRap));
 157     double tCSsp = cos((*aPhiS) - (*aPhiP));
 158     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
 159     double tGamma = 1.0/sqrt(1-mBWVt*mBWVt);
 160     double tPU  = tMt*tGamma*tCHay-(*aPt)*mBWVt*tGamma*tCSsp - aMiu;
 161     if (fabs(floor(aSpin) - aSpin)< 0.01)
 162       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
 163     else
 164       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
 165   }
 166   else if (sModel == 6) { // Blast-Wave with Vt
 167     double tMt  = TMath::Hypot(aMass, (*aPt));
 168     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
 169     double tCHay = cosh((*aAlfaP)-(*aRap));
 170     double tCSsp = cos((*aPhiS) - (*aPhiP));
 171     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
 172     double tGamma = 1.0/sqrt(1-mBWVt*mBWVt);
 173     double tPU  = tMt*tGamma*tCHay-(*aPt)*mBWVt*tGamma*tCSsp - aMiu;
 174     if (fabs(floor(aSpin) - aSpin)< 0.01)
 175       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
 176     else
 177       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
 178   }
 179   else if (sModel == 4) { // Blast-Wave with linear Vt profile
 180     double tMt  = TMath::Hypot(aMass, (*aPt));
 181     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
 182     double tCHay = cosh((*aAlfaP)-(*aRap));
 183     double tCSsp = cos((*aPhiS) - (*aPhiP));
 184     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
 185     double tVt = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
 186     double tGamma = 1.0/sqrt(1-tVt*tVt);
 187     double tPU  = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
 188     if (fabs(floor(aSpin) - aSpin)< 0.01)
 189       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
 190     else
 191       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
 192   }
 193   else if (sModel == 7) {
 194     // Blast-Wave with linear Vt profile
 195     // and delay in particle emission point - formation time
 196     double tMt  = TMath::Hypot(aMass, (*aPt));
 197     double tPre = (mTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
 198     double tCHay = cosh((*aAlfaP)-(*aRap));
 199     double tCSsp = cos((*aPhiS) - (*aPhiP));
 200     double tBra = tMt*tCHay - mBWA*(*aPt)*tCSsp;
 201     double tVt = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
 202     double tGamma = 1.0/sqrt(1-tVt*tVt);
 203     double tPU  = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
 204     if (fabs(floor(aSpin) - aSpin)< 0.01)
 205       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
 206     else
 207       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
 208   }
 209   else if (sModel == 8) {
 210     // Blast-Wave with linear Vt profile
 211     // and delay in emission time - exponential decay
 212     double tTau   = (*aTime);
 213     double tMt    = TMath::Hypot(aMass, (*aPt));
 214     double tPre   = (tTau + mBWA*(*aRho))*(*aPt)*(*aRho)/kTwoPi3;
 215     double tCHay  = cosh((*aAlfaP)-(*aRap));
 216     double tCSsp  = cos((*aPhiS) - (*aPhiP));
 217     double tBra   = tMt*tCHay - mBWA*(*aPt)*tCSsp;
 218     double tVt    = ((*aRho)/mRhoMax)/(mBWVt+((*aRho)/mRhoMax));
 219     double tGamma = 1.0/sqrt(1-tVt*tVt);
 220     double tPU    = tMt*tGamma*tCHay-(*aPt)*tVt*tGamma*tCSsp - aMiu;
 221     if (fabs(floor(aSpin) - aSpin)< 0.01)
 222       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcBE(tPU/mTemp);
 223     else
 224       tFpod = 1/((1-tZet)*(1-tZet))*tPre*tBra*CalcFD(tPU/mTemp);
 225   }
 226 /* -=[ Testing with Mathematica ]=-*/
 227 /*if ((sModel == 10)||(sModel == 11)) {
 228   cout.precision(20);
 229   cout << endl;
 230   cout << "{phi  = " << (*aPhiS)  << ","  << endl;
 231   cout << " zeta = " <<   tZeta   << ","  << endl;
 232   cout << " z    = " <<   tZet    << ","  << endl;
 233   cout << " phip = " << (*aPhiP)  << ","  << endl;
 234   cout << " ap   = " << (*aAlfaP) << ","  << endl;
 235   cout << " Y    = " << (*aRap)   << ","  << endl;
 236   cout << " m    = " <<   aMass   << ","  << endl;
 237   cout << " s    = " <<   aSpin   << ","  << endl;
 238   cout << " mu   = " <<   aMiu    << "};" << endl;
 239   cout << endl;
 240   cout.precision(6);
 241   cout << "TF0  = " <<   tTemp  << endl;
 242   cout << "t    = " << (*aTime) << endl;
 243   cout << "rho  = " << (*aRho)  << endl;
 244   cout << "pT   = " << (*aPt)   << endl;
 245   cout << "dpT  = " <<   tdPt   << endl;
 246   cout << "mT   = " <<   tMt    << endl;
 247   cout << endl;
 248   cout << "aHS     = " << taHS     << endl;
 249   cout << "vHS     = " << tvHS     << endl;
 250   cout << "dHS     = " << tdHS     << endl;
 251   cout << "DpdHS   = " << tDpdHS   << endl;
 252   cout << "DzdHS   = " << tDzdHS   << endl;
 253   cout << "PU      = " << tPU      << endl;
 254   cout << "FC      = " << tFC      << endl;
 255   cout << "part1   = " << (tPU - aMiu)/tTemp << endl;
 256   cout << "part2BE = " << CalcBE( (tPU - aMiu)/tTemp ) << endl;
 257   cout << "part2FD = " << CalcFD( (tPU - aMiu)/tTemp ) << endl;
 258   cout << "Fpod    = " << tFpod    << endl;
 259   getchar();
 260 }*/
 261
 262   return tFpod;
 263 }
 264
 265 double Integrator::GetMiu(double aIzo, double aBar, double aStr)
 266 {
 267   return (mMiu_i*aIzo+mMiu_s*aStr+mMiu_b*aBar);
 268 }
 269
 270 double Integrator::CalcFunPodCalk(double aMass, double aIzo, double aBar, double aStr, double aSpin)
 271 {
 272   //  int proc = mNPart/100;
 273   double tFpod;
 274   double tMax = 0.0;
 275   double tMiu = GetMiu(aIzo, aBar, aStr);
 276   double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime;
 277   double tPt;
 278
 279   for (int tIpart=0; tIpart<mNPart; tIpart++)
 280     {
 281       tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
 282       if (tFpod>tMax) tMax = tFpod;
 283     }
 284   return tMax;
 285 }
 286
 287 double Integrator::Integrate(double aMass, double aIzo, double aBar, double aStr, double aSpin)
 288 {
 289   double tFpod=0.0;
 290   double tFtest=0.0;
 291   double tMiu = GetMiu(aIzo, aBar, aStr);
 292   double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime;
 293   double tPt;
 294
 295   for (int tIpart=0; tIpart<mNPart; tIpart++)
 296     {
 297       tFpod = Calka(aMass,tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,aSpin);
 298
 299       tFtest += tFpod;
 300     }
 301
 302   double tCalka;
 303   if (sModel == 5)
 304     tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*(mTauf-mTau0)* tFtest / (mNPart*1.0);
 305 /*MCH begin*/
 306   else if (sModel == 10)
 307     tCalka = mRapRange*mAlfaRange*(1.0)*2.0*TMath::Pi()*2.0*TMath::Pi()*0.5*TMath::Pi()*tFtest / (mNPart*1.0);
 308   else if (sModel == 11)
 309     tCalka = mRapRange*(1.0)*2.0*TMath::Pi()*2.0*TMath::Pi()*0.5*TMath::Pi()*tFtest / (mNPart*1.0);
 310 /*MCH end*/
 311   else
 312     tCalka = mRapRange*(1.0)*mAlfaRange*4*TMath::Pi()*TMath::Pi()*mRhoMax* tFtest / (mNPart*1.0);
 313   return tCalka;
 314 }
 315
 316 void
 317 Integrator::Generate(ParticleType *aPartType, int aPartCount, Particle ***oParticles)
 318 {
 319   int tIter = 0;
 320   double tFpod;
 321   double tFtest;
 322   double tMiu = GetMiu(1.0*aPartType->GetI3(), 1.0*aPartType->GetBarionN(), 1.0*aPartType->GetStrangeness());
 323   double tFMax;
 324   double tSpin = aPartType->GetSpin();
 325   double tRap, tPhiP, tAlfaP, tRho, tPhiS, tTime;
 326   double tPt;
 327
 328   PRINT_DEBUG_3("Gen for " << (aPartType->GetName()) << " i3 b s " << 1.0*aPartType->GetI3() << " " << 1.0*aPartType->GetBarionN() << " " << 1.0*aPartType->GetStrangeness() << " " << tMiu);
 329
 330   tFMax = aPartType->GetFMax();
 331
 332   (*oParticles) = (Particle **) malloc(sizeof(Particle *) * aPartCount);
 333   Particle *tBuf;
 334
 335   while (tIter<aPartCount)
 336     {
 337       tFpod = Calka(aPartType->GetMass(),tMiu,&tRap,&tPt,&tPhiP,&tAlfaP,&tRho,&tPhiS,&tTime,tSpin);
 338       tFtest = mRandom->Rndm()*tFMax;
 339       if (tFtest<tFpod)
 340         {
 341           if ((sModel == 0) || (sModel == 3)) // Single freeze-out
 342             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, TMath::Hypot(mTau, tRho), aPartType);
 343 /*MCH begin*/
 344           else if ((sModel == 10) || (sModel == 11)) {
 345             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTime/cosh(tAlfaP), aPartType);
 346           }
 347 /*MCH end*/
 348           else if ((sModel == 1) || (sModel == 2) || (sModel == 4)) { // Blast-wave
 349             double tTau = mTau + mBWA * tRho;
 350             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTau, aPartType);
 351           }
 352           else if (sModel == 7) { // Blast-wave
 353             double tTau = mTau + mBWA * tRho;
 354             double px = tPt*TMath::Cos(tPhiP);
 355             double py = tPt*TMath::Sin(tPhiP);
 356             double tMt = TMath::Hypot(aPartType->GetMass(),tPt);
 357             double pz = tMt*TMath::SinH(tRap);
 358
 359             double rx = tRho*TMath::Cos(tPhiS);
 360             double ry = tRho*TMath::Sin(tPhiS);
 361
 362             double rz = tTau*TMath::SinH(tAlfaP);
 363             double rt = tTau*TMath::CosH(tAlfaP);
 364             double dt = -mBWDelay * TMath::Log(mRandom->Rndm());
 365             rt += dt;
 366             double en = sqrt(tMt*tMt + pz*pz);
 367             rx += dt * px/en;
 368             ry += dt * py/en;
 369             rz += dt * pz/en;
 370
 371             tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
 372           }
 373           else if (sModel == 8) { // Blast-wave
 374             double tTau = tTime + mBWA * tRho;
 375             tBuf = new Particle(tRap, tPt, tPhiP, tAlfaP, tRho, tPhiS, tTau, aPartType);
 376           }
 377           else if (sModel == 6) { // Blast-wave
 378             double tTau = mTau + mBWA * tRho;
 379             double px = tPt*TMath::Cos(tPhiP);
 380             double py = tPt*TMath::Sin(tPhiP);
 381             double tMt = TMath::Hypot(aPartType->GetMass(),tPt);
 382             double pz = tMt*TMath::SinH(tRap);
 383
 384             double rx = tRho*TMath::Cos(tPhiS);
 385             double ry = tRho*TMath::Sin(tPhiS);
 386
 387             double rz = tTau*TMath::SinH(tAlfaP);
 388             double rt = tTau*TMath::CosH(tAlfaP);
 389             rt += -mBWDelay * TMath::Log(mRandom->Rndm());
 390
 391             tBuf = new Particle(aPartType, px, py, pz, rx, ry, rz, rt);
 392           }
 393           else if (sModel == 5) {
 394             double tTau = sqrt(tTime*tTime - tAlfaP*tAlfaP);
 395             double tAlfa = 0.5*log((tTime+tAlfaP)/(tTime-tAlfaP));
 396             tBuf = new Particle(tRap, tPt, tPhiP, tAlfa, tRho, tPhiS, tTau, aPartType);
 397           }
 398           (*oParticles)[tIter] = tBuf;
 399           tIter++;
 400         }
 401     }
 402 }
 403
 404 void
 405 Integrator::Randomize()
 406 {
 407   TDatime tDat;
 408
 409   //  mRandom->SetSeed2(tDat.Get(), (tDat.Get() % 11) * 7 + (tDat.Get() / 7));
 410   mRandom->SetSeed(tDat.Get());
 411 }
 412
 413 void
 414 Integrator::ReadParameters()
 415 {
 416   STR tModel;
 417   STR tTable;
 418
 419   // First read the model parameters
 420   try {
 421
 422     mRhoMax = atof(sRPInstance->getPar("RhoMax").Data()) / kFmToGev;
 423     mTau    = atof(sRPInstance->getPar("Tau").Data()) / kFmToGev;
 424     mTemp   = atof(sRPInstance->getPar("Temperature").Data());
 425     mMiu_i  = atof(sRPInstance->getPar("MiuI").Data());
 426     mMiu_s  = atof(sRPInstance->getPar("MiuS").Data());
 427     mMiu_b  = atof(sRPInstance->getPar("MiuB").Data());
 428   }
 429   catch (STR tError) {
 430     PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
 431     PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
 432     exit(0);
 433   }
 434
 435   // Read additional parameters for BW
 436   if ((sModel == 1) || (sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
 437     try {
 438       mBWA  = atof(sRPInstance->getPar("BWA").Data());
 439     }
 440     catch (STR tError) {
 441       // Using default value of 0
 442       mBWA = 0.0;
 443     }
 444   }
 445
 446   // Read additional parameters for MBWVt
 447   if ((sModel == 2) || (sModel == 4) || (sModel == 7) || (sModel == 8)) {
 448     try {
 449       mBWVt  = atof(sRPInstance->getPar("BWVt").Data());
 450     }
 451     catch (STR tError) {
 452       PRINT_DEBUG_1("Integrator::ReadParameters (BW Vt part) - Caught exception " << tError);
 453       PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
 454       exit(0);
 455     }
 456   }
 457
 458   // Read additional parameters for MBWVt
 459   if ((sModel ==6) || (sModel == 7) || (sModel == 8))
 460     try {
 461       mBWDelay  = atof(sRPInstance->getPar("BWDelay").Data()) / kFmToGev;
 462     }
 463     catch (STR tError) {
 464       PRINT_DEBUG_1("Integrator::ReadParameters (BW Vt Delay part) - Caught exception " << tError);
 465       PRINT_MESSAGE("Did not find one of the neccessary model parameters.");
 466       exit(0);
 467     }
 468
 469
 470   // Then the integration range
 471   try {
 472     mAlfaRange = atof(sRPInstance->getPar("AlphaRange").Data());
 473     mRapRange  = atof(sRPInstance->getPar("RapidityRange").Data());
 474   }
 475   catch (STR tError) {
 476     PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
 477     PRINT_MESSAGE("Did not find one of the neccessary integration ranges.");
 478     exit(0);
 479   }
 480
 481   // Read hydro-specific parameters
 482   if (sModel == 5) {
 483     try {
 484       mTauf   = atof(sRPInstance->getPar("TauF").Data());
 485       mTau0   = atof(sRPInstance->getPar("Tau0").Data());
 486       mLambda = atof(sRPInstance->getPar("Lambda").Data());
 487       mBN     = atof(sRPInstance->getPar("BN").Data());
 488       mAlfa   = atof(sRPInstance->getPar("Alfa").Data());
 489     }
 490     catch (STR tError) {
 491       PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
 492       PRINT_MESSAGE("Did not find one of the neccessary hydro parameters.");
 493       exit(0);
 494     }
 495   }
 496 }
 497
 498 char *
 499 Integrator::ParameterHash()
 500 {
 501   char *tBuf;
 502
 503   tBuf = (char *) malloc(sizeof(char *) * (15 * 3 + 3));
 504
 505   if (sModel == 0) {
 506     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);
 507   }
 508 /*MCH begin*/
 509   else if ((sModel == 10) || (sModel == 11)) {
 510     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);
 511   }
 512 /*MCH end*/
 513   else if ((sModel == 2) || (sModel == 4))  {
 514     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);
 515   }
 516   else if ((sModel == 6) || (sModel == 7) || (sModel == 8)) {
 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%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);
 518   }
 519   else if (sModel == 5)  {
 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%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);
 521   }
 522
 523   return tBuf;
 524 }
 525
 526 void
 527 Integrator::ReadMultInteg(ParticleDB *aDB)
 528 {
 529   // Make or read table with propabilities
 530   int tK;
 531   char *tHash;
 532   char  tIntName[100];
 533   char  tMultName[100];
 534   ifstream *tFileIn = NULL;
 535   ofstream *tFileOut = NULL;
 536
 537   tHash = ParameterHash();
 538
 539   strcpy(tIntName, "fintegrandmax_");
 540   strcat(tIntName, tHash);
 541   strcat(tIntName, ".txt");
 542
 543   tFileIn = new ifstream(tIntName);
 544   if ((tFileIn) && (tFileIn->is_open())) {
 545     PRINT_MESSAGE("Reading Max Integrand values from " << tIntName);
 546
 547     char tPart[255];
 548     float tRFunMax;
 549     std::string tStrPart;
 550
 551     while (!tFileIn->eof()) {
 552       (*tFileIn) >> tPart >> tRFunMax;
 553
 554       tStrPart = tPart;
 555       aDB->GetParticleType(tStrPart)->SetFMax(tRFunMax);
 556     }
 557   }
 558   else {
 559     float tRFunMax;
 560
 561     PRINT_MESSAGE("Max Integrand file " << tIntName << " not found. Generating...");
 562
 563     tFileOut = new ofstream(tIntName);
 564
 565     for(tK=0; tK<aDB->GetParticleTypeCount(); tK++)
 566       {
 567         tRFunMax = CalcFunPodCalk(aDB->GetParticleType(tK)->GetMass(),
 568                                aDB->GetParticleType(tK)->GetI3(),
 569                                aDB->GetParticleType(tK)->GetBarionN()*1.0,
 570                                aDB->GetParticleType(tK)->GetStrangeness()*1.0,
 571                                aDB->GetParticleType(tK)->GetSpin());
 572
 573         (*tFileOut) << aDB->GetParticleType(tK)->GetName() << " "
 574                 << tRFunMax <<endl;
 575         aDB->GetParticleType(tK)->SetFMax(tRFunMax);
 576         PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tK)->GetNumber()<<":"
 577                       <<aDB->GetParticleType(tK)->GetName()
 578                       <<" done");
 579       }
 580     tFileOut->close();
 581   }
 582
 583   // Calculate or read multiplicities
 584   strcpy(tMultName, "fmultiplicity_");
 585   strcat(tMultName, tHash);
 586   strcat(tMultName, ".txt");
 587
 588   tFileIn = new ifstream(tMultName);
 589   if ((tFileIn) && (tFileIn->is_open())) {
 590     PRINT_MESSAGE("Reading Multiplicities from " << tMultName);
 591   }
 592   else {
 593     PRINT_MESSAGE("Multiplicities file " << tMultName << " not found. Generating...");
 594
 595     tFileOut = new ofstream(tMultName);
 596
 597     for (int tPart=0; tPart<aDB->GetParticleTypeCount(); tPart++) {
 598       double tMult = Integrate(aDB->GetParticleType(tPart)->GetMass(),
 599                                aDB->GetParticleType(tPart)->GetI3(),
 600                                aDB->GetParticleType(tPart)->GetBarionN()*1.0,
 601                                aDB->GetParticleType(tPart)->GetStrangeness()*1.0,
 602                                aDB->GetParticleType(tPart)->GetSpin());
 603       (*tFileOut) << (aDB->GetParticleType(tPart)->GetName()) << " " << tMult*TMath::Abs(aDB->GetParticleType(tPart)->GetSpin()*2+1) << endl;
 604
 605       PRINT_DEBUG_1("particle "<<aDB->GetParticleType(tPart)->GetNumber()<<":"
 606                     <<aDB->GetParticleType(tPart)->GetName()
 607                     <<" done");
 608     }
 609     tFileOut->close();
 610   }
 611
 612   free (tHash);
 613 }
 614