]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG3/hfe/AliHFEpostAnalysis.cxx
Update of the HFE package and of the macro to run on the
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpostAnalysis.cxx
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
37 ClassImp(AliHFEpostAnalysis)
38
39 //____________________________________________________________
40 AliHFEpostAnalysis::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 //____________________________________________________________
54 AliHFEpostAnalysis::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 //____________________________________________________________
68 AliHFEpostAnalysis& 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 //____________________________________________________________
82 AliHFEpostAnalysis::~AliHFEpostAnalysis(){
83   //
84   // Do not delete objects where we are not owner
85   //
86   if(fResults) delete fResults;
87 }
88
89 //____________________________________________________________
90 Int_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("trackContainer"));
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 //____________________________________________________________
121 void 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 //____________________________________________________________
135 void 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(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
156   // Prepare Canvas
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 //____________________________________________________________
177 TH1 *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
211   hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
212   hEff->GetYaxis()->SetTitle("Efficiency");
213   hEff->SetStats(kFALSE);
214   hEff->SetLineColor(kBlack);
215   hEff->SetLineWidth(1);
216   hEff->SetMarkerStyle(22);
217   hEff->SetMarkerColor(mycolor);
218
219   delete hNom; delete hDenom;
220   return hEff;
221 }
222
223 //____________________________________________________________
224 void 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");
263   effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
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");
276   effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCInAcceptance);
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 //____________________________________________________________
293 void AliHFEpostAnalysis::DrawPIDperformance(){
294   //
295   // Plotting Ratio histograms
296   // + All electrons / all candidates (Purity for Electrons)
297   // + All signal electrons / all electrons (Purity for signals)
298   //
299
300   if(!fPIDperformance) return;
301   // Make projection
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);
306
307   signalPurity[0] = CreateHistoPIDperformance(1, 0);
308   signalPurity[1] = CreateHistoPIDperformance(1, 1);
309   signalPurity[2] = CreateHistoPIDperformance(1, 2);
310
311   fakeContamination[0] = CreateHistoPIDperformance(2, 0);
312   fakeContamination[1] = CreateHistoPIDperformance(2, 1);
313   fakeContamination[2] = CreateHistoPIDperformance(2, 2);
314
315   // Draw output
316   TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600);
317   const char *chargename[3] = {"All Charges", "Negative Charge", "Positive Charge"};
318   cRatios->Divide(2,2);
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   }
333 }
334
335 //____________________________________________________________
336 TH1 *AliHFEpostAnalysis::CreateHistoPIDperformance(Int_t mode, Int_t charge){
337   //
338   // Make Histograms for PID performance plots
339   //
340   fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
341   fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
342
343   TH1 *hNom = NULL, *hDenom = NULL;
344   char hname[256], htitle[256], cname[256];
345   Color_t mycolor = kBlack;
346   if(charge) fPIDperformance->GetAxis(3)->SetRange(charge, charge);
347   // Normalisation by all candidates - no restriction in axis 4 - only for mode == 1 
348   if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(2,3);
349   hDenom = fPIDperformance->Projection(0);
350   hDenom->Sumw2();
351   hDenom->SetName("hDenom");
352   if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins() + 1);
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
391   fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
392   if(charge) fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
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 //____________________________________________________________
413 void AliHFEpostAnalysis::DrawCutEfficiency(Bool_t MC, Int_t source){
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
430   if(MC){
431     if(source > -1 && source < 4){
432       AliInfo(Form("Setting source to %d", source))
433       for(Int_t istep = 0; istep < fEfficiencyContainer->GetNStep(); istep++) fEfficiencyContainer->GetAxis(4, istep)->SetRange(source + 1, source + 1);
434     }
435   }
436   AliCFEffGrid effcalc("cutEfficiency", "Cut step efficiency calculation", *fEfficiencyContainer);
437   Bool_t first = kTRUE;
438   // Calculate efficiency for MC Steps
439   if(MC){
440     //
441     // Draw plots related to MC values
442     //
443     effcalc.CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
444     hTemp = effcalc.Project(0);
445     hTemp->SetName("hEff1");
446     hTemp->SetTitle("Cut Step Efficiency");
447     hTemp->SetMarkerColor(color[0]);
448     hTemp->SetMarkerStyle(marker[0]);
449     hTemp->GetXaxis()->SetTitle("P / GeV/c");
450     hTemp->GetYaxis()->SetTitle("Efficiency");
451     hTemp->GetYaxis()->SetRangeUser(0., 2.);
452     hTemp->SetStats(kFALSE);
453     hTemp->Draw("ep");
454     leg->AddEntry(hTemp, names[0], "p");
455     effcalc.CalculateEfficiency(AliHFEcuts::kStepRecNoCut + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
456     hTemp = effcalc.Project(0);
457     hTemp->SetName("hEff2");
458     hTemp->SetTitle("Cut Step Efficiency");
459     hTemp->SetMarkerColor(color[1]);
460     hTemp->SetMarkerStyle(marker[1]);
461     hTemp->GetXaxis()->SetTitle("P / GeV/c");
462     hTemp->GetYaxis()->SetTitle("Efficiency");
463     hTemp->GetYaxis()->SetRangeUser(0., 2.);
464     hTemp->SetStats(kFALSE);
465     hTemp->Draw("epsame");
466     leg->AddEntry(hTemp, names[1], "p");
467     first = kFALSE;
468   }
469   for(Int_t istep = AliHFEcuts::kStepRecKineITSTPC; istep <= AliHFEcuts::kStepPID; istep++){
470     effcalc.CalculateEfficiency(istep + 2*AliHFEcuts::kNcutStepsESDtrack, istep-1 + 2*AliHFEcuts::kNcutStepsESDtrack); 
471     hTemp = effcalc.Project(0);
472     hTemp->SetName(Form("hEff%d", istep));
473     hTemp->SetTitle("Cut Step Efficiency");
474     hTemp->SetMarkerColor(color[istep-2]);
475     hTemp->SetMarkerStyle(marker[istep-2]);
476     hTemp->SetStats(kFALSE);
477     hTemp->Draw(first ? "ep" : "epsame");
478     hTemp->GetXaxis()->SetTitle("P / GeV/c");
479     hTemp->GetYaxis()->SetTitle("Efficiency");
480     hTemp->GetYaxis()->SetRangeUser(0., 2.);
481     leg->AddEntry(hTemp, names[istep-2], "p");
482     first = kFALSE;
483   }
484   leg->Draw();
485   // undo axis restriction
486   if(MC && source > -1 && source < 4){
487     for(Int_t istep = 0; istep < fEfficiencyContainer->GetNStep(); istep++) fEfficiencyContainer->GetAxis(4, istep)->SetRange(0, 4);
488   }
489 }