]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythia.cxx
Fix to take into account the fact that Pythia6 from ROOT
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
index 3c19f929cdee75f5b666f651dae0aac0c489a583..daaa22261a46f363ed7cf3b07771d7b6432567b2 100644 (file)
 
 ClassImp(AliGenPythia)
 
-AliGenPythia::AliGenPythia()
-                 :AliGenMC()
+
+AliGenPythia::AliGenPythia():
+    AliGenMC(),
+    fProcess(kPyCharm),          
+    fStrucFunc(kCTEQ5L), 
+    fEnergyCMS(5500.),
+    fKineBias(0.),
+    fTrials(0),
+    fTrialsRun(0),
+    fQ(0.),
+    fX1(0.),
+    fX2(0.),
+    fEventTime(0.),
+    fInteractionRate(0.),
+    fTimeWindow(0.),
+    fCurSubEvent(0),
+    fEventsTime(0),
+    fNev(0),
+    fFlavorSelect(0),
+    fXsection(0.),
+    fPythia(0),
+    fPtHardMin(0.),
+    fPtHardMax(1.e4),
+    fYHardMin(-1.e10),
+    fYHardMax(1.e10),
+    fGinit(1),
+    fGfinal(1),
+    fHadronisation(1),
+    fNpartons(0),
+    fReadFromFile(0),
+    fQuench(0),
+    fPtKick(1.),
+    fFullEvent(kTRUE),
+    fDecayer(new AliDecayerPythia()),
+    fDebugEventFirst(-1),
+    fDebugEventLast(-1),
+    fEtMinJet(0.),      
+    fEtMaxJet(1.e4),      
+    fEtaMinJet(-20.),     
+    fEtaMaxJet(20.),     
+    fPhiMinJet(0.),     
+    fPhiMaxJet(2.* TMath::Pi()),     
+    fJetReconstruction(kCell),
+    fEtaMinGamma(-20.),      
+    fEtaMaxGamma(20.),      
+    fPhiMinGamma(0.),      
+    fPhiMaxGamma(2. * TMath::Pi()),      
+    fPycellEtaMax(2.),     
+    fPycellNEta(274),       
+    fPycellNPhi(432),       
+    fPycellThreshold(0.),  
+    fPycellEtSeed(4.),     
+    fPycellMinEtJet(10.),  
+    fPycellMaxRadius(1.), 
+    fStackFillOpt(kFlavorSelection),   
+    fFeedDownOpt(kTRUE),    
+    fFragmentation(kTRUE),
+    fSetNuclei(kFALSE),
+    fNewMIS(kFALSE),   
+    fHFoff(kFALSE),    
+    fTriggerParticle(0),
+    fTriggerEta(0.9),     
+    fCountMode(kCountAll),      
+    fHeader(0),  
+    fRL(0),      
+    fFileName(0),
+    fFragPhotonInCalo(kFALSE),
+    fPi0InCalo(kFALSE) ,
+    fCheckEMCAL(kFALSE),
+    fCheckPHOS(kFALSE),
+    fFragPhotonOrPi0MinPt(0), 
+    fPHOSMinPhi(219.),
+    fPHOSMaxPhi(321.),
+    fPHOSEta(0.13),
+    fEMCALMinPhi(79.),
+    fEMCALMaxPhi(191.),
+    fEMCALEta(0.71)
+
 {
 // Default Constructor
-  fParticles = 0;
-  fPythia    = 0;
-  fHeader = 0;
-  fReadFromFile = 0;
-  fEventTime = 0.;
-  fInteractionRate = 0.;
-  fTimeWindow = 0.;
-  fEventsTime = 0;
-  fCurSubEvent = 0;
-  fDecayer   = new AliDecayerPythia();
-  SetEventListRange();
-  SetJetPhiRange();
-  SetJetEtaRange();
-  SetJetEtRange();
-  SetGammaPhiRange();
-  SetGammaEtaRange();
-  SetPtKick();
-  SetQuench();
-  SetHadronisation();  
-  SetTriggerParticle();
   SetNuclei(0,0);
-  fSetNuclei = kFALSE;
-  fNewMIS    = kFALSE;
-  fHFoff     = kFALSE;
   if (!AliPythiaRndm::GetPythiaRandom()) 
-    AliPythiaRndm::SetPythiaRandom(GetRandom());
+      AliPythiaRndm::SetPythiaRandom(GetRandom());
 }
 
 AliGenPythia::AliGenPythia(Int_t npart)
-                 :AliGenMC(npart)
+    :AliGenMC(npart),
+     fProcess(kPyCharm),          
+     fStrucFunc(kCTEQ5L), 
+     fEnergyCMS(5500.),
+     fKineBias(0.),
+     fTrials(0),
+     fTrialsRun(0),
+     fQ(0.),
+     fX1(0.),
+     fX2(0.),
+     fEventTime(0.),
+     fInteractionRate(0.),
+     fTimeWindow(0.),
+     fCurSubEvent(0),
+     fEventsTime(0),
+     fNev(0),
+     fFlavorSelect(0),
+     fXsection(0.),
+     fPythia(0),
+     fPtHardMin(0.),
+     fPtHardMax(1.e4),
+     fYHardMin(-1.e10),
+     fYHardMax(1.e10),
+     fGinit(kTRUE),
+     fGfinal(kTRUE),
+     fHadronisation(kTRUE),
+     fNpartons(0),
+     fReadFromFile(kFALSE),
+     fQuench(kFALSE),
+     fPtKick(1.),
+     fFullEvent(kTRUE),
+     fDecayer(new AliDecayerPythia()),
+     fDebugEventFirst(-1),
+     fDebugEventLast(-1),
+     fEtMinJet(0.),      
+     fEtMaxJet(1.e4),      
+     fEtaMinJet(-20.),     
+     fEtaMaxJet(20.),     
+     fPhiMinJet(0.),     
+     fPhiMaxJet(2.* TMath::Pi()),     
+     fJetReconstruction(kCell),
+     fEtaMinGamma(-20.),      
+     fEtaMaxGamma(20.),      
+     fPhiMinGamma(0.),      
+     fPhiMaxGamma(2. * TMath::Pi()),      
+     fPycellEtaMax(2.),     
+     fPycellNEta(274),       
+     fPycellNPhi(432),       
+     fPycellThreshold(0.),  
+     fPycellEtSeed(4.),     
+     fPycellMinEtJet(10.),  
+     fPycellMaxRadius(1.), 
+     fStackFillOpt(kFlavorSelection),   
+     fFeedDownOpt(kTRUE),    
+     fFragmentation(kTRUE),
+     fSetNuclei(kFALSE),
+     fNewMIS(kFALSE),   
+     fHFoff(kFALSE),    
+     fTriggerParticle(0),
+     fTriggerEta(0.9),     
+     fCountMode(kCountAll),      
+     fHeader(0),  
+     fRL(0),      
+     fFileName(0),
+     fFragPhotonInCalo(kFALSE),
+     fPi0InCalo(kFALSE) ,
+     fCheckEMCAL(kFALSE),
+     fCheckPHOS(kFALSE),
+     fFragPhotonOrPi0MinPt(0),
+     fPHOSMinPhi(219.),
+     fPHOSMaxPhi(321.),
+     fPHOSEta(0.13),
+     fEMCALMinPhi(79.),
+     fEMCALMaxPhi(191.),
+     fEMCALEta(0.71)
 {
 // default charm production at 5. 5 TeV
 // semimuonic decay
@@ -86,59 +213,13 @@ AliGenPythia::AliGenPythia(Int_t npart)
 //
     fName = "Pythia";
     fTitle= "Particle Generator using PYTHIA";
-    fXsection  = 0.;
-    fReadFromFile = 0;
-    fEventTime = 0.;
-    fInteractionRate = 0.;
-    fTimeWindow = 0.;
-    fEventsTime = 0;
-    fCurSubEvent = 0;
-    SetProcess();
-    SetStrucFunc();
     SetForceDecay();
-    SetPtHard();
-    SetYHard();
-    SetEnergyCMS();
-    fDecayer = new AliDecayerPythia();
     // Set random number generator 
     if (!AliPythiaRndm::GetPythiaRandom()) 
       AliPythiaRndm::SetPythiaRandom(GetRandom());
-    fFlavorSelect   = 0;
-    // Produced particles  
     fParticles = new TClonesArray("TParticle",1000);
-    fHeader = 0;
-    SetEventListRange();
-    SetJetPhiRange();
-    SetJetEtaRange();
-    SetJetEtRange();
-    SetGammaPhiRange();
-    SetGammaEtaRange();
-    SetJetReconstructionMode();
-    SetQuench();
-    SetHadronisation();
-    SetPtKick();
-    SetTriggerParticle();
     SetNuclei(0,0);
-    // Options determining what to keep in the stack (Heavy flavour generation)
-    fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
-    fFeedDownOpt = kTRUE;             // allow feed down from higher family
-    // Fragmentation on/off
-    fFragmentation = kTRUE;
-    // Default counting mode
-    fCountMode = kCountAll;
-    // Pycel
-    SetPycellParameters();
-    fSetNuclei = kFALSE;
-    fNewMIS    = kFALSE;
-    fHFoff     = kFALSE;
-}
-
-AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
-    :AliGenMC(Pythia)
-{
-// copy constructor
-    Pythia.Copy(*this);
-}
+ }
 
 AliGenPythia::~AliGenPythia()
 {
@@ -348,6 +429,7 @@ void AliGenPythia::Init()
     case kPyMbMSEL1:
     case kPyJets:
     case kPyDirectGamma:
+    case kPyLhwgMb:    
        break;
     case kPyW:
     case kPyZ:
@@ -392,6 +474,18 @@ void AliGenPythia::Init()
     if (fQuench) {
        fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
     }
+    fPythia->SetPARJ(200, 0.0);
+
+    if (fQuench == 3) {
+       // Nestor's change of the splittings
+       fPythia->SetPARJ(200, 0.8);
+       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
+    }
 }
 
 void AliGenPythia::Generate()
@@ -452,7 +546,10 @@ void AliGenPythia::Generate()
            fPythia->Quench();
        } else if (fQuench == 2){
            fPythia->Pyquen(208., 0, 0.);
+       } else if (fQuench == 3) {
+           // Quenching is via multiplicative correction of the splittings
        }
+       
 //
 // Switch hadronisation on
 //
@@ -730,6 +827,36 @@ Int_t  AliGenPythia::GenerateMB()
        }
     }
 
+    // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
+    if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
+
+      Bool_t ok = kFALSE;
+
+      Int_t pdg  = 0; 
+      if (fFragPhotonInCalo) pdg = 22   ; // Photon
+      else if (fPi0InCalo) pdg = 111 ; // Pi0
+
+      for (i=0; i< np; i++) {
+       TParticle* iparticle = (TParticle *) fParticles->At(i);
+       if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg && 
+          iparticle->Pt() > fFragPhotonOrPi0MinPt){
+         Int_t imother = iparticle->GetFirstMother() - 1;
+         TParticle* pmother = (TParticle *) fParticles->At(imother);
+         if(pdg == 111 || 
+            (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
+           {
+             Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+             Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax         
+             if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
+                (fCheckPHOS    && IsInPHOS(phi,eta)) )
+               ok =kTRUE;
+           }
+       }
+      }
+      if(!ok)
+       return 0;
+    }
+    
     if (fTriggerParticle) {
        Bool_t triggered = kFALSE;
        for (i = 0; i < np; i++) {
@@ -1000,7 +1127,7 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 //
        GetJets(njets, ntrig, jets);
        
-       if (ntrig) triggered = kTRUE;
+       if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
 //
     } else {
        Int_t ij = 0;
@@ -1055,14 +1182,6 @@ Bool_t AliGenPythia::CheckKinematicsOnChild(){
     return checking;
 }
 
-         
-AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
-{
-// Assignment operator
-    rhs.Copy(*this);
-    return *this;
-}
-
 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
 {
 //
@@ -1216,6 +1335,31 @@ void AliGenPythia::GetSubEventTime()
   return;
 }
 
+Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
+{
+  // Is particle in EMCAL acceptance? 
+  // phi in degrees, etamin=-etamax
+  if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi && 
+     eta < fEMCALEta  ) 
+    return kTRUE;
+  else 
+    return kFALSE;
+}
+
+Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
+{
+  // Is particle in PHOS acceptance? 
+  // Acceptance slightly larger considered.
+  // phi in degrees, etamin=-etamax
+  if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi && 
+     eta < fPHOSEta  ) 
+    return kTRUE;
+  else 
+    return kFALSE;
+}
+
+
+
 #ifdef never
 void AliGenPythia::Streamer(TBuffer &R__b)
 {