updating makefile to include evgen headers to compile systematics selector
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaAnalysisESDSelector.cxx
1 /* $Id$ */
2
3 #include "AlidNdEtaAnalysisESDSelector.h"
4
5 #include <TStyle.h>
6 #include <TSystem.h>
7 #include <TCanvas.h>
8 #include <TVector3.h>
9 #include <TChain.h>
10 #include <TFile.h>
11
12 #include <AliLog.h>
13 #include <AliESDVertex.h>
14 #include <AliESD.h>
15
16 #include "esdTrackCuts/AliESDtrackCuts.h"
17 #include "dNdEta/dNdEtaAnalysis.h"
18 #include "dNdEta/AlidNdEtaCorrection.h"
19 #include "AliPWG0Helper.h"
20
21 ClassImp(AlidNdEtaAnalysisESDSelector)
22
23 AlidNdEtaAnalysisESDSelector::AlidNdEtaAnalysisESDSelector() :
24   AliSelector(),
25   fdNdEtaAnalysisMBVtx(0),
26   fdNdEtaAnalysisMB(0),
27   fdNdEtaAnalysis(0),
28   fEsdTrackCuts(0),
29   fdNdEtaCorrection(0)
30   {
31   //
32   // Constructor. Initialization of pointers
33   //
34
35   AliLog::SetClassDebugLevel("AlidNdEtaAnalysisESDSelector", AliLog::kDebug);
36 }
37
38 AlidNdEtaAnalysisESDSelector::~AlidNdEtaAnalysisESDSelector()
39 {
40   //
41   // Destructor
42   //
43
44   // histograms are in the output list and deleted when the output
45   // list is deleted by the TSelector dtor
46 }
47
48 void AlidNdEtaAnalysisESDSelector::Begin(TTree* tree)
49 {
50   // Begin function
51
52   ReadUserObjects(tree);
53 }
54
55 void AlidNdEtaAnalysisESDSelector::ReadUserObjects(TTree* tree)
56 {
57   // read the user objects, called from slavebegin and begin
58
59   if (!fEsdTrackCuts && fInput)
60     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fInput->FindObject("AliESDtrackCuts"));
61
62   if (!fEsdTrackCuts && tree)
63     fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (tree->GetUserInfo()->FindObject("AliESDtrackCuts"));
64
65   if (!fEsdTrackCuts)
66      AliDebug(AliLog::kError, "ERROR: Could not read EsdTrackCuts from input list.");
67
68
69   if (!fdNdEtaCorrection && fInput)
70     fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (fInput->FindObject("dndeta_correction"));
71
72   if (!fdNdEtaCorrection && tree)
73     fdNdEtaCorrection = dynamic_cast<AlidNdEtaCorrection*> (tree->GetUserInfo()->FindObject("dndeta_correction"));
74
75   if (!fdNdEtaCorrection)
76      AliDebug(AliLog::kError, "ERROR: Could not read dndeta_correction from input list.");
77 }
78
79 void AlidNdEtaAnalysisESDSelector::SlaveBegin(TTree* tree)
80 {
81   // The SlaveBegin() function is called after the Begin() function.
82   // When running with PROOF SlaveBegin() is called on each slave server.
83   // The tree argument is deprecated (on PROOF 0 is passed).
84
85   AliSelector::SlaveBegin(tree);
86
87   ReadUserObjects(tree);
88
89   fdNdEtaAnalysisMBVtx = new dNdEtaAnalysis("dndeta_mbvtx", "dndeta_mbvtx");
90   fdNdEtaAnalysisMB = new dNdEtaAnalysis("dndeta_mb", "dndeta_mb");
91   fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta");
92 }
93
94 void AlidNdEtaAnalysisESDSelector::Init(TTree* tree)
95 {
96   // read the user objects
97
98   AliSelector::Init(tree);
99
100   // Enable only the needed branches
101   if (tree)
102   {
103     tree->SetBranchStatus("*", 0);
104     tree->SetBranchStatus("fTriggerMask", 1);
105     tree->SetBranchStatus("fSPDVertex*", 1);
106
107     AliESDtrackCuts::EnableNeededBranches(tree);
108   }
109 }
110
111 Bool_t AlidNdEtaAnalysisESDSelector::Process(Long64_t entry)
112 {
113   // The Process() function is called for each entry in the tree (or possibly
114   // keyed object in the case of PROOF) to be processed. The entry argument
115   // specifies which entry in the currently loaded tree is to be processed.
116   // It can be passed to either TTree::GetEntry() or TBranch::GetEntry()
117   // to read either all or the required parts of the data. When processing
118   // keyed objects with PROOF, the object is already loaded and is available
119   // via the fObject pointer.
120   //
121   // This function should contain the "body" of the analysis. It can contain
122   // simple or elaborate selection criteria, run algorithms on the data
123   // of the event and typically fill histograms.
124
125   // WARNING when a selector is used with a TChain, you must use
126   //  the pointer to the current TTree to call GetEntry(entry).
127   //  The entry is always the local entry number in the current tree.
128   //  Assuming that fTree is the pointer to the TChain being processed,
129   //  use fTree->GetTree()->GetEntry(entry).
130
131   if (AliSelector::Process(entry) == kFALSE)
132     return kFALSE;
133
134   // Check prerequisites
135   if (!fESD)
136   {
137     AliDebug(AliLog::kError, "ESD branch not available");
138     return kFALSE;
139   }
140
141   if (!fEsdTrackCuts)
142   {
143     AliDebug(AliLog::kError, "fESDTrackCuts not available");
144     return kFALSE;
145   }
146
147   if (!fdNdEtaCorrection)
148   {
149     AliDebug(AliLog::kError, "fdNdEtaCorrection not available");
150     return kFALSE;
151   }
152
153   if (AliPWG0Helper::IsEventTriggered(fESD) == kFALSE)
154   {
155     AliDebug(AliLog::kDebug+1, Form("Skipping event %d because it was not triggered", (Int_t) entry));
156     return kTRUE;
157   }
158
159   if (AliPWG0Helper::IsVertexReconstructed(fESD) == kFALSE)
160   {
161     AliDebug(AliLog::kDebug+1, Form("Skipping event %d because its vertex was not reconstructed", (Int_t) entry));
162     return kTRUE;
163   }
164
165   // ########################################################
166   // get the EDS vertex
167   const AliESDVertex* vtxESD = fESD->GetVertex();
168   Double_t vtx[3];
169   vtxESD->GetXYZ(vtx);
170
171   // get number of "good" tracks
172   TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
173   Int_t nGoodTracks = list->GetEntries();
174
175   Float_t vertexRecoCorr = fdNdEtaCorrection->GetVertexRecoCorrection(vtx[2], nGoodTracks);
176   if (vertexRecoCorr <= 0)
177   {
178     AliDebug(AliLog::kError, Form("INFO: Skipping event because vertexRecoCorr is <= 0 (%f)", vertexRecoCorr));
179     delete list;
180     return kTRUE;
181   }
182
183   Float_t triggerCorr = fdNdEtaCorrection->GetTriggerCorrection(vtx[2], nGoodTracks);
184   if (triggerCorr <= 0)
185   {
186     AliDebug(AliLog::kError, Form("INFO: Skipping event because triggerCorr is <= 0 (%f)", triggerCorr));
187     delete list;
188     return kTRUE;
189   }
190
191   // loop over esd tracks
192   for (Int_t t=0; t<nGoodTracks; t++)
193   {
194     AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(t));
195     if (!esdTrack)
196     {
197       AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", t));
198       continue;
199     }
200
201     Double_t p[3];
202     esdTrack->GetConstrainedPxPyPz(p); // ### TODO should be okay because we have a vertex, however GetInnerPxPyPy / GetOuterPxPyPy also exist
203     TVector3 vector(p);
204
205     Float_t theta = vector.Theta();
206     Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
207     Float_t pt = vector.Pt();
208
209     Float_t track2particleCorr = fdNdEtaCorrection->GetTrack2ParticleCorrection(vtx[2], eta, pt);
210
211     Float_t weight = track2particleCorr * vertexRecoCorr * triggerCorr;
212     if (weight <= 0)
213     {
214       AliDebug(AliLog::kError, Form("INFO: Skipping track because weight is <= 0 (track %d, weight %f) (vtx %f, eta %f, pt %f)", t, weight, vtx[2], eta, pt));
215       continue;
216     }
217
218     fdNdEtaAnalysisMBVtx->FillTrack(vtx[2], eta, pt, track2particleCorr);
219     fdNdEtaAnalysisMB->FillTrack(vtx[2], eta, pt, track2particleCorr * vertexRecoCorr);
220     fdNdEtaAnalysis->FillTrack(vtx[2], eta, pt, weight);
221   } // end of track loop
222
223   delete list;
224   list = 0;
225
226   // for event count per vertex
227   fdNdEtaAnalysisMBVtx->FillEvent(vtx[2], 1);
228   fdNdEtaAnalysisMB->FillEvent(vtx[2], vertexRecoCorr);
229   fdNdEtaAnalysis->FillEvent(vtx[2], vertexRecoCorr * triggerCorr);
230
231   return kTRUE;
232 }
233
234 void AlidNdEtaAnalysisESDSelector::SlaveTerminate()
235 {
236   // The SlaveTerminate() function is called after all entries or objects
237   // have been processed. When running with PROOF SlaveTerminate() is called
238   // on each slave server.
239
240   AliSelector::SlaveTerminate();
241
242   // Add the histograms to the output on each slave server
243   if (!fOutput)
244   {
245     AliDebug(AliLog::kError, Form("ERROR: Output list not initialized."));
246     return;
247   }
248
249   // Add the objects to the output list and set them to 0, so that the destructor does not delete them.
250
251   fOutput->Add(fdNdEtaAnalysis);
252   fdNdEtaAnalysis = 0;
253
254   fOutput->Add(fdNdEtaAnalysisMB);
255   fdNdEtaAnalysisMB = 0;
256
257   fOutput->Add(fdNdEtaAnalysisMBVtx);
258   fdNdEtaAnalysisMBVtx = 0;
259 }
260
261 void AlidNdEtaAnalysisESDSelector::Terminate()
262 {
263   // The Terminate() function is the last function to be called during
264   // a query. It always runs on the client, it can be used to present
265   // the results graphically or save the results to file.
266
267   AliSelector::Terminate();
268
269   fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
270   fdNdEtaAnalysisMB = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta_mb"));
271   fdNdEtaAnalysisMBVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta_mbvtx"));
272
273   if (!fdNdEtaAnalysis || !fdNdEtaAnalysisMB || !fdNdEtaAnalysisMBVtx)
274   {
275     AliDebug(AliLog::kError, Form("ERROR: Histograms not available %p %p %p", (void*) fdNdEtaAnalysis, (void*) fdNdEtaAnalysisMB, (void*) fdNdEtaAnalysisMBVtx));
276     return;
277   }
278
279   if (fdNdEtaAnalysis)
280     fdNdEtaAnalysis->Finish(fdNdEtaCorrection, 0.3);
281
282   if (fdNdEtaAnalysisMB)
283     fdNdEtaAnalysisMB->Finish(fdNdEtaCorrection, 0.3);
284
285   if (fdNdEtaAnalysisMBVtx)
286     fdNdEtaAnalysisMBVtx->Finish(fdNdEtaCorrection, 0.3);
287
288   TFile* fout = new TFile("analysis_esd.root","RECREATE");
289
290   if (fdNdEtaAnalysis)
291     fdNdEtaAnalysis->SaveHistograms();
292
293   if (fdNdEtaAnalysisMB)
294     fdNdEtaAnalysisMB->SaveHistograms();
295
296   if (fdNdEtaAnalysisMBVtx)
297     fdNdEtaAnalysisMBVtx->SaveHistograms();
298
299   if (fEsdTrackCuts)
300     fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
301
302   if (fdNdEtaCorrection)
303     fdNdEtaCorrection->SaveHistograms();
304
305   fout->Write();
306   fout->Close();
307 }