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   gAlice->SetRunLoader(rl);
 119   rl->MakeHeader();
 120   rl->MakeStack();
 121   AliStack *stack = rl->Stack();
 122   fGen->SetStack(stack);
 123   fGen->Init();
 124
 125   if (!(InputEvent()->FindListObject(fTracksName))) {
 126     fOutTracks = new TClonesArray("AliPicoTrack", 1000);
 127     fOutTracks->SetName(fTracksName);
 128     InputEvent()->AddObject(fOutTracks);
 129     fNTracks = 0;
 130   }
 131
 132   if (!(InputEvent()->FindListObject(fPartonInfoName))) {
 133   fStackPartonInfo = new AliStackPartonInfo("PartonsInfo");
 134   fStackPartonInfo->SetName(fPartonInfoName);
 135   InputEvent()->AddObject(fStackPartonInfo);
 136   }
 137   return kTRUE;
 138 }
 139
 140 //________________________________________________________________________
 141 void AliJetEmbeddingFromGenTask::Run()
 142 {
 143   // Embed particles.
 144
 145   if (fCopyArray)
 146     CopyTracks();
 147
 148   AliStack *stack = fGen->GetStack();
 149   stack->Reset();
 150   fGen->Generate();
 151   const Int_t nprim = stack->GetNprimary();
 152   // reject if partons are missing from stack for some reason
 153   if(nprim < 8) return;
 154   TParticle *part6 = stack->Particle(6);
 155   TParticle *part7 = stack->Particle(7);
 156
 157   fStackPartonInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
 158   fStackPartonInfo->SetPartonPt6(part6->Pt());
 159   fStackPartonInfo->SetPartonEta6(part6->Eta());
 160   fStackPartonInfo->SetPartonPhi6(part6->Phi());
 161
 162   fStackPartonInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
 163   fStackPartonInfo->SetPartonPt7(part7->Pt());
 164   fStackPartonInfo->SetPartonEta7(part7->Eta());
 165   fStackPartonInfo->SetPartonPhi7(part7->Phi());
 166
 167   for (Int_t i=0;i<nprim;++i) {
 168     if (!stack->IsPhysicalPrimary(i))
 169       continue;
 170     TParticle *part = stack->Particle(i);
 171     TParticlePDG *pdg = part->GetPDG(1);
 172     if (!pdg)
 173       continue;
 174     Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
 175     if (fChargedOnly && c==0)  continue;
 176     Double_t pt = part->Pt();
 177     Double_t eta = part->Eta();
 178     Double_t phi = part->Phi();
 179     if (eta<fEtaMin)
 180       continue;
 181     if (eta>fEtaMax)
 182       continue;
 183     if (phi<fPhiMin)
 184       continue;
 185     if (phi>fPhiMax)
 186       continue;
 187     if (pt<fPtMin)
 188       continue;
 189     if (pt>fPtMax)
 190       continue;
 191     Double_t mass = part->GetMass();
 192     if(fMassless) mass = 0.;
 193     fHistPt->Fill(pt);
 194     fHistEtaPhi->Fill(eta,phi);
 195     AddTrack(pt, eta, phi,0,0,0,0,0,0,c,mass);
 196   }
 197
 198   FillPythiaHistograms();
 199 }
 200
 201 //________________________________________________________________________
 202 void AliJetEmbeddingFromGenTask::FillPythiaHistograms() {
 203   //Get PYTHIA info: pt-hard, x-section, trials
 204
 205   if (!fQAhistos)
 206     return;
 207
 208   AliRunLoader *rl = AliRunLoader::Instance();
 209   AliGenPythiaEventHeader *genPH = dynamic_cast<AliGenPythiaEventHeader*>(rl->GetHeader()->GenEventHeader());
 210   if(genPH) {
 211     Float_t xsec = genPH->GetXsection();
 212     Int_t trials = genPH->Trials();
 213     Float_t pthard = genPH->GetPtHard();
 214
 215     fHistXsection->Fill(0.5,xsec);
 216     fHistTrials->Fill(0.5,trials);
 217     fHistPtHard->Fill(pthard);
 218   }
 219 }