Optionally the log term of Rossi parameterization for mult.scattering can be
[u/mrichter/AliRoot.git] / TTherminator / AliGenTherminator.cxx
index 149f0b4..4f42251 100644 (file)
@@ -31,6 +31,7 @@ AliGenTherminator::AliGenTherminator():
   fEventNumber(0),
   fFileName(""),
   fFreezeOutModel(""),
+  fFOHSlocation(""),
   fTemperature(0.1656),
   fMiuI(-0.0009),
   fMiuS(0.0),
@@ -44,7 +45,15 @@ AliGenTherminator::AliGenTherminator():
   fBWDelay(0.0)
 {
   // Default constructor
-  fParticles = new TClonesArray("TParticle",20000);    
+  fFOHSlocation = "";
+
+  fEnergyCMS = 5500.;
+  fAProjectile = 208;
+  fZProjectile = 82;
+  fProjectile = "A";
+  fATarget = 208;
+  fZTarget = 82;
+  fTarget = "A";
 }
 AliGenTherminator::AliGenTherminator(Int_t npart):
   AliGenMC(npart),
@@ -52,6 +61,7 @@ AliGenTherminator::AliGenTherminator(Int_t npart):
   fEventNumber(0),
   fFileName(""),
   fFreezeOutModel(""),
+  fFOHSlocation(""),
   fTemperature(0.1656),
   fMiuI(-0.0009),
   fMiuS(0.0),
@@ -65,8 +75,14 @@ AliGenTherminator::AliGenTherminator(Int_t npart):
   fBWDelay(0.0)
 {
   // Constructor specifying the size of the particle table
-  fParticles = new TClonesArray("TParticle",20000);    
   fNprimaries = 0;
+  fEnergyCMS = 5500.;
+  fAProjectile = 208;
+  fZProjectile = 82;
+  fProjectile = "A";
+  fATarget = 208;
+  fZTarget = 82;
+  fTarget = "A";
 }
 
 AliGenTherminator::~AliGenTherminator()
@@ -97,56 +113,131 @@ void AliGenTherminator::Generate()
 
   ((TTherminator *) fMCEvGen)->GenerateEvent();
   AliWarning("Generated");
-  ((TTherminator *) fMCEvGen)->ImportParticles(fParticles);
+  ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles);
 
-  Int_t np = fParticles->GetEntriesFast();
+  Int_t np = fParticles.GetEntriesFast();
   AliWarning(Form("Imported %d particles", np));
 
+  Int_t *idsOnStack;
+  idsOnStack = new Int_t[np];
+
   TParticle *iparticle;
   Double_t evrot = gRandom->Rndm()*TMath::Pi();
   
   for (int i = 0; i < np; i++) {
-    iparticle = (TParticle *) fParticles->At(i);
+    iparticle = (TParticle *) fParticles.At(i);
     Bool_t  hasMother   = (iparticle->GetFirstMother()     >=0);
     Bool_t  hasDaughter = (iparticle->GetFirstDaughter()   >=0);
     
-    kf   = iparticle->GetPdgCode();
-    ks   = iparticle->GetStatusCode();
-    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]);
-    
-    Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
-    Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
-    origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
-    origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
-    origin[2] = origin0[2]+iparticle->Vz();
-             
-    imo = -1;
-    TParticle* mother = 0;
-    if (hasMother) {
-      imo = iparticle->GetFirstMother();
-      mother = (TParticle *) fParticles->At(imo);
-    } // if has mother   
-    Bool_t tFlag = (!hasDaughter);
-
-    printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
-    PushTrack(tFlag,imo,kf,
-             p[0],p[1],p[2],energy,
-             origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000,
-             polar[0],polar[1],polar[2],
-             hasMother ? kPDecay:kPNoProcess,nt);
-    
-    fNprimaries++;
-    KeepTrack(nt);
+    if (hasDaughter) {
+      // This particle has decayed
+      // It will not be tracked
+      // Add it only once with coorduinates not
+      // smeared with primary vertex position
+
+      kf   = iparticle->GetPdgCode();
+      ks   = iparticle->GetStatusCode();
+      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]);
+      
+      Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
+      Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
+      origin[0] = vrho*TMath::Cos(vphi + evrot);
+      origin[1] = vrho*TMath::Sin(vphi + evrot);
+      origin[2] = iparticle->Vz();
+      
+      imo = -1;
+      //      TParticle* mother = 0;
+      if (hasMother) {
+       imo = iparticle->GetFirstMother();
+       //      mother = (TParticle *) fParticles.At(imo);
+      } // if has mother   
+      Bool_t tFlag = (!hasDaughter);
+      
+      //      printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
+      PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
+               p[0],p[1],p[2],energy,
+               origin[0],origin[1],origin[2],iparticle->T(),
+               polar[0],polar[1],polar[2],
+               hasMother ? kPDecay:kPNoProcess,nt);
+      idsOnStack[i] = nt;
+      fNprimaries++;
+      KeepTrack(nt);
+    }
+    else {
+      // This is a final state particle
+      // It will be tracked
+      // Add it TWICE to the stack !!!
+      // First time with event-wide coordicates (for femtoscopy) - 
+      //   this one will not be tracked
+      // Second time with event-wide ccordiantes and vertex smearing
+      //   this one will be tracked
+      
+      kf   = iparticle->GetPdgCode();
+      ks   = iparticle->GetStatusCode();
+      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]);
+      
+      Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx());
+      Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy());
+      origin[0] = vrho*TMath::Cos(vphi + evrot);
+      origin[1] = vrho*TMath::Sin(vphi + evrot);
+      origin[2] = iparticle->Vz();
+      
+      imo = -1;
+      //      TParticle* mother = 0;
+      if (hasMother) {
+       imo = iparticle->GetFirstMother();
+       //      mother = (TParticle *) fParticles.At(imo);
+      } // if has mother   
+      Bool_t tFlag = (hasDaughter);
+      
+//       printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo);
+//       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo);
+      PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf,
+               p[0],p[1],p[2],energy,
+               origin[0],origin[1],origin[2],iparticle->T(),
+               polar[0],polar[1],polar[2],
+               hasMother ? kPDecay:kPNoProcess,nt);
+      idsOnStack[i] = nt;
+      fNprimaries++;
+      KeepTrack(nt);
+
+      origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot);
+      origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot);
+      origin[2] = origin0[2]+iparticle->Vz();
+
+      imo = nt;
+      //      mother = (TParticle *) fParticles.At(nt);
+      tFlag = (!hasDaughter);
+      
+//       printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo);
+      PushTrack(tFlag,imo,kf,
+               p[0],p[1],p[2],energy,
+               origin[0],origin[1],origin[2],iparticle->T(),
+               polar[0],polar[1],polar[2],
+               hasMother ? kPDecay:kPNoProcess,nt);
+      fNprimaries++;
+      KeepTrack(nt);
+    }
   }
 
+
+
   SetHighWaterMark(fNprimaries);
 
   TArrayF eventVertex;
@@ -205,7 +296,9 @@ void AliGenTherminator::Generate()
   AddHeader(header);
   fCollisionGeometry = (AliGenHijingEventHeader*)  header;
 
-  gAlice->SetGenEventHeader(header); 
+  delete [] idsOnStack;
+
+  //  gAlice->SetGenEventHeader(header); 
 }
 
 void AliGenTherminator::Init()
@@ -319,6 +412,8 @@ void AliGenTherminator::SetModel(const char *model)
 {
   // Set the freeze-out model to use
   fFreezeOutModel = model;
+  AliWarning(Form("Selected model %s", fFreezeOutModel.Data()));
+  AliWarning(Form("FOHSLocation is %s", fFOHSlocation.Data()));
 }
 
 void AliGenTherminator::SetLhyquidSet(const char *set)
@@ -327,27 +422,112 @@ void AliGenTherminator::SetLhyquidSet(const char *set)
   const char *aroot = gSystem->Getenv("ALICE_ROOT");
   if (strstr(set, "LHC500C0005")) {
     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
-    AliWarning(Form("  Pb-Pb collisions, centrality 0-5%"));
+    AliWarning(Form("  Pb-Pb collisions, centrality 0-5 percent"));
     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt"));
-    fFOHSlocation = aroot;
-    fFOHSlocation += "/TTherminator/data/LHC500C0005";
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC500C0005");
+  }
+  if (strstr(set, "LHC500C0510")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, centrality 5-10 percent"));
+    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0510/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC500C0510");
+  }
+  if (strstr(set, "LHC500C1020")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, centrality 10-20 percent"));
+    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C1020/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC500C1020");
   }
   else if (strstr(set, "LHC500C2030")) {
     AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
-    AliWarning(Form("  Pb-Pb collisions, centrality 20-30%"));
+    AliWarning(Form("  Pb-Pb collisions, centrality 20-30 percent"));
     AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
     AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
     AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt"));
-    fFOHSlocation = aroot;
-    fFOHSlocation += "/TTherminator/data/LHC500C2030";
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC500C2030");
+  }
+  else if (strstr(set, "LHC500C3040")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, centrality 30-40 percent"));
+    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C3040/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC500C3040");
+  }
+  else if (strstr(set, "LHC500C4050")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, centrality 40-50 percent"));
+    AliWarning(Form("  initial temperature at tau=1 fm in the center Ti=500 MeV"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC500C4050/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC500C4050");
+  }
+  else if (strstr(set, "LHC276TC0005")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 00-05 percent"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0005/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC276TC0005");
+  }
+  else if (strstr(set, "LHC276TC0510")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 05-10 percent"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0510/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC276TC0510");
+  }
+  else if (strstr(set, "LHC276TC1020")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 10-20 percent"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC1020/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC276TC1020");
+  }
+  else if (strstr(set, "LHC276TC2030")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 20-30 percent"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC2030/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC276TC2030");
+  }
+  else if (strstr(set, "LHC276TC3040")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 30-40 percent"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC3040/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC276TC3040");
+  }
+  else if (strstr(set, "LHC276TC4050")) {
+    AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface"));
+    AliWarning(Form("  Pb-Pb collisions, 2.76TeV, centrality 40-50 percent"));
+    AliWarning(Form("  freeze-out criteria Tf=145 MeV"));
+    AliWarning(Form("  for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC4050/FO.txt"));
+    fFOHSlocation.Append(aroot);
+    fFOHSlocation.Append("/TTherminator/data/LHC276TC4050");
   }
   else {
     AliWarning(Form("Did not find Lhyquid set %s", set));
     AliWarning(Form("Reverting to default: current directory"));
-    fFOHSlocation = "";
+    fFOHSlocation += "";
   }
+
 }
 
 void AliGenTherminator::SetLhyquidInputDir(const char *inputdir)