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