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