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   TH1 *hNom = NULL, *hDenom = NULL;
 341   char hname[256], htitle[256], cname[256];
 342   Color_t mycolor = kBlack;
 343   if(charge) fSignalToBackgroundMC->GetAxis(3)->SetRange(charge, charge);
 344   // Normalisation by all candidates - no restriction in axis 4 - only for mode == 1
 345   if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(2,2);
 346   hDenom = fPIDperformance->Projection(0);
 347   hDenom->Sumw2();
 348   hDenom->SetName("hDenom");
 349   if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetLast() + 1);
 350   // Nominator need a restriction in the 4th axis
 351   switch(mode){
 352     case 0: // Electron purity
 353       fPIDperformance->GetAxis(4)->SetRange(2,3);
 354       sprintf(hname, "electronPurity");
 355       sprintf(htitle, "Electron Purity");
 356       break;
 357     case 1: // Signal purity
 358       fPIDperformance->GetAxis(4)->SetRange(3,3);   // here signal not divided into charm and beauty
 359       sprintf(hname, "signalPurity");
 360       sprintf(htitle, "Signal Purity");
 361       break;
 362     case 2: // Fake contamination
 363       fPIDperformance->GetAxis(4)->SetRange(1,1);
 364       sprintf(hname, "fakeContamination");
 365       sprintf(htitle, "Contamination of misidentified hadrons");
 366       break;
 367     default: break;
 368   }
 369   switch(charge){
 370     case 0:
 371       sprintf(cname, "All");
 372       mycolor = kBlue;
 373       break;
 374     case 1:
 375       sprintf(cname, "Neg");
 376       mycolor = kBlack;
 377       break;
 378     case 2:
 379       sprintf(cname, "Pos");
 380       mycolor = kRed;
 381       break;
 382   }
 383   sprintf(hname, "%s%s", hname, cname);
 384   hNom = fPIDperformance->Projection(0);
 385   hNom->Sumw2();
 386   hNom->SetName("hNom");
 387   // Reset axis
 388   fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetLast() + 1);
 389   if(charge) fSignalToBackgroundMC->GetAxis(3)->SetRange(0, fSignalToBackgroundMC->GetAxis(3)->GetLast() + 1);
 390
 391   // Create Efficiency histogram
 392   TH1 *hEff = dynamic_cast<TH1D *>(hNom->Clone());
 393   hEff->SetName(hname);
 394   hEff->SetTitle(htitle);
 395   hEff->Divide(hDenom);
 396   hEff->Scale(100.);
 397   hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
 398   hEff->GetYaxis()->SetTitle(mode < 2 ? "Purity / %" : "Contamination / %");
 399   hEff->GetYaxis()->SetRangeUser(0., 100.);
 400   hEff->SetStats(kFALSE);
 401   hEff->SetLineColor(kBlack);
 402   hEff->SetLineWidth(1);
 403   hEff->SetMarkerColor(mycolor);
 404   hEff->SetMarkerStyle(22);
 405   delete hNom; delete hDenom;
 406   return hEff;
 407 }
 408
 409 //____________________________________________________________
 410 void AliHFEpostAnalysis::DrawCutEfficiency(){
 411   //
 412   // Calculate efficiency for each cut step
 413   // Starting from MC steps
 414   //
 415   TCanvas *output = new TCanvas("effCut", "Cut Step efficiency", 800, 600);
 416   output->cd();
 417   TLegend *leg = new TLegend(0.6, 0.7, 0.89, 0.89);
 418   leg->SetHeader("Cut Step:");
 419   leg->SetBorderSize(0);
 420   leg->SetFillColor(kWhite);
 421   leg->SetFillStyle(0);
 422   const char* names[AliHFEcuts::kNcutStepsTrack - 1] = {"Acceptance", "No Cuts", "Rec Tracks" "Rec Kine(TPC/ITS)", "Primary", "HFE ITS", "HFE TRD", "PID", };
 423   Color_t color[AliHFEcuts::kNcutStepsTrack - 1] = {kRed, kGreen, kMagenta, kBlack, kOrange, kTeal, kViolet, kBlue};
 424   Marker_t marker[AliHFEcuts::kNcutStepsTrack - 1] = {24, 25, 26, 27, 28, 29, 30};
 425   TH1 *hTemp = NULL;
 426   AliCFEffGrid effcalc("cutEfficiency", "Cut step efficiency calculation", *fEfficiencyContainer);
 427   // Calculate efficiency for MC Steps
 428   effcalc.CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
 429   hTemp = effcalc.Project(0);
 430   hTemp->SetName("hEff1");
 431   hTemp->SetMarkerColor(color[0]);
 432   hTemp->SetMarkerStyle(marker[0]);
 433   hTemp->GetXaxis()->SetTitle("P / GeV/c");
 434   hTemp->GetYaxis()->SetTitle("Efficiency");
 435   hTemp->GetYaxis()->SetRangeUser(0., 2.);
 436   hTemp->SetStats(kFALSE);
 437   hTemp->Draw("ep");
 438   leg->AddEntry(hTemp, names[0], "p");
 439   effcalc.CalculateEfficiency(AliHFEcuts::kStepRecNoCut + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
 440   hTemp = effcalc.Project(0);
 441   hTemp->SetName("hEff2");
 442   hTemp->SetMarkerColor(color[1]);
 443   hTemp->SetMarkerStyle(marker[1]);
 444   hTemp->GetXaxis()->SetTitle("P / GeV/c");
 445   hTemp->GetYaxis()->SetTitle("Efficiency");
 446   hTemp->GetYaxis()->SetRangeUser(0., 2.);
 447   hTemp->SetStats(kFALSE);
 448   hTemp->Draw("epsame");
 449   leg->AddEntry(hTemp, names[1], "p");
 450   for(Int_t istep = AliHFEcuts::kStepRecKineITSTPC; istep <= AliHFEcuts::kStepPID; istep++){
 451     effcalc.CalculateEfficiency(istep + 2*AliHFEcuts::kNcutStepsESDtrack, istep-1 + 2*AliHFEcuts::kNcutStepsESDtrack);
 452     hTemp = effcalc.Project(0);
 453     hTemp->SetName(Form("hEff%d", istep));
 454     hTemp->SetMarkerColor(color[istep-2]);
 455     hTemp->SetMarkerStyle(marker[istep-2]);
 456     hTemp->SetStats(kFALSE);
 457     hTemp->Draw("epsame");
 458     hTemp->GetXaxis()->SetTitle("P / GeV/c");
 459     hTemp->GetYaxis()->SetTitle("Efficiency");
 460     hTemp->GetYaxis()->SetRangeUser(0., 2.);
 461     leg->AddEntry(hTemp, names[istep-2], "p");
 462   }
 463   leg->Draw();
 464 }