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
  32 ClassImp(AliJetEmbeddingFromGenTask)
  33
  34 //________________________________________________________________________
  35 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask() :
  36   AliJetModelBaseTask("AliJetEmbeddingFromGenTask"),
  37   fGen(0),
  38   fMassless(kFALSE),
  39   fHistPt(0),
  40   fHistTrials(0),
  41   fHistXsection(0),
  42   fHistPtHard(0)
  43 {
  44   // Default constructor.
  45   SetSuffix("EmbeddedFromGen");
  46 }
  47
  48 //________________________________________________________________________
  49 AliJetEmbeddingFromGenTask::AliJetEmbeddingFromGenTask(const char *name, Bool_t drawqa) :
  50   AliJetModelBaseTask(name,drawqa),
  51   fGen(0),
  52   fMassless(kFALSE),
  53   fHistPt(0),
  54   fHistTrials(0),
  55   fHistXsection(0),
  56   fHistPtHard(0)
  57 {
  58   // Standard constructor.
  59   SetSuffix("EmbeddedFromGen");
  60 }
  61
  62 //________________________________________________________________________
  63 AliJetEmbeddingFromGenTask::~AliJetEmbeddingFromGenTask()
  64 {
  65   // Destructor
  66 }
  67
  68 //________________________________________________________________________
  69 void AliJetEmbeddingFromGenTask::UserCreateOutputObjects()
  70 {
  71   // Create user output.
  72
  73   if (!fQAhistos)
  74     return;
  75
  76   AliJetModelBaseTask::UserCreateOutputObjects();
  77
  78   fHistPt = new TH1F("fHistpt","fHistPt;#it{p}_{T};N",100,0.,100.);
  79   fOutput->Add(fHistPt);
  80
  81   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
  82   fHistTrials->GetYaxis()->SetTitle("trials");
  83   fOutput->Add(fHistTrials);
  84
  85   fHistXsection = new TProfile("fHistXsection", "fHistXsection", 1, 0, 1);
  86   fHistXsection->GetYaxis()->SetTitle("xsection");
  87   fOutput->Add(fHistXsection);
  88
  89   fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", 500, 0., 500.);
  90   fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
  91   fHistPtHard->GetYaxis()->SetTitle("counts");
  92   fOutput->Add(fHistPtHard);
  93
  94   PostData(1, fOutput);
  95 }
  96
  97 //________________________________________________________________________
  98 Bool_t AliJetEmbeddingFromGenTask::ExecOnce()
  99 {
 100   // Exec only once.
 101
 102   if (!gAlice) {
 103     new AliRun("gAlice","The ALICE Off-line Simulation Framework");
 104     delete gRandom;
 105     gRandom = new TRandom3(0);
 106   }
 107
 108   TFolder *folder = new TFolder(GetName(),GetName());
 109   AliRunLoader *rl = new AliRunLoader(folder);
 110   rl->MakeHeader();
 111   rl->MakeStack();
 112   AliStack *stack = rl->Stack();
 113   fGen->SetStack(stack);
 114   fGen->Init();
 115
 116   if (!(InputEvent()->FindListObject(fTracksName))) {
 117     fOutTracks = new TClonesArray("AliPicoTrack", 1000);
 118     fOutTracks->SetName(fTracksName);
 119     InputEvent()->AddObject(fOutTracks);
 120     fNTracks = 0;
 121   }
 122
 123   return kTRUE;
 124 }
 125
 126 //________________________________________________________________________
 127 void AliJetEmbeddingFromGenTask::Run()
 128 {
 129   // Embed particles.
 130
 131   if (fCopyArray)
 132     CopyTracks();
 133
 134   AliStack *stack = fGen->GetStack();
 135   stack->Reset();
 136   fGen->Generate();
 137   const Int_t nprim = stack->GetNprimary();
 138   for (Int_t i=0;i<nprim;++i) {
 139     if (!stack->IsPhysicalPrimary(i))
 140       continue;
 141     TParticle *part = stack->Particle(i);
 142     TParticlePDG *pdg = part->GetPDG(1);
 143     if (!pdg)
 144       continue;
 145     Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
 146     if (c==0)
 147       continue;
 148     Double_t pt = part->Pt();
 149     Double_t eta = part->Eta();
 150     Double_t phi = part->Phi();
 151     if (eta<fEtaMin)
 152       continue;
 153     if (eta>fEtaMax)
 154       continue;
 155     if (phi<fPhiMin)
 156       continue;
 157     if (phi>fPhiMax)
 158       continue;
 159     if (pt<fPtMin)
 160       continue;
 161     if (pt>fPtMax)
 162       continue;
 163     Double_t mass = part->GetMass();
 164     if(fMassless) mass = 0.;
 165     fHistPt->Fill(pt);
 166     AddTrack(pt, eta, phi,0,0,0,0,0,0,c,mass);
 167   }
 168
 169   FillPythiaHistograms();
 170 }
 171
 172 //________________________________________________________________________
 173 void AliJetEmbeddingFromGenTask::FillPythiaHistograms() {
 174   //Get PYTHIA info: pt-hard, x-section, trials
 175
 176   if (!fQAhistos)
 177     return;
 178
 179   AliRunLoader *rl = AliRunLoader::Instance();
 180   AliGenPythiaEventHeader *genPH = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
 181   if(genPH) {
 182     Float_t xsec = genPH->GetXsection();
 183     Int_t trials = genPH->Trials();
 184     Float_t pthard = genPH->GetPtHard();
 185
 186     fHistXsection->Fill(0.5,xsec);
 187     fHistTrials->Fill(0.5,trials);
 188     fHistPtHard->Fill(pthard);
 189   }
 190 }