Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDpidCuts.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: The ALICE Off-line Project.                                    *\r
5  * Contributors are mentioned in the code where appropriate.              *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 #include <TCanvas.h>\r
16 #include <TClass.h>\r
17 #include <TCollection.h>\r
18 #include <TDirectory.h>\r
19 #include <TH1F.h>\r
20 #include <TH1I.h>\r
21 #include <TH2I.h>\r
22 #include <TMath.h>\r
23 #include <TIterator.h>\r
24 #include <TString.h>\r
25 #include <TList.h>\r
26 \r
27 #include "AliAnalysisManager.h"\r
28 #include "AliInputEventHandler.h"\r
29 #include "AliESDtrack.h"\r
30 #include "AliESDEvent.h"\r
31 #include "AliLog.h"\r
32 #include "AliPIDResponse.h"\r
33 \r
34 #include "AliESDpidCuts.h"\r
35 \r
36 ClassImp(AliESDpidCuts)\r
37 \r
38 const Int_t AliESDpidCuts::kNcuts = 3;\r
39 \r
40 //_____________________________________________________________________\r
41 AliESDpidCuts::AliESDpidCuts(const Char_t *name, const Char_t *title):\r
42     AliAnalysisCuts(name, title)\r
43   , fPIDresponse(NULL)\r
44   , fTPCsigmaCutRequired(0)\r
45   , fTOFsigmaCutRequired(0)\r
46   , fCutTPCclusterRatio(0.)\r
47   , fMinMomentumTOF(0.5)\r
48   , fHcutStatistics(NULL)\r
49   , fHcutCorrelation(NULL)\r
50 {\r
51   //\r
52   // Default constructor\r
53   //\r
54   \r
55   memset(fCutTPCnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);\r
56   memset(fCutTOFnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);\r
57 \r
58   memset(fHclusterRatio, 0, sizeof(TH1F *) * 2);\r
59   memset(fHnSigmaTPC, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);\r
60   memset(fHnSigmaTOF, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);\r
61 }\r
62 \r
63 //_____________________________________________________________________\r
64 AliESDpidCuts::AliESDpidCuts(const AliESDpidCuts &ref):\r
65     AliAnalysisCuts(ref)\r
66   , fPIDresponse(NULL)\r
67   , fTPCsigmaCutRequired(ref.fTPCsigmaCutRequired)\r
68   , fTOFsigmaCutRequired(ref.fTOFsigmaCutRequired)\r
69   , fCutTPCclusterRatio(ref.fCutTPCclusterRatio)\r
70   , fMinMomentumTOF(ref.fMinMomentumTOF)\r
71   , fHcutStatistics(NULL)\r
72   , fHcutCorrelation(NULL)\r
73 {\r
74   //\r
75   // Copy constructor\r
76   //\r
77   fPIDresponse = ref.fPIDresponse;\r
78   memcpy(fCutTPCnSigma, ref.fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
79   memcpy(fCutTOFnSigma, ref.fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
80   \r
81   if(ref.fHcutStatistics) fHcutStatistics = dynamic_cast<TH1I *>(ref.fHcutStatistics->Clone());\r
82   if(ref.fHcutCorrelation) fHcutCorrelation = dynamic_cast<TH2I *>(ref.fHcutCorrelation->Clone());\r
83   for(Int_t imode = 0; imode < 2; imode++){\r
84     if(ref.fHclusterRatio[imode]) fHclusterRatio[imode] = dynamic_cast<TH1F *>(ref.fHclusterRatio[imode]->Clone());\r
85     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
86       if(fHnSigmaTPC[ispec][imode]) fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());\r
87       if(fHnSigmaTOF[ispec][imode]) fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());\r
88     }\r
89   }\r
90 }\r
91 \r
92 //_____________________________________________________________________\r
93 AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){\r
94   //\r
95   // Assignment operator\r
96   //\r
97   if(this != &ref)\r
98     ref.Copy(*this);\r
99   return *this;\r
100 }\r
101 \r
102 //_____________________________________________________________________\r
103 AliESDpidCuts::~AliESDpidCuts(){\r
104   //\r
105   // Destructor\r
106   //\r
107 \r
108   delete fHcutCorrelation;\r
109   for(Int_t imode = 0; imode < 2; imode++){\r
110     delete fHclusterRatio[imode];\r
111     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
112       delete fHnSigmaTPC[ispec][imode];\r
113       delete fHnSigmaTOF[ispec][imode];\r
114     }\r
115   }\r
116 }\r
117 \r
118 //_____________________________________________________________________\r
119 void AliESDpidCuts::Init(){\r
120   //\r
121   // Init function, get PID response from the Analysis Manager\r
122   //\r
123   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
124   if(mgr){\r
125     AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());\r
126     if(handler){\r
127       fPIDresponse = handler->GetPIDResponse();\r
128     }\r
129   }\r
130 }\r
131 \r
132 //_____________________________________________________________________\r
133 Bool_t AliESDpidCuts::IsSelected(TObject *obj){\r
134   //\r
135   // Select Track\r
136   // \r
137   AliESDtrack * trk = dynamic_cast<AliESDtrack*>(obj);\r
138   if(!trk){\r
139     AliError("Provided object is not AliESDtrack!");\r
140     return kFALSE;\r
141   }\r
142   const AliESDEvent* evt = trk->GetESDEvent();\r
143   if(!evt){\r
144     AliError("No AliESDEvent!");\r
145     return kFALSE;\r
146   }\r
147   return AcceptTrack(trk, evt);\r
148 }\r
149 \r
150 //_____________________________________________________________________\r
151 void AliESDpidCuts::Copy(TObject &c) const {\r
152   //\r
153   // Copy function\r
154   //\r
155   AliESDpidCuts &target = dynamic_cast<AliESDpidCuts &>(c);\r
156 \r
157   target.fPIDresponse = fPIDresponse;\r
158 \r
159   target.fCutTPCclusterRatio = fCutTPCclusterRatio;\r
160   target.fMinMomentumTOF = fMinMomentumTOF;\r
161 \r
162   target.fTPCsigmaCutRequired = fTPCsigmaCutRequired;\r
163   target.fTOFsigmaCutRequired = fTOFsigmaCutRequired;\r
164   \r
165   if(fHcutStatistics) target.fHcutStatistics = dynamic_cast<TH1I *>(fHcutStatistics->Clone());\r
166   if(fHcutCorrelation) target.fHcutCorrelation = dynamic_cast<TH2I *>(fHcutCorrelation->Clone());\r
167   for(Int_t imode = 0; imode < 2; imode++){\r
168     if(fHclusterRatio[imode]) target.fHclusterRatio[imode] = dynamic_cast<TH1F *>(fHclusterRatio[imode]->Clone());\r
169     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
170       if(fHnSigmaTPC[ispec][imode]) target.fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());\r
171       if(fHnSigmaTOF[ispec][imode]) target.fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTOF[ispec][imode]->Clone());\r
172     }\r
173   }\r
174  \r
175   memcpy(target.fCutTPCnSigma, fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
176   memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);\r
177  \r
178   AliESDpidCuts::Copy(c);\r
179 }\r
180 \r
181 //_____________________________________________________________________\r
182 Long64_t AliESDpidCuts::Merge(TCollection *coll){\r
183   //\r
184   // Merge Cut objects\r
185   //\r
186   if(!coll) return 0;\r
187   if(coll->IsEmpty())   return 1;\r
188   if(!HasHistograms())  return 0;\r
189   \r
190   TIterator *iter = coll->MakeIterator();\r
191   TObject *o = NULL; \r
192   AliESDpidCuts *ref = NULL;\r
193   Int_t counter = 0;\r
194   while((o = iter->Next())){\r
195     ref = dynamic_cast<AliESDpidCuts *>(o);\r
196     if(!ref) continue;\r
197     if(!ref->HasHistograms()) continue;\r
198 \r
199     fHcutStatistics->Add(ref->fHcutStatistics);\r
200     fHcutCorrelation->Add(ref->fHcutCorrelation);\r
201     for(Int_t imode = 0; imode < 2; imode++){\r
202       fHclusterRatio[imode]->Add(ref->fHclusterRatio[imode]);\r
203       for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
204         fHnSigmaTPC[ispec][imode]->Add(ref->fHnSigmaTPC[ispec][imode]);\r
205         fHnSigmaTOF[ispec][imode]->Add(ref->fHnSigmaTOF[ispec][imode]);\r
206       }\r
207     }\r
208     ++counter;\r
209   }\r
210   return ++counter;\r
211 }\r
212 \r
213 //_____________________________________________________________________\r
214 void AliESDpidCuts::DefineHistograms(Color_t color){\r
215   //\r
216   // Swich on QA and create the histograms\r
217   //\r
218   SetBit(kHasHistograms, kTRUE);\r
219   fHcutStatistics = new TH1I("fHcutStatistics", "Cut Statistics", kNcuts, 0, kNcuts);\r
220   fHcutStatistics->SetLineColor(color);\r
221   fHcutCorrelation = new TH2I("fHcutCorrelation", "Cut Correlation", kNcuts, 0, kNcuts, kNcuts, 0, kNcuts);\r
222   TString cutname[kNcuts] = {"TPCclusterRatio", "TPC sigma", "TOF sigma"};\r
223   for(Int_t icut = 0; icut < kNcuts; icut++){\r
224     fHcutStatistics->GetXaxis()->SetBinLabel(fHcutStatistics->GetXaxis()->GetFirst() + icut, cutname[icut].Data());\r
225     fHcutCorrelation->GetXaxis()->SetBinLabel(fHcutCorrelation->GetXaxis()->GetFirst() + icut, cutname[icut].Data());\r
226     fHcutCorrelation->GetYaxis()->SetBinLabel(fHcutCorrelation->GetYaxis()->GetFirst() + icut, cutname[icut].Data());\r
227   }\r
228   Char_t hname[256], htitle[256];\r
229   for(Int_t imode = 0; imode < 2; imode++){\r
230     snprintf(hname, 256, "fHclusterRatio%s", imode ? "After" : "Before");\r
231     snprintf(htitle, 256, "TPC cluster Ratio %s cuts;Ratio;Entries", imode ? "after" : "before");\r
232     fHclusterRatio[imode] = new TH1F(hname, htitle, 20, 0., 1.);\r
233     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
234       snprintf(hname, 256, "fHnSigma%sTPC%s", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
235       snprintf(htitle, 256, "TPC sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
236       fHnSigmaTPC[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);\r
237       snprintf(hname, 256, "fHnSigma%sTOF%s", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
238       snprintf(htitle, 256, "TOF sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");\r
239       fHnSigmaTOF[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);\r
240     }\r
241   }\r
242 }\r
243 \r
244 //_____________________________________________________________________\r
245 Bool_t AliESDpidCuts::AcceptTrack(const AliESDtrack *track, const AliESDEvent *event){\r
246   //\r
247   // Check whether the tracks survived the cuts\r
248   //\r
249   if(!fPIDresponse){\r
250     Init();\r
251     if (!fPIDresponse)\r
252       AliError("PID Response not available");\r
253     return 0;\r
254   }\r
255   enum{\r
256     kCutClusterRatioTPC,\r
257     kCutNsigmaTPC,\r
258     kCutNsigmaTOF\r
259   };\r
260   Long64_t cutRequired=0, cutFullfiled = 0;\r
261   if(fTOFsigmaCutRequired && event == 0)  {\r
262     AliError("No event pointer. Need event pointer for T0 for TOF cut");\r
263     return (0);\r
264   }\r
265   Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;\r
266   if(fCutTPCclusterRatio > 0.){\r
267     SETBIT(cutRequired, kCutClusterRatioTPC);\r
268     if(clusterRatio >= fCutTPCclusterRatio) \r
269       SETBIT(cutFullfiled, kCutClusterRatioTPC);\r
270   }\r
271   // check TPC nSigma cut\r
272   Float_t nsigmaTPC[AliPID::kSPECIES], nsigma;   // need all sigmas for QA plotting\r
273   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
274     nsigmaTPC[ispec] = nsigma = fPIDresponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);\r
275     if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;\r
276     SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required\r
277     if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC);    // Fullfiled for at least one species\r
278   }\r
279   // check TOF nSigma cut\r
280   Float_t nsigmaTOF[AliPID::kSPECIES];    // see above\r
281   Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set\r
282   Double_t times[AliPID::kSPECIES];\r
283   track->GetIntegratedTimes(times);\r
284   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
285     \r
286     if(hasTOFpid && event) nsigmaTOF[ispec] = nsigma = fPIDresponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);\r
287     if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;\r
288     SETBIT(cutRequired, kCutNsigmaTOF);\r
289     if(track->GetOuterParam()->P() >= fMinMomentumTOF){\r
290       if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);\r
291     }\r
292   }\r
293 \r
294   // Fill Histograms before cuts\r
295   if(HasHistograms()){\r
296     fHclusterRatio[0]->Fill(clusterRatio);\r
297     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
298       fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);\r
299       if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);\r
300     }\r
301   }\r
302   if(cutRequired != cutFullfiled){\r
303     // Fill cut statistics\r
304     if(HasHistograms()){\r
305       for(Int_t icut = 0; icut < kNcuts; icut++){\r
306         if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){\r
307           // cut not fullfiled\r
308           fHcutStatistics->Fill(icut);\r
309           for(Int_t jcut = 0; jcut <= icut; jcut++)\r
310             if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);\r
311         }\r
312       }\r
313     }\r
314     return kFALSE;    // At least one cut is not fullfiled\r
315   }\r
316 \r
317   // Fill Histograms after cuts\r
318   if(HasHistograms()){\r
319     fHclusterRatio[1]->Fill(clusterRatio);\r
320     for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
321       fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);\r
322       if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);\r
323     }\r
324   }\r
325 \r
326   return kTRUE;\r
327 }\r
328 \r
329 //_____________________________________________________________________\r
330 void AliESDpidCuts::SaveHistograms(const Char_t * location){\r
331   //\r
332   // Save the histograms to a file\r
333   //\r
334   if(!HasHistograms()){\r
335     AliError("Histograms not on - Exiting");\r
336     return;\r
337   }\r
338   if(!location) location = GetName();\r
339   gDirectory->mkdir(location);\r
340   gDirectory->cd(location);\r
341   fHcutStatistics->Write();\r
342   fHcutCorrelation->Write();\r
343 \r
344   gDirectory->mkdir("before_cuts");\r
345   gDirectory->mkdir("after_cuts");\r
346 \r
347   gDirectory->cd("before_cuts");\r
348   fHclusterRatio[0]->Write();\r
349   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
350     fHnSigmaTPC[ispec][0]->Write();\r
351     fHnSigmaTOF[ispec][0]->Write();\r
352   }\r
353 \r
354   gDirectory->cd("../after_cuts");\r
355   fHclusterRatio[1]->Write();\r
356   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
357     fHnSigmaTPC[ispec][1]->Write();\r
358     fHnSigmaTOF[ispec][1]->Write();\r
359   }\r
360 \r
361   gDirectory->cd("..");\r
362 }\r
363 \r
364 //_____________________________________________________________________\r
365 void AliESDpidCuts::DrawHistograms(){\r
366   //\r
367   // Draw the Histograms\r
368   //\r
369   TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);\r
370   stat->cd();\r
371   fHcutStatistics->SetStats(kFALSE);\r
372   fHcutStatistics->Draw();\r
373   stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));\r
374 \r
375   TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);\r
376   correl->cd();\r
377   fHcutCorrelation->SetStats(kFALSE);\r
378   fHcutCorrelation->Draw("colz");\r
379   correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));\r
380 \r
381   TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);\r
382   cRatio->cd();\r
383   fHclusterRatio[0]->SetLineColor(kRed);\r
384   fHclusterRatio[0]->SetStats(kFALSE);\r
385   fHclusterRatio[0]->Draw();\r
386   fHclusterRatio[1]->SetLineColor(kBlue);\r
387   fHclusterRatio[1]->SetStats(kFALSE);\r
388   fHclusterRatio[1]->Draw("same");\r
389   cRatio->SaveAs(Form("%s_%s.gif",  GetName(), cRatio->GetName()));\r
390 \r
391   TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);\r
392   cNsigTPC->Divide(3,2);\r
393   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
394     cNsigTPC->cd(ispec + 1);\r
395     fHnSigmaTPC[ispec][0]->SetLineColor(kRed);\r
396     fHnSigmaTPC[ispec][0]->SetStats(kFALSE);\r
397     fHnSigmaTPC[ispec][0]->Draw();\r
398     fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);\r
399     fHnSigmaTPC[ispec][1]->SetStats(kFALSE);\r
400     fHnSigmaTPC[ispec][1]->Draw("same");\r
401   }\r
402   cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));\r
403 \r
404   TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);\r
405   cNsigTOF->Divide(3,2);\r
406   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){\r
407     cNsigTOF->cd(ispec + 1);\r
408     fHnSigmaTOF[ispec][0]->SetLineColor(kRed);\r
409     fHnSigmaTOF[ispec][0]->SetStats(kFALSE);\r
410     fHnSigmaTOF[ispec][0]->Draw();\r
411     fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);\r
412     fHnSigmaTOF[ispec][1]->SetStats(kFALSE);\r
413     fHnSigmaTOF[ispec][1]->Draw("same");\r
414   }\r
415   cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));\r
416 }\r
417 \r