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