Completing the setup for TTherminator model.

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

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index b1ef23a6d9132779a2e245714d6bec421e576b45..649e123b03023e6ef1001b6e09582b5ecc154cd0 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),
@@ -335,7 +341,8 @@ void AliGenPythia::Init()
 //  pt - kick
     if (fPtKick > 0.) {
        fPythia->SetMSTP(91,1);
-       fPythia->SetPARP(91,fPtKick);
+       fPythia->SetPARP(91,fPtKick);   
+       fPythia->SetPARP(93, 4. * fPtKick);
     } else {
        fPythia->SetMSTP(91,0);
     }
@@ -348,18 +355,23 @@ void AliGenPythia::Init()
     } else {
        fRL = 0x0;
     }
+ //
+    fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
+    //  Forward Paramters to the AliPythia object
+    fDecayer->SetForceDecay(fForceDecay);    
 // 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();
     }
- //
-    fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
-    //  Forward Paramters to the AliPythia object
-    fDecayer->SetForceDecay(fForceDecay);    
+
     fDecayer->Init();
 
+
 //  Parent and Children Selection
     switch (fProcess) 
     {
@@ -864,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++) {
@@ -1106,6 +1153,7 @@ void AliGenPythia::MakeHeader()
 //  Pass header
 //
     AddHeader(fHeader);
+    fHeader = 0x0;
 }
 
 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
@@ -1363,6 +1411,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