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 }