]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythia.cxx
warning removal
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
index 43f03ec80e4b1484dfb0238afa64c0ba7d6d594f..6999a14204f02aad9745eab7db8572dd9eb7bd49 100644 (file)
@@ -52,6 +52,7 @@ ClassImp(AliGenPythia)
 AliGenPythia::AliGenPythia():
     AliGenMC(),
     fProcess(kPyCharm),          
+    fItune(-1),
     fStrucFunc(kCTEQ5L), 
     fKineBias(0.),
     fTrials(0),
@@ -75,11 +76,13 @@ AliGenPythia::AliGenPythia():
     fGinit(1),
     fGfinal(1),
     fHadronisation(1),
+    fPatchOmegaDalitz(0), 
     fNpartons(0),
     fReadFromFile(0),
     fQuench(0),
     fQhat(0.),
     fLength(0.),
+    fImpact(0.),
     fPtKick(1.),
     fFullEvent(kTRUE),
     fDecayer(new AliDecayerPythia()),
@@ -114,18 +117,21 @@ AliGenPythia::AliGenPythia():
     fTriggerEta(0.9),     
     fTriggerMultiplicity(0),
     fTriggerMultiplicityEta(0),
+    fTriggerMultiplicityPtMin(0),
     fCountMode(kCountAll),      
     fHeader(0),  
     fRL(0),      
-    fFileName(0),
+    fkFileName(0),
     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),
@@ -136,7 +142,6 @@ AliGenPythia::AliGenPythia():
 {
 // Default Constructor
   fEnergyCMS = 5500.;
-  SetNuclei(0,0);
   if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
 }
@@ -144,6 +149,7 @@ AliGenPythia::AliGenPythia():
 AliGenPythia::AliGenPythia(Int_t npart)
     :AliGenMC(npart),
      fProcess(kPyCharm),          
+     fItune(-1),
      fStrucFunc(kCTEQ5L), 
      fKineBias(0.),
      fTrials(0),
@@ -167,11 +173,13 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fGinit(kTRUE),
      fGfinal(kTRUE),
      fHadronisation(kTRUE),
+     fPatchOmegaDalitz(0), 
      fNpartons(0),
      fReadFromFile(kFALSE),
      fQuench(kFALSE),
      fQhat(0.),
      fLength(0.),
+     fImpact(0.),
      fPtKick(1.),
      fFullEvent(kTRUE),
      fDecayer(new AliDecayerPythia()),
@@ -206,18 +214,21 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fTriggerEta(0.9),     
      fTriggerMultiplicity(0),
      fTriggerMultiplicityEta(0),
+     fTriggerMultiplicityPtMin(0),
      fCountMode(kCountAll),      
      fHeader(0),  
      fRL(0),      
-     fFileName(0),
+     fkFileName(0),
      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),
@@ -236,7 +247,6 @@ AliGenPythia::AliGenPythia(Int_t npart)
     // Set random number generator 
     if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
-    SetNuclei(0,0);
  }
 
 AliGenPythia::~AliGenPythia()
@@ -360,14 +370,14 @@ void AliGenPythia::Init()
 
 
     if (fReadFromFile) {
-       fRL  =  AliRunLoader::Open(fFileName, "Partons");
+       fRL  =  AliRunLoader::Open(fkFileName, "Partons");
        fRL->LoadKinematics();
        fRL->LoadHeader();
     } else {
        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  
@@ -423,7 +433,12 @@ void AliGenPythia::Init()
        fParentSelect[0] =   431;
        fFlavorSelect    =   4; 
        break;
+    case kPyLambdacppMNR:
+       fParentSelect[0] =  4122;
+       fFlavorSelect    =   4; 
+       break;      
     case kPyBeauty:
+    case kPyBeautyJets:
     case kPyBeautyPbPbMNR:
     case kPyBeautypPbMNR:
     case kPyBeautyppMNR:
@@ -452,6 +467,7 @@ void AliGenPythia::Init()
        fParentSelect[0] = 443;
        break;
     case kPyMbDefault:
+    case kPyMbAtlasTuneMC09:
     case kPyMb:
     case kPyMbWithDirectPhoton:
     case kPyMbNonDiffr:
@@ -528,15 +544,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 +588,9 @@ 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;
+           fPythia->Qpygin0();
        } 
 //
 // Either produce new event or read partons from file
@@ -617,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");
-       Boost();
+       
+       if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
+       if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
+       
 //
 //
 //
@@ -648,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 &&
@@ -656,7 +679,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 +909,7 @@ 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) {
+     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
        TParticle* jet1 = (TParticle *) fParticles.At(6);
        TParticle* jet2 = (TParticle *) fParticles.At(7);
        if (!CheckTrigger(jet1, jet2)) {
@@ -901,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);
@@ -923,7 +947,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) {
       Int_t multiplicity = 0;
@@ -933,17 +982,15 @@ Int_t  AliGenPythia::GenerateMB()
        Int_t statusCode = iparticle->GetStatusCode();
        
        // Initial state particle
-       if (statusCode > 20)
+       if (statusCode != 1)
          continue;
-       
-       // skip quarks and gluons
-       Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
-       if (pdgCode <= 10 || pdgCode == 21)
-         continue;
-       
+       // eta cut
        if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
          continue;
-       
+       // pt cut
+       if (iparticle->Pt() < fTriggerMultiplicityPtMin) 
+           continue;
+
        TParticlePDG* pdgPart = iparticle->GetPDG();
        if (pdgPart && pdgPart->Charge() == 0)
          continue;
@@ -955,8 +1002,7 @@ Int_t  AliGenPythia::GenerateMB()
        delete [] pParent;
        return 0;
       }
-      
-      Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
+      Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
     }    
     
      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
@@ -1013,7 +1059,8 @@ 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) {
+
+    if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
       TParticle *partCheck;
       TParticle *mother;
       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
@@ -1059,6 +1106,7 @@ Int_t  AliGenPythia::GenerateMB()
          fProcess == kPyZ ||
          fProcess == kPyMbDefault ||
          fProcess == kPyMb ||
+         fProcess == kPyMbAtlasTuneMC09 ||
          fProcess == kPyMbWithDirectPhoton ||
          fProcess == kPyMbNonDiffr)  
         && (fCutOnChild == 1) ) {
@@ -1077,7 +1125,7 @@ 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)) {
+           ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
            nc++;
            if (ks == 1) trackIt = 1;
            Int_t ipa = iparticle->GetFirstMother()-1;
@@ -1179,14 +1227,15 @@ void AliGenPythia::MakeHeader()
 //
 // Event Vertex 
     fHeader->SetPrimaryVertex(fVertex);
-    
+    fHeader->SetInteractionTime(fEventTime);
 //
 // Number of primaries
     fHeader->SetNProduced(fNprimaries);
 //
 // 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 || fProcess == kPyBeautyppMNRwmi)
     {
        Int_t ntrig, njet;
        Float_t jets[4][10];
@@ -1241,7 +1290,7 @@ void AliGenPythia::MakeHeader()
            AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
            xp = xy[0];
            yp = xy[1];
-           ((AliGenPythiaEventHeader*) fHeader)->SetInMediumLength(2. * i0i1[1] / i0i1[0]);
+           ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
        }
        
            ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
@@ -1272,7 +1321,7 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
     pdg[1] = jet2->GetPdgCode();    
     Bool_t   triggered = kFALSE;
 
-    if (fProcess == kPyJets) {
+    if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi) {
        Int_t njets = 0;
        Int_t ntrig = 0;
        Float_t jets[4][10];
@@ -1582,8 +1631,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();
@@ -1591,21 +1640,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 
   
@@ -1617,46 +1666,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
-
-