]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythia.cxx
ProcessOutput update: deal with both TFileMerged merged data (histos inside TList...
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
index 6523bc601b4fa5a7efb3dff789335648b4d550aa..79827a074ad630e77d9c8ae0873f3c01c769233d 100644 (file)
@@ -76,6 +76,7 @@ AliGenPythia::AliGenPythia():
     fGinit(1),
     fGfinal(1),
     fHadronisation(1),
+    fPatchOmegaDalitz(0), 
     fNpartons(0),
     fReadFromFile(0),
     fQuench(0),
@@ -120,7 +121,7 @@ AliGenPythia::AliGenPythia():
     fCountMode(kCountAll),      
     fHeader(0),  
     fRL(0),      
-    fFileName(0),
+    fkFileName(0),
     fFragPhotonInCalo(kFALSE),
     fPi0InCalo(kFALSE) ,
     fPhotonInCalo(kFALSE),
@@ -172,6 +173,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fGinit(kTRUE),
      fGfinal(kTRUE),
      fHadronisation(kTRUE),
+     fPatchOmegaDalitz(0), 
      fNpartons(0),
      fReadFromFile(kFALSE),
      fQuench(kFALSE),
@@ -216,7 +218,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fCountMode(kCountAll),      
      fHeader(0),  
      fRL(0),      
-     fFileName(0),
+     fkFileName(0),
      fFragPhotonInCalo(kFALSE),
      fPi0InCalo(kFALSE) ,
      fPhotonInCalo(kFALSE),
@@ -368,7 +370,7 @@ void AliGenPythia::Init()
 
 
     if (fReadFromFile) {
-       fRL  =  AliRunLoader::Open(fFileName, "Partons");
+       fRL  =  AliRunLoader::Open(fkFileName, "Partons");
        fRL->LoadKinematics();
        fRL->LoadHeader();
     } else {
@@ -431,6 +433,10 @@ void AliGenPythia::Init()
        fParentSelect[0] =   431;
        fFlavorSelect    =   4; 
        break;
+    case kPyLambdacppMNR:
+       fParentSelect[0] =  4122;
+       fFlavorSelect    =   4; 
+       break;      
     case kPyBeauty:
     case kPyBeautyJets:
     case kPyBeautyPbPbMNR:
@@ -550,9 +556,9 @@ void AliGenPythia::Generate()
     if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
     fDecayer->ForceDecay();
 
-    Float_t polar[3]   =   {0,0,0};
-    Float_t origin[3]  =   {0,0,0};
-    Float_t p[4];
+    Double_t polar[3]   =   {0,0,0};
+    Double_t origin[3]  =   {0,0,0};
+    Double_t p[4];
 //  converts from mm/c to s
     const Float_t kconv=0.001/2.999792458e8;
 //
@@ -584,6 +590,7 @@ void AliGenPythia::Generate()
            AliFastGlauber::Instance()->GetRandomBHard(bimp);
            fPythia->SetPARJ(197, bimp);
            fImpact = bimp;
+           fPythia->Qpygin0();
        } 
 //
 // Either produce new event or read partons from file
@@ -622,11 +629,21 @@ void AliGenPythia::Generate()
 //
 // .. and perform hadronisation
 //     printf("Calling hadronisation %d\n", fPythia->GetN());
-           fPythia->Pyexec();  
+
+           if (fPatchOmegaDalitz) {
+             fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+             fPythia->Pyexec();
+             fPythia->DalitzDecays();
+             fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+           } 
+           fPythia->Pyexec();
        }
        fTrials++;
        fPythia->ImportParticles(&fParticles,"All");
+       
        if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
+       if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
+       
 //
 //
 //
@@ -653,6 +670,7 @@ void AliGenPythia::Generate()
        Int_t nTkbles = 0;   // Trackable particles
        if (fProcess != kPyMbDefault && 
            fProcess != kPyMb && 
+           fProcess != kPyMbAtlasTuneMC09 && 
            fProcess != kPyMbWithDirectPhoton && 
            fProcess != kPyJets && 
            fProcess != kPyDirectGamma &&
@@ -877,9 +895,9 @@ Int_t  AliGenPythia::GenerateMB()
 //
     Int_t i, kf, nt, iparent;
     Int_t nc = 0;
-    Float_t p[4];
-    Float_t polar[3]   =   {0,0,0};
-    Float_t origin[3]  =   {0,0,0};
+    Double_t p[4];
+    Double_t polar[3]   =   {0,0,0};
+    Double_t origin[3]  =   {0,0,0};
 //  converts from mm/c to s
     const Float_t kconv=0.001/2.999792458e8;
     
@@ -907,7 +925,7 @@ Int_t  AliGenPythia::GenerateMB()
 
       Int_t pdg  = 0; 
       if (fFragPhotonInCalo) pdg = 22   ; // Photon
-      else if (fPi0InCalo) pdg = 111 ;    // Pi0
+      else if (fPi0InCalo)   pdg = 111 ;    // Pi0
 
       for (i=0; i< np; i++) {
        TParticle* iparticle = (TParticle *) fParticles.At(i);
@@ -926,8 +944,10 @@ Int_t  AliGenPythia::GenerateMB()
            }
        }
       }
-      if(!ok)
-       return 0;
+      if(!ok) {
+         delete[] pParent;
+         return 0;
+      }
     }
 
     // Select beauty jets with electron in EMCAL
@@ -951,11 +971,13 @@ Int_t  AliGenPythia::GenerateMB()
            ok =kTRUE;
        }
       }
-      if(!ok)
+      if(!ok) {
+       delete[] pParent;
        return 0;
+      }
+      
       AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
     }
-    
     // Check for minimum multiplicity
     if (fTriggerMultiplicity > 0) {
       Int_t multiplicity = 0;
@@ -985,7 +1007,6 @@ Int_t  AliGenPythia::GenerateMB()
        delete [] pParent;
        return 0;
       }
-      
       Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
     }    
     
@@ -1019,8 +1040,10 @@ Int_t  AliGenPythia::GenerateMB()
       if(!okd && iphcand != -1) // execute rotation in phi 
           RotatePhi(iphcand,okd);
       
-      if(!okd)
-       return 0;
+      if(!okd) {
+         delete [] pParent;
+         return 0;
+      }
     }
     
     if (fTriggerParticle) {
@@ -1049,6 +1072,7 @@ Int_t  AliGenPythia::GenerateMB()
       TParticle *mother;
       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
       Bool_t  theChild=kFALSE;
+      Bool_t  triggered=kFALSE;
       Float_t y;  
       Int_t   pdg,mpdg,mpdgUpperFamily;
       for(i=0; i<np; i++) {
@@ -1060,6 +1084,9 @@ Int_t  AliGenPythia::GenerateMB()
                             (partCheck->Energy()-partCheck->Pz()+1.e-13));
          if(y>fYMin && y<fYMax) inYcut=kTRUE;
        }
+       if(fTriggerParticle) {
+         if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
+       }
        if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
          Int_t mi = partCheck->GetFirstMother() - 1;
          if(mi<0) continue;
@@ -1082,6 +1109,10 @@ Int_t  AliGenPythia::GenerateMB()
        delete[] pParent;
        return 0;       
       }
+      if(fTriggerParticle && !triggered) { // particle requested is not produced
+       delete[] pParent;
+       return 0;       
+      }
 
     }
 
@@ -1090,6 +1121,7 @@ Int_t  AliGenPythia::GenerateMB()
          fProcess == kPyZ ||
          fProcess == kPyMbDefault ||
          fProcess == kPyMb ||
+         fProcess == kPyMbAtlasTuneMC09 ||
          fProcess == kPyMbWithDirectPhoton ||
          fProcess == kPyMbNonDiffr)  
         && (fCutOnChild == 1) ) {
@@ -1210,7 +1242,7 @@ void AliGenPythia::MakeHeader()
 //
 // Event Vertex 
     fHeader->SetPrimaryVertex(fVertex);
-    
+    fHeader->SetInteractionTime(fEventTime);
 //
 // Number of primaries
     fHeader->SetNProduced(fNprimaries);
@@ -1250,8 +1282,10 @@ void AliGenPythia::MakeHeader()
 // Store quenching parameters
 //
     if (fQuench){
-       Double_t z[4];
-       Double_t xp, yp;
+       Double_t z[4] = {0.};
+       Double_t xp = 0.;
+       Double_t yp = 0.;
+       
        if (fQuench == 1) {
            // Pythia::Quench()
            fPythia->GetQuenchingParameters(xp, yp, z);
@@ -1289,7 +1323,7 @@ void AliGenPythia::MakeHeader()
     fHeader = 0x0;
 }
 
-Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
+Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
 {
 // Check the kinematic trigger condition
 //
@@ -1422,7 +1456,7 @@ void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
   }
 }
 
-void  AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
+void  AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
 {
   //
   // Load event into Pythia Common Block
@@ -1441,17 +1475,19 @@ void  AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
   
   for (Int_t part = 0; part < npart; part++) {
     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
+    if (!mPart) continue;
+    
     Int_t kf     =  mPart->GetPdgCode();
     Int_t ks     =  mPart->GetStatusCode();
     Int_t idf    =  mPart->GetFirstDaughter();
     Int_t idl    =  mPart->GetLastDaughter();
     
     if (reHadr) {
-           if (ks == 11 || ks == 12) {
-        ks  -= 10;
-        idf  = -1;
-        idl  = -1;
-           }
+       if (ks == 11 || ks == 12) {
+           ks  -= 10;
+           idf  = -1;
+           idl  = -1;
+       }
     }
     
     Float_t px = mPart->Px();
@@ -1574,7 +1610,7 @@ void AliGenPythia::GetSubEventTime()
   return;
 }
 
-Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
+Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
 {
   // Is particle in EMCAL acceptance? 
   // phi in degrees, etamin=-etamax
@@ -1585,7 +1621,7 @@ Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
     return kFALSE;
 }
 
-Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
+Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta) const
 {
   // Is particle in PHOS acceptance? 
   // Acceptance slightly larger considered.
@@ -1614,8 +1650,8 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
   //loop for all particles and produce the phi rotation
   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
   Double_t oldphi, newphi;
-  Double_t newVx, newVy, R, Vz, time; 
-  Double_t newPx, newPy, pt, Pz, e;
+  Double_t newVx, newVy, r, vZ, time; 
+  Double_t newPx, newPy, pt, pz, e;
   for(Int_t i=0; i< np; i++) {
       TParticle* iparticle = (TParticle *) fParticles.At(i);
       oldphi = iparticle->Phi();
@@ -1623,21 +1659,21 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle 
       if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
       
-      R = iparticle->R();
-      newVx = R*TMath::Cos(newphi);
-      newVy = R*TMath::Sin(newphi);
-      Vz = iparticle->Vz(); // don't transform
+      r = iparticle->R();
+      newVx = r * TMath::Cos(newphi);
+      newVy = r * TMath::Sin(newphi);
+      vZ   = iparticle->Vz(); // don't transform
       time = iparticle->T(); // don't transform
       
       pt = iparticle->Pt();
-      newPx = pt*TMath::Cos(newphi);
-      newPy = pt*TMath::Sin(newphi);
-      Pz = iparticle->Pz(); // don't transform
-      e = iparticle->Energy(); // don't transform
+      newPx = pt * TMath::Cos(newphi);
+      newPy = pt * TMath::Sin(newphi);
+      pz = iparticle->Pz(); // don't transform
+      e  = iparticle->Energy(); // don't transform
       
       // apply rotation 
-      iparticle->SetProductionVertex(newVx, newVy, Vz, time);
-      iparticle->SetMomentum(newPx, newPy, Pz, e);
+      iparticle->SetProductionVertex(newVx, newVy, vZ, time);
+      iparticle->SetMomentum(newPx, newPy, pz, e);
       
   } //end particle loop 
   
@@ -1649,46 +1685,4 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
 }
 
 
-#ifdef never
-void AliGenPythia::Streamer(TBuffer &R__b)
-{
-   // Stream an object of class AliGenPythia.
-
-   if (R__b.IsReading()) {
-      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
-      AliGenerator::Streamer(R__b);
-      R__b >> (Int_t&)fProcess;
-      R__b >> (Int_t&)fStrucFunc;
-      R__b >> (Int_t&)fForceDecay;
-      R__b >> fEnergyCMS;
-      R__b >> fKineBias;
-      R__b >> fTrials;
-      fParentSelect.Streamer(R__b);
-      fChildSelect.Streamer(R__b);
-      R__b >> fXsection;
-//      (AliPythia::Instance())->Streamer(R__b);
-      R__b >> fPtHardMin;
-      R__b >> fPtHardMax;
-//      if (fDecayer) fDecayer->Streamer(R__b);
-   } else {
-      R__b.WriteVersion(AliGenPythia::IsA());
-      AliGenerator::Streamer(R__b);
-      R__b << (Int_t)fProcess;
-      R__b << (Int_t)fStrucFunc;
-      R__b << (Int_t)fForceDecay;
-      R__b << fEnergyCMS;
-      R__b << fKineBias;
-      R__b << fTrials;
-      fParentSelect.Streamer(R__b);
-      fChildSelect.Streamer(R__b);
-      R__b << fXsection;
-//      R__b << fPythia;
-      R__b << fPtHardMin;
-      R__b << fPtHardMax;
-      //     fDecayer->Streamer(R__b);
-   }
-}
-#endif
-
-