8abc3c69f9d3eb24f91eec8e9d3b75506ebb61de
[u/mrichter/AliRoot.git] / 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   // check if the event was triggered
138   Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD);
139
140   fdNdEtaCorrection->IncreaseEventCount();
141   if (eventTriggered)
142     fdNdEtaCorrection->IncreaseTriggeredEventCount();
143
144   // get the MC vertex
145   AliGenEventHeader* genHeader = header->GenEventHeader();
146
147   TArrayF vtxMC(3);
148   genHeader->PrimaryVertex(vtxMC);
149
150   // loop over mc particles
151   Int_t nPrim  = stack->GetNprimary();
152
153   for (Int_t iMc = 0; iMc < nPrim; ++iMc)
154   {
155     AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
156
157     TParticle* particle = stack->Particle(iMc);
158
159     if (!particle)
160     {
161       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
162       continue;
163     }
164
165     if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
166       continue;
167
168     Float_t eta = particle->Eta();
169     Float_t pt = particle->Pt();
170
171     if (vertexReconstructed)
172       fdNdEtaCorrection->FillParticle(vtxMC[2], eta, pt);
173
174     fdNdEtaCorrection->FillParticleAllEvents(eta, pt);
175     if (eventTriggered)
176       fdNdEtaCorrection->FillParticleWhenEventTriggered(eta, pt);
177   }// end of mc particle
178
179   // ########################################################
180   // loop over esd tracks
181   Int_t nTracks = fESD->GetNumberOfTracks();
182
183   // count the number of "good" tracks for vertex reconstruction efficiency
184   // TODO change to number of ITS clusters or similar
185   Int_t nGoodTracks = 0;
186   for (Int_t t=0; t<nTracks; t++)
187   {
188     AliDebug(AliLog::kDebug+1, Form("ESD Loop: Processing track %d.", t));
189
190     AliESDtrack* esdTrack = fESD->GetTrack(t);
191
192     // cut the esd track?
193     if (!fEsdTrackCuts->AcceptTrack(esdTrack))
194       continue;
195
196     nGoodTracks++;
197
198     // using the properties of the mc particle
199     Int_t label = TMath::Abs(esdTrack->GetLabel());
200     if (label == 0)
201     {
202       AliDebug(AliLog::kWarning, Form("WARNING: cannot find corresponding mc part for track %d.", t));
203       continue;
204     }
205
206     TParticle* particle = stack->Particle(label);
207     if (!particle)
208     {
209       AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (track loop).", label));
210       continue;
211     }
212
213     if (vertexReconstructed)
214       fdNdEtaCorrection->FillParticleWhenMeasuredTrack(vtxMC[2], particle->Eta(), particle->Pt());
215   } // end of track loop
216
217   fdNdEtaCorrection->FillEvent(vtxMC[2], nGoodTracks);
218   if (vertexReconstructed)
219     fdNdEtaCorrection->FillEventWithReconstructedVertex(vtxMC[2], nGoodTracks);
220
221   return kTRUE;
222 }
223
224 void AlidNdEtaCorrectionSelector::SlaveTerminate()
225 {
226   // The SlaveTerminate() function is called after all entries or objects
227   // have been processed. When running with PROOF SlaveTerminate() is called
228   // on each slave server.
229
230   AliSelectorRL::SlaveTerminate();
231
232   // Add the histograms to the output on each slave server
233   if (!fOutput)
234   {
235     AliDebug(AliLog::kError, "ERROR: Output list not initialized");
236     return;
237   }
238
239   fOutput->Add(fdNdEtaCorrection);
240 }
241
242 void AlidNdEtaCorrectionSelector::Terminate()
243 {
244   // The Terminate() function is the last function to be called during
245   // a query. It always runs on the client, it can be used to present
246   // the results graphically or save the results to file.
247
248   AliSelectorRL::Terminate();
249
250   fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fOutput->FindObject("dndeta_correction"));
251
252   fdNdEtaCorrection->Finish();
253
254   TFile* fout = new TFile("correction_map.root","RECREATE");
255
256   fEsdTrackCuts->SaveHistograms("esd_track_cuts");
257   fdNdEtaCorrection->SaveHistograms();
258
259   fout->Write();
260   fout->Close();
261
262   fdNdEtaCorrection->DrawHistograms();
263 }