fixing warnings
[u/mrichter/AliRoot.git] / 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 "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 }