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
  20 #include "dNdEta/dNdEtaAnalysis.h"
  21 #include "AliPWG0Helper.h"
  22
  23
  24 ClassImp(AlidNdEtaAnalysisMCSelector)
  25
  26 AlidNdEtaAnalysisMCSelector::AlidNdEtaAnalysisMCSelector() :
  27   AliSelectorRL(),
  28   fdNdEtaAnalysis(0),
  29   fVertex(0),
  30   fPartEta(0),
  31   fEvents(0)
  32 {
  33   //
  34   // Constructor. Initialization of pointers
  35   //
  36 }
  37
  38 AlidNdEtaAnalysisMCSelector::~AlidNdEtaAnalysisMCSelector()
  39 {
  40   //
  41   // Destructor
  42   //
  43 }
  44
  45 void AlidNdEtaAnalysisMCSelector::SlaveBegin(TTree * tree)
  46 {
  47   // The SlaveBegin() function is called after the Begin() function.
  48   // When running with PROOF SlaveBegin() is called on each slave server.
  49   // The tree argument is deprecated (on PROOF 0 is passed).
  50
  51   AliSelectorRL::SlaveBegin(tree);
  52
  53   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
  54 }
  55
  56 void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
  57 {
  58   AliSelectorRL::Init(tree);
  59
  60   tree->SetBranchStatus("ESD", 0);
  61
  62   fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
  63   fPartEta = new TH1F("dndeta_check", "dndeta_check", 120, -6, 6);
  64   fPartEta->Sumw2();
  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   TTree* particleTree = GetKinematics();
  75   if (!particleTree)
  76   {
  77     AliDebug(AliLog::kError, "Kinematics 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   particleTree->SetBranchStatus("*", 0);
  95   particleTree->SetBranchStatus("fDaughter[2]", 1);
  96   particleTree->SetBranchStatus("fPdgCode", 1);
  97   particleTree->SetBranchStatus("fPx", 1);
  98   particleTree->SetBranchStatus("fPy", 1);
  99   particleTree->SetBranchStatus("fPz", 1);
 100   particleTree->SetBranchStatus("fVx", 1);
 101   particleTree->SetBranchStatus("fVy", 1);
 102   particleTree->SetBranchStatus("fVz", 1);
 103
 104   TParticle* particle = 0;
 105   particleTree->SetBranchAddress("Particles", &particle);
 106
 107   Int_t nPrim  = header->GetNprimary();
 108   Int_t nTotal = header->GetNtrack();
 109
 110   for (Int_t i_mc = nTotal - nPrim; i_mc < nTotal; ++i_mc)
 111   {
 112     particleTree->GetEntry(i_mc);
 113
 114     if (!particle)
 115       continue;
 116
 117     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 118       continue;
 119
 120     AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", i_mc, particle->GetUniqueID()));
 121
 122     fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
 123     fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 124
 125     fPartEta->Fill(particle->Eta());
 126   }
 127   fdNdEtaAnalysis->FillEvent(vtxMC[2], 1);
 128
 129   ++fEvents;
 130
 131   return kTRUE;
 132 }
 133
 134 void AlidNdEtaAnalysisMCSelector::SlaveTerminate()
 135 {
 136   // The SlaveTerminate() function is called after all entries or objects
 137   // have been processed. When running with PROOF SlaveTerminate() is called
 138   // on each slave server.
 139
 140   AliSelectorRL::SlaveTerminate();
 141
 142   // Add the histograms to the output on each slave server
 143   if (!fOutput)
 144   {
 145     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 146     return;
 147   }
 148
 149   fOutput->Add(fdNdEtaAnalysis);
 150 }
 151
 152 void AlidNdEtaAnalysisMCSelector::Terminate()
 153 {
 154   //
 155
 156   AliSelectorRL::Terminate();
 157
 158   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
 159
 160   if (!fdNdEtaAnalysis)
 161   {
 162     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p", (void*) fdNdEtaAnalysis));
 163     return;
 164   }
 165
 166   TFile* fout = new TFile("analysis_mc.root","RECREATE");
 167
 168   fdNdEtaAnalysis->SaveHistograms();
 169
 170   fout->Write();
 171   fout->Close();
 172
 173   fPartEta->Scale(1.0/fEvents);
 174   fPartEta->Scale(1.0/fPartEta->GetBinWidth(1));
 175
 176   TCanvas* canvas = new TCanvas("control", "control", 900, 450);
 177   canvas->Divide(2, 1);
 178
 179   canvas->cd(1);
 180   fVertex->Draw();
 181
 182   canvas->cd(2);
 183   fPartEta->Draw();
 184 }