PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AlidNdEtaCorrectionSelector.h"
   4
   5 #include <TStyle.h>
   6 #include <TSystem.h>
   7 #include <TCanvas.h>
   8 #include <TParticle.h>
   9 #include <TParticlePDG.h>
  10
  11 #include <TChain.h>
  12 #include <TSelector.h>
  13 #include <TFile.h>
  14
  15 #include <AliLog.h>
  16 #include <AliGenEventHeader.h>
  17 #include <AliTracker.h>
  18 #include <AliHeader.h>
  19 #include <AliESDVertex.h>
  20 #include <AliESD.h>
  21 #include <AliESDtrack.h>
  22 #include <AliRunLoader.h>
  23 #include <AliStack.h>
  24
  25 #include "esdTrackCuts/AliESDtrackCuts.h"
  26 #include "dNdEta/AlidNdEtaCorrection.h"
  27 #include "AliPWG0Helper.h"
  28
  29 ClassImp(AlidNdEtaCorrectionSelector)
  30
  31 AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector() :
  32   AliSelectorRL(),
  33   fEsdTrackCuts(0),
  34   fdNdEtaCorrection(0)
  35 {
  36   //
  37   // Constructor. Initialization of pointers
  38   //
  39 }
  40
  41 AlidNdEtaCorrectionSelector::~AlidNdEtaCorrectionSelector()
  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 AlidNdEtaCorrectionSelector::Begin(TTree * tree)
  52 {
  53   // The Begin() function is called at the start of the query.
  54   // When running with PROOF Begin() is only called on the client.
  55   // The tree argument is deprecated (on PROOF 0 is passed).
  56
  57   AliSelectorRL::Begin(tree);
  58 }
  59
  60 void AlidNdEtaCorrectionSelector::SlaveBegin(TTree * tree)
  61 {
  62   // The SlaveBegin() function is called after the Begin() function.
  63   // When running with PROOF SlaveBegin() is called on each slave server.
  64   // The tree argument is deprecated (on PROOF 0 is passed).
  65
  66   AliSelectorRL::SlaveBegin(tree);
  67
  68   fdNdEtaCorrection = new AlidNdEtaCorrection();
  69
  70   if (fTree)
  71     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fTree->GetUserInfo()->FindObject("AliESDtrackCuts"));
  72
  73   if (!fEsdTrackCuts)
  74     AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from user info");
  75 }
  76
  77 Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
  78 {
  79   // The Process() function is called for each entry in the tree (or possibly
  80   // keyed object in the case of PROOF) to be processed. The entry argument
  81   // specifies which entry in the currently loaded tree is to be processed.
  82   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
  83   // to read either all or the required parts of the data. When processing
  84   // keyed objects with PROOF, the object is already loaded and is available
  85   // via the fObject pointer.
  86   //
  87   // This function should contain the "body" of the analysis. It can contain
  88   // simple or elaborate selection criteria, run algorithms on the data
  89   // of the event and typically fill histograms.
  90
  91   // WARNING when a selector is used with a TChain, you must use
  92   //  the pointer to the current TTree to call GetEntry(entry).
  93   //  The entry is always the local entry number in the current tree.
  94   //  Assuming that fTree is the pointer to the TChain being processed,
  95   //  use fTree->GetTree()->GetEntry(entry).
  96
  97   if (AliSelectorRL::Process(entry) == kFALSE)
  98     return kFALSE;
  99
 100   // check prerequesites
 101   if (!fESD)
 102   {
 103     AliDebug(AliLog::kError, "ESD branch not available");
 104     return kFALSE;
 105   }
 106
 107   AliHeader* header = GetHeader();
 108   if (!header)
 109   {
 110     AliDebug(AliLog::kError, "Header not available");
 111     return kFALSE;
 112   }
 113
 114   AliRunLoader* runLoader = GetRunLoader();
 115   if (!runLoader)
 116   {
 117     AliDebug(AliLog::kError, "RunLoader not available");
 118     return kFALSE;
 119   }
 120
 121   runLoader->LoadKinematics();
 122   AliStack* stack = runLoader->Stack();
 123   if (!stack)
 124   {
 125     AliDebug(AliLog::kError, "Stack not available");
 126     return kFALSE;
 127   }
 128
 129   if (!fEsdTrackCuts)
 130   {
 131     AliDebug(AliLog::kError, "fESDTrackCuts not available");
 132     return kFALSE;
 133   }
 134
 135   Bool_t vertexReconstructed = AliPWG0Helper::IsVertexReconstructed(fESD);
 136
 137   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
 138
 139   fdNdEtaCorrection->IncreaseEventCount();
 140   if (eventTriggered)
 141     fdNdEtaCorrection->IncreaseTriggeredEventCount();
 142
 143   // get the MC vertex
 144   AliGenEventHeader* genHeader = header->GenEventHeader();
 145
 146   TArrayF vtxMC(3);
 147   genHeader->PrimaryVertex(vtxMC);
 148
 149   // loop over mc particles
 150   Int_t nPrim  = stack->GetNprimary();
 151
 152   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 153   {
 154     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 155
 156     TParticle* particle = stack->Particle(iMc);
 157
 158     if (!particle)
 159     {
 160       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 161       continue;
 162     }
 163
 164     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
 165       continue;
 166
 167     Float_t eta = particle->Eta();
 168     Float_t pt = particle->Pt();
 169
 170     if (vertexReconstructed)
 171       fdNdEtaCorrection->FillParticle(vtxMC[2], eta, pt);
 172
 173     fdNdEtaCorrection->FillParticleAllEvents(eta, pt);
 174     if (eventTriggered)
 175       fdNdEtaCorrection->FillParticleWhenEventTriggered(eta, pt);
 176   }// end of mc particle
 177
 178   // ########################################################
 179   // loop over esd tracks
 180   Int_t nTracks = fESD->GetNumberOfTracks();
 181
 182   // count the number of "good" tracks for vertex reconstruction efficiency
 183   // TODO change to number of ITS clusters or similar
 184   Int_t nGoodTracks = 0;
 185   for (Int_t t=0; t<nTracks; t++)
 186   {
 187     AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t));
 188
 189     AliESDtrack* esdTrack = fESD->GetTrack(t);
 190
 191     // cut the esd track?
 192     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
 193       continue;
 194
 195     nGoodTracks++;
 196
 197     // using the properties of the mc particle
 198     Int_t label = TMath::Abs(esdTrack->GetLabel());
 199     if (label == 0)
 200     {
 201       AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
 202       continue;
 203     }
 204
 205     TParticle* particle = stack->Particle(label);
 206     if (!particle)
 207     {
 208       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
 209       continue;
 210     }
 211
 212     if (vertexReconstructed)
 213       fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
 214   } // end of track loop
 215
 216
 217   if (eventTriggered)
 218   {
 219     fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
 220     if (vertexReconstructed)
 221       fdNdEtaCorrection->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
 222   }
 223
 224   return kTRUE;
 225 }
 226
 227 void AlidNdEtaCorrectionSelector::SlaveTerminate()
 228 {
 229   // The SlaveTerminate() function is called after all entries or objects
 230   // have been processed. When running with PROOF SlaveTerminate() is called
 231   // on each slave server.
 232
 233   AliSelectorRL::SlaveTerminate();
 234
 235   // Add the histograms to the output on each slave server
 236   if (!fOutput)
 237   {
 238     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
 239     return;
 240   }
 241
 242   fOutput->Add(fdNdEtaCorrection);
 243 }
 244
 245 void AlidNdEtaCorrectionSelector::Terminate()
 246 {
 247   // The Terminate() function is the last function to be called during
 248   // a query. It always runs on the client, it can be used to present
 249   // the results graphically or save the results to file.
 250
 251   AliSelectorRL::Terminate();
 252
 253   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
 254
 255   fdNdEtaCorrection->Finish();
 256
 257   TFile* fout = new TFile("correction_map.root","RECREATE");
 258
 259   fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 260   fdNdEtaCorrection->SaveHistograms();
 261
 262   fout->Write();
 263   fout->Close();
 264
 265   fdNdEtaCorrection->DrawHistograms();
 266 }