Disable caching for async prefetching
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDpidCuts.cxx
CommitLineData
a65a7e70 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
36ClassImp(AliESDpidCuts)
37
38const Int_t AliESDpidCuts::kNcuts = 3;
39
40//_____________________________________________________________________
41AliESDpidCuts::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//_____________________________________________________________________
64AliESDpidCuts::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//_____________________________________________________________________
93AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){
94 //
95 // Assignment operator
96 //
97 if(this != &ref)
98 ref.Copy(*this);
99 return *this;
100}
101
102//_____________________________________________________________________
103AliESDpidCuts::~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//_____________________________________________________________________
119void 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//_____________________________________________________________________
133Bool_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//_____________________________________________________________________
151void 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//_____________________________________________________________________
182Long64_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//_____________________________________________________________________
214void 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//_____________________________________________________________________
245Bool_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//_____________________________________________________________________
330void 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//_____________________________________________________________________
365void 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