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