PWG0/dNdEta/AlidNdEtaVertexRecEffSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaVertexRecEffSelector.h"
   4
   5 #include <TCanvas.h>
   6 #include <TH1F.h>
   7 #include <TTree.h>
   8 #include <TParticle.h>
   9 #include <TFile.h>
  10
  11 #include <AliLog.h>
  12 #include <AliGenEventHeader.h>
  13 #include <AliHeader.h>
  14 #include <AliESD.h>
  15 #include <AliESDVertex.h>
  16
  17 #include "AliPWG0Helper.h"
  18
  19 //
  20 // This class plots the vertex reconstruction efficiency
  21 // If a vertex was reconstructed is decided by the function CheckVertex()
  22 // In any case the *generated* multiplicity is filled into the histogram
  23 //
  24
  25 ClassImp(AlidNdEtaVertexRecEffSelector)
  26
  27 const Float_t AlidNdEtaVertexRecEffSelector::fkEtaRange = 0.9;
  28
  29 AlidNdEtaVertexRecEffSelector::AlidNdEtaVertexRecEffSelector() :
  30   AliSelectorRL(),
  31   fdNGen(0),
  32   fdNRec(0),
  33   fVtxGen(0),
  34   fVtxRec(0)
  35 {
  36   //
  37   // Constructor. Initialization of pointers
  38   //
  39 }
  40
  41 AlidNdEtaVertexRecEffSelector::~AlidNdEtaVertexRecEffSelector()
  42 {
  43   //
  44   // Destructor
  45   //
  46
  47   // histograms are in the output list and deleted when the output
  48   // list is deleted by the TSelector dtor
  49 }
  50
  51 void AlidNdEtaVertexRecEffSelector::SlaveBegin(TTree * tree)
  52 {
  53   // initializes the histograms
  54
  55   AliSelectorRL::SlaveBegin(tree);
  56
  57   fdNGen = new TH1F("dNGen", "dNGen", 90, 0, 50);
  58   fdNRec = dynamic_cast<TH1F*>(fdNGen->Clone("dNRec"));
  59
  60   fdNGen->SetTitle("Generated Events;dN_{Gen};Count");
  61   fdNRec->SetTitle("Events with reconstructed vertex;dN_{Gen};Count");
  62
  63   fVtxGen = new TH1F("VtxGen", "VtxGen", 200, -20, 20);
  64   fVtxRec = dynamic_cast<TH1F*>(fVtxGen->Clone("VtxRec"));
  65 }
  66
  67 Bool_t AlidNdEtaVertexRecEffSelector::CheckVertex()
  68 {
  69   //
  70   // check if the vertex has been reconstructed well enough
  71   //
  72   if (!fESD)
  73     return kFALSE;
  74
  75   const AliESDVertex* vtxESD = fESD->GetVertex();
  76
  77   // the vertex should be reconstructed
  78   if (strcmp(vtxESD->GetName(),"default")==0)
  79     return kFALSE;
  80
  81   Double_t vtxRes[3];
  82   vtxRes[0] = vtxESD->GetXRes();
  83   vtxRes[1] = vtxESD->GetYRes();
  84   vtxRes[2] = vtxESD->GetZRes();
  85
  86   // the resolution should be reasonable???
  87   if (vtxRes[2]==0 || vtxRes[2]>0.1)
  88     return kFALSE;
  89
  90   return kTRUE;
  91 }
  92
  93 Bool_t AlidNdEtaVertexRecEffSelector::Process(Long64_t entry)
  94 {
  95   //
  96   // fills fdNGen and fdNRec
  97   //
  98
  99   if (AliSelectorRL::Process(entry) == kFALSE)
 100     return kFALSE;
 101
 102   // check prerequisites
 103   if (!fESD)
 104   {
 105     AliDebug(AliLog::kError, "ESD branch not available");
 106     return kFALSE;
 107   }
 108
 109   AliHeader* header = GetHeader();
 110   if (!header)
 111   {
 112     AliDebug(AliLog::kError, "Header not available");
 113     return kFALSE;
 114   }
 115
 116
 117   // loop over mc particles
 118   TTree* particleTree = GetKinematics();
 119   TParticle* particle = 0;
 120   particleTree->SetBranchAddress("Particles", &particle);
 121
 122   Int_t nPrim  = header->GetNprimary();
 123   Int_t nTotal = header->GetNtrack();
 124
 125   Int_t n = 0;
 126
 127   for (Int_t iMc = nTotal - nPrim; iMc < nTotal; ++iMc)
 128   {
 129     particleTree->GetEntry(iMc);
 130
 131     if (!particle)
 132       continue;
 133
 134     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 135       continue;
 136
 137     if (TMath::Abs(particle->Eta()) < fkEtaRange)
 138       ++n;
 139   }// end of mc particle
 140
 141   Float_t dN = (Float_t) n / (fkEtaRange*2);
 142
 143   fdNGen->Fill(dN);
 144
 145   // check vertex reconstruction
 146   if (CheckVertex() != kFALSE)
 147     fdNRec->Fill(dN);
 148
 149   AliGenEventHeader* genHeader = header->GenEventHeader();
 150
 151   TArrayF vtxMC(3);
 152   genHeader->PrimaryVertex(vtxMC);
 153
 154   fVtxGen->Fill(vtxMC[2]);
 155
 156   if (CheckVertex() != kFALSE)
 157     fVtxRec->Fill(vtxMC[2]);
 158
 159   return kTRUE;
 160 }
 161
 162 void AlidNdEtaVertexRecEffSelector::SlaveTerminate()
 163 {
 164   // The SlaveTerminate() function is called after all entries or objects
 165   // have been processed. When running with PROOF SlaveTerminate() is called
 166   // on each slave server.
 167
 168   AliSelectorRL::SlaveTerminate();
 169
 170   // Add the histograms to the output on each slave server
 171   if (!fOutput)
 172   {
 173     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
 174     return;
 175   }
 176
 177   fOutput->Add(fdNGen);
 178   fOutput->Add(fdNRec);
 179   fOutput->Add(fVtxGen);
 180   fOutput->Add(fVtxRec);
 181 }
 182
 183 void AlidNdEtaVertexRecEffSelector::Terminate()
 184 {
 185   // The Terminate() function is the last function to be called during
 186   // a query. It always runs on the client, it can be used to present
 187   // the results graphically or save the results to file.
 188
 189   AliSelectorRL::Terminate();
 190
 191   if (!fOutput)
 192   {
 193     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
 194     return;
 195   }
 196
 197   fdNGen = dynamic_cast<TH1F*> (fOutput->FindObject("dNGen"));
 198   fdNRec = dynamic_cast<TH1F*> (fOutput->FindObject("dNRec"));
 199   fVtxGen = dynamic_cast<TH1F*> (fOutput->FindObject("VtxGen"));
 200   fVtxRec = dynamic_cast<TH1F*> (fOutput->FindObject("VtxRec"));
 201   if (!fdNGen || !fdNRec || !fVtxGen || !fVtxRec)
 202   {
 203     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p %p %p", (void*) fdNGen, (void*) fdNRec, (void*) fVtxGen, (void*) fVtxRec));
 204     return;
 205   }
 206
 207   TFile* fout = new TFile("vertexRecEff.root","RECREATE");
 208
 209   fdNGen->Write();
 210   fdNRec->Write();
 211
 212   fVtxGen->Write();
 213   fVtxRec->Write();
 214
 215   fout->Write();
 216   fout->Close();
 217
 218   TCanvas* canvas = new TCanvas("dN", "dN", 900, 450);
 219   canvas->Divide(2, 1);
 220
 221   canvas->cd(1);
 222   fdNGen->Draw();
 223   fdNRec->SetLineColor(kRed);
 224   fdNRec->Draw("SAME");
 225
 226   canvas->cd(2);
 227   fVtxGen->Draw();
 228   fVtxRec->SetLineColor(kRed);
 229   fVtxRec->Draw("SAME");
 230 }