Comments for Pytune added.

[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index e46828981db4342665e2d815f8bdf02e1ec44bc3..947bf64ba7c7f9b85f52484c4652873caed3ffa0 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -52,6 +52,7 @@ ClassImp(AliGenPythia)
 AliGenPythia::AliGenPythia():
     AliGenMC(),
     fProcess(kPyCharm),          
+    fItune(-1),
     fStrucFunc(kCTEQ5L), 
     fKineBias(0.),
     fTrials(0),
@@ -80,6 +81,7 @@ AliGenPythia::AliGenPythia():
     fQuench(0),
     fQhat(0.),
     fLength(0.),
+    fImpact(0.),
     fPtKick(1.),
     fFullEvent(kTRUE),
     fDecayer(new AliDecayerPythia()),
@@ -121,11 +123,13 @@ AliGenPythia::AliGenPythia():
     fFragPhotonInCalo(kFALSE),
     fPi0InCalo(kFALSE) ,
     fPhotonInCalo(kFALSE),
+    fEleInEMCAL(kFALSE),
     fCheckEMCAL(kFALSE),
     fCheckPHOS(kFALSE),
     fCheckPHOSeta(kFALSE),
     fFragPhotonOrPi0MinPt(0), 
     fPhotonMinPt(0), 
+    fElectronMinPt(0), 
     fPHOSMinPhi(219.),
     fPHOSMaxPhi(321.),
     fPHOSEta(0.13),
@@ -144,6 +148,7 @@ AliGenPythia::AliGenPythia():
 AliGenPythia::AliGenPythia(Int_t npart)
     :AliGenMC(npart),
      fProcess(kPyCharm),          
+     fItune(-1),
      fStrucFunc(kCTEQ5L), 
      fKineBias(0.),
      fTrials(0),
@@ -172,6 +177,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fQuench(kFALSE),
      fQhat(0.),
      fLength(0.),
+     fImpact(0.),
      fPtKick(1.),
      fFullEvent(kTRUE),
      fDecayer(new AliDecayerPythia()),
@@ -213,11 +219,13 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fFragPhotonInCalo(kFALSE),
      fPi0InCalo(kFALSE) ,
      fPhotonInCalo(kFALSE),
+     fEleInEMCAL(kFALSE),
      fCheckEMCAL(kFALSE),
      fCheckPHOS(kFALSE),
      fCheckPHOSeta(kFALSE),
      fFragPhotonOrPi0MinPt(0),
      fPhotonMinPt(0),
+     fElectronMinPt(0),
      fPHOSMinPhi(219.),
      fPHOSMaxPhi(321.),
      fPHOSEta(0.13),
@@ -367,7 +375,7 @@ void AliGenPythia::Init()
        fRL = 0x0;
     }
  //
-    fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
+    fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
     //  Forward Paramters to the AliPythia object
     fDecayer->SetForceDecay(fForceDecay);    
 // Switch off Heavy Flavors on request  
@@ -424,6 +432,7 @@ void AliGenPythia::Init()
        fFlavorSelect    =   4; 
        break;
     case kPyBeauty:
+    case kPyBeautyJets:
     case kPyBeautyPbPbMNR:
     case kPyBeautypPbMNR:
     case kPyBeautyppMNR:
@@ -528,15 +537,9 @@ void AliGenPythia::Init()
        fPythia->SetPARJ(200, 1.);
        fPythia->SetPARJ(198, fQhat);
        fPythia->SetPARJ(199, fLength);
-       
-       fPythia->SetMSTJ(41, 1);  // QCD radiation only
        fPythia->SetMSTJ(42, 2);  // angular ordering
        fPythia->SetMSTJ(44, 2);  // option to run alpha_s
-       //fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
-       //fPythia->SetMSTJ(50, 0);  // No coherence in first branching
        fPythia->SetPARJ(82, 1.); // Cut off for parton showers
-       //    MSTJ(41) must NOT be 11 or 12, as then FSR may go through PYPTFS
-       //   (kt-ordered cascade) in which medium effects have not been introduced.
     }
 }
 
@@ -578,7 +581,8 @@ void AliGenPythia::Generate()
            Double_t bimp;
            // Quenching comes through medium-modified splitting functions.
            AliFastGlauber::Instance()->GetRandomBHard(bimp);
-           fPythia->SetPARJ(197,bimp);
+           fPythia->SetPARJ(197, bimp);
+           fImpact = bimp;
        } 
 //
 // Either produce new event or read partons from file
@@ -656,7 +660,8 @@ void AliGenPythia::Generate()
            fProcess != kPyW && 
            fProcess != kPyZ &&
            fProcess != kPyCharmppMNRwmi && 
-           fProcess != kPyBeautyppMNRwmi) {
+           fProcess != kPyBeautyppMNRwmi &&
+           fProcess != kPyBeautyJets) {
            
            for (i = 0; i < np; i++) {
                TParticle* iparticle = (TParticle *) fParticles.At(i);
@@ -885,7 +890,10 @@ 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) {
+    //
+    //TO BE CHECKED:  Should we require this for Beauty Jets?
+    //
+    if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets) {
        TParticle* jet1 = (TParticle *) fParticles.At(6);
        TParticle* jet2 = (TParticle *) fParticles.At(7);
        if (!CheckTrigger(jet1, jet2)) {
@@ -923,6 +931,32 @@ Int_t  AliGenPythia::GenerateMB()
       if(!ok)
        return 0;
     }
+
+    // Select beauty jets with electron in EMCAL
+    if (fProcess == kPyBeautyJets && fEleInEMCAL) {
+
+      Bool_t ok = kFALSE;
+
+      Int_t pdg  = 11; //electron
+
+      Float_t pt  = 0.;
+      Float_t eta = 0.;
+      Float_t phi = 0.;
+      for (i=0; i< np; i++) {
+       TParticle* iparticle = (TParticle *) fParticles.At(i);
+       if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg && 
+          iparticle->Pt() > fElectronMinPt){
+         pt = iparticle->Pt();
+         phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+         eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax    
+         if(IsInEMCAL(phi,eta))
+           ok =kTRUE;
+       }
+      }
+      if(!ok)
+       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) {
@@ -1013,7 +1047,10 @@ Int_t  AliGenPythia::GenerateMB()
 
     // Check if there is a ccbar or bbbar pair with at least one of the two
     // in fYMin < y < fYMax
-    if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
+    //
+    // TO BE CHECKED:  Should we require this for beauty jets?
+    //
+    if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
       TParticle *partCheck;
       TParticle *mother;
       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
@@ -1077,7 +1114,10 @@ Int_t  AliGenPythia::GenerateMB()
        Int_t km = iparticle->GetFirstMother();
        if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
            (ks != 1) ||
-           (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
+           //
+           //TO BE CHECKED:  Should we require this for beauty jets?
+           //
+           ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
            nc++;
            if (ks == 1) trackIt = 1;
            Int_t ipa = iparticle->GetFirstMother()-1;
@@ -1186,7 +1226,8 @@ void AliGenPythia::MakeHeader()
 //
 // Jets that have triggered
 
-    if (fProcess == kPyJets || fProcess == kPyDirectGamma)
+    //Need to store jets for b-jet studies too!
+    if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets)
     {
        Int_t ntrig, njet;
        Float_t jets[4][10];
@@ -1236,9 +1277,12 @@ void AliGenPythia::MakeHeader()
        } else if (fQuench == 4) {
            // QPythia
            Double_t xy[2];
-           AliFastGlauber::Instance()->GetLastXY(xy);
+           Double_t i0i1[2];
+           AliFastGlauber::Instance()->GetSavedXY(xy);
+           AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
            xp = xy[0];
            yp = xy[1];
+           ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
        }
        
            ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
@@ -1269,7 +1313,10 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
     pdg[1] = jet2->GetPdgCode();    
     Bool_t   triggered = kFALSE;
 
-    if (fProcess == kPyJets) {
+    //
+    //TO BE CHECKED: If we call this method for kPyBeautyJets, we need it here
+    //
+    if (fProcess == kPyJets || fProcess == kPyBeautyJets) {
        Int_t njets = 0;
        Int_t ntrig = 0;
        Float_t jets[4][10];