extra check on pointer
[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;
95 fEfficiencyContainer = dynamic_cast<AliCFContainer *>(input->FindObject("container"));
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
142 fSignalToBackgroundMC->GetAxis(4)->SetRange(2,2);
143 // For Signal electrons we project axis 3 to everything > 0
144 fSignalToBackgroundMC->GetAxis(3)->SetRange(2,3);
145 TH1 *hSignal = fSignalToBackgroundMC->Projection(0);
146 hSignal->SetName("hSignal");
147 hSignal->SetTitle("Signal Electrons");
148 // For Background studies project axis 3 to bin 0
149 fSignalToBackgroundMC->GetAxis(3)->SetRange(1,1);
150 TH1 *hBackground = fSignalToBackgroundMC->Projection(0);
151 hBackground->SetName("hBackground");
152 hBackground->SetTitle("Background Electrons");
153 // For All electrons we undo the projection of axis 3
154 fSignalToBackgroundMC->GetAxis(3)->SetRange(0, fSignalToBackgroundMC->GetAxis(3)->GetNbins());
155 TH1 *hAll = fSignalToBackgroundMC->Projection(0);
156 hAll->SetName("hAll");
157 hAll->SetTitle("All Electrons");
158 // Prepare Canvas
159 TCanvas *cEff = new TCanvas("cEff", "MC Sig/Backg studies", 800, 400);
160 cEff->Divide(2);
161 // Plot Efficiency
162 TH1 *hEff = (TH1 *)hSignal->Clone();
163 hEff->Divide(hAll);
164 hEff->SetName("sigEff");
165 hEff->SetTitle("Signal/(Signal+Background)");
166 hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
167 hEff->GetYaxis()->SetTitle("Efficiency");
168 hEff->SetStats(kFALSE);
169 hEff->SetMarkerStyle(22);
170 hEff->SetMarkerColor(kBlue);
171 cEff->cd(1);
172 hEff->Draw("p");
173 // Plot Signal/Background
174 TH1 *hSB = (TH1 *)hSignal->Clone();
175 hSB->Divide(hBackground);
176 hSB->SetName("sigEff");
177 hSB->SetTitle("Signal/Background");
178 hSB->GetXaxis()->SetTitle("p_{T} / GeV/c");
179 hSB->GetYaxis()->SetTitle("Signal/Background");
180 hSB->SetStats(kFALSE);
181 hSB->SetMarkerStyle(22);
182 hSB->SetMarkerColor(kBlue);
183 cEff->cd(2);
184 hSB->Draw("p");
185
186 // Undo projections
187 fSignalToBackgroundMC->GetAxis(4)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetNbins());
188}
189
190//____________________________________________________________
191void AliHFEpostAnalysis::DrawEfficiency(){
192 //
193 // Draw the Efficiency
194 // We show:
195 // + InAcceptance / Generated
196 // + Signal / Generated
197 // + Selected / Generated
198 // + Selected / InAcceptance (Reconstructible)
199 //
200 TCanvas *cEff = new TCanvas("cEff", "Efficiency", 800, 600);
201 cEff->Divide(2,2);
202 if(!fEfficiencyContainer) return;
203 AliCFEffGrid *effCalc = new AliCFEffGrid("effCalc", "Efficiency Calculation Grid", *fEfficiencyContainer);
204 effCalc->CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
205 TH1 *effReconstructibleP = effCalc->Project(0);
206 effReconstructibleP->SetName("effReconstructibleP");
207 effReconstructibleP->SetTitle("Efficiency of reconstructible tracks");
208 effReconstructibleP->GetXaxis()->SetTitle("p_{T} / GeV/c");
209 effReconstructibleP->GetYaxis()->SetTitle("Efficiency");
210 effReconstructibleP->GetYaxis()->SetRangeUser(0.,1.);
211 effReconstructibleP->SetMarkerStyle(22);
212 effReconstructibleP->SetMarkerColor(kBlue);
213 effReconstructibleP->SetLineColor(kBlack);
214 effReconstructibleP->SetStats(kFALSE);
215 cEff->cd(1);
216 effReconstructibleP->Draw("e");
217 effCalc->CalculateEfficiency(AliHFEcuts::kStepMCsignal, AliHFEcuts::kStepMCGenerated);
218 TH1 *effSignal = effCalc->Project(0);
219 effSignal->SetName("effSignal");
220 effSignal->SetTitle("Efficiency of Signal Electrons");
221 effSignal->GetXaxis()->SetTitle("p_{T} / GeV/c");
222 effSignal->GetYaxis()->SetTitle("Efficiency");
223 effSignal->GetYaxis()->SetRangeUser(0., 1.);
224 effSignal->SetMarkerStyle(22);
225 effSignal->SetMarkerColor(kBlue);
226 effSignal->SetLineColor(kBlack);
227 effSignal->SetStats(kFALSE);
228 cEff->cd(2);
229 effSignal->Draw("e");
230 effCalc->CalculateEfficiency(AliHFEcuts::kStepHFEcutsTRD + 1, AliHFEcuts::kStepMCGenerated);
231 TH1 *effPIDP = effCalc->Project(0);
232 effPIDP->SetName("effPIDP");
233 effPIDP->SetTitle("Efficiency of selected tracks");
234 effPIDP->GetXaxis()->SetTitle("p_{T} / GeV/c");
235 effPIDP->GetYaxis()->SetTitle("Efficiency");
236 effPIDP->GetYaxis()->SetRangeUser(0.,1.);
237 effPIDP->SetMarkerStyle(22);
238 effPIDP->SetMarkerColor(kBlue);
239 effPIDP->SetLineColor(kBlack);
240 effPIDP->SetStats(kFALSE);
241 cEff->cd(3);
242 effPIDP->Draw("e");
243 effCalc->CalculateEfficiency(AliHFEcuts::kStepHFEcutsTRD + 1, AliHFEcuts::kStepMCInAcceptance);
244 TH1 *effPIDAcc = effCalc->Project(0);
245 effPIDAcc->SetName("effPIDAcc");
246 effPIDAcc->SetTitle("Efficiency of selected tracks in acceptance");
247 effPIDAcc->GetXaxis()->SetTitle("p_{T} / GeV/c");
248 effPIDAcc->GetYaxis()->SetTitle("Efficiency");
249 effPIDAcc->GetYaxis()->SetRangeUser(0.,1.);
250 effPIDAcc->SetMarkerStyle(22);
251 effPIDAcc->SetMarkerColor(kBlue);
252 effPIDAcc->SetLineColor(kBlack);
253 effPIDAcc->SetStats(kFALSE);
254 cEff->cd(4);
255 effPIDAcc->Draw("e");
256 delete effCalc;
257}
258
259//____________________________________________________________
260void AliHFEpostAnalysis::DrawPIDperformance(){
261 //
262 // Plotting Ratio histograms
263 // + All electrons / all candidates (Purity for Electrons)
264 // + All signal electrons / all electrons (Purity for signals)
265 // For this the following pt-histograms have to be projected from the THnSparse
266 // + All Electron candidates
267 // + All Real electrons
268 // + All Signal Electrons
269 // + All misidentified electrons
270 // Additionally we draw Efficiency histograms:
271 // + Reconstructible Electrons
272 // + Signal Electrons
273 // + Selected Electrons
274 // + Selected Electrons in Acceptance
275 //
276
277 if(!fPIDperformance) return;
278 // Make projection
279 // always project to pt dimension
280 // get the histograms under our control
281 TH1 *allCandidates = NULL, *allElectrons = NULL, *allSignals = NULL, *allFakes = NULL;
282 allCandidates = fPIDperformance->Projection(0);
283 allCandidates->SetName("hAllCandidates");
284 allCandidates->SetTitle("All Candidates");
285 allCandidates->Sumw2();
286 Int_t firstDim3 = fPIDperformance->GetAxis(3)->GetFirst(), lastDim3 = fPIDperformance->GetAxis(3)->GetLast();
287 fPIDperformance->GetAxis(3)->SetRange(firstDim3 + 1, lastDim3);
288 allElectrons = fPIDperformance->Projection(0);
289 allElectrons->Sumw2();
290 allElectrons->SetName("hAllElectrons");
291 allElectrons->SetTitle("All Electrons");
292 Int_t firstDim4 = fPIDperformance->GetAxis(4)->GetFirst(), lastDim4 = fPIDperformance->GetAxis(4)->GetLast();
293 fPIDperformance->GetAxis(4)->SetRange(firstDim4 + 1, lastDim4);
294 allSignals = fPIDperformance->Projection(0);
295 allSignals->Sumw2();
296 allSignals->SetName("hAllSignals");
297 allSignals->SetTitle("All Signal Electrons");
298 fPIDperformance->GetAxis(4)->SetRange(firstDim4, lastDim4); // Reset 4th axis
299 fPIDperformance->GetAxis(3)->SetRange(firstDim3, firstDim3); // Select fakes
300 allFakes = fPIDperformance->Projection(0);
301 allFakes->Sumw2();
302 allFakes->SetName("hAllFakes");
303 allFakes->SetTitle("All Fakes");
304 fPIDperformance->GetAxis(3)->SetRange(firstDim3, lastDim3); // Reset also 3rd axis
305
306 // Make Ratios
307 TH1D *electronPurity = dynamic_cast<TH1D *>(allElectrons->Clone());
308 electronPurity->Divide(allCandidates);
309 electronPurity->SetName("electronPurity");
310 electronPurity->SetTitle("Electron Purity");
311 electronPurity->GetXaxis()->SetTitle("p_{T} / GeV/c");
312 electronPurity->GetYaxis()->SetTitle("Purity / %");
313 TH1D *signalPurity = dynamic_cast<TH1D *>(allSignals->Clone());
314 signalPurity->Divide(allElectrons);
315 signalPurity->SetName("signalPurity");
316 signalPurity->SetTitle("Purity of Electrons coming from Heavy flavours");
317 signalPurity->GetXaxis()->SetTitle("p_{T} / GeV/c");
318 signalPurity->GetYaxis()->SetTitle("Purity / %");
319 TH1D *fakeContamination = dynamic_cast<TH1D *>(allFakes->Clone());
320 fakeContamination->Divide(allCandidates);
321 fakeContamination->SetName("fakeContamination");
322 fakeContamination->SetTitle("Contamination of misidentified hadrons");
323 fakeContamination->GetXaxis()->SetTitle("p_{T} / GeV/c");
324 fakeContamination->GetYaxis()->SetTitle("Purity / %");
325
326 // Graphics settings
327 const Int_t nHistos = 7;
328 TH1 *hptr[7] = {allCandidates, allElectrons, allSignals, allFakes, electronPurity, signalPurity, fakeContamination}, *histo;
329 for(Int_t ihist = 0; ihist < nHistos; ihist++){
330 (histo = hptr[ihist])->SetStats(kFALSE);
331 histo->SetLineColor(kBlack);
332 histo->SetMarkerColor(kBlue);
333 histo->SetMarkerStyle(22);
334 }
335
336 // Make percent output
337 electronPurity->Scale(100);
338 signalPurity->Scale(100);
339 fakeContamination->Scale(100);
340
341 // Draw output
342 TCanvas *cSpectra = new TCanvas("cSpectra","pt Spectra", 800, 600);
343 cSpectra->Divide(2,2);
344 TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600);
345 cRatios->Divide(2,2);
346 cSpectra->cd(1);
347 gPad->SetLogy();
348 allCandidates->GetXaxis()->SetTitle("p_{T} / GeV/c");
349 allCandidates->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
350 allCandidates->Draw("e");
351 cSpectra->cd(2);
352 gPad->SetLogy();
353 allElectrons->GetXaxis()->SetTitle("p_{T} / GeV/c");
354 allElectrons->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
355 allElectrons->Draw("e");
356 cSpectra->cd(3);
357 gPad->SetLogy();
358 allSignals->GetXaxis()->SetTitle("p_{T} / GeV/c");
359 allSignals->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
360 allSignals->Draw("e");
361 cSpectra->cd(4);
362 gPad->SetLogy();
363 allFakes->GetXaxis()->SetTitle("p_{T} / GeV/c");
364 allFakes->GetYaxis()->SetTitle("dN/dp_{T} 1/GeV/c");
365 allFakes->Draw("e");
366 cRatios->cd(1);
367 electronPurity->Draw("e");
368 cRatios->cd(2);
369 signalPurity->Draw("e");
370 cRatios->cd(3);
371 fakeContamination->Draw("e");
372}
373