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("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 //____________________________________________________________
 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(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 //____________________________________________________________
 191 void 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 //____________________________________________________________
 260 void 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