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