SetSeed implementation
[u/mrichter/AliRoot.git] / EVGEN / AliGenCorrHF.cxx
index 8e9c11b..4e5ba35 100644 (file)
@@ -21,9 +21,9 @@
 // Is a generalisation of AliGenParam class for correlated pairs of hadrons.
 // In this version quark pairs and fragmentation functions are obtained from
 // ~2.10^6 Pythia6.214 events generated with kCharmppMNRwmi & kBeautyppMNRwmi, 
-// CTEQ5L PDF and Pt_hard = 2.76 GeV/c for p-p collisions at 7, 10 and 14 TeV,
+// CTEQ5L PDF and Pt_hard = 2.76 GeV/c for p-p collisions at 2.76, 7, 8, 10 and 14 TeV,
 // and with kCharmppMNR (Pt_hard = 2.10 GeV/c) & kBeautyppMNR (Pt_hard = 2.75 GeV/c), 
-// CTEQ4L PDF for Pb-Pb at 3.94 TeV, for p-Pb & Pb-p at 8.8 TeV. 
+// CTEQ4L PDF for Pb-Pb at 2.76 and 3.94 TeV, for p-Pb & Pb-p at 5 and 8.8 TeV. 
 // Decays are performed by Pythia.
 // Author: S. Grigoryan, LPC Clermont-Fd & YerPhI, Smbat.Grigoryan@cern.ch
 // July 07: added quarks in the stack (B. Vulpescu)
 // Oct 09: added energy choice between 7, 10, 14 TeV (for p-p), 4 TeV (for Pb-Pb),
 // 9 TeV (for p-Pb) and -9 TeV (for Pb-p) (S. Grigoryan)
 // April 10: removed "static" from definition of some variables (B. Vulpescu)
-// May 11: Added Flag for transportation Background While using SetForceDecay() function (L. Manceau)
-// June 11 modifications allowing the setting of cuts on Children added 
-//(L. Manceau)
-// Quarks, hadrons and decayed particles are loaded in the stack outside the 
-// main loop on heavy hadrons. 
-// The particles are loaded only when a pair containing
-// two heavy hadrons given children wich statify cut conditions 
-// are tagged in the main loop
+// May 11: added Flag for transportation of background particles while using 
+// SetForceDecay() function (L. Manceau)
+// June 11: added modifications allowing the setting of cuts on HF-hadron children.
+// Quarks, hadrons and decay particles are loaded in the stack outside the loop
+// of HF-hadrons, when the cuts on their children are satisfied (L. Manceau)
+// Oct 11: added Pb-Pb at 2.76 TeV (S. Grigoryan)
+// June 12: added p-Pb & Pb-p at 5 TeV (S. Grigoryan)
+// April 13: added p-p at 2.76 and 8 TeV (S. Grigoryan)
 // 
 //-------------------------------------------------------------------------
 // How it works (for the given flavor and p-p energy):
@@ -74,7 +74,6 @@
 //
 // add the following typical lines in Config.C
 /*
-  if (!strcmp(option,"corr")) {
     // An example for correlated charm or beauty hadron pair production at 14 TeV
 
     // AliGenCorrHF *gener = new AliGenCorrHF(1, 4, 14);  // for charm, 1 pair per event
     gener->SetOrigin(0,0,0);          //vertex position    
     gener->SetSigma(0,0,0);           //Sigma in (X,Y,Z) (cm) on IP position
     gener->SetForceDecay(kSemiMuonic);
-    gener->SetTrackingFlag(0);
+    gener->SetSelectAll(kTRUE);      //Force the transport of all particles. 
+                                     //Necessary while using a different 
+                                     //option than kAll for SetForceDecay
+    gener->SetTrackingFlag(1);       //1: Decay during transport, 
+                                     //0: No Decay during transport
     gener->Init();
-}
 */
-// and in aliroot do e.g. gAlice->Run(10,"Config.C") to produce 10 events.
 // One can include AliGenCorrHF in an AliGenCocktail generator.
 //--------------------------------------------------------------------------
 
@@ -166,12 +167,20 @@ AliGenCorrHF::AliGenCorrHF(Int_t npart, Int_t idquark, Int_t energy):
     if (fQuark == 5) {
       if (fEnergy == 7)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPP7PythiaMNRwmi.root";
+      else if (fEnergy == 8)
+           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPP8PythiaMNRwmi.root";
       else if (fEnergy == 10)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPP10PythiaMNRwmi.root";
       else if (fEnergy == 14)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPP14PythiaMNRwmi.root";
+      else if (fEnergy == 2)
+           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPP276PythiaMNRwmi.root";
+      else if (fEnergy == 3)
+          fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPbPb276PythiaMNR.root";
       else if (fEnergy == 4)
           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPbPb394PythiaMNR.root";
+      else if (fEnergy == 5 || fEnergy == -5)
+          fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPPb5PythiaMNR.root";
       else if (fEnergy == 9 || fEnergy == -9)
           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPPb88PythiaMNR.root";
       else fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/BeautyPbPb394PythiaMNR.root";
@@ -180,12 +189,20 @@ AliGenCorrHF::AliGenCorrHF(Int_t npart, Int_t idquark, Int_t energy):
       fQuark = 4;
       if (fEnergy == 7)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPP7PythiaMNRwmi.root";
+      else if (fEnergy == 8)
+           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPP8PythiaMNRwmi.root";
       else if (fEnergy == 10)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPP10PythiaMNRwmi.root";
       else if (fEnergy == 14)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPP14PythiaMNRwmi.root";
+      else if (fEnergy == 2)
+           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPP276PythiaMNRwmi.root";
+      else if (fEnergy == 3)
+           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPbPb276PythiaMNR.root";
       else if (fEnergy == 4)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPbPb394PythiaMNR.root";
+      else if (fEnergy == 5 || fEnergy == -5)
+           fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPPb5PythiaMNR.root";
       else if (fEnergy == 9 || fEnergy == -9)
            fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPPb88PythiaMNR.root";
       else fFileName = "$ALICE_ROOT/EVGEN/dataCorrHF/CharmPbPb394PythiaMNR.root";
@@ -256,7 +273,7 @@ void AliGenCorrHF::Init()
 
 // particle decay related initialization
 
-    if (gMC) fDecayer = gMC->GetDecayer();
+    if (TVirtualMC::GetMC()) fDecayer = TVirtualMC::GetMC()->GetDecayer();
     fDecayer->SetForceDecay(fForceDecay);
     fDecayer->Init();
 
@@ -317,13 +334,17 @@ void AliGenCorrHF::Generate()
   TDatabasePDG *pDataBase = TDatabasePDG::Instance();
 
   // Calculating vertex position per event
-  for (i=0;i<2;i++){
-    for (j=0;j<3;j++) origin0[i][j]=fOrigin[j];
-    if (fVertexSmear==kPerEvent) {
-      Vertex();
+  if (fVertexSmear==kPerEvent) {
+    Vertex();
+    for (i=0;i<2;i++){
       for (j=0;j<3;j++) origin0[i][j]=fVertex[j];
     }
   }
+  else {
+    for (i=0;i<2;i++){
+      for (j=0;j<3;j++) origin0[i][j]=fOrigin[j];
+    }
+  }
   
   ipa  = 0;
   ipa1 = 0;
@@ -338,7 +359,8 @@ void AliGenCorrHF::Generate()
     
     GetHadronPair(fFile, fQuark, yq[0], yq[1], ptq[0], ptq[1], ihadron[0], ihadron[1], plh[0], plh[1], pth[0], pth[1]);
     
-    if (fEnergy == 9 || fEnergy == -9) {      // boost particles from c.m.s. to ALICE lab frame
+    // Boost particles from c.m.s. to ALICE lab frame for p-Pb & Pb-p collisions
+    if (fEnergy == 5 || fEnergy == -5 || fEnergy == 9 || fEnergy == -9) {
       Double_t dyBoost = 0.47;
       Double_t beta  = TMath::TanH(dyBoost);
       Double_t gamma = 1./TMath::Sqrt((1.-beta)*(1.+beta));
@@ -347,7 +369,7 @@ void AliGenCorrHF::Generate()
       yq[1] += dyBoost;
       plh[0] = gb * TMath::Sqrt(plh[0]*plh[0] + pth[0]*pth[0]) + gamma * plh[0];
       plh[1] = gb * TMath::Sqrt(plh[1]*plh[1] + pth[1]*pth[1]) + gamma * plh[1];
-      if (fEnergy == 9) {
+      if (fEnergy == 5 || fEnergy == 9) {
        yq[0] *= -1;
        yq[1] *= -1;
        plh[0] *= -1;
@@ -587,8 +609,8 @@ void AliGenCorrHF::Generate()
   delete particleshad2;
   delete particleslist;
  
-  delete pSelected;
-  delete trackIt;
+  delete[] pSelected;
+  delete[] trackIt;
 }
 //____________________________________________________________________________________
 void AliGenCorrHF::IpCharm(TH2F *hProbHH, Int_t &pdg3, Int_t &pdg4)
@@ -798,10 +820,11 @@ void AliGenCorrHF::LoadTracks(Int_t iquark, Float_t *pq,
        TParticle* iparticle = (TParticle *) particles->At(i);
        Int_t kf  = iparticle->GetPdgCode();
        Int_t jpa = iparticle->GetFirstMother()-1;
-
-       och[0] = origin0[0]+iparticle->Vx()/10;
-       och[1] = origin0[1]+iparticle->Vy()/10;
-       och[2] = origin0[2]+iparticle->Vz()/10;
+        Int_t ksc  = iparticle->GetStatusCode();
+       // RS: note, the conversion mm->cm is done now in the decayer. The time is ignored here!
+       och[0] = origin0[0]+iparticle->Vx();
+       och[1] = origin0[1]+iparticle->Vy();
+       och[2] = origin0[2]+iparticle->Vz();
        pc[0]  = iparticle->Px();
        pc[1]  = iparticle->Py();
        pc[2]  = iparticle->Pz();
@@ -814,7 +837,7 @@ void AliGenCorrHF::LoadTracks(Int_t iquark, Float_t *pq,
        
        PushTrack(fTrackIt*trackIt[i], iparent, kf,
                  pc, och, polar,
-                 0, kPDecay, nt, wgtch);
+                 0, kPDecay, nt, wgtch,ksc);
        pParent[i] = nt;
        KeepTrack(nt); 
        fNprimaries++;
@@ -826,4 +849,3 @@ void AliGenCorrHF::LoadTracks(Int_t iquark, Float_t *pq,
  
   return;
 }
-