1 /**************************************************************************
\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
\r
4 * Author: The ALICE Off-line Project. *
\r
5 * Contributors are mentioned in the code where appropriate. *
\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
17 #include <TCollection.h>
\r
18 #include <TDirectory.h>
\r
23 #include <TIterator.h>
\r
24 #include <TString.h>
\r
27 #include "AliAnalysisManager.h"
\r
28 #include "AliInputEventHandler.h"
\r
29 #include "AliESDtrack.h"
\r
30 #include "AliESDEvent.h"
\r
32 #include "AliPIDResponse.h"
\r
34 #include "AliESDpidCuts.h"
\r
36 ClassImp(AliESDpidCuts)
\r
38 const Int_t AliESDpidCuts::kNcuts = 3;
\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
52 // Default constructor
\r
55 memset(fCutTPCnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);
\r
56 memset(fCutTOFnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);
\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
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
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
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
92 //_____________________________________________________________________
\r
93 AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){
\r
95 // Assignment operator
\r
102 //_____________________________________________________________________
\r
103 AliESDpidCuts::~AliESDpidCuts(){
\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
118 //_____________________________________________________________________
\r
119 void AliESDpidCuts::Init(){
\r
121 // Init function, get PID response from the Analysis Manager
\r
123 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
\r
125 AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
\r
127 fPIDresponse = handler->GetPIDResponse();
\r
132 //_____________________________________________________________________
\r
133 Bool_t AliESDpidCuts::IsSelected(TObject *obj){
\r
137 AliESDtrack * trk = dynamic_cast<AliESDtrack*>(obj);
\r
139 AliError("Provided object is not AliESDtrack!");
\r
142 const AliESDEvent* evt = trk->GetESDEvent();
\r
144 AliError("No AliESDEvent!");
\r
147 return AcceptTrack(trk, evt);
\r
150 //_____________________________________________________________________
\r
151 void AliESDpidCuts::Copy(TObject &c) const {
\r
155 AliESDpidCuts &target = dynamic_cast<AliESDpidCuts &>(c);
\r
157 target.fPIDresponse = fPIDresponse;
\r
159 target.fCutTPCclusterRatio = fCutTPCclusterRatio;
\r
160 target.fMinMomentumTOF = fMinMomentumTOF;
\r
162 target.fTPCsigmaCutRequired = fTPCsigmaCutRequired;
\r
163 target.fTOFsigmaCutRequired = fTOFsigmaCutRequired;
\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
175 memcpy(target.fCutTPCnSigma, fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
\r
176 memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
\r
178 AliESDpidCuts::Copy(c);
\r
181 //_____________________________________________________________________
\r
182 Long64_t AliESDpidCuts::Merge(TCollection *coll){
\r
184 // Merge Cut objects
\r
186 if(!coll) return 0;
\r
187 if(coll->IsEmpty()) return 1;
\r
188 if(!HasHistograms()) return 0;
\r
190 TIterator *iter = coll->MakeIterator();
\r
191 TObject *o = NULL;
\r
192 AliESDpidCuts *ref = NULL;
\r
194 while((o = iter->Next())){
\r
195 ref = dynamic_cast<AliESDpidCuts *>(o);
\r
197 if(!ref->HasHistograms()) continue;
\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
213 //_____________________________________________________________________
\r
214 void AliESDpidCuts::DefineHistograms(Color_t color){
\r
216 // Swich on QA and create the histograms
\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
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
244 //_____________________________________________________________________
\r
245 Bool_t AliESDpidCuts::AcceptTrack(const AliESDtrack *track, const AliESDEvent *event){
\r
247 // Check whether the tracks survived the cuts
\r
250 AliError("PID Response not available");
\r
254 kCutClusterRatioTPC,
\r
258 Long64_t cutRequired=0, cutFullfiled = 0;
\r
259 if(fTOFsigmaCutRequired && event == 0) {
\r
260 AliError("No event pointer. Need event pointer for T0 for TOF cut");
\r
263 Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
\r
264 if(fCutTPCclusterRatio > 0.){
\r
265 SETBIT(cutRequired, kCutClusterRatioTPC);
\r
266 if(clusterRatio >= fCutTPCclusterRatio)
\r
267 SETBIT(cutFullfiled, kCutClusterRatioTPC);
\r
269 // check TPC nSigma cut
\r
270 Float_t nsigmaTPC[AliPID::kSPECIES], nsigma; // need all sigmas for QA plotting
\r
271 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
272 nsigmaTPC[ispec] = nsigma = fPIDresponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
\r
273 if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;
\r
274 SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required
\r
275 if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC); // Fullfiled for at least one species
\r
277 // check TOF nSigma cut
\r
278 Float_t nsigmaTOF[AliPID::kSPECIES]; // see above
\r
279 Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set
\r
280 Double_t times[AliPID::kSPECIES];
\r
281 track->GetIntegratedTimes(times);
\r
282 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
284 if(hasTOFpid && event) nsigmaTOF[ispec] = nsigma = fPIDresponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);
\r
285 if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;
\r
286 SETBIT(cutRequired, kCutNsigmaTOF);
\r
287 if(track->GetOuterParam()->P() >= fMinMomentumTOF){
\r
288 if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);
\r
292 // Fill Histograms before cuts
\r
293 if(HasHistograms()){
\r
294 fHclusterRatio[0]->Fill(clusterRatio);
\r
295 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
296 fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);
\r
297 if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);
\r
300 if(cutRequired != cutFullfiled){
\r
301 // Fill cut statistics
\r
302 if(HasHistograms()){
\r
303 for(Int_t icut = 0; icut < kNcuts; icut++){
\r
304 if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){
\r
305 // cut not fullfiled
\r
306 fHcutStatistics->Fill(icut);
\r
307 for(Int_t jcut = 0; jcut <= icut; jcut++)
\r
308 if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);
\r
312 return kFALSE; // At least one cut is not fullfiled
\r
315 // Fill Histograms after cuts
\r
316 if(HasHistograms()){
\r
317 fHclusterRatio[1]->Fill(clusterRatio);
\r
318 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
319 fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);
\r
320 if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);
\r
327 //_____________________________________________________________________
\r
328 void AliESDpidCuts::SaveHistograms(const Char_t * location){
\r
330 // Save the histograms to a file
\r
332 if(!HasHistograms()){
\r
333 AliError("Histograms not on - Exiting");
\r
336 if(!location) location = GetName();
\r
337 gDirectory->mkdir(location);
\r
338 gDirectory->cd(location);
\r
339 fHcutStatistics->Write();
\r
340 fHcutCorrelation->Write();
\r
342 gDirectory->mkdir("before_cuts");
\r
343 gDirectory->mkdir("after_cuts");
\r
345 gDirectory->cd("before_cuts");
\r
346 fHclusterRatio[0]->Write();
\r
347 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
348 fHnSigmaTPC[ispec][0]->Write();
\r
349 fHnSigmaTOF[ispec][0]->Write();
\r
352 gDirectory->cd("../after_cuts");
\r
353 fHclusterRatio[1]->Write();
\r
354 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
355 fHnSigmaTPC[ispec][1]->Write();
\r
356 fHnSigmaTOF[ispec][1]->Write();
\r
359 gDirectory->cd("..");
\r
362 //_____________________________________________________________________
\r
363 void AliESDpidCuts::DrawHistograms(){
\r
365 // Draw the Histograms
\r
367 TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);
\r
369 fHcutStatistics->SetStats(kFALSE);
\r
370 fHcutStatistics->Draw();
\r
371 stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));
\r
373 TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);
\r
375 fHcutCorrelation->SetStats(kFALSE);
\r
376 fHcutCorrelation->Draw("colz");
\r
377 correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));
\r
379 TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);
\r
381 fHclusterRatio[0]->SetLineColor(kRed);
\r
382 fHclusterRatio[0]->SetStats(kFALSE);
\r
383 fHclusterRatio[0]->Draw();
\r
384 fHclusterRatio[1]->SetLineColor(kBlue);
\r
385 fHclusterRatio[1]->SetStats(kFALSE);
\r
386 fHclusterRatio[1]->Draw("same");
\r
387 cRatio->SaveAs(Form("%s_%s.gif", GetName(), cRatio->GetName()));
\r
389 TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);
\r
390 cNsigTPC->Divide(3,2);
\r
391 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
392 cNsigTPC->cd(ispec + 1);
\r
393 fHnSigmaTPC[ispec][0]->SetLineColor(kRed);
\r
394 fHnSigmaTPC[ispec][0]->SetStats(kFALSE);
\r
395 fHnSigmaTPC[ispec][0]->Draw();
\r
396 fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);
\r
397 fHnSigmaTPC[ispec][1]->SetStats(kFALSE);
\r
398 fHnSigmaTPC[ispec][1]->Draw("same");
\r
400 cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));
\r
402 TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);
\r
403 cNsigTOF->Divide(3,2);
\r
404 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
\r
405 cNsigTOF->cd(ispec + 1);
\r
406 fHnSigmaTOF[ispec][0]->SetLineColor(kRed);
\r
407 fHnSigmaTOF[ispec][0]->SetStats(kFALSE);
\r
408 fHnSigmaTOF[ispec][0]->Draw();
\r
409 fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);
\r
410 fHnSigmaTOF[ispec][1]->SetStats(kFALSE);
\r
411 fHnSigmaTOF[ispec][1]->Draw("same");
\r
413 cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));
\r