PWG0/dNdEta/AlidNdEtaAnalysisMCSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaAnalysisMCSelector.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TParticle.h>
   9 #include <TParticlePDG.h>
  10 #include <TVector3.h>
  11 #include <TH1F.h>
  12 #include <TH3F.h>
  13 #include <TTree.h>
  14 #include <TFile.h>
  15
  16 #include <AliLog.h>
  17 #include <AliGenEventHeader.h>
  18 #include <AliHeader.h>
  19 #include <AliStack.h>
  20
  21 #include "dNdEta/dNdEtaAnalysis.h"
  22 #include "AliPWG0Helper.h"
  23
  24
  25 ClassImp(AlidNdEtaAnalysisMCSelector)
  26
  27 AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
  28   AliSelectorRL(),
  29   fdNdEtaAnalysis(0),
  30   fVertex(0),
  31   fPartEta(0),
  32   fEvents(0)
  33 {
  34   //
  35   // Constructor. Initialization of pointers
  36   //
  37 }
  38
  39 AlidNdEtaAnalysisMCSelector::~AlidNdEtaAnalysisMCSelector()
  40 {
  41   //
  42   // Destructor
  43   //
  44 }
  45
  46 void AlidNdEtaAnalysisMCSelector::SlaveBegin(TTree * tree)
  47 {
  48   // The SlaveBegin() function is called after the Begin() function.
  49   // When running with PROOF SlaveBegin() is called on each slave server.
  50   // The tree argument is deprecated (on PROOF 0 is passed).
  51
  52   AliSelectorRL::SlaveBegin(tree);
  53
  54   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
  55   fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
  56   fPartEta = new TH1F("dndeta_check", "dndeta_check", 120, -6, 6);
  57   fPartEta->Sumw2();
  58 }
  59
  60 void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
  61 {
  62   AliSelectorRL::Init(tree);
  63
  64   tree->SetBranchStatus("ESD", 0);
  65 }
  66
  67 Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
  68 {
  69   // fill the dNdEtaAnalysis class from the monte carlo
  70
  71   if (AliSelectorRL::Process(entry) == kFALSE)
  72     return kFALSE;
  73
  74   AliStack* stack = GetStack();
  75   if (!stack)
  76   {
  77     AliDebug(AliLog::kError, "Stack not available");
  78     return kFALSE;
  79   }
  80
  81   AliHeader* header = GetHeader();
  82   if (!header)
  83   {
  84     AliDebug(AliLog::kError, "Header not available");
  85     return kFALSE;
  86   }
  87
  88   // get the MC vertex
  89   AliGenEventHeader* genHeader = header->GenEventHeader();
  90
  91   TArrayF vtxMC(3);
  92   genHeader->PrimaryVertex(vtxMC);
  93
  94   // loop over mc particles
  95   Int_t nPrim  = stack->GetNprimary();
  96
  97   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
  98   {
  99     TParticle* particle = stack->Particle(iMc);
 100
 101     if (!particle)
 102       continue;
 103
 104     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 105       continue;
 106
 107     AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
 108
 109     fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
 110     fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 111
 112     fPartEta->Fill(particle->Eta());
 113   }
 114   fdNdEtaAnalysis->FillEvent(vtxMC[2], 1);
 115
 116   ++fEvents;
 117
 118   return kTRUE;
 119 }
 120
 121 void AlidNdEtaAnalysisMCSelector::SlaveTerminate()
 122 {
 123   // The SlaveTerminate() function is called after all entries or objects
 124   // have been processed. When running with PROOF SlaveTerminate() is called
 125   // on each slave server.
 126
 127   AliSelectorRL::SlaveTerminate();
 128
 129   // Add the histograms to the output on each slave server
 130   if (!fOutput)
 131   {
 132     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 133     return;
 134   }
 135
 136   fOutput->Add(fdNdEtaAnalysis);
 137 }
 138
 139 void AlidNdEtaAnalysisMCSelector::Terminate()
 140 {
 141   //
 142
 143   AliSelectorRL::Terminate();
 144
 145   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
 146
 147   if (!fdNdEtaAnalysis)
 148   {
 149     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fdNdEtaAnalysis));
 150     return;
 151   }
 152
 153   fdNdEtaAnalysis->Finish(0, -1);
 154
 155   TFile* fout = new TFile("analysis_mc.root","RECREATE");
 156
 157   fdNdEtaAnalysis->SaveHistograms();
 158
 159   fout->Write();
 160   fout->Close();
 161
 162   if (fPartEta)
 163   {
 164     fPartEta->Scale(1.0/fEvents);
 165     fPartEta->Scale(1.0/fPartEta->GetBinWidth(1));
 166
 167     TCanvas* canvas = new TCanvas("control", "control", 900, 450);
 168     canvas->Divide(2, 1);
 169
 170     canvas->cd(1);
 171     fVertex->Draw();
 172
 173     canvas->cd(2);
 174     fPartEta->Draw();
 175   }
 176 }