PWG0/dNdEta/AlidNdEtaCorrectionSelector.cxx

   1 #include "AlidNdEtaCorrectionSelector.h"
   2
   3 #include <TStyle.h>
   4 #include <TSystem.h>
   5 #include <TCanvas.h>
   6 #include <TParticle.h>
   7 #include <TParticlePDG.h>
   8
   9 #include <AliLog.h>
  10 #include <AliGenEventHeader.h>
  11 #include <AliTracker.h>
  12
  13 #include "../esdTrackCuts/AliESDtrackCuts.h"
  14 #include "dNdEtaCorrection.h"
  15
  16 ClassImp(AlidNdEtaCorrectionSelector)
  17
  18 AlidNdEtaCorrectionSelector::AlidNdEtaCorrectionSelector(TTree *) :
  19   AliSelector(),
  20   fEsdTrackCuts(0),
  21   fdNdEtaCorrection(0)
  22 {
  23   //
  24   // Constructor. Initialization of pointers
  25   //
  26 }
  27
  28 AlidNdEtaCorrectionSelector::~AlidNdEtaCorrectionSelector()
  29 {
  30   //
  31   // Destructor
  32   //
  33
  34   // histograms are in the output list and deleted when the output
  35   // list is deleted by the TSelector dtor
  36 }
  37
  38 void AlidNdEtaCorrectionSelector::Begin(TTree * tree)
  39 {
  40   // The Begin() function is called at the start of the query.
  41   // When running with PROOF Begin() is only called on the client.
  42   // The tree argument is deprecated (on PROOF 0 is passed).
  43
  44   AliSelector::Begin(tree);
  45 }
  46
  47 void AlidNdEtaCorrectionSelector::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   AliSelector::SlaveBegin(tree);
  54
  55   fdNdEtaCorrection = new dNdEtaCorrection();
  56
  57   if (fChain)
  58     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fChain->GetUserInfo()->FindObject("AliESDtrackCuts"));
  59
  60   if (!fEsdTrackCuts)
  61     printf("ERROR: Could not read EsdTrackCuts from user info\n");
  62
  63   AliLog::SetClassDebugLevel("ESDtrackQualityCuts",1);
  64 }
  65
  66 Bool_t AlidNdEtaCorrectionSelector::Notify()
  67 {
  68   // The Notify() function is called when a new file is opened. This
  69   // can be either for a new TTree in a TChain or when when a new TTree
  70   // is started when using PROOF. Typically here the branch pointers
  71   // will be retrieved. It is normaly not necessary to make changes
  72   // to the generated code, but the routine can be extended by the
  73   // user if needed.
  74
  75   if (AliSelector::Notify() == kFALSE)
  76     return kFALSE;
  77
  78   return kTRUE;
  79 }
  80
  81 Bool_t AlidNdEtaCorrectionSelector::IsPrimary(const TParticle* aParticle, Int_t aTotalPrimaries)
  82 {
  83   //
  84   // Returns if the given particle is a primary particle
  85   // This function or a equivalent should be available in some common place of AliRoot
  86   //
  87
  88   // if the particle has a daughter primary, we do not want to count it
  89   if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries)
  90     return kFALSE;
  91
  92   Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode());
  93
  94   // skip quarks and gluon
  95   if (pdgCode > 10 && pdgCode != 21)
  96     return kTRUE;
  97
  98   return kFALSE;
  99 }
 100
 101 Bool_t AlidNdEtaCorrectionSelector::Process(Long64_t entry)
 102 {
 103   // The Process() function is called for each entry in the tree (or possibly
 104   // keyed object in the case of PROOF) to be processed. The entry argument
 105   // specifies which entry in the currently loaded tree is to be processed.
 106   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
 107   // to read either all or the required parts of the data. When processing
 108   // keyed objects with PROOF, the object is already loaded and is available
 109   // via the fObject pointer.
 110   //
 111   // This function should contain the "body" of the analysis. It can contain
 112   // simple or elaborate selection criteria, run algorithms on the data
 113   // of the event and typically fill histograms.
 114
 115   // WARNING when a selector is used with a TChain, you must use
 116   //  the pointer to the current TTree to call GetEntry(entry).
 117   //  The entry is always the local entry number in the current tree.
 118   //  Assuming that fChain is the pointer to the TChain being processed,
 119   //  use fChain->GetTree()->GetEntry(entry).
 120
 121   if (AliSelector::Process(entry) == kFALSE)
 122     return kFALSE;
 123
 124   if (!fESD || !fHeader)
 125     return kFALSE;
 126
 127   // ########################################################
 128   // get the EDS vertex
 129   const AliESDVertex* vtxESD = fESD->GetVertex();
 130
 131   // the vertex should be reconstructed
 132   if (strcmp(vtxESD->GetName(),"default")==0)
 133     return kTRUE;
 134
 135   Double_t vtx_res[3];
 136   vtx_res[0] = vtxESD->GetXRes();
 137   vtx_res[1] = vtxESD->GetYRes();
 138   vtx_res[2] = vtxESD->GetZRes();
 139
 140   // the resolution should be reasonable???
 141   if (vtx_res[2]==0 || vtx_res[2]>0.1)
 142     return kTRUE;
 143
 144   // ########################################################
 145   // get the MC vertex
 146   AliGenEventHeader* genHeader = fHeader->GenEventHeader();
 147
 148   TArrayF vtxMC(3);
 149   genHeader->PrimaryVertex(vtxMC);
 150
 151   // ########################################################
 152   // loop over mc particles
 153   TTree* particleTree = GetKinematics();
 154   TParticle* particle = 0;
 155   particleTree->SetBranchAddress("Particles", &particle);
 156
 157   Int_t nPrim  = fHeader->GetNprimary();
 158   Int_t nTotal = fHeader->GetNtrack();
 159
 160   for (Int_t i_mc = nTotal - nPrim; i_mc < nTotal; ++i_mc)
 161   {
 162     particleTree->GetEntry(i_mc);
 163
 164     if (!particle)
 165       continue;
 166
 167     if (strcmp(particle->GetName(),"XXX") == 0)
 168     {
 169       printf("WARNING: There is a particle named XXX (%d).\n", i_mc);
 170       continue;
 171     }
 172
 173     TParticlePDG* pdgPart = particle->GetPDG();
 174
 175     if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0)
 176     {
 177       printf("WARNING: There is a particle with an unknown particle class (%d pdg code %d).\n", i_mc, particle->GetPdgCode());
 178       continue;
 179     }
 180
 181     if (IsPrimary(particle, nPrim) == kFALSE)
 182       continue;
 183
 184     if (pdgPart->Charge() == 0)
 185       continue;
 186
 187     fdNdEtaCorrection->FillGene(vtxMC[2], particle->Eta());
 188
 189   }// end of mc particle
 190
 191   // ########################################################
 192   // loop over esd tracks
 193   Int_t nTracks = fESD->GetNumberOfTracks();
 194   for (Int_t t=0; t<nTracks; t++)
 195   {
 196     AliESDtrack* esdTrack = fESD->GetTrack(t);
 197
 198     // cut the esd track?
 199     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
 200       continue;
 201
 202     // using the eta of the mc particle
 203     Int_t label = TMath::Abs(esdTrack->GetLabel());
 204     if (label == 0)
 205     {
 206       printf("WARNING: cannot find corresponding mc part for track %d.", t);
 207       continue;
 208     }
 209     particleTree->GetEntry(nTotal - nPrim + label);
 210
 211     fdNdEtaCorrection->FillMeas(vtxMC[2], particle->Eta());
 212
 213   } // end of track loop
 214
 215   return kTRUE;
 216 }
 217
 218 void AlidNdEtaCorrectionSelector::SlaveTerminate()
 219 {
 220   // The SlaveTerminate() function is called after all entries or objects
 221   // have been processed. When running with PROOF SlaveTerminate() is called
 222   // on each slave server.
 223
 224   AliSelector::SlaveTerminate();
 225
 226   // Add the histograms to the output on each slave server
 227   if (!fOutput)
 228   {
 229     printf("ERROR: Output list not initialized\n");
 230     return;
 231   }
 232
 233   fOutput->Add(fdNdEtaCorrection->GetGeneratedHistogram());
 234   fOutput->Add(fdNdEtaCorrection->GetMeasuredHistogram());
 235 }
 236
 237 void AlidNdEtaCorrectionSelector::Terminate()
 238 {
 239   // The Terminate() function is the last function to be called during
 240   // a query. It always runs on the client, it can be used to present
 241   // the results graphically or save the results to file.
 242
 243   AliSelector::Terminate();
 244
 245   fdNdEtaCorrectionFinal = new dNdEtaCorrection();
 246   TH2F* measuredHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("etaVsVtx_meas"));
 247   TH2F* generatedHistogram = dynamic_cast<TH2F*> (fOutput->FindObject("etaVsVtx_gene"));
 248   if (!measuredHistogram || !generatedHistogram)
 249   {
 250     printf("ERROR: Histograms not available %p %p\n", (void*) generatedHistogram, (void*) measuredHistogram);
 251     return;
 252   }
 253   fdNdEtaCorrectionFinal->SetGeneratedHistogram(generatedHistogram);
 254   fdNdEtaCorrectionFinal->SetMeasuredHistogram(measuredHistogram);
 255
 256   fdNdEtaCorrectionFinal->Finish();
 257
 258   TFile* fout = new TFile("correction_map.root","RECREATE");
 259
 260   fEsdTrackCuts->SaveHistograms("esd_track_cuts");
 261   fdNdEtaCorrectionFinal->SaveHistograms();
 262
 263   fout->Write();
 264   fout->Close();
 265
 266   fdNdEtaCorrectionFinal->DrawHistograms();
 267 }