Introduction of jet-by-jet correction for detector effects

[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliJetEmbeddingFromGenTask.cxx

diff --git a/PWGJE/EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx
index 5e5492580d579b5436dd9dc0ecbb1d9a6e705a69..29fa272c66decab9302e17d37a0ec384a7c28348 100644 (file)
--- a/PWGJE/EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx
+++ b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx
@@ -24,17 +24,21 @@
 #include "AliRun.h"
 #include "AliRunLoader.h"
 #include "AliStack.h"
-#include "AliStack.h"
 #include "AliVCluster.h"
 #include "AliVEvent.h"
 #include "AliGenPythiaEventHeader.h"
-
+#include "AliPythiaInfo.h"
+#include "AliPythiaRndm.h"
 ClassImp(AliJetEmbeddingFromGenTask)
 
 //________________________________________________________________________
 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask() : 
   AliJetModelBaseTask("AliJetEmbeddingFromGenTask"),
   fGen(0),
+  fMassless(kFALSE),
+  fChargedOnly(kFALSE),
+  fHistPt(0),
+  fHistEtaPhi(0),
   fHistTrials(0),
   fHistXsection(0),
   fHistPtHard(0)
@@ -47,6 +51,10 @@ AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask() :
 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask(const char *name, Bool_t drawqa) :
   AliJetModelBaseTask(name,drawqa),
   fGen(0),
+  fMassless(kFALSE),
+  fChargedOnly(kFALSE),
+  fHistPt(0),
+  fHistEtaPhi(0),
   fHistTrials(0),
   fHistXsection(0),
   fHistPtHard(0)
@@ -71,6 +79,12 @@ void AliJetEmbeddingFromGenTask::UserCreateOutputObjects()
 
   AliJetModelBaseTask::UserCreateOutputObjects();
 
+  fHistPt = new TH1F("fHistpt","fHistPt;#it{p}_{T};N",100,0.,100.);
+  fOutput->Add(fHistPt);
+
+  fHistEtaPhi = new TH2F("fHistEtapHI","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
+  fOutput->Add(fHistEtaPhi);
+
   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
   fHistTrials->GetYaxis()->SetTitle("trials");
   fOutput->Add(fHistTrials);
@@ -100,6 +114,7 @@ Bool_t AliJetEmbeddingFromGenTask::ExecOnce()
 
   TFolder *folder = new TFolder(GetName(),GetName());
   AliRunLoader *rl = new AliRunLoader(folder);
+  gAlice->SetRunLoader(rl);
   rl->MakeHeader();
   rl->MakeStack();
   AliStack *stack = rl->Stack();
@@ -112,7 +127,14 @@ Bool_t AliJetEmbeddingFromGenTask::ExecOnce()
     InputEvent()->AddObject(fOutTracks);
     fNTracks = 0;
   }
-  
+
+  if(!fPythiaInfoName.IsNull()) {
+    if (!(InputEvent()->FindListObject(fPythiaInfoName))) {
+      fPythiaInfo = new AliPythiaInfo("PythiaInfo");
+      fPythiaInfo->SetName(fPythiaInfoName);
+      InputEvent()->AddObject(fPythiaInfo);
+    }
+  }
   return kTRUE;
 }
 
@@ -123,11 +145,29 @@ void AliJetEmbeddingFromGenTask::Run()
 
   if (fCopyArray) 
     CopyTracks();
-
+  AliPythiaRndm::SetPythiaRandom(new TRandom3());
+  AliPythiaRndm::GetPythiaRandom()->SetSeed(clock()+gSystem->GetPid());
   AliStack *stack = fGen->GetStack();
   stack->Reset();
   fGen->Generate();
   const Int_t nprim = stack->GetNprimary();
+  // reject if partons are missing from stack for some reason
+  if(nprim < 8) return;
+  if(fPythiaInfo) {
+    TParticle *part6 = stack->Particle(6);
+    TParticle *part7 = stack->Particle(7);
+    
+    fPythiaInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
+    fPythiaInfo->SetPartonPt6(part6->Pt());
+    fPythiaInfo->SetPartonEta6(part6->Eta());
+    fPythiaInfo->SetPartonPhi6(part6->Phi());
+    
+    fPythiaInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
+    fPythiaInfo->SetPartonPt7(part7->Pt());
+    fPythiaInfo->SetPartonEta7(part7->Eta());
+    fPythiaInfo->SetPartonPhi7(part7->Phi());
+  }
+
   for (Int_t i=0;i<nprim;++i) {
     if (!stack->IsPhysicalPrimary(i))
       continue;
@@ -136,8 +176,7 @@ void AliJetEmbeddingFromGenTask::Run()
     if (!pdg) 
       continue;
     Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
-    if (c==0) 
-      continue;
+    if (fChargedOnly && c==0)  continue;
     Double_t pt = part->Pt();
     Double_t eta = part->Eta();
     Double_t phi = part->Phi();
@@ -153,7 +192,11 @@ void AliJetEmbeddingFromGenTask::Run()
       continue;
     if (pt>fPtMax)
       continue;
-    AddTrack(pt, eta, phi,0,0,0,0,0,0,c,part->GetMass());
+    Double_t mass = part->GetMass();
+    if(fMassless) mass = 0.;
+    fHistPt->Fill(pt);
+    fHistEtaPhi->Fill(eta,phi);
+    AddTrack(pt, eta, phi,0,0,0,0,0,0,c,mass);
   }
 
   FillPythiaHistograms();
@@ -169,11 +212,11 @@ void AliJetEmbeddingFromGenTask::FillPythiaHistograms() {
   AliRunLoader *rl = AliRunLoader::Instance();
   AliGenPythiaEventHeader *genPH = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
   if(genPH) {
-    // Printf("found pythia event header. pThard: %f Trials: %d xsec: %f",genPH->GetPtHard(),genPH->Trials(),genPH->GetXsection());
     Float_t xsec = genPH->GetXsection();
     Int_t trials = genPH->Trials();
     Float_t pthard = genPH->GetPtHard();
-
+    Float_t ptWeight=genPH->EventWeight(); 
+    if(fPythiaInfo) fPythiaInfo->SetPythiaEventWeight(ptWeight);
     fHistXsection->Fill(0.5,xsec);
     fHistTrials->Fill(0.5,trials);
     fHistPtHard->Fill(pthard);