Correct time for primary ionisation steps.

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

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 043b69663c536470187771faaeec706e498832e2..af5924524d73975f3e96b1c66a860df4143004f1 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -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),
@@ -870,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++) {
@@ -1112,6 +1153,7 @@ void AliGenPythia::MakeHeader()
 //  Pass header
 //
     AddHeader(fHeader);
+    fHeader = 0x0;
 }
 
 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
@@ -1369,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