applying fix for bug https://savannah.cern.ch/bugs/index.php?71490
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpostAnalysis.cxx
CommitLineData
d2af20c5 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//
16// Post analysis code
17// Drawing nice pictures containing
18// - Efficiency
19// - Signal/Background
20// - PID Performance
21// More Post Analysis code will be added in time
22//
23// Autor: Markus Fasel
24//
25#include <TAxis.h>
26#include <TCanvas.h>
27#include <TFile.h>
28#include <TH1D.h>
29#include <TLegend.h>
30#include <TList.h>
31#include "AliCFContainer.h"
32#include "AliCFEffGrid.h"
33
34#include "AliHFEcuts.h"
35#include "AliHFEpostAnalysis.h"
36
37ClassImp(AliHFEpostAnalysis)
38
39//____________________________________________________________
40AliHFEpostAnalysis::AliHFEpostAnalysis():
41 TObject(),
42 fResults(NULL),
43 fAnalysisObjects(0),
44 fEfficiencyContainer(NULL),
45 fPIDperformance(NULL),
46 fSignalToBackgroundMC(NULL)
47{
48 //
49 // Default Constructor
50 //
51}
52
53//____________________________________________________________
54AliHFEpostAnalysis::AliHFEpostAnalysis(const AliHFEpostAnalysis &ref):
55 TObject(ref),
56 fResults(ref.fResults),
57 fAnalysisObjects(ref.fAnalysisObjects),
58 fEfficiencyContainer(ref.fEfficiencyContainer),
59 fPIDperformance(ref.fPIDperformance),
60 fSignalToBackgroundMC(ref.fSignalToBackgroundMC)
61{
62 //
63 // Copy Constructor
64 //
65}
66
67//____________________________________________________________
68AliHFEpostAnalysis& AliHFEpostAnalysis::operator=(const AliHFEpostAnalysis &ref){
69 //
70 // Assignment Operator
71 //
72 TObject::operator=(ref);
73 fResults = ref.fResults;
74 fAnalysisObjects = ref.fAnalysisObjects;
75 fPIDperformance = ref.fPIDperformance;
76 fSignalToBackgroundMC = ref.fSignalToBackgroundMC;
77
78 return *this;
79}
80
81//____________________________________________________________
82AliHFEpostAnalysis::~AliHFEpostAnalysis(){
83 //
84 // Do not delete objects where we are not owner
85 //
86 if(fResults) delete fResults;
87}
88
89//____________________________________________________________
90Int_t AliHFEpostAnalysis::SetResults(TList *input){
91 //
92 // Publish the results to the post analysis
93 //
94 Int_t nFound = 0;
70da6c5a 95 fEfficiencyContainer = dynamic_cast<AliCFContainer *>(input->FindObject("trackContainer"));
d2af20c5 96 if(!fEfficiencyContainer){
97 AliError("Efficiency Correction Framework Container not found in the list of outputs");
98 } else {
99 SETBIT(fAnalysisObjects, kCFC);
100 nFound++;
101 }
102 fPIDperformance = dynamic_cast<THnSparseF *>(input->FindObject("PIDperformance"));
103 if(!fPIDperformance){
104 AliError("Histogram fPIDperformance not found in the List of Outputs");
105 } else {
106 SETBIT(fAnalysisObjects, kPIDperf);
107 nFound++;
108 }
109 fSignalToBackgroundMC = dynamic_cast<THnSparseF *>(input->FindObject("SignalToBackgroundMC"));
110 if(!fSignalToBackgroundMC){
111 AliError("Histogram fSignalToBackgroundMC not found in the list of outputs");
112 } else {
113 SETBIT(fAnalysisObjects, kSigBackg);
114 nFound++;
115 }
116 AliInfo(Form("Found %d analysis objects", nFound));
117 return nFound;
118}
119
120//____________________________________________________________
121void AliHFEpostAnalysis::StoreOutput(const char *filename){
122 //
123 // Save the results produced in a rootfile
124 //
125 if(fResults){
126 TFile *outfile = new TFile(filename, "RECREATE");
127 outfile->cd();
128 fResults->Write("HFEresults", TObject::kSingleKey);
129 outfile->Close();
130 delete outfile;
131 }
132}
133
134//____________________________________________________________
135void AliHFEpostAnalysis::DrawMCSignal2Background(){
136 //
137 // Draw the MC signal/background plots
138 //
139 if(!fSignalToBackgroundMC) return;
140
141 // First Select everything within the first ITS Layer
70da6c5a 142 fSignalToBackgroundMC->GetAxis(5)->SetRange(2,2);
143 TH1 *hEff[3], *hSB[3];
144 // Select for different charge
145 hEff[0] = CreateHistoSignalToBackgroundMC(0, 0);
146 hEff[1] = CreateHistoSignalToBackgroundMC(0, 1);
147 hEff[2] = CreateHistoSignalToBackgroundMC(0, 2);
148
149 hSB[0] = CreateHistoSignalToBackgroundMC(1, 0);
150 hSB[1] = CreateHistoSignalToBackgroundMC(1, 1);
151 hSB[2] = CreateHistoSignalToBackgroundMC(1, 2);
152
153 // Undo projections
154 fSignalToBackgroundMC->GetAxis(5)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetNbins());
155
d2af20c5 156 // Prepare Canvas
70da6c5a 157 TCanvas *cMCSB = new TCanvas("cMCSB", "MC Sig/Backg studies", 800, 400);
158 cMCSB->Divide(2);
159 TLegend *leg;
160 TH1 **sample[2] = {&hEff[0], &hSB[0]};
161 const char *chargename[3] = {"All Charges", "Negative Charge", "Positive Charge"};
162 for(Int_t isample = 0; isample < 2; isample++){
163 leg = new TLegend(0.7, 0.1, 0.89, 0.3);
164 leg->SetBorderSize(1);
165 leg->SetFillColor(kWhite);
166 cMCSB->cd(isample + 1);
167 for(Int_t icharge = 0; icharge < 3; icharge++){
168 sample[isample][icharge]->Draw(icharge > 0? "psame" : "p");
169 leg->AddEntry(sample[isample][icharge], chargename[icharge], "p");
170 }
171 leg->Draw();
172 gPad->Update();
173 }
174}
175
176//____________________________________________________________
177TH1 *AliHFEpostAnalysis::CreateHistoSignalToBackgroundMC(Int_t mode, Int_t charge){
178 //
179 // Make Efficiency / SB histograms for different charges
180 //
181 TH1 *hNom = NULL, *hDenom = NULL;
182 if(charge) fSignalToBackgroundMC->GetAxis(3)->SetRange(charge, charge);
183 // For Signal electrons we project axis 4 to everything > 0
184 fSignalToBackgroundMC->GetAxis(4)->SetRange(2,3);
185 hNom = fSignalToBackgroundMC->Projection(0);
186 hNom->SetName("nom");
187 fSignalToBackgroundMC->GetAxis(4)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetLast() + 1);
188 if(mode == 1) fSignalToBackgroundMC->GetAxis(4)->SetRange(1,1);
189 hDenom = fSignalToBackgroundMC->Projection(0);
190 hDenom->SetName("denom");
191 if(mode == 1) fSignalToBackgroundMC->GetAxis(4)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetLast() + 1);
192 if(charge) fSignalToBackgroundMC->GetAxis(3)->SetRange(0, fSignalToBackgroundMC->GetAxis(3)->GetLast() + 1);
193
194 TH1 *hEff = dynamic_cast<TH1D *>(hNom->Clone());
195 char hname[256];
196 sprintf(hname, mode ? "sigToBack" : "sigEff");
197 char cname[256];
198 Color_t mycolor = kBlack;
199 switch(charge){
200 case 0: mycolor = kBlue; sprintf(cname, "All"); break;
201 case 1: mycolor = kBlack; sprintf(cname, "Neg"); break;
202 case 2: mycolor = kRed; sprintf(cname, "Pos"); break;
203 default: break;
204 }
205 sprintf(hname, "%s%s", hname, cname);
206 hEff->SetName(hname);
207 hEff->SetTitle(mode ? "Signal/Background" : "Signal/(Signal+Background)");
208 hEff->Divide(hDenom);
209
210 // Make nice plots
d2af20c5 211 hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
212 hEff->GetYaxis()->SetTitle("Efficiency");
213 hEff->SetStats(kFALSE);
70da6c5a 214 hEff->SetLineColor(kBlack);
215 hEff->SetLineWidth(1);
d2af20c5 216 hEff->SetMarkerStyle(22);
70da6c5a 217 hEff->SetMarkerColor(mycolor);
d2af20c5 218
70da6c5a 219 delete hNom; delete hDenom;
220 return hEff;
d2af20c5 221}
222
223//____________________________________________________________
224void AliHFEpostAnalysis::DrawEfficiency(){
225 //
226 // Draw the Efficiency
227 // We show:
228 // + InAcceptance / Generated
229 // + Signal / Generated
230 // + Selected / Generated
231 // + Selected / InAcceptance (Reconstructible)
232 //
233 TCanvas *cEff = new TCanvas("cEff", "Efficiency", 800, 600);
234 cEff->Divide(2,2);
235 if(!fEfficiencyContainer) return;
236 AliCFEffGrid *effCalc = new AliCFEffGrid("effCalc", "Efficiency Calculation Grid", *fEfficiencyContainer);
237 effCalc->CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
238 TH1 *effReconstructibleP = effCalc->Project(0);
239 effReconstructibleP->SetName("effReconstructibleP");
240 effReconstructibleP->SetTitle("Efficiency of reconstructible tracks");
241 effReconstructibleP->GetXaxis()->SetTitle("p_{T} / GeV/c");
242 effReconstructibleP->GetYaxis()->SetTitle("Efficiency");
243 effReconstructibleP->GetYaxis()->SetRangeUser(0.,1.);
244 effReconstructibleP->SetMarkerStyle(22);
245 effReconstructibleP->SetMarkerColor(kBlue);
246 effReconstructibleP->SetLineColor(kBlack);
247 effReconstructibleP->SetStats(kFALSE);
248 cEff->cd(1);
249 effReconstructibleP->Draw("e");
250 effCalc->CalculateEfficiency(AliHFEcuts::kStepMCsignal, AliHFEcuts::kStepMCGenerated);
251 TH1 *effSignal = effCalc->Project(0);
252 effSignal->SetName("effSignal");
253 effSignal->SetTitle("Efficiency of Signal Electrons");
254 effSignal->GetXaxis()->SetTitle("p_{T} / GeV/c");
255 effSignal->GetYaxis()->SetTitle("Efficiency");
256 effSignal->GetYaxis()->SetRangeUser(0., 1.);
257 effSignal->SetMarkerStyle(22);
258 effSignal->SetMarkerColor(kBlue);
259 effSignal->SetLineColor(kBlack);
260 effSignal->SetStats(kFALSE);
261 cEff->cd(2);
262 effSignal->Draw("e");
70da6c5a 263 effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
d2af20c5 264 TH1 *effPIDP = effCalc->Project(0);
265 effPIDP->SetName("effPIDP");
266 effPIDP->SetTitle("Efficiency of selected tracks");
267 effPIDP->GetXaxis()->SetTitle("p_{T} / GeV/c");
268 effPIDP->GetYaxis()->SetTitle("Efficiency");
269 effPIDP->GetYaxis()->SetRangeUser(0.,1.);
270 effPIDP->SetMarkerStyle(22);
271 effPIDP->SetMarkerColor(kBlue);
272 effPIDP->SetLineColor(kBlack);
273 effPIDP->SetStats(kFALSE);
274 cEff->cd(3);
275 effPIDP->Draw("e");
70da6c5a 276 effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCInAcceptance);
d2af20c5 277 TH1 *effPIDAcc = effCalc->Project(0);
278 effPIDAcc->SetName("effPIDAcc");
279 effPIDAcc->SetTitle("Efficiency of selected tracks in acceptance");
280 effPIDAcc->GetXaxis()->SetTitle("p_{T} / GeV/c");
281 effPIDAcc->GetYaxis()->SetTitle("Efficiency");
282 effPIDAcc->GetYaxis()->SetRangeUser(0.,1.);
283 effPIDAcc->SetMarkerStyle(22);
284 effPIDAcc->SetMarkerColor(kBlue);
285 effPIDAcc->SetLineColor(kBlack);
286 effPIDAcc->SetStats(kFALSE);
287 cEff->cd(4);
288 effPIDAcc->Draw("e");
289 delete effCalc;
290}
291
292//____________________________________________________________
293void AliHFEpostAnalysis::DrawPIDperformance(){
294 //
295 // Plotting Ratio histograms
296 // + All electrons / all candidates (Purity for Electrons)
297 // + All signal electrons / all electrons (Purity for signals)
d2af20c5 298 //
299
300 if(!fPIDperformance) return;
301 // Make projection
70da6c5a 302 TH1 *electronPurity[3], *signalPurity[3], *fakeContamination[3];
303 electronPurity[0] = CreateHistoPIDperformance(0, 0);
304 electronPurity[1] = CreateHistoPIDperformance(0, 1);
305 electronPurity[2] = CreateHistoPIDperformance(0, 2);
d2af20c5 306
70da6c5a 307 signalPurity[0] = CreateHistoPIDperformance(1, 0);
308 signalPurity[1] = CreateHistoPIDperformance(1, 1);
309 signalPurity[2] = CreateHistoPIDperformance(1, 2);
d2af20c5 310
70da6c5a 311 fakeContamination[0] = CreateHistoPIDperformance(2, 0);
312 fakeContamination[1] = CreateHistoPIDperformance(2, 1);
313 fakeContamination[2] = CreateHistoPIDperformance(2, 2);
d2af20c5 314
315 // Draw output
d2af20c5 316 TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600);
70da6c5a 317 const char *chargename[3] = {"All Charges", "Negative Charge", "Positive Charge"};
d2af20c5 318 cRatios->Divide(2,2);
70da6c5a 319 TH1 **sample[3] = {&electronPurity[0], &signalPurity[0], &fakeContamination[0]};
320 TLegend *leg;
321 for(Int_t isample = 0; isample < 3; isample++){
322 cRatios->cd(isample + 1);
323 leg = new TLegend(0.7, 0.1, 0.89, 0.3);
324 leg->SetBorderSize(1);
325 leg->SetFillColor(kWhite);
326 for(Int_t icharge = 0; icharge < 3; icharge++){
327 leg->AddEntry(sample[isample][icharge], chargename[icharge], "p");
328 sample[isample][icharge]->Draw(icharge > 0 ? "esame" : "e");
329 }
330 leg->Draw();
331 gPad->Update();
332 }
d2af20c5 333}
334
70da6c5a 335//____________________________________________________________
336TH1 *AliHFEpostAnalysis::CreateHistoPIDperformance(Int_t mode, Int_t charge){
337 //
338 // Make Histograms for PID performance plots
339 //
faee3b18 340 fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
341 fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
342
70da6c5a 343 TH1 *hNom = NULL, *hDenom = NULL;
344 char hname[256], htitle[256], cname[256];
345 Color_t mycolor = kBlack;
faee3b18 346 if(charge) fPIDperformance->GetAxis(3)->SetRange(charge, charge);
70da6c5a 347 // Normalisation by all candidates - no restriction in axis 4 - only for mode == 1
faee3b18 348 if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(2,3);
70da6c5a 349 hDenom = fPIDperformance->Projection(0);
350 hDenom->Sumw2();
351 hDenom->SetName("hDenom");
faee3b18 352 if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins() + 1);
70da6c5a 353 // Nominator need a restriction in the 4th axis
354 switch(mode){
355 case 0: // Electron purity
356 fPIDperformance->GetAxis(4)->SetRange(2,3);
357 sprintf(hname, "electronPurity");
358 sprintf(htitle, "Electron Purity");
359 break;
360 case 1: // Signal purity
361 fPIDperformance->GetAxis(4)->SetRange(3,3); // here signal not divided into charm and beauty
362 sprintf(hname, "signalPurity");
363 sprintf(htitle, "Signal Purity");
364 break;
365 case 2: // Fake contamination
366 fPIDperformance->GetAxis(4)->SetRange(1,1);
367 sprintf(hname, "fakeContamination");
368 sprintf(htitle, "Contamination of misidentified hadrons");
369 break;
370 default: break;
371 }
372 switch(charge){
373 case 0:
374 sprintf(cname, "All");
375 mycolor = kBlue;
376 break;
377 case 1:
378 sprintf(cname, "Neg");
379 mycolor = kBlack;
380 break;
381 case 2:
382 sprintf(cname, "Pos");
383 mycolor = kRed;
384 break;
385 }
386 sprintf(hname, "%s%s", hname, cname);
387 hNom = fPIDperformance->Projection(0);
388 hNom->Sumw2();
389 hNom->SetName("hNom");
390 // Reset axis
faee3b18 391 fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
392 if(charge) fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
70da6c5a 393
394 // Create Efficiency histogram
395 TH1 *hEff = dynamic_cast<TH1D *>(hNom->Clone());
396 hEff->SetName(hname);
397 hEff->SetTitle(htitle);
398 hEff->Divide(hDenom);
399 hEff->Scale(100.);
400 hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
401 hEff->GetYaxis()->SetTitle(mode < 2 ? "Purity / %" : "Contamination / %");
402 hEff->GetYaxis()->SetRangeUser(0., 100.);
403 hEff->SetStats(kFALSE);
404 hEff->SetLineColor(kBlack);
405 hEff->SetLineWidth(1);
406 hEff->SetMarkerColor(mycolor);
407 hEff->SetMarkerStyle(22);
408 delete hNom; delete hDenom;
409 return hEff;
410}
411
412//____________________________________________________________
413void AliHFEpostAnalysis::DrawCutEfficiency(){
414 //
415 // Calculate efficiency for each cut step
416 // Starting from MC steps
417 //
418 TCanvas *output = new TCanvas("effCut", "Cut Step efficiency", 800, 600);
419 output->cd();
420 TLegend *leg = new TLegend(0.6, 0.7, 0.89, 0.89);
421 leg->SetHeader("Cut Step:");
422 leg->SetBorderSize(0);
423 leg->SetFillColor(kWhite);
424 leg->SetFillStyle(0);
425 const char* names[AliHFEcuts::kNcutStepsTrack - 1] = {"Acceptance", "No Cuts", "Rec Tracks" "Rec Kine(TPC/ITS)", "Primary", "HFE ITS", "HFE TRD", "PID", };
426 Color_t color[AliHFEcuts::kNcutStepsTrack - 1] = {kRed, kGreen, kMagenta, kBlack, kOrange, kTeal, kViolet, kBlue};
427 Marker_t marker[AliHFEcuts::kNcutStepsTrack - 1] = {24, 25, 26, 27, 28, 29, 30};
428 TH1 *hTemp = NULL;
429 AliCFEffGrid effcalc("cutEfficiency", "Cut step efficiency calculation", *fEfficiencyContainer);
430 // Calculate efficiency for MC Steps
431 effcalc.CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
432 hTemp = effcalc.Project(0);
433 hTemp->SetName("hEff1");
434 hTemp->SetMarkerColor(color[0]);
435 hTemp->SetMarkerStyle(marker[0]);
436 hTemp->GetXaxis()->SetTitle("P / GeV/c");
437 hTemp->GetYaxis()->SetTitle("Efficiency");
438 hTemp->GetYaxis()->SetRangeUser(0., 2.);
439 hTemp->SetStats(kFALSE);
440 hTemp->Draw("ep");
441 leg->AddEntry(hTemp, names[0], "p");
442 effcalc.CalculateEfficiency(AliHFEcuts::kStepRecNoCut + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
443 hTemp = effcalc.Project(0);
444 hTemp->SetName("hEff2");
445 hTemp->SetMarkerColor(color[1]);
446 hTemp->SetMarkerStyle(marker[1]);
447 hTemp->GetXaxis()->SetTitle("P / GeV/c");
448 hTemp->GetYaxis()->SetTitle("Efficiency");
449 hTemp->GetYaxis()->SetRangeUser(0., 2.);
450 hTemp->SetStats(kFALSE);
451 hTemp->Draw("epsame");
452 leg->AddEntry(hTemp, names[1], "p");
453 for(Int_t istep = AliHFEcuts::kStepRecKineITSTPC; istep <= AliHFEcuts::kStepPID; istep++){
454 effcalc.CalculateEfficiency(istep + 2*AliHFEcuts::kNcutStepsESDtrack, istep-1 + 2*AliHFEcuts::kNcutStepsESDtrack);
455 hTemp = effcalc.Project(0);
456 hTemp->SetName(Form("hEff%d", istep));
457 hTemp->SetMarkerColor(color[istep-2]);
458 hTemp->SetMarkerStyle(marker[istep-2]);
459 hTemp->SetStats(kFALSE);
460 hTemp->Draw("epsame");
461 hTemp->GetXaxis()->SetTitle("P / GeV/c");
462 hTemp->GetYaxis()->SetTitle("Efficiency");
463 hTemp->GetYaxis()->SetRangeUser(0., 2.);
464 leg->AddEntry(hTemp, names[istep-2], "p");
465 }
466 leg->Draw();
467}