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   fdNdEtaAnalysisTr(0),
  31   fdNdEtaAnalysisTrVtx(0),
  32   fVertex(0),
  33   fPartPt(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   fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr");
  57   fdNdEtaAnalysisTrVtx = new dNdEtaAnalysis("dndetaTrVtx", "dndetaTrVtx");
  58   fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
  59   for (Int_t i=0; i<3; ++i)
  60   {
  61     fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 120, -6, 6);
  62     fPartEta[i]->Sumw2();
  63   }
  64   fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
  65   fPartPt->Sumw2();
  66 }
  67
  68 void AlidNdEtaAnalysisMCSelector::Init(TTree *tree)
  69 {
  70   AliSelectorRL::Init(tree);
  71
  72   tree->SetBranchStatus("*", 0);
  73   tree->SetBranchStatus("fTriggerMask", 1);
  74   tree->SetBranchStatus("fSPDVertex*", 1);
  75 }
  76
  77 Bool_t AlidNdEtaAnalysisMCSelector::Process(Long64_t entry)
  78 {
  79   // fill the dNdEtaAnalysis class from the monte carlo
  80
  81   if (AliSelectorRL::Process(entry) == kFALSE)
  82     return kFALSE;
  83
  84   // Check prerequisites
  85   if (!fESD)
  86   {
  87     AliDebug(AliLog::kError, "ESD branch not available");
  88     return kFALSE;
  89   }
  90
  91   AliStack* stack = GetStack();
  92   if (!stack)
  93   {
  94     AliDebug(AliLog::kError, "Stack not available");
  95     return kFALSE;
  96   }
  97
  98   AliHeader* header = GetHeader();
  99   if (!header)
 100   {
 101     AliDebug(AliLog::kError, "Header not available");
 102     return kFALSE;
 103   }
 104
 105   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 106   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 107
 108   // get the MC vertex
 109   AliGenEventHeader* genHeader = header->GenEventHeader();
 110
 111   TArrayF vtxMC(3);
 112   genHeader->PrimaryVertex(vtxMC);
 113
 114   // loop over mc particles
 115   Int_t nPrim  = stack->GetNprimary();
 116
 117   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 118   {
 119     TParticle* particle = stack->Particle(iMc);
 120
 121     if (!particle)
 122       continue;
 123
 124     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 125       continue;
 126
 127     AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
 128
 129     fdNdEtaAnalysis->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
 130     fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 131
 132     if (eventTriggered)
 133     {
 134       fdNdEtaAnalysisTr->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
 135       if (vertexReconstructed)
 136         fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], particle->Eta(), particle->Pt(), 1);
 137     }
 138
 139     if (TMath::Abs(vtxMC[2]) < 10)
 140     {
 141       fPartEta[0]->Fill(particle->Eta());
 142
 143       if (vtxMC[2] < 0)
 144         fPartEta[1]->Fill(particle->Eta());
 145       else
 146         fPartEta[2]->Fill(particle->Eta());
 147     }
 148
 149     if (TMath::Abs(particle->Eta()) < 0.8)
 150       fPartPt->Fill(particle->Pt());
 151   }
 152   fdNdEtaAnalysis->FillEvent(vtxMC[2], 1);
 153   if (eventTriggered)
 154   {
 155     fdNdEtaAnalysisTr->FillEvent(vtxMC[2], 1);
 156     if (vertexReconstructed)
 157       fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], 1);
 158   }
 159
 160   return kTRUE;
 161 }
 162
 163 void AlidNdEtaAnalysisMCSelector::SlaveTerminate()
 164 {
 165   // The SlaveTerminate() function is called after all entries or objects
 166   // have been processed. When running with PROOF SlaveTerminate() is called
 167   // on each slave server.
 168
 169   AliSelectorRL::SlaveTerminate();
 170
 171   // Add the histograms to the output on each slave server
 172   if (!fOutput)
 173   {
 174     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 175     return;
 176   }
 177
 178   fOutput->Add(fdNdEtaAnalysis);
 179   fOutput->Add(fdNdEtaAnalysisTr);
 180   fOutput->Add(fdNdEtaAnalysisTrVtx);
 181   fOutput->Add(fPartPt);
 182   for (Int_t i=0; i<3; ++i)
 183     fOutput->Add(fPartEta[i]);
 184 }
 185
 186 void AlidNdEtaAnalysisMCSelector::Terminate()
 187 {
 188   //
 189
 190   AliSelectorRL::Terminate();
 191
 192   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
 193   fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
 194   fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
 195   fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
 196   for (Int_t i=0; i<3; ++i)
 197     fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
 198
 199   if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
 200   {
 201     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p", (void*) fdNdEtaAnalysis, (void*) fPartPt));
 202     return;
 203   }
 204
 205   fdNdEtaAnalysis->Finish(0, -1);
 206   fdNdEtaAnalysisTr->Finish(0, -1);
 207   fdNdEtaAnalysisTrVtx->Finish(0, -1);
 208
 209   Int_t events = (Int_t) fdNdEtaAnalysis->GetVtxHistogram()->Integral();
 210   fPartPt->Scale(1.0/events);
 211   fPartPt->Scale(1.0/fPartPt->GetBinWidth(1));
 212
 213   if (fPartEta[0])
 214   {
 215     TH1* vtxHist = fdNdEtaAnalysis->GetVtxHistogram();
 216     Int_t events1 = (Int_t) vtxHist->Integral(vtxHist->GetXaxis()->FindBin(-9.9), vtxHist->GetXaxis()->FindBin(-0.1));
 217     Int_t events2 = (Int_t) vtxHist->Integral(vtxHist->GetXaxis()->FindBin(0.1), vtxHist->GetXaxis()->FindBin(9.9));
 218
 219     fPartEta[0]->Scale(1.0 / (events1 + events2));
 220     fPartEta[1]->Scale(1.0 / events1);
 221     fPartEta[2]->Scale(1.0 / events2);
 222
 223     for (Int_t i=0; i<3; ++i)
 224       fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
 225
 226     new TCanvas("control", "control", 500, 500);
 227     for (Int_t i=0; i<3; ++i)
 228     {
 229       fPartEta[i]->SetLineColor(i+1);
 230       fPartEta[i]->Draw((i==0) ? "" : "SAME");
 231     }
 232   }
 233
 234   TFile* fout = new TFile("analysis_mc.root","RECREATE");
 235
 236   fdNdEtaAnalysis->SaveHistograms();
 237   fdNdEtaAnalysisTr->SaveHistograms();
 238   fdNdEtaAnalysisTrVtx->SaveHistograms();
 239   fPartPt->Write();
 240
 241   fout->Write();
 242   fout->Close();
 243
 244   if (fPartPt)
 245   {
 246     new TCanvas("control2", "control2", 500, 500);
 247     fPartPt->DrawCopy();
 248   }
 249 }