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