Add reaction plane to the header. Update the Lhyquid hypersurface parameterization...

diff --git a/TTherminator/AliGenTherminator.cxx b/TTherminator/AliGenTherminator.cxx
index dffe99c..0aedb11 100644 (file)
--- a/TTherminator/AliGenTherminator.cxx
+++ b/TTherminator/AliGenTherminator.cxx
@@ -14,6 +14,7 @@
 #include "AliConst.h"
 #include "AliDecayer.h"
 #include "AliGenEventHeader.h"
+#include "AliGenHijingEventHeader.h"
 #include "AliGenTherminator.h"
 #include "AliLog.h"
 #include "AliRun.h"
@@ -100,6 +101,7 @@ void AliGenTherminator::Generate()
   AliWarning(Form("Imported %d particles", np));
 
   TParticle *iparticle;
+  Double_t evrot = gRandom->Rndm()*TMath::Pi();
   
   for (int i = 0; i < np; i++) {
     iparticle = (TParticle *) fParticles->At(i);
@@ -108,8 +110,12 @@ void AliGenTherminator::Generate()
     
     kf   = iparticle->GetPdgCode();
     ks   = iparticle->GetStatusCode();
-    p[0] = iparticle->Px();
-    p[1] = iparticle->Py();
+    Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px());
+    Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py());
+    p[0] = arho*TMath::Cos(aphi + evrot);
+    p[1] = arho*TMath::Sin(aphi + evrot);
+//     p[0] = iparticle->Px();
+//     p[1] = iparticle->Py();
     p[2] = iparticle->Pz();
     mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass();
     energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
@@ -144,11 +150,56 @@ void AliGenTherminator::Generate()
   eventVertex[1] = origin0[1];
   eventVertex[2] = origin0[2];
 
+// Builds the event header, to be called after each event
+  AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator");
+
   // Header
-  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
+  //  AliGenEventHeader* header = new AliGenEventHeader("Therminator");
   // Event Vertex
-  header->SetPrimaryVertex(eventVertex);
-  header->SetNProduced(fNprimaries);
+//   header->SetPrimaryVertex(eventVertex);
+//   header->SetNProduced(fNprimaries);
+
+  ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries);
+  ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex);
+  ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0);
+  ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0);
+  ((AliGenHijingEventHeader*) header)->SetHardScatters(0);
+  ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0);
+  ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0);
+  ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0);
+  ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot);
+
+
+// 4-momentum vectors of the triggered jets.
+//
+// Before final state gluon radiation.
+//     TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), 
+//                                           fHijing->GetHINT1(22),
+//                                           fHijing->GetHINT1(23),
+//                                           fHijing->GetHINT1(24));
+
+//     TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), 
+//                                           fHijing->GetHINT1(32),
+//                                           fHijing->GetHINT1(33),
+//                                           fHijing->GetHINT1(34));
+// // After final state gluon radiation.
+//     TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), 
+//                                           fHijing->GetHINT1(27),
+//                                           fHijing->GetHINT1(28),
+//                                           fHijing->GetHINT1(29));
+
+//     TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), 
+//                                           fHijing->GetHINT1(37),
+//                                           fHijing->GetHINT1(38),
+//                                           fHijing->GetHINT1(39));
+//     ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
+// Bookkeeping for kinematic bias
+//     ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
+// Event Vertex
+  header->SetPrimaryVertex(fVertex);
+  AddHeader(header);
+  fCollisionGeometry = (AliGenHijingEventHeader*)  header;
+
   gAlice->SetGenEventHeader(header); 
 }
 

diff --git a/TTherminator/AliGenTherminator.h b/TTherminator/AliGenTherminator.h
index 48d025d..037f12d 100644 (file)
--- a/TTherminator/AliGenTherminator.h
+++ b/TTherminator/AliGenTherminator.h
@@ -43,9 +43,9 @@ class AliGenTherminator : public AliGenMC
   void SetBWDelay(Double_t bwd); 
 
   void SetModel(const char *model);
+  void ReadShareParticleTable();
 
  protected:
-  void     ReadShareParticleTable();
   void     CreateTherminatorInputFile();
 
   Int_t    fNt;                // CurrentTrack;

diff --git a/TTherminator/Therminator/Hypersurface.cxx b/TTherminator/Therminator/Hypersurface.cxx
index 056df9c..9b3b4eb 100644 (file)
--- a/TTherminator/Therminator/Hypersurface.cxx
+++ b/TTherminator/Therminator/Hypersurface.cxx
@@ -24,6 +24,9 @@ Hypersurface::Hypersurface(void) {
 // freeze-out temperature
     HSFile->getline(buff,100);
     TFO=atof(buff);
+// hydro initial time
+    HSFile->getline(buff,100);
+    tau0=atof(buff);
 
     HSFile->close();
     dp = (fp-ip)/(Np-1);

diff --git a/TTherminator/Therminator/Hypersurface.h b/TTherminator/Therminator/Hypersurface.h
index 5981f6b..0dba0a2 100644 (file)
--- a/TTherminator/Therminator/Hypersurface.h
+++ b/TTherminator/Therminator/Hypersurface.h
@@ -38,6 +38,7 @@ class Hypersurface {
     ~Hypersurface(void);
     // variables
     double TFO;
+    double tau0;
     // functions
     double fahs  (double p, double z);
     double fvhs  (double p, double z);

diff --git a/TTherminator/Therminator/Integrator.cxx b/TTherminator/Therminator/Integrator.cxx
index 50edc25..cdd7c30 100644 (file)
--- a/TTherminator/Therminator/Integrator.cxx
+++ b/TTherminator/Therminator/Integrator.cxx
@@ -103,22 +103,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 +131,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*(