Clean-up in includes.

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

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 8b2a360d474230f045c7f3e17e49b4f6207202fe..d8238ce621d391878dd4dd68a1dead3cc2983b2f 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -42,7 +42,7 @@
 #include "AliStack.h"
 #include "AliRunLoader.h"
 #include "AliMC.h"
-#include "pyquenCommon.h"
+#include "PyquenCommon.h"
 
 ClassImp(AliGenPythia)
 
@@ -114,9 +114,12 @@ AliGenPythia::AliGenPythia():
     fFileName(0),
     fFragPhotonInCalo(kFALSE),
     fPi0InCalo(kFALSE) ,
+    fPhotonInCalo(kFALSE),
     fCheckEMCAL(kFALSE),
     fCheckPHOS(kFALSE),
+    fCheckPHOSeta(kFALSE),
     fFragPhotonOrPi0MinPt(0), 
+    fPhotonMinPt(0), 
     fPHOSMinPhi(219.),
     fPHOSMaxPhi(321.),
     fPHOSEta(0.13),
@@ -198,9 +201,12 @@ AliGenPythia::AliGenPythia(Int_t npart)
      fFileName(0),
      fFragPhotonInCalo(kFALSE),
      fPi0InCalo(kFALSE) ,
+     fPhotonInCalo(kFALSE),
      fCheckEMCAL(kFALSE),
      fCheckPHOS(kFALSE),
+     fCheckPHOSeta(kFALSE),
      fFragPhotonOrPi0MinPt(0),
+     fPhotonMinPt(0),
      fPHOSMinPhi(219.),
      fPHOSMaxPhi(321.),
      fPHOSEta(0.13),
@@ -353,14 +359,18 @@ void AliGenPythia::Init()
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
     //  Forward Paramters to the AliPythia object
     fDecayer->SetForceDecay(fForceDecay);    
-    fDecayer->Init();
 // Switch off Heavy Flavors on request  
     if (fHFoff) {
+       // Maximum number of quark flavours used in pdf 
        fPythia->SetMSTP(58, 3);
+       // Maximum number of flavors that can be used in showers
        fPythia->SetMSTJ(45, 3);        
-       for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
+       // Switch off g->QQbar splitting in decay table
+       ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
     }
 
+    fDecayer->Init();
+
 
 //  Parent and Children Selection
     switch (fProcess) 
@@ -866,6 +876,41 @@ Int_t  AliGenPythia::GenerateMB()
        return 0;
     }
     
+    
+     // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
+    if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
+
+      Bool_t okd = kFALSE;
+
+      Int_t pdg  = 22; 
+      Int_t iphcand = -1;
+      for (i=0; i< np; i++) {
+        TParticle* iparticle = (TParticle *) fParticles->At(i);
+        Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+        Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax 
+        
+        if(iparticle->GetStatusCode() == 1 
+           && iparticle->GetPdgCode() == pdg   
+           && iparticle->Pt() > fPhotonMinPt    
+           && eta > fPHOSEta){                 
+           
+           // first check if the photon is in PHOS phi
+           if(IsInPHOS(phi,eta)){ 
+               okd = kTRUE;
+               break;
+           } 
+           if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
+            
+        }
+      }
+      
+      if(!okd && iphcand != -1) // execute rotation in phi 
+          RotatePhi(iphcand,okd);
+      
+      if(!okd)
+       return 0;
+    }
+    
     if (fTriggerParticle) {
        Bool_t triggered = kFALSE;
        for (i = 0; i < np; i++) {
@@ -1365,6 +1410,56 @@ Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
     return kFALSE;
 }
 
+void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
+{
+  //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi 
+  Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
+  Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
+  Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
+  
+  //calculate deltaphi
+  TParticle* ph = (TParticle *) fParticles->At(iphcand);
+  Double_t phphi = ph->Phi();
+  Double_t deltaphi = phiPHOS - phphi;
+
+  
+  
+  //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;
+  for(Int_t i=0; i< np; i++) {
+      TParticle* iparticle = (TParticle *) fParticles->At(i);
+      oldphi = iparticle->Phi();
+      newphi = oldphi + deltaphi;
+      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
+      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
+      
+      // apply rotation 
+      iparticle->SetProductionVertex(newVx, newVy, Vz, time);
+      iparticle->SetMomentum(newPx, newPy, Pz, e);
+      
+  } //end particle loop 
+  
+   // now let's check that we put correctly the candidate photon in PHOS
+   Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
+   Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax 
+   if(IsInPHOS(phi,eta)) 
+      okdd = kTRUE;
+}
 
 
 #ifdef never