]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythia.cxx
Fix axis label
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
index e490e6f0c8244f9d0c07146c385a66632eca0b7d..f2a4023d2c33c55de511ebc6ffd65d7550268307 100644 (file)
@@ -46,7 +46,6 @@
 #include "AliMC.h"
 #include "AliLog.h"
 #include "PyquenCommon.h"
-#include "AliLog.h"
 
 ClassImp(AliGenPythia)
 
@@ -71,18 +70,20 @@ AliGenPythia::AliGenPythia():
     fFlavorSelect(0),
     fXsection(0.),
     fPythia(0),
+    fWeightPower(0.),
     fPtHardMin(0.),
     fPtHardMax(1.e4),
     fYHardMin(-1.e10),
     fYHardMax(1.e10),
     fGinit(1),
     fGfinal(1),
+    fCRoff(0),
     fHadronisation(1),
     fPatchOmegaDalitz(0), 
+    fDecayerExodus(0),
     fNpartons(0),
     fReadFromFile(0),
     fReadLHEF(0),
-    fPartonShower(1),
     fQuench(0),
     fQhat(0.),
     fLength(0.),
@@ -119,15 +120,21 @@ AliGenPythia::AliGenPythia():
     fFeedDownOpt(kTRUE),    
     fFragmentation(kTRUE),
     fSetNuclei(kFALSE),
+    fUseNuclearPDF(kFALSE),
+    fUseLorentzBoost(kTRUE),
     fNewMIS(kFALSE),   
     fHFoff(kFALSE),    
     fNucPdf(0),
     fTriggerParticle(0),
     fTriggerEta(0.9),     
+    fTriggerY(999.),     
+    fTriggerEtaMin(0.9),     
     fTriggerMinPt(-1),  
     fTriggerMaxPt(1000),  
     fTriggerMultiplicity(0),
     fTriggerMultiplicityEta(0),
+    fTriggerMultiplicityEtaMin(0),
+    fTriggerMultiplicityEtaMax(0),
     fTriggerMultiplicityPtMin(0),
     fCountMode(kCountAll),      
     fHeader(0),  
@@ -186,18 +193,20 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fFlavorSelect(0),
      fXsection(0.),
      fPythia(0),
+     fWeightPower(0.),
      fPtHardMin(0.),
      fPtHardMax(1.e4),
      fYHardMin(-1.e10),
      fYHardMax(1.e10),
      fGinit(kTRUE),
      fGfinal(kTRUE),
+     fCRoff(kFALSE),
      fHadronisation(kTRUE),
-     fPatchOmegaDalitz(0), 
+     fPatchOmegaDalitz(0),
+     fDecayerExodus(0), 
      fNpartons(0),
      fReadFromFile(kFALSE),
      fReadLHEF(0),
-     fPartonShower(1),
      fQuench(kFALSE),
      fQhat(0.),
      fLength(0.),
@@ -234,15 +243,21 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fFeedDownOpt(kTRUE),    
      fFragmentation(kTRUE),
      fSetNuclei(kFALSE),
+     fUseNuclearPDF(kFALSE),
+     fUseLorentzBoost(kTRUE),
      fNewMIS(kFALSE),   
      fHFoff(kFALSE),    
      fNucPdf(0),
      fTriggerParticle(0),
      fTriggerEta(0.9), 
+     fTriggerY(999.), 
+     fTriggerEtaMin(0.9),     
      fTriggerMinPt(-1),  
      fTriggerMaxPt(1000),      
      fTriggerMultiplicity(0),
      fTriggerMultiplicityEta(0),
+     fTriggerMultiplicityEtaMin(0),
+     fTriggerMultiplicityEtaMax(0),
      fTriggerMultiplicityPtMin(0),
      fCountMode(kCountAll),      
      fHeader(0),  
@@ -378,16 +393,17 @@ void AliGenPythia::Init()
 //
     fParentWeight=1./Float_t(fNpart);
 //
-
-
+    if (fWeightPower != 0)
+      fPythia->SetWeightPower(fWeightPower);
     fPythia->SetCKIN(3,fPtHardMin);
     fPythia->SetCKIN(4,fPtHardMax);
     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);
@@ -400,6 +416,8 @@ void AliGenPythia::Init()
     fPythia->SetMSTP(61,fGinit);
 //  final state radiation
     fPythia->SetMSTP(71,fGfinal);
+    //color reconnection strength
+    if(fCRoff==1)fPythia->SetMSTP(95,0);
 //  pt - kick
     if (fPtKick > 0.) {
        fPythia->SetMSTP(91,1);
@@ -521,8 +539,10 @@ void AliGenPythia::Init()
     case kPyDirectGamma:
     case kPyLhwgMb:    
        break;
+    case kPyWPWHG:
     case kPyW:
     case kPyZ:
+    case kPyZgamma:
     case kPyMBRSingleDiffraction:
     case kPyMBRDoubleDiffraction:
     case kPyMBRCentralDiffraction:
@@ -556,14 +576,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);
@@ -598,6 +623,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()
@@ -686,8 +716,29 @@ void AliGenPythia::Generate()
              fPythia->DalitzDecays();
              fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
            } 
-           fPythia->Pyexec();
-       }
+         
+      else  if (fDecayerExodus) {
+
+        fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
+        fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
+        fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
+        fPythia->Pyexec();
+        fPythia->OmegaDalitz();
+        fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
+        fPythia->PizeroDalitz();
+        fPythia->PhiDalitz();
+        fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
+        fPythia->EtaDalitz();
+        fPythia->EtaprimeDalitz();
+        fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
+        fPythia->RhoDirect();
+        fPythia->OmegaDirect();
+        fPythia->PhiDirect();
+        fPythia->JPsiDirect();
+      }
+    
+      fPythia->Pyexec();
+  }
        fTrials++;
        fPythia->ImportParticles(&fParticles,"All");
        
@@ -728,12 +779,14 @@ void AliGenPythia::Generate()
            fProcess != kPyMbMSEL1     &&
            fProcess != kPyW && 
            fProcess != kPyZ &&
+            fProcess != kPyZgamma &&
            fProcess != kPyCharmppMNRwmi && 
            fProcess != kPyBeautyppMNRwmi &&
            fProcess != kPyBeautyJets &&
-     fProcess != kPyJetsPWHG &&
-     fProcess != kPyCharmPWHG &&
-     fProcess != kPyBeautyPWHG) {
+            fProcess != kPyWPWHG &&
+            fProcess != kPyJetsPWHG &&
+            fProcess != kPyCharmPWHG &&
+            fProcess != kPyBeautyPWHG) {
            
            for (i = 0; i < np; i++) {
                TParticle* iparticle = (TParticle *) fParticles.At(i);
@@ -1010,6 +1063,10 @@ Int_t  AliGenPythia::GenerateMB()
        // eta cut
        if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
          continue;
+       //multiplicity check for a given eta range
+       if ((fTriggerMultiplicityEtaMin != fTriggerMultiplicityEtaMax) && 
+           (iparticle->Eta() < fTriggerMultiplicityEtaMin || iparticle->Eta() > fTriggerMultiplicityEtaMax))
+         continue;
        // pt cut
        if (iparticle->Pt() < fTriggerMultiplicityPtMin) 
            continue;
@@ -1036,7 +1093,12 @@ Int_t  AliGenPythia::GenerateMB()
            kf = CheckPDGCode(iparticle->GetPdgCode());
            if (kf != fTriggerParticle) continue;
            if (iparticle->Pt() == 0.) continue;
-           if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+           if (TMath::Abs(iparticle->Y()) > fTriggerY) continue;
+           if (fTriggerEtaMin == fTriggerEta) {
+             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+           } else {
+             if (iparticle->Eta() > fTriggerEta || iparticle->Eta() < fTriggerEtaMin) continue;
+           }
            if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
            triggered = kTRUE;
            break;
@@ -1058,8 +1120,8 @@ Int_t  AliGenPythia::GenerateMB()
       Bool_t  theChild=kFALSE;
       Bool_t  triggered=kFALSE;
       Float_t y;  
-      Int_t   pdg,mpdg,mpdgUpperFamily;
-      for(i=0; i<np; i++) {
+      Int_t   pdg, mpdg, mpdgUpperFamily;
+      for(i = 0; i < np; i++) {
        partCheck = (TParticle*)fParticles.At(i);
        pdg = partCheck->GetPdgCode();  
        if(TMath::Abs(pdg) == fFlavorSelect) { // quark  
@@ -1101,8 +1163,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 ||
@@ -1114,7 +1179,7 @@ Int_t  AliGenPythia::GenerateMB()
        return 0;
       }
     }
-  
+    
 
     for (i = 0; i < np; i++) {
        Int_t trackIt = 0;
@@ -1199,10 +1264,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;
 }
 
 
@@ -1239,6 +1305,9 @@ 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!
@@ -1306,13 +1375,15 @@ 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));
-                               
+    ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7)*fPythia->GetPARI(10));
+    // PARI(7) is 1 or -1, for weighted generation with accept/reject, e.g. POWHEG
+    // PARI(10) is a weight associated with reweighted generation, using Pyevwt
 //
 //  Pass header
 //