Added AliRsnAction, AliRsnPIDRange

[u/mrichter/AliRoot.git] / TTherminator / Therminator / Integrator.cxx

diff --git a/TTherminator/Therminator/Integrator.cxx b/TTherminator/Therminator/Integrator.cxx
index 50edc25c90d05ab4463ceee22b89394babc9a45a..3b412ab0346d849075695cdc13e937793446aa90 100644 (file)
--- a/TTherminator/Therminator/Integrator.cxx
+++ b/TTherminator/Therminator/Integrator.cxx
@@ -40,7 +40,10 @@ Integrator::Integrator(int aNpart)
   kFmToGev  = 0.197326960277;                          /*MCH updated: kFmToGev  = 0.197;*/
   ReadParameters();
 
-  PRINT_MESSAGE("Hash for these parameters is: " << ParameterHash());
+  char *tHash;
+  tHash = ParameterHash();
+
+  PRINT_MESSAGE("Hash for these parameters is: " << tHash);
   
   mNPart = aNpart;
   kTwoPi2 = TMath::Pi()*TMath::Pi()*2*2;               /*MCH*/
@@ -49,16 +52,20 @@ Integrator::Integrator(int aNpart)
   //  mRandom->SetSeed2(41321, 8457);
   mRandom->SetSeed(41321);
 
-  mFOHS = new Hypersurface();                          /*MCH*/
+  mFOHS = new Hypersurface(mFOHSlocation.Data());                              /*MCH*/
+
+  free (tHash);
 }
 
 double Integrator::CalcBE(double aX)
 {
+  if (aX>200.0) return 0.0;
   return 1.0/(TMath::Exp(aX)-1);
 }
 
 double Integrator::CalcFD(double aX)
 {
+  if (aX>200.0) return 0.0;
   return 1.0/(TMath::Exp(aX)+1);
 }
 
@@ -103,22 +110,24 @@ double Integrator::Calka(double aMass, double aMiu,
     double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over (*aPhiS) [GeV^-1]
     double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over tZeta    [GeV^-1]
     double tTemp  = mFOHS->TFO;                                        // freeze-out temparature
+    double ttau0  = mFOHS->tau0/kFmToGev;                      // tau 0
     double tMt    = TMath::Hypot(aMass, (*aPt));               // transverse mass
     (*aRho)       = tdHS*cos(tZeta);                           // rho
     double tdPt   = 1.0/((1-tZet)*(1-tZet));                   // dPt
-    (*aTime)      = tdHS*sin(tZeta)*cosh(*aAlfaP);             // t
+    (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);     // t
 
     double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt*cosh((*aRap)-(*aAlfaP))-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
-    double tFC    = tdHS*tdHS*sin(tZeta)*(
+    double tFC    = tdHS*(ttau0+tdHS*sin(tZeta))*(
                      tdHS  *cos(tZeta)*( tMt*sin(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*cos(tZeta)*cos((*aPhiS)-(*aPhiP)))+
                      tDzdHS*cos(tZeta)*(-tMt*cos(tZeta)*cosh((*aRap)-(*aAlfaP))+(*aPt)*sin(tZeta)*cos((*aPhiS)-(*aPhiP)))+
                      tDpdHS*(*aPt)*sin((*aPhiS)-(*aPhiP))
                    );
-    if(tFC < 0.0) tFC = 0.0;
+   if(tFC < 0.0) tFC = 0.0;
     if (fabs(floor(aSpin)-aSpin) < 0.01)
       tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcBE((tPU-aMiu)/tTemp);
     else
       tFpod =1.0/kTwoPi3 * tdPt*(*aPt) * tFC * CalcFD((tPU-aMiu)/tTemp);
+
   }
   else if (sModel == 11) { // Lhyquid2D
     (*aAlfaP)     = (*aRap);                                   // dirac delta (Y-ap)
@@ -129,10 +138,11 @@ double Integrator::Calka(double aMass, double aMiu,
     double tDpdHS = mFOHS->fDpdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over (*aPhiS) [GeV^-1]
     double tDzdHS = mFOHS->fDzdhs((*aPhiS),tZeta)/kFmToGev;    // distance derivative over tZeta    [GeV^-1]
     double tTemp  = mFOHS->TFO;                                        // freeze-out temparature
+    double ttau0  = mFOHS->tau0/kFmToGev;                      // tau 0
     double tMt    = TMath::Hypot(aMass, (*aPt));               // transverse mass
     (*aRho)       = tdHS*cos(tZeta);                           // rho
     double tdPt   = 1.0/((1-tZet)*(1-tZet));                   // dPt
-    (*aTime)      = tdHS*sin(tZeta)*cosh(*aAlfaP);             // t
+    (*aTime)      = (ttau0+tdHS*sin(tZeta))*cosh(*aAlfaP);     // t
 
     double tPU    = 1.0/sqrt(1.0-tvHS*tvHS)*(tMt-(*aPt)*tvHS*cos((*aPhiS)-(*aPhiP)+taHS));
     double tFC    = tdHS*(
@@ -327,7 +337,7 @@ Integrator::Generate(ParticleType *aPartType, int aPartCount, Particle ***oParti
   tFMax = aPartType->GetFMax();
 
   (*oParticles) = (Particle **) malloc(sizeof(Particle *) * aPartCount);
-  Particle *tBuf;
+  Particle *tBuf=0;
   
   while (tIter<aPartCount)
     {
@@ -422,6 +432,7 @@ Integrator::ReadParameters()
     mMiu_i  = atof(sRPInstance->getPar("MiuI").Data());
     mMiu_s  = atof(sRPInstance->getPar("MiuS").Data());
     mMiu_b  = atof(sRPInstance->getPar("MiuB").Data());
+    mFOHSlocation = sRPInstance->getPar("FOHSLocation");
   }
   catch (STR tError) {
     PRINT_DEBUG_1("Integrator::ReadParameters - Caught exception " << tError);
@@ -526,18 +537,28 @@ Integrator::ReadMultInteg(ParticleDB *aDB)
   // Make or read table with propabilities
   int tK;
   char *tHash;
-  char  tIntName[100];
-  char  tMultName[100];
+  char  tIntName[1000];
+  char  tMultName[1000];
   ifstream *tFileIn = NULL;
   ofstream *tFileOut = NULL;
   
   tHash = ParameterHash();
   
-  strcpy(tIntName, "fintegrandmax_");
-  strcat(tIntName, tHash);
-  strcat(tIntName, ".txt");
+  if (mFOHSlocation != "") {
+    strcpy(tIntName, mFOHSlocation.Data());
+    strcat(tIntName, "/");
+    strcat(tIntName, "fintegrandmax_");
+    strcat(tIntName, tHash);
+    strcat(tIntName, ".txt");
+    tFileIn = new ifstream(tIntName);
+  }
+  else if (!((tFileIn) && (tFileIn->is_open()))) {
+    strcpy(tIntName, "fintegrandmax_");
+    strcat(tIntName, tHash);
+    strcat(tIntName, ".txt");
+    tFileIn = new ifstream(tIntName);
+  }
   
-  tFileIn = new ifstream(tIntName);
   if ((tFileIn) && (tFileIn->is_open())) {
     PRINT_MESSAGE("Reading Max Integrand values from " << tIntName);
 
@@ -578,9 +599,20 @@ Integrator::ReadMultInteg(ParticleDB *aDB)
   }
   
   // Calculate or read multiplicities
-  strcpy(tMultName, "fmultiplicity_");
-  strcat(tMultName, tHash);
-  strcat(tMultName, ".txt");
+  if (mFOHSlocation != "") {
+    strcpy(tMultName, mFOHSlocation.Data());
+    strcat(tMultName, "/");
+    strcat(tMultName, "fmultiplicity_");
+    strcat(tMultName, tHash);
+    strcat(tMultName, ".txt");
+    tFileIn = new ifstream(tMultName);
+  }
+  else if (!((tFileIn) && (tFileIn->is_open()))) {
+    strcpy(tMultName, "fmultiplicity_");
+    strcat(tMultName, tHash);
+    strcat(tMultName, ".txt");
+    tFileIn = new ifstream(tMultName);
+  }
 
   tFileIn = new ifstream(tMultName);
   if ((tFileIn) && (tFileIn->is_open())) {