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