]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythia.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
making FillQA consistent with AOD
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 8df77430d2d9b9ff2736c2d218db464cb8cfa2d0..f50efaaf2d7eb8e5d72edeeecedb3c5e314bc1a6 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -46,7 +46,6 @@
 #include "AliMC.h"
 #include "AliLog.h"
 #include "PyquenCommon.h"
-#include "AliLog.h"
 
 ClassImp(AliGenPythia)
 
@@ -81,6 +80,7 @@ AliGenPythia::AliGenPythia():
     fPatchOmegaDalitz(0), 
     fNpartons(0),
     fReadFromFile(0),
+    fReadLHEF(0),
     fQuench(0),
     fQhat(0.),
     fLength(0.),
@@ -117,6 +117,8 @@ AliGenPythia::AliGenPythia():
     fFeedDownOpt(kTRUE),    
     fFragmentation(kTRUE),
     fSetNuclei(kFALSE),
+    fUseNuclearPDF(kFALSE),
+    fUseLorentzBoost(kTRUE),
     fNewMIS(kFALSE),   
     fHFoff(kFALSE),    
     fNucPdf(0),
@@ -131,6 +133,7 @@ AliGenPythia::AliGenPythia():
     fHeader(0),  
     fRL(0),      
     fkFileName(0),
+    fkNameLHEF(0),
     fFragPhotonInCalo(kFALSE),
     fHadronInCalo(kFALSE) ,
     fPi0InCalo(kFALSE) ,
@@ -193,6 +196,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fPatchOmegaDalitz(0), 
      fNpartons(0),
      fReadFromFile(kFALSE),
+     fReadLHEF(0),
      fQuench(kFALSE),
      fQhat(0.),
      fLength(0.),
@@ -229,6 +233,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fFeedDownOpt(kTRUE),    
      fFragmentation(kTRUE),
      fSetNuclei(kFALSE),
+     fUseNuclearPDF(kFALSE),
+     fUseLorentzBoost(kTRUE),
      fNewMIS(kFALSE),   
      fHFoff(kFALSE),    
      fNucPdf(0),
@@ -243,6 +249,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fHeader(0),  
      fRL(0),      
      fkFileName(0),
+     fkNameLHEF(0),
      fFragPhotonInCalo(kFALSE),
      fHadronInCalo(kFALSE) ,
      fPi0InCalo(kFALSE) ,
@@ -379,7 +386,10 @@ void AliGenPythia::Init()
     fPythia->SetCKIN(7,fYHardMin);
     fPythia->SetCKIN(8,fYHardMax);
     
-    if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);  
+    if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
+  
+    if(fUseNuclearPDF)
+      fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
     // Fragmentation?
     if (fFragmentation) {
       fPythia->SetMSTP(111,1);
@@ -401,7 +411,8 @@ void AliGenPythia::Init()
        fPythia->SetMSTP(91,0);
     }
 
-
+          if (fReadLHEF) fPythia->OpenFortranFile(97, const_cast<char*>(fkNameLHEF));
+       
     if (fReadFromFile) {
        fRL  =  AliRunLoader::Open(fkFileName, "Partons");
        fRL->LoadKinematics();
@@ -439,6 +450,7 @@ void AliGenPythia::Init()
     case kPyCharmpPbMNR:
     case kPyCharmppMNR:
     case kPyCharmppMNRwmi:
+    case kPyCharmPWHG:
        fParentSelect[0] =   411;
        fParentSelect[1] =   421;
        fParentSelect[2] =   431;
@@ -476,6 +488,7 @@ void AliGenPythia::Init()
     case kPyBeautypPbMNR:
     case kPyBeautyppMNR:
     case kPyBeautyppMNRwmi:
+    case kPyBeautyPWHG:
        fParentSelect[0]=  511;
        fParentSelect[1]=  521;
        fParentSelect[2]=  531;
@@ -506,11 +519,17 @@ void AliGenPythia::Init()
     case kPyMbNonDiffr:
     case kPyMbMSEL1:
     case kPyJets:
+    case kPyJetsPWHG:
     case kPyDirectGamma:
     case kPyLhwgMb:    
        break;
+    case kPyWPWHG:
     case kPyW:
     case kPyZ:
+    case kPyZgamma:
+    case kPyMBRSingleDiffraction:
+    case kPyMBRDoubleDiffraction:
+    case kPyMBRCentralDiffraction:
         break;
     }
 //
@@ -541,14 +560,19 @@ void AliGenPythia::Init()
 //
 //
     AliGenMC::Init();
+
+    // Reset Lorentz boost if demanded
+    if(!fUseLorentzBoost) {
+      fDyBoost = 0;
+      Warning("Init","Demand to discard Lorentz boost.\n");
+    }
 //
 //
 //  
     if (fSetNuclei) {
-       fDyBoost = 0;
-       Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
+      fDyBoost = 0;
+      Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
     }
-    
     fPythia->SetPARJ(200, 0.0);
     fPythia->SetPARJ(199, 0.0);
     fPythia->SetPARJ(198, 0.0);
@@ -583,6 +607,11 @@ void AliGenPythia::Init()
        fPythia->SetMSTJ(44, 2);  // option to run alpha_s
        fPythia->SetPARJ(82, 1.); // Cut off for parton showers
     }
+  
+  if ( AliLog::GetDebugLevel("","AliGenPythia") >= 1 ) {
+    fPythia->Pystat(4);
+    fPythia->Pystat(5);
+  }
 }
 
 void AliGenPythia::Generate()
@@ -713,9 +742,14 @@ void AliGenPythia::Generate()
            fProcess != kPyMbMSEL1     &&
            fProcess != kPyW && 
            fProcess != kPyZ &&
+      fProcess != kPyZgamma &&
            fProcess != kPyCharmppMNRwmi && 
            fProcess != kPyBeautyppMNRwmi &&
-           fProcess != kPyBeautyJets) {
+           fProcess != kPyBeautyJets &&
+     fProcess != kPyWPWHG &&
+     fProcess != kPyJetsPWHG &&
+     fProcess != kPyCharmPWHG &&
+     fProcess != kPyBeautyPWHG) {
            
            for (i = 0; i < np; i++) {
                TParticle* iparticle = (TParticle *) fParticles.At(i);
@@ -944,10 +978,12 @@ Int_t  AliGenPythia::GenerateMB()
 
     Int_t* pParent = new Int_t[np];
     for (i=0; i< np; i++) pParent[i] = -1;
-     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
-       TParticle* jet1 = (TParticle *) fParticles.At(6);
+    if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
+       && fEtMinJet > 0.) {
+       TParticle* jet1 = (TParticle *) fParticles.At(6);
        TParticle* jet2 = (TParticle *) fParticles.At(7);
-       if (!CheckTrigger(jet1, jet2)) {
+       
+       if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
          delete [] pParent;
          return 0;
        }
@@ -1081,8 +1117,11 @@ Int_t  AliGenPythia::GenerateMB()
     }
 
     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
-    if ( (fProcess == kPyW ||
+    if ( (
+    fProcess == kPyWPWHG ||
+    fProcess == kPyW ||
          fProcess == kPyZ ||
+    fProcess == kPyZgamma ||
          fProcess == kPyMbDefault ||
          fProcess == kPyMb ||
          fProcess == kPyMbAtlasTuneMC09 ||
@@ -1104,7 +1143,7 @@ Int_t  AliGenPythia::GenerateMB()
        Int_t km = iparticle->GetFirstMother();
        if (
            (((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
-           ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)
+           ((fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) && ks == 21 && km == 0 && i>1)
            ) 
          {
             nc++;
@@ -1179,10 +1218,11 @@ void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
 {
 // Treat protons as inside nuclei with mass numbers a1 and a2  
 
-    fAProjectile = a1;
-    fATarget     = a2;
-    fNucPdf      = pdfset;  // 0 EKS98 9 EPS09LO 19 EPS09NLO
-    fSetNuclei   = kTRUE;
+    fAProjectile   = a1;
+    fATarget       = a2;
+    fNucPdf        = pdfset;  // 0 EKS98 9 EPS09LO 19 EPS09NLO
+    fUseNuclearPDF = kTRUE;
+    fSetNuclei     = kTRUE;
 }
 
 
@@ -1219,10 +1259,13 @@ void AliGenPythia::MakeHeader()
 // Number of primaries
     fHeader->SetNProduced(fNprimaries);
 //
+// Event weight
+    fHeader->SetEventWeight(fPythia->GetVINT(97));
+//
 // Jets that have triggered
 
     //Need to store jets for b-jet studies too!
-    if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
+    if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
     {
        Int_t ntrig, njet;
        Float_t jets[4][10];
@@ -1286,8 +1329,14 @@ void AliGenPythia::MakeHeader()
            ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
     }
 //
-// Store pt^hard 
+// Store pt^hard and cross-section
     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
+    ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1)); 
+               
+//
+// Store Event Weight
+    ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7));
+                               
 //
 //  Pass header
 //
@@ -1310,7 +1359,7 @@ Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
     pdg[1] = jet2->GetPdgCode();    
     Bool_t   triggered = kFALSE;
 
-    if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi) {
+    if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) {
        Int_t njets = 0;
        Int_t ntrig = 0;
        Float_t jets[4][10];
@@ -1791,7 +1840,9 @@ Bool_t AliGenPythia::CheckDiffraction()
        TParticle *  part = (TParticle *) fParticles.At(iPart);
        Double_t E= part->Energy();
        Double_t P= part->P();
-       M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
+       Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
+       if(M2<0)  return kFALSE;
+       M= TMath::Sqrt(M2);
      }
 
      Double_t Mmin, Mmax, wSD, wDD, wND;