PWGJE/EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx

   1 // $Id$
   2 //
   3 // Jet embedding task.
   4 //
   5 // Author: S.Aiola, C.Loizides
   6
   7 #include "AliJetEmbeddingFromGenTask.h"
   8
   9 #include <TClonesArray.h>
  10 #include <TFolder.h>
  11 #include <TLorentzVector.h>
  12 #include <TParticle.h>
  13 #include <TParticlePDG.h>
  14 #include <TRandom3.h>
  15 #include <TProfile.h>
  16 #include "AliAnalysisManager.h"
  17 #include "AliEMCALDigit.h"
  18 #include "AliEMCALGeometry.h"
  19 #include "AliEMCALRecPoint.h"
  20 #include "AliGenerator.h"
  21 #include "AliHeader.h"
  22 #include "AliLog.h"
  23 #include "AliPicoTrack.h"
  24 #include "AliRun.h"
  25 #include "AliRunLoader.h"
  26 #include "AliStack.h"
  27 #include "AliStack.h"
  28 #include "AliVCluster.h"
  29 #include "AliVEvent.h"
  30 #include "AliGenPythiaEventHeader.h"
  31 #include "AliStackPartonInfo.h"
  32
  33 ClassImp(AliJetEmbeddingFromGenTask)
  34
  35 //________________________________________________________________________
  36 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask() :
  37   AliJetModelBaseTask("AliJetEmbeddingFromGenTask"),
  38   fGen(0),
  39   fMassless(kFALSE),
  40   fChargedOnly(kFALSE),
  41   fHistPt(0),
  42   fHistEtaPhi(0),
  43   fHistTrials(0),
  44   fHistXsection(0),
  45   fHistPtHard(0)
  46 {
  47   // Default constructor.
  48   SetSuffix("EmbeddedFromGen");
  49 }
  50
  51 //________________________________________________________________________
  52 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask(const char *name, Bool_t drawqa) :
  53   AliJetModelBaseTask(name,drawqa),
  54   fGen(0),
  55   fMassless(kFALSE),
  56   fChargedOnly(kFALSE),
  57   fHistPt(0),
  58   fHistEtaPhi(0),
  59   fHistTrials(0),
  60   fHistXsection(0),
  61   fHistPtHard(0)
  62 {
  63   // Standard constructor.
  64   SetSuffix("EmbeddedFromGen");
  65 }
  66
  67 //________________________________________________________________________
  68 AliJetEmbeddingFromGenTask::~AliJetEmbeddingFromGenTask()
  69 {
  70   // Destructor
  71 }
  72
  73 //________________________________________________________________________
  74 void AliJetEmbeddingFromGenTask::UserCreateOutputObjects()
  75 {
  76   // Create user output.
  77
  78   if (!fQAhistos)
  79     return;
  80
  81   AliJetModelBaseTask::UserCreateOutputObjects();
  82
  83   fHistPt = new TH1F("fHistpt","fHistPt;#it{p}_{T};N",100,0.,100.);
  84   fOutput->Add(fHistPt);
  85
  86   fHistEtaPhi = new TH2F("fHistEtapHI","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
  87   fOutput->Add(fHistEtaPhi);
  88
  89   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
  90   fHistTrials->GetYaxis()->SetTitle("trials");
  91   fOutput->Add(fHistTrials);
  92
  93   fHistXsection = new TProfile("fHistXsection", "fHistXsection", 1, 0, 1);
  94   fHistXsection->GetYaxis()->SetTitle("xsection");
  95   fOutput->Add(fHistXsection);
  96
  97   fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", 500, 0., 500.);
  98   fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
  99   fHistPtHard->GetYaxis()->SetTitle("counts");
 100   fOutput->Add(fHistPtHard);
 101
 102   PostData(1, fOutput);
 103 }
 104
 105 //________________________________________________________________________
 106 Bool_t AliJetEmbeddingFromGenTask::ExecOnce()
 107 {
 108   // Exec only once.
 109
 110   if (!gAlice) {
 111     new AliRun("gAlice","The ALICE Off-line Simulation Framework");
 112     delete gRandom;
 113     gRandom = new TRandom3(0);
 114   }
 115
 116   TFolder *folder = new TFolder(GetName(),GetName());
 117   AliRunLoader *rl = new AliRunLoader(folder);
 118   rl->MakeHeader();
 119   rl->MakeStack();
 120   AliStack *stack = rl->Stack();
 121   fGen->SetStack(stack);
 122   fGen->Init();
 123
 124   if (!(InputEvent()->FindListObject(fTracksName))) {
 125     fOutTracks = new TClonesArray("AliPicoTrack", 1000);
 126     fOutTracks->SetName(fTracksName);
 127     InputEvent()->AddObject(fOutTracks);
 128     fNTracks = 0;
 129   }
 130
 131   if (!(InputEvent()->FindListObject(fPartonInfoName))) {
 132   fStackPartonInfo = new AliStackPartonInfo("PartonsInfo");
 133   fStackPartonInfo->SetName(fPartonInfoName);
 134   InputEvent()->AddObject(fStackPartonInfo);
 135   }
 136   return kTRUE;
 137 }
 138
 139 //________________________________________________________________________
 140 void AliJetEmbeddingFromGenTask::Run()
 141 {
 142   // Embed particles.
 143
 144   if (fCopyArray)
 145     CopyTracks();
 146
 147   AliStack *stack = fGen->GetStack();
 148   stack->Reset();
 149   fGen->Generate();
 150   const Int_t nprim = stack->GetNprimary();
 151   TParticle *part6 = stack->Particle(6);
 152   TParticle *part7 = stack->Particle(7);
 153
 154   fStackPartonInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
 155   fStackPartonInfo->SetPartonPt6(part6->Pt());
 156   fStackPartonInfo->SetPartonEta6(part6->Eta());
 157   fStackPartonInfo->SetPartonPhi6(part6->Phi());
 158
 159   fStackPartonInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
 160   fStackPartonInfo->SetPartonPt7(part7->Pt());
 161   fStackPartonInfo->SetPartonEta7(part7->Eta());
 162   fStackPartonInfo->SetPartonPhi7(part7->Phi());
 163
 164   for (Int_t i=0;i<nprim;++i) {
 165     if (!stack->IsPhysicalPrimary(i))
 166       continue;
 167     TParticle *part = stack->Particle(i);
 168     TParticlePDG *pdg = part->GetPDG(1);
 169     if (!pdg)
 170       continue;
 171     Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
 172     if (fChargedOnly && c==0)  continue;
 173     Double_t pt = part->Pt();
 174     Double_t eta = part->Eta();
 175     Double_t phi = part->Phi();
 176     if (eta<fEtaMin)
 177       continue;
 178     if (eta>fEtaMax)
 179       continue;
 180     if (phi<fPhiMin)
 181       continue;
 182     if (phi>fPhiMax)
 183       continue;
 184     if (pt<fPtMin)
 185       continue;
 186     if (pt>fPtMax)
 187       continue;
 188     Double_t mass = part->GetMass();
 189     if(fMassless) mass = 0.;
 190     fHistPt->Fill(pt);
 191     fHistEtaPhi->Fill(eta,phi);
 192     AddTrack(pt, eta, phi,0,0,0,0,0,0,c,mass);
 193   }
 194
 195   FillPythiaHistograms();
 196 }
 197
 198 //________________________________________________________________________
 199 void AliJetEmbeddingFromGenTask::FillPythiaHistograms() {
 200   //Get PYTHIA info: pt-hard, x-section, trials
 201
 202   if (!fQAhistos)
 203     return;
 204
 205   AliRunLoader *rl = AliRunLoader::Instance();
 206   AliGenPythiaEventHeader *genPH = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
 207   if(genPH) {
 208     Float_t xsec = genPH->GetXsection();
 209     Int_t trials = genPH->Trials();
 210     Float_t pthard = genPH->GetPtHard();
 211
 212     fHistXsection->Fill(0.5,xsec);
 213     fHistTrials->Fill(0.5,trials);
 214     fHistPtHard->Fill(pthard);
 215   }
 216 }