PWG0/esdTrackCuts/AliTestESDtrackCutsSelector.cxx

   1 /* $Id$ */
   2
   3 #include "AliTestESDtrackCutsSelector.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 <TChain.h>
  12 #include <TFile.h>
  13 #include <TH1F.h>
  14 #include <TH3F.h>
  15
  16 #include <TSelector.h>
  17 #include <TFile.h>
  18
  19 #include <AliLog.h>
  20 #include <AliESD.h>
  21 #include <AliRunLoader.h>
  22 #include <AliStack.h>
  23
  24 #include "esdTrackCuts/AliESDtrackCuts.h"
  25 #include "AliPWG0Helper.h"
  26
  27
  28 ClassImp(AliTestESDtrackCutsSelector)
  29
  30 AliTestESDtrackCutsSelector::AliTestESDtrackCutsSelector() :
  31   AliSelectorRL(),
  32   fEsdTrackCutsAll(0),
  33   fEsdTrackCutsNoVtx(0),
  34   fEsdTrackCutsPri(0),
  35   fEsdTrackCutsSec(0),
  36   fEsdTrackCutsPlusZ(0),
  37   fEsdTrackCutsMinusZ(0),
  38   fEsdTrackCutsPos(0),
  39   fEsdTrackCutsNeg(0),
  40   fPIDAfterCutNoVtx(0),
  41   fPIDAfterCutAll(0),
  42   fVertex(0)
  43 {
  44   //
  45   // Constructor. Initialization of pointers
  46   //
  47 }
  48
  49 AliTestESDtrackCutsSelector::~AliTestESDtrackCutsSelector()
  50 {
  51   //
  52   // Destructor
  53   //
  54
  55   // histograms are in the output list and deleted when the output
  56   // list is deleted by the TSelector dtor
  57 }
  58
  59 void AliTestESDtrackCutsSelector::Begin(TTree* tree)
  60 {
  61   // Begin function
  62
  63   AliSelectorRL::Begin(tree);
  64
  65   ReadUserObjects(tree);
  66 }
  67
  68 void AliTestESDtrackCutsSelector::ReadUserObjects(TTree* tree)
  69 {
  70   // read the user objects, called from slavebegin and begin
  71
  72   // only do it once
  73   if (fEsdTrackCutsAll)
  74     return;
  75
  76   if (!fEsdTrackCutsAll && fInput)
  77     fEsdTrackCutsAll = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("esdTrackCutsAll")->Clone());
  78
  79   if (!fEsdTrackCutsAll && tree)
  80     fEsdTrackCutsAll = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("esdTrackCutsAll"));
  81
  82   if (!fEsdTrackCutsAll)
  83      AliDebug(AliLog::kError, "ERROR: Could not read fEsdTrackCutsAll from input list.");
  84
  85   fEsdTrackCutsNoVtx =    dynamic_cast<AliESDtrackCuts*> (fEsdTrackCutsAll->Clone("fEsdTrackCutsNoVtx"));
  86   fEsdTrackCutsNoVtx->SetRequireSigmaToVertex(kFALSE);
  87
  88   fEsdTrackCutsPri =    dynamic_cast<AliESDtrackCuts*> (fEsdTrackCutsAll->Clone("fEsdTrackCutsPri"));
  89   fEsdTrackCutsSec =    dynamic_cast<AliESDtrackCuts*> (fEsdTrackCutsAll->Clone("fEsdTrackCutsSec"));
  90   fEsdTrackCutsPlusZ =  dynamic_cast<AliESDtrackCuts*> (fEsdTrackCutsAll->Clone("fEsdTrackCutsPlusZ"));
  91   fEsdTrackCutsMinusZ = dynamic_cast<AliESDtrackCuts*> (fEsdTrackCutsAll->Clone("fEsdTrackCutsMinusZ"));
  92   fEsdTrackCutsPos =    dynamic_cast<AliESDtrackCuts*> (fEsdTrackCutsAll->Clone("fEsdTrackCutsPos"));
  93   fEsdTrackCutsNeg =    dynamic_cast<AliESDtrackCuts*> (fEsdTrackCutsAll->Clone("fEsdTrackCutsNeg"));
  94 }
  95
  96 void AliTestESDtrackCutsSelector::SlaveBegin(TTree* tree)
  97 {
  98   // The SlaveBegin() function is called after the Begin() function.
  99   // When running with PROOF SlaveBegin() is called on each slave server.
 100   // The tree argument is deprecated (on PROOF 0 is passed).
 101
 102   AliSelectorRL::SlaveBegin(tree);
 103
 104   ReadUserObjects(tree);
 105
 106   fPIDAfterCutNoVtx = new TH1F("fPIDAfterCutNoVtx", "fPIDAfterCutNoVtx", 5001, -2500.5, 2500.5);
 107   fPIDAfterCutAll   = new TH1F("fPIDAfterCutAll", "fPIDAfterCutAll", 5001, -2500.5, 2500.5);
 108
 109   fVertex = new TH3F("fVertex", "fVertex", 100, -10, 10, 100, -10, 10, 100, -10, 10);
 110 }
 111
 112 void AliTestESDtrackCutsSelector::Init(TTree* tree)
 113 {
 114   // read the user objects
 115
 116   AliSelectorRL::Init(tree);
 117
 118   // Enable only the needed branches
 119   if (tree)
 120   {
 121     tree->SetBranchStatus("*", 0);
 122     tree->SetBranchStatus("fTriggerMask", 1);
 123     tree->SetBranchStatus("fSPDVertex*", 1);
 124     tree->SetBranchStatus("fTracks.fLabel", 1);
 125
 126     AliESDtrackCuts::EnableNeededBranches(tree);
 127   }
 128 }
 129
 130 Bool_t AliTestESDtrackCutsSelector::Process(Long64_t entry)
 131 {
 132   // The Process() function is called for each entry in the tree (or possibly
 133   // keyed object in the case of PROOF) to be processed. The entry argument
 134   // specifies which entry in the currently loaded tree is to be processed.
 135   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 136   // to read either all or the required parts of the data. When processing
 137   // keyed objects with PROOF, the object is already loaded and is available
 138   // via the fObject pointer.
 139   //
 140   // This function should contain the "body" of the analysis. It can contain
 141   // simple or elaborate selection criteria, run algorithms on the data
 142   // of the event and typically fill histograms.
 143
 144   // WARNING when a selector is used with a TChain, you must use
 145   //  the pointer to the current TTree to call GetEntry(entry).
 146   //  The entry is always the local entry number in the current tree.
 147   //  Assuming that fTree is the pointer to the TChain being processed,
 148   //  use fTree->GetTree()->GetEntry(entry).
 149
 150   if (AliSelectorRL::Process(entry) == kFALSE)
 151     return kFALSE;
 152
 153   // Check prerequisites
 154   if (!fESD) {
 155     AliDebug(AliLog::kError, "ESD branch not available");
 156     return kFALSE;
 157   }
 158
 159   if (!AliPWG0Helper::IsVertexReconstructed(fESD)) {
 160     AliDebug(AliLog::kDebug+5, "Vertex is not reconstructed");
 161     return kFALSE;
 162   }
 163
 164   // check if the esd track cut objects are there
 165   if (!fEsdTrackCutsAll || !fEsdTrackCutsPri || !fEsdTrackCutsSec || !fEsdTrackCutsPlusZ || !fEsdTrackCutsMinusZ || !fEsdTrackCutsPos || !fEsdTrackCutsNeg) {
 166     AliDebug(AliLog::kError, "fEsdTrackCutsXXX not available");
 167     return kFALSE;
 168   }
 169
 170   // get particle stack
 171   AliStack* stack = GetStack();
 172   if (!stack) {
 173     AliDebug(AliLog::kError, "Stack not available");
 174     return kFALSE;
 175   }
 176   Int_t nPrim  = stack->GetNprimary();
 177
 178   // ########################################################
 179   // loop over esd tracks
 180   Int_t nTracks = fESD->GetNumberOfTracks();
 181
 182   // count the number of "good" tracks as parameter for vertex reconstruction efficiency
 183   for (Int_t t=0; t<nTracks; t++) {
 184     AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t));
 185
 186     AliESDtrack* esdTrack = fESD->GetTrack(t);
 187
 188     Bool_t passed = fEsdTrackCutsAll->AcceptTrack(esdTrack);
 189
 190     // using the properties of the mc particle
 191     Int_t label = TMath::Abs(esdTrack->GetLabel());
 192     if (label == 0) {
 193       AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
 194       continue;
 195     }
 196     TParticle* particle = stack->Particle(label);
 197     if (!particle) {
 198       AliDebug(AliLog::kError, Form("UNEXPECTED: part with label %d not found in stack (track loop).", label));
 199       continue;
 200     }
 201
 202     if (label < nPrim)
 203       fEsdTrackCutsPri->AcceptTrack(esdTrack);
 204     else
 205     {
 206       fEsdTrackCutsSec->AcceptTrack(esdTrack);
 207       if (passed)
 208       {
 209         fPIDAfterCutAll->Fill(particle->GetPdgCode());
 210         fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
 211       }
 212
 213       if (fEsdTrackCutsNoVtx->AcceptTrack(esdTrack))
 214         fPIDAfterCutNoVtx->Fill(particle->GetPdgCode());
 215     }
 216
 217     TParticlePDG* pdgPart = particle->GetPDG();
 218     if (pdgPart)
 219     {
 220       if (pdgPart->Charge() > 0)
 221         fEsdTrackCutsPos->AcceptTrack(esdTrack);
 222       else if (pdgPart->Charge() < 0)
 223         fEsdTrackCutsNeg->AcceptTrack(esdTrack);
 224     }
 225
 226     if (particle->Eta() < 0)
 227       fEsdTrackCutsPlusZ->AcceptTrack(esdTrack);
 228     else
 229       fEsdTrackCutsMinusZ->AcceptTrack(esdTrack);
 230   }
 231
 232   return kTRUE;
 233 }
 234
 235 void AliTestESDtrackCutsSelector::SlaveTerminate()
 236 {
 237   // The SlaveTerminate() function is called after all entries or objects
 238   // have been processed. When running with PROOF SlaveTerminate() is called
 239   // on each slave server.
 240
 241   AliSelectorRL::SlaveTerminate();
 242
 243   // Add the histograms to the output on each slave server
 244   if (!fOutput)
 245   {
 246     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
 247     return;
 248   }
 249
 250   fOutput->Add(fEsdTrackCutsAll);
 251   fOutput->Add(fEsdTrackCutsNoVtx);
 252   fOutput->Add(fEsdTrackCutsPri);
 253   fOutput->Add(fEsdTrackCutsSec);
 254   fOutput->Add(fEsdTrackCutsPlusZ);
 255   fOutput->Add(fEsdTrackCutsMinusZ);
 256   fOutput->Add(fEsdTrackCutsPos);
 257   fOutput->Add(fEsdTrackCutsNeg);
 258   fOutput->Add(fPIDAfterCutNoVtx);
 259   fOutput->Add(fPIDAfterCutAll);
 260   fOutput->Add(fVertex);
 261 }
 262
 263 void AliTestESDtrackCutsSelector::Terminate()
 264 {
 265   // The Terminate() function is the last function to be called during
 266   // a query. It always runs on the client, it can be used to present
 267   // the results graphically or save the results to file.
 268
 269   AliSelectorRL::Terminate();
 270
 271   fEsdTrackCutsAll = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("esdTrackCutsAll"));
 272   fEsdTrackCutsNoVtx = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsNoVtx"));
 273   fEsdTrackCutsPri = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPri"));
 274   fEsdTrackCutsSec = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsSec"));
 275   fEsdTrackCutsPlusZ = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPlusZ"));
 276   fEsdTrackCutsMinusZ = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsMinusZ"));
 277   fEsdTrackCutsPos = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsPos"));
 278   fEsdTrackCutsNeg = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCutsNeg"));
 279   fPIDAfterCutNoVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDAfterCutNoVtx"));
 280   fPIDAfterCutAll = dynamic_cast<TH1F*> (fOutput->FindObject("fPIDAfterCutAll"));
 281   fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("fVertex"));
 282
 283   // check if the esd track cut objects are there
 284   if (!fEsdTrackCutsAll || !fEsdTrackCutsPri || !fEsdTrackCutsSec || !fEsdTrackCutsPlusZ || !fEsdTrackCutsMinusZ || !fEsdTrackCutsPos || !fEsdTrackCutsNeg) {
 285     AliDebug(AliLog::kError, Form("fEsdTrackCutsXXX not available %p %p %p %p %p %p %p", fEsdTrackCutsAll, fEsdTrackCutsPri, fEsdTrackCutsSec, fEsdTrackCutsPlusZ, fEsdTrackCutsMinusZ, fEsdTrackCutsPos, fEsdTrackCutsNeg));
 286     return;
 287   }
 288
 289   TFile* file = TFile::Open("trackCuts.root", "RECREATE");
 290
 291   fEsdTrackCutsAll->SaveHistograms();
 292   fEsdTrackCutsNoVtx->SaveHistograms();
 293   fEsdTrackCutsPri->SaveHistograms();
 294   fEsdTrackCutsSec->SaveHistograms();
 295   fEsdTrackCutsPlusZ->SaveHistograms();
 296   fEsdTrackCutsMinusZ->SaveHistograms();
 297   fEsdTrackCutsPos->SaveHistograms();
 298   fEsdTrackCutsNeg->SaveHistograms();
 299   fPIDAfterCutNoVtx->Write();
 300   fPIDAfterCutAll->Write();
 301   fVertex->Write();
 302
 303   file->Close();
 304
 305         fEsdTrackCutsAll->DrawHistograms();
 306 }