Added histos in a new output container (Chiara)
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDpidCuts.cxx
CommitLineData
7e1f5ba5 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>
9eeae5d5 11#include <TList.h>
7e1f5ba5 12
13#include "AliESDtrack.h"
10d100d4 14#include "AliESDEvent.h"
7e1f5ba5 15#include "AliLog.h"
10d100d4 16#include "AliESDpid.h"
7e1f5ba5 17
18#include "AliESDpidCuts.h"
19
20ClassImp(AliESDpidCuts)
21
22const Int_t AliESDpidCuts::kNcuts = 3;
23
24//_____________________________________________________________________
25AliESDpidCuts::AliESDpidCuts(const Char_t *name, const Char_t *title):
26 AliAnalysisCuts(name, title)
10d100d4 27 , fESDpid(NULL)
f61ab382 28 , fTPCsigmaCutRequired(0)
29 , fTOFsigmaCutRequired(0)
7e1f5ba5 30 , fCutTPCclusterRatio(0.)
31 , fMinMomentumTOF(0.5)
32 , fHcutStatistics(NULL)
33 , fHcutCorrelation(NULL)
34{
35 //
36 // Default constructor
37 //
38
10d100d4 39 fESDpid = new AliESDpid;
f61ab382 40 memset(fCutTPCnSigma, 0, sizeof(Float_t) * 2);
41 memset(fCutTOFnSigma, 0, sizeof(Float_t) * 2);
7e1f5ba5 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//_____________________________________________________________________
49AliESDpidCuts::AliESDpidCuts(const AliESDpidCuts &ref):
50 AliAnalysisCuts(ref)
10d100d4 51 , fESDpid(NULL)
f61ab382 52 , fTPCsigmaCutRequired(ref.fTPCsigmaCutRequired)
53 , fTOFsigmaCutRequired(ref.fTOFsigmaCutRequired)
7e1f5ba5 54 , fCutTPCclusterRatio(ref.fCutTPCclusterRatio)
55 , fMinMomentumTOF(ref.fMinMomentumTOF)
56 , fHcutStatistics(NULL)
57 , fHcutCorrelation(NULL)
58{
59 //
60 // Copy constructor
61 //
10d100d4 62 fESDpid = new AliESDpid(*ref.fESDpid);
f61ab382 63 memcpy(fCutTPCnSigma, ref.fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
64 memcpy(fCutTOFnSigma, ref.fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
7e1f5ba5 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//_____________________________________________________________________
78AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){
79 //
80 // Assignment operator
81 //
82 if(this != &ref)
83 ref.Copy(*this);
84 return *this;
85}
86
87//_____________________________________________________________________
88AliESDpidCuts::~AliESDpidCuts(){
89 //
90 // Destructor
91 //
10d100d4 92 delete fESDpid;
7e1f5ba5 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//_____________________________________________________________________
95621324 106Bool_t AliESDpidCuts::IsSelected(TObject *obj){
7e1f5ba5 107 //
108 // Select Track
109 //
95621324 110 AliESDtrack * trk = dynamic_cast<AliESDtrack*>(obj);
9eeae5d5 111 if(!trk){
112 AliError("Provided object is not AliESDtrack!");
7e1f5ba5 113 return kFALSE;
114 }
95621324 115 AliESDEvent* evt = trk->GetESDEvent();
9eeae5d5 116 if(!evt){
117 AliError("No AliESDEvent!");
118 return kFALSE;
119 }
95621324 120 return AcceptTrack(trk, evt);
7e1f5ba5 121}
122
123//_____________________________________________________________________
124void AliESDpidCuts::Copy(TObject &c) const {
125 //
126 // Copy function
127 //
128 AliESDpidCuts &target = dynamic_cast<AliESDpidCuts &>(c);
129
10d100d4 130 target.fESDpid = new AliESDpid(*fESDpid);
7e1f5ba5 131
132 target.fCutTPCclusterRatio = fCutTPCclusterRatio;
133 target.fMinMomentumTOF = fMinMomentumTOF;
f61ab382 134
135 target.fTPCsigmaCutRequired = fTPCsigmaCutRequired;
136 target.fTOFsigmaCutRequired = fTOFsigmaCutRequired;
7e1f5ba5 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
f61ab382 148 memcpy(target.fCutTPCnSigma, fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
7e1f5ba5 149 memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES);
150
151 AliESDpidCuts::Copy(c);
152}
153
154//_____________________________________________________________________
155Long64_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//_____________________________________________________________________
187void 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//_____________________________________________________________________
10d100d4 218Bool_t AliESDpidCuts::AcceptTrack(const AliESDtrack *track, const AliESDEvent *event){
7e1f5ba5 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;
10d100d4 228 if(fTOFsigmaCutRequired && event == 0)
229 AliError("No event pointer. Need event pointer for T0 for TOF cut");
7e1f5ba5 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
f61ab382 237 Float_t nsigmaTPC[AliPID::kSPECIES], nsigma; // need all sigmas for QA plotting
7e1f5ba5 238 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
10d100d4 239 nsigmaTPC[ispec] = nsigma = fESDpid->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
f61ab382 240 if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;
7e1f5ba5 241 SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required
f61ab382 242 if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC); // Fullfiled for at least one species
7e1f5ba5 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
10d100d4 247 Double_t times[AliPID::kSPECIES];
248 track->GetIntegratedTimes(times);
7e1f5ba5 249 for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
10d100d4 250
251 if(hasTOFpid) nsigmaTOF[ispec] = nsigma = fESDpid->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec, event->GetT0());
f61ab382 252 if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;
7e1f5ba5 253 SETBIT(cutRequired, kCutNsigmaTOF);
254 if(track->GetOuterParam()->P() >= fMinMomentumTOF){
f61ab382 255 if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);
7e1f5ba5 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//_____________________________________________________________________
295void 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//_____________________________________________________________________
330void 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