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