[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpostAnalysis.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/
//
// Post analysis code
// Drawing nice pictures containing
//  - Efficiency
//  - Signal/Background
//  - PID Performance
//  More Post Analysis code will be added in time
//
//  Autor: Markus Fasel
//
#include <TAxis.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TH1D.h>
#include <TLegend.h>
#include <TList.h>
#include "AliCFContainer.h"
#include "AliCFEffGrid.h"

#include "AliHFEcuts.h"
#include "AliHFEpostAnalysis.h"

ClassImp(AliHFEpostAnalysis)

//____________________________________________________________
AliHFEpostAnalysis::AliHFEpostAnalysis():
  TObject(),
  fResults(NULL),
  fAnalysisObjects(0),
  fEfficiencyContainer(NULL),
  fPIDperformance(NULL),
  fSignalToBackgroundMC(NULL)
{
  //
  // Default Constructor
  //
}

//____________________________________________________________
AliHFEpostAnalysis::AliHFEpostAnalysis(const AliHFEpostAnalysis &ref):
  TObject(ref),
  fResults(ref.fResults),
  fAnalysisObjects(ref.fAnalysisObjects),
  fEfficiencyContainer(ref.fEfficiencyContainer),
  fPIDperformance(ref.fPIDperformance),
  fSignalToBackgroundMC(ref.fSignalToBackgroundMC)
{
  //
  // Copy Constructor
  //
}

//____________________________________________________________
AliHFEpostAnalysis& AliHFEpostAnalysis::operator=(const AliHFEpostAnalysis &ref){
  //
  // Assignment Operator
  //
  TObject::operator=(ref);
  fResults = ref.fResults;
  fAnalysisObjects = ref.fAnalysisObjects;
  fPIDperformance = ref.fPIDperformance;
  fSignalToBackgroundMC = ref.fSignalToBackgroundMC;

  return *this;
}

//____________________________________________________________
AliHFEpostAnalysis::~AliHFEpostAnalysis(){
  //
  // Do not delete objects where we are not owner
  //
  if(fResults) delete fResults;
}

//____________________________________________________________
Int_t AliHFEpostAnalysis::SetResults(TList *input){
  //
  // Publish the results to the post analysis
  //
  Int_t nFound = 0;
  fEfficiencyContainer = dynamic_cast<AliCFContainer *>(input->FindObject("trackContainer"));
  if(!fEfficiencyContainer){
    AliError("Efficiency Correction Framework Container not found in the list of outputs");
  } else {
    SETBIT(fAnalysisObjects, kCFC);
    nFound++;
  }
  fPIDperformance = dynamic_cast<THnSparseF *>(input->FindObject("PIDperformance"));
  if(!fPIDperformance){
    AliError("Histogram fPIDperformance not found in the List of Outputs");
  } else {
    SETBIT(fAnalysisObjects, kPIDperf);
    nFound++;
  }
  fSignalToBackgroundMC = dynamic_cast<THnSparseF *>(input->FindObject("SignalToBackgroundMC"));
  if(!fSignalToBackgroundMC){
    AliError("Histogram fSignalToBackgroundMC not found in the list of outputs");
  } else {
    SETBIT(fAnalysisObjects, kSigBackg);
    nFound++;
  }
  AliInfo(Form("Found %d analysis objects", nFound));
  return nFound;
}

//____________________________________________________________
void AliHFEpostAnalysis::StoreOutput(const char *filename){
  //
  // Save the results produced in a rootfile
  //
  if(fResults){
    TFile *outfile = new TFile(filename, "RECREATE");
    outfile->cd();
    fResults->Write("HFEresults", TObject::kSingleKey);
    outfile->Close();
    delete outfile;
  }
}

//____________________________________________________________
void AliHFEpostAnalysis::DrawMCSignal2Background(){
  //
  // Draw the MC signal/background plots
  //
  if(!fSignalToBackgroundMC) return;

  // First Select everything within the first ITS Layer
  fSignalToBackgroundMC->GetAxis(5)->SetRange(2,2);
  TH1 *hEff[3], *hSB[3];
  // Select for different charge
  hEff[0] = CreateHistoSignalToBackgroundMC(0, 0);
  hEff[1] = CreateHistoSignalToBackgroundMC(0, 1);
  hEff[2] = CreateHistoSignalToBackgroundMC(0, 2);

  hSB[0] = CreateHistoSignalToBackgroundMC(1, 0);
  hSB[1] = CreateHistoSignalToBackgroundMC(1, 1);
  hSB[2] = CreateHistoSignalToBackgroundMC(1, 2);
 
  // Undo projections
  fSignalToBackgroundMC->GetAxis(5)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetNbins());

  // Prepare Canvas
  TCanvas *cMCSB = new TCanvas("cMCSB", "MC Sig/Backg studies", 800, 400);
  cMCSB->Divide(2);
  TLegend *leg;
  TH1 **sample[2] = {&hEff[0], &hSB[0]};
  const char *chargename[3] = {"All Charges", "Negative Charge", "Positive Charge"};
  for(Int_t isample = 0; isample < 2; isample++){
    leg = new TLegend(0.7, 0.1, 0.89, 0.3);
    leg->SetBorderSize(1); 
    leg->SetFillColor(kWhite);
    cMCSB->cd(isample + 1);
    for(Int_t icharge = 0; icharge < 3; icharge++){
      sample[isample][icharge]->Draw(icharge > 0?  "psame" : "p");
      leg->AddEntry(sample[isample][icharge], chargename[icharge], "p");
    }
    leg->Draw();
    gPad->Update();
  }
}

//____________________________________________________________
TH1 *AliHFEpostAnalysis::CreateHistoSignalToBackgroundMC(Int_t mode, Int_t charge){ 
  //
  // Make Efficiency / SB histograms for different charges
  //
  TH1 *hNom = NULL, *hDenom = NULL;
  if(charge) fSignalToBackgroundMC->GetAxis(3)->SetRange(charge, charge);
  // For Signal electrons we project axis 4 to everything > 0
  fSignalToBackgroundMC->GetAxis(4)->SetRange(2,3);
  hNom = fSignalToBackgroundMC->Projection(0);
  hNom->SetName("nom");
  fSignalToBackgroundMC->GetAxis(4)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetLast() + 1);
  if(mode == 1) fSignalToBackgroundMC->GetAxis(4)->SetRange(1,1);
  hDenom =  fSignalToBackgroundMC->Projection(0);
  hDenom->SetName("denom");
  if(mode == 1)  fSignalToBackgroundMC->GetAxis(4)->SetRange(0, fSignalToBackgroundMC->GetAxis(4)->GetLast() + 1);
  if(charge) fSignalToBackgroundMC->GetAxis(3)->SetRange(0, fSignalToBackgroundMC->GetAxis(3)->GetLast() + 1);

  TH1 *hEff = dynamic_cast<TH1D *>(hNom->Clone());
  char hname[256];
  sprintf(hname, mode ? "sigToBack" : "sigEff");
  char cname[256];
  Color_t mycolor = kBlack;
  switch(charge){
    case 0: mycolor = kBlue; sprintf(cname, "All"); break;
    case 1: mycolor = kBlack; sprintf(cname, "Neg"); break;
    case 2: mycolor = kRed; sprintf(cname, "Pos"); break;
    default: break;
  }
  sprintf(hname, "%s%s", hname, cname);
  hEff->SetName(hname);
  hEff->SetTitle(mode ? "Signal/Background" : "Signal/(Signal+Background)");
  hEff->Divide(hDenom);

  // Make nice plots
  hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
  hEff->GetYaxis()->SetTitle("Efficiency");
  hEff->SetStats(kFALSE);
  hEff->SetLineColor(kBlack);
  hEff->SetLineWidth(1);
  hEff->SetMarkerStyle(22);
  hEff->SetMarkerColor(mycolor);

  delete hNom; delete hDenom;
  return hEff;
}

//____________________________________________________________
void AliHFEpostAnalysis::DrawEfficiency(){
  //
  // Draw the Efficiency
  // We show: 
  // + InAcceptance / Generated
  // + Signal / Generated
  // + Selected / Generated
  // + Selected / InAcceptance (Reconstructible)
  //
  TCanvas *cEff = new TCanvas("cEff", "Efficiency", 800, 600);
  cEff->Divide(2,2);
  if(!fEfficiencyContainer) return;
  AliCFEffGrid *effCalc = new AliCFEffGrid("effCalc", "Efficiency Calculation Grid", *fEfficiencyContainer);
  effCalc->CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
  TH1 *effReconstructibleP = effCalc->Project(0);
  effReconstructibleP->SetName("effReconstructibleP");
  effReconstructibleP->SetTitle("Efficiency of reconstructible tracks");
  effReconstructibleP->GetXaxis()->SetTitle("p_{T} / GeV/c");
  effReconstructibleP->GetYaxis()->SetTitle("Efficiency");
  effReconstructibleP->GetYaxis()->SetRangeUser(0.,1.);
  effReconstructibleP->SetMarkerStyle(22);
  effReconstructibleP->SetMarkerColor(kBlue);
  effReconstructibleP->SetLineColor(kBlack);
  effReconstructibleP->SetStats(kFALSE);
  cEff->cd(1);
  effReconstructibleP->Draw("e");
  effCalc->CalculateEfficiency(AliHFEcuts::kStepMCsignal, AliHFEcuts::kStepMCGenerated);
  TH1 *effSignal = effCalc->Project(0);
  effSignal->SetName("effSignal");
  effSignal->SetTitle("Efficiency of Signal Electrons");
  effSignal->GetXaxis()->SetTitle("p_{T} / GeV/c");
  effSignal->GetYaxis()->SetTitle("Efficiency");
  effSignal->GetYaxis()->SetRangeUser(0., 1.);
  effSignal->SetMarkerStyle(22);
  effSignal->SetMarkerColor(kBlue);
  effSignal->SetLineColor(kBlack);
  effSignal->SetStats(kFALSE);
  cEff->cd(2);
  effSignal->Draw("e");
  effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
  TH1 *effPIDP = effCalc->Project(0);
  effPIDP->SetName("effPIDP");
  effPIDP->SetTitle("Efficiency of selected tracks");
  effPIDP->GetXaxis()->SetTitle("p_{T} / GeV/c");
  effPIDP->GetYaxis()->SetTitle("Efficiency");
  effPIDP->GetYaxis()->SetRangeUser(0.,1.);
  effPIDP->SetMarkerStyle(22);
  effPIDP->SetMarkerColor(kBlue);
  effPIDP->SetLineColor(kBlack);
  effPIDP->SetStats(kFALSE);
  cEff->cd(3);
  effPIDP->Draw("e");
  effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCInAcceptance);
  TH1 *effPIDAcc = effCalc->Project(0);
  effPIDAcc->SetName("effPIDAcc");
  effPIDAcc->SetTitle("Efficiency of selected tracks in acceptance");
  effPIDAcc->GetXaxis()->SetTitle("p_{T} / GeV/c");
  effPIDAcc->GetYaxis()->SetTitle("Efficiency");
  effPIDAcc->GetYaxis()->SetRangeUser(0.,1.);
  effPIDAcc->SetMarkerStyle(22);
  effPIDAcc->SetMarkerColor(kBlue);
  effPIDAcc->SetLineColor(kBlack);
  effPIDAcc->SetStats(kFALSE);
  cEff->cd(4);
  effPIDAcc->Draw("e");
  delete effCalc;
}

//____________________________________________________________
void AliHFEpostAnalysis::DrawPIDperformance(){
  //
  // Plotting Ratio histograms
  // + All electrons / all candidates (Purity for Electrons)
  // + All signal electrons / all electrons (Purity for signals)
  //

  if(!fPIDperformance) return;
  // Make projection
  TH1 *electronPurity[3], *signalPurity[3], *fakeContamination[3];
  electronPurity[0] = CreateHistoPIDperformance(0, 0);
  electronPurity[1] = CreateHistoPIDperformance(0, 1);
  electronPurity[2] = CreateHistoPIDperformance(0, 2);

  signalPurity[0] = CreateHistoPIDperformance(1, 0);
  signalPurity[1] = CreateHistoPIDperformance(1, 1);
  signalPurity[2] = CreateHistoPIDperformance(1, 2);

  fakeContamination[0] = CreateHistoPIDperformance(2, 0);
  fakeContamination[1] = CreateHistoPIDperformance(2, 1);
  fakeContamination[2] = CreateHistoPIDperformance(2, 2);

  // Draw output
  TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600);
  const char *chargename[3] = {"All Charges", "Negative Charge", "Positive Charge"};
  cRatios->Divide(2,2);
  TH1 **sample[3] = {&electronPurity[0], &signalPurity[0], &fakeContamination[0]};
  TLegend *leg;
  for(Int_t isample = 0; isample < 3; isample++){
    cRatios->cd(isample + 1);
    leg = new TLegend(0.7, 0.1, 0.89, 0.3);
    leg->SetBorderSize(1);
    leg->SetFillColor(kWhite);
    for(Int_t icharge = 0; icharge < 3; icharge++){
      leg->AddEntry(sample[isample][icharge], chargename[icharge], "p");
      sample[isample][icharge]->Draw(icharge > 0 ? "esame" : "e");
    }
    leg->Draw();
    gPad->Update();
  }
}

//____________________________________________________________
TH1 *AliHFEpostAnalysis::CreateHistoPIDperformance(Int_t mode, Int_t charge){
  //
  // Make Histograms for PID performance plots
  //
  fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
  fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);

  TH1 *hNom = NULL, *hDenom = NULL;
  char hname[256], htitle[256], cname[256];
  Color_t mycolor = kBlack;
  if(charge) fPIDperformance->GetAxis(3)->SetRange(charge, charge);
  // Normalisation by all candidates - no restriction in axis 4 - only for mode == 1 
  if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(2,3);
  hDenom = fPIDperformance->Projection(0);
  hDenom->Sumw2();
  hDenom->SetName("hDenom");
  if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins() + 1);
  // Nominator need a restriction in the 4th axis
  switch(mode){
    case 0: // Electron purity
      fPIDperformance->GetAxis(4)->SetRange(2,3);
      sprintf(hname, "electronPurity");
      sprintf(htitle, "Electron Purity");
      break;
    case 1: // Signal purity
      fPIDperformance->GetAxis(4)->SetRange(3,3);   // here signal not divided into charm and beauty
      sprintf(hname, "signalPurity");
      sprintf(htitle, "Signal Purity");
      break;
    case 2: // Fake contamination
      fPIDperformance->GetAxis(4)->SetRange(1,1);
      sprintf(hname, "fakeContamination");
      sprintf(htitle, "Contamination of misidentified hadrons");
      break;
    default: break;
  }
  switch(charge){
    case 0: 
      sprintf(cname, "All"); 
      mycolor = kBlue;
      break;
    case 1: 
      sprintf(cname, "Neg"); 
      mycolor = kBlack;
      break;
    case 2: 
      sprintf(cname, "Pos"); 
      mycolor = kRed;
      break;
  }
  sprintf(hname, "%s%s", hname, cname);
  hNom = fPIDperformance->Projection(0);
  hNom->Sumw2();
  hNom->SetName("hNom");
  // Reset axis
  fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
  if(charge) fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);

  // Create Efficiency histogram
  TH1 *hEff = dynamic_cast<TH1D *>(hNom->Clone());
  hEff->SetName(hname);
  hEff->SetTitle(htitle);
  hEff->Divide(hDenom);
  hEff->Scale(100.);
  hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
  hEff->GetYaxis()->SetTitle(mode < 2 ? "Purity / %" : "Contamination / %");
  hEff->GetYaxis()->SetRangeUser(0., 100.);
  hEff->SetStats(kFALSE);
  hEff->SetLineColor(kBlack);
  hEff->SetLineWidth(1);
  hEff->SetMarkerColor(mycolor);
  hEff->SetMarkerStyle(22);
  delete hNom; delete hDenom;
  return hEff;
}

//____________________________________________________________
void AliHFEpostAnalysis::DrawCutEfficiency(){
  //
  // Calculate efficiency for each cut step
  // Starting from MC steps 
  //
  TCanvas *output = new TCanvas("effCut", "Cut Step efficiency", 800, 600);
  output->cd();
  TLegend *leg = new TLegend(0.6, 0.7, 0.89, 0.89);
  leg->SetHeader("Cut Step:");
  leg->SetBorderSize(0);
  leg->SetFillColor(kWhite);
  leg->SetFillStyle(0);
  const char* names[AliHFEcuts::kNcutStepsTrack - 1] = {"Acceptance", "No Cuts", "Rec Tracks" "Rec Kine(TPC/ITS)", "Primary", "HFE ITS", "HFE TRD", "PID", };
  Color_t color[AliHFEcuts::kNcutStepsTrack - 1] = {kRed, kGreen, kMagenta, kBlack, kOrange, kTeal, kViolet, kBlue};
  Marker_t marker[AliHFEcuts::kNcutStepsTrack - 1] = {24, 25, 26, 27, 28, 29, 30}; 
  TH1 *hTemp = NULL;
  AliCFEffGrid effcalc("cutEfficiency", "Cut step efficiency calculation", *fEfficiencyContainer);
  // Calculate efficiency for MC Steps
  effcalc.CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
  hTemp = effcalc.Project(0);
  hTemp->SetName("hEff1");
  hTemp->SetMarkerColor(color[0]);
  hTemp->SetMarkerStyle(marker[0]);
  hTemp->GetXaxis()->SetTitle("P / GeV/c");
  hTemp->GetYaxis()->SetTitle("Efficiency");
  hTemp->GetYaxis()->SetRangeUser(0., 2.);
  hTemp->SetStats(kFALSE);
  hTemp->Draw("ep");
  leg->AddEntry(hTemp, names[0], "p");
  effcalc.CalculateEfficiency(AliHFEcuts::kStepRecNoCut + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
  hTemp = effcalc.Project(0);
  hTemp->SetName("hEff2");
  hTemp->SetMarkerColor(color[1]);
  hTemp->SetMarkerStyle(marker[1]);
  hTemp->GetXaxis()->SetTitle("P / GeV/c");
  hTemp->GetYaxis()->SetTitle("Efficiency");
  hTemp->GetYaxis()->SetRangeUser(0., 2.);
  hTemp->SetStats(kFALSE);
  hTemp->Draw("epsame");
  leg->AddEntry(hTemp, names[1], "p");
  for(Int_t istep = AliHFEcuts::kStepRecKineITSTPC; istep <= AliHFEcuts::kStepPID; istep++){
    effcalc.CalculateEfficiency(istep + 2*AliHFEcuts::kNcutStepsESDtrack, istep-1 + 2*AliHFEcuts::kNcutStepsESDtrack); 
    hTemp = effcalc.Project(0);
    hTemp->SetName(Form("hEff%d", istep));
    hTemp->SetMarkerColor(color[istep-2]);
    hTemp->SetMarkerStyle(marker[istep-2]);
    hTemp->SetStats(kFALSE);
    hTemp->Draw("epsame");
    hTemp->GetXaxis()->SetTitle("P / GeV/c");
    hTemp->GetYaxis()->SetTitle("Efficiency");
    hTemp->GetYaxis()->SetRangeUser(0., 2.);
    leg->AddEntry(hTemp, names[istep-2], "p");
  }
  leg->Draw();
}
Commit	Line	Data
d2af20c5	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;
70da6c5a	95	fEfficiencyContainer = dynamic_cast<AliCFContainer *>(input->FindObject("trackContainer"));
d2af20c5	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
70da6c5a	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
d2af20c5	156	// Prepare Canvas
70da6c5a	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
d2af20c5	211	hEff->GetXaxis()->SetTitle("p_{T} / GeV/c");
	212	hEff->GetYaxis()->SetTitle("Efficiency");
	213	hEff->SetStats(kFALSE);
70da6c5a	214	hEff->SetLineColor(kBlack);
70da6c5a	215	hEff->SetLineWidth(1);
d2af20c5	216	hEff->SetMarkerStyle(22);
70da6c5a	217	hEff->SetMarkerColor(mycolor);
d2af20c5	218
70da6c5a	219	delete hNom; delete hDenom;
70da6c5a	220	return hEff;
d2af20c5	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");
70da6c5a	263	effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
d2af20c5	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");
70da6c5a	276	effCalc->CalculateEfficiency(AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCInAcceptance);
d2af20c5	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)
d2af20c5	298	//
	299
	300	if(!fPIDperformance) return;
	301	// Make projection
70da6c5a	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);
d2af20c5	306
70da6c5a	307	signalPurity[0] = CreateHistoPIDperformance(1, 0);
	308	signalPurity[1] = CreateHistoPIDperformance(1, 1);
	309	signalPurity[2] = CreateHistoPIDperformance(1, 2);
d2af20c5	310
70da6c5a	311	fakeContamination[0] = CreateHistoPIDperformance(2, 0);
	312	fakeContamination[1] = CreateHistoPIDperformance(2, 1);
	313	fakeContamination[2] = CreateHistoPIDperformance(2, 2);
d2af20c5	314
d2af20c5	315	// Draw output
d2af20c5	316	TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600);
70da6c5a	317	const char *chargename[3] = {"All Charges", "Negative Charge", "Positive Charge"};
d2af20c5	318	cRatios->Divide(2,2);
70da6c5a	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	}
d2af20c5	333	}
d2af20c5	334
70da6c5a	335	//____________________________________________________________
	336	TH1 *AliHFEpostAnalysis::CreateHistoPIDperformance(Int_t mode, Int_t charge){
	337	//
	338	// Make Histograms for PID performance plots
	339	//
faee3b18	340	fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
	341	fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
	342
70da6c5a	343	TH1 hNom = NULL, hDenom = NULL;
	344	char hname[256], htitle[256], cname[256];
	345	Color_t mycolor = kBlack;
faee3b18	346	if(charge) fPIDperformance->GetAxis(3)->SetRange(charge, charge);
70da6c5a	347	// Normalisation by all candidates - no restriction in axis 4 - only for mode == 1
faee3b18	348	if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(2,3);
70da6c5a	349	hDenom = fPIDperformance->Projection(0);
	350	hDenom->Sumw2();
	351	hDenom->SetName("hDenom");
faee3b18	352	if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins() + 1);
70da6c5a	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
faee3b18	391	fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
faee3b18	392	if(charge) fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
70da6c5a	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(){
	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	AliCFEffGrid effcalc("cutEfficiency", "Cut step efficiency calculation", *fEfficiencyContainer);
	430	// Calculate efficiency for MC Steps
	431	effcalc.CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
	432	hTemp = effcalc.Project(0);
	433	hTemp->SetName("hEff1");
	434	hTemp->SetMarkerColor(color[0]);
	435	hTemp->SetMarkerStyle(marker[0]);
	436	hTemp->GetXaxis()->SetTitle("P / GeV/c");
	437	hTemp->GetYaxis()->SetTitle("Efficiency");
	438	hTemp->GetYaxis()->SetRangeUser(0., 2.);
	439	hTemp->SetStats(kFALSE);
	440	hTemp->Draw("ep");
	441	leg->AddEntry(hTemp, names[0], "p");
	442	effcalc.CalculateEfficiency(AliHFEcuts::kStepRecNoCut + 2*AliHFEcuts::kNcutStepsESDtrack, AliHFEcuts::kStepMCGenerated);
	443	hTemp = effcalc.Project(0);
	444	hTemp->SetName("hEff2");
	445	hTemp->SetMarkerColor(color[1]);
	446	hTemp->SetMarkerStyle(marker[1]);
	447	hTemp->GetXaxis()->SetTitle("P / GeV/c");
	448	hTemp->GetYaxis()->SetTitle("Efficiency");
	449	hTemp->GetYaxis()->SetRangeUser(0., 2.);
	450	hTemp->SetStats(kFALSE);
	451	hTemp->Draw("epsame");
	452	leg->AddEntry(hTemp, names[1], "p");
	453	for(Int_t istep = AliHFEcuts::kStepRecKineITSTPC; istep <= AliHFEcuts::kStepPID; istep++){
	454	effcalc.CalculateEfficiency(istep + 2AliHFEcuts::kNcutStepsESDtrack, istep-1 + 2AliHFEcuts::kNcutStepsESDtrack);
	455	hTemp = effcalc.Project(0);
	456	hTemp->SetName(Form("hEff%d", istep));
457	hTemp->SetMarkerColor(color[istep-2]);
458	hTemp->SetMarkerStyle(marker[istep-2]);
459	hTemp->SetStats(kFALSE);
460	hTemp->Draw("epsame");
461	hTemp->GetXaxis()->SetTitle("P / GeV/c");
462	hTemp->GetYaxis()->SetTitle("Efficiency");
463	hTemp->GetYaxis()->SetRangeUser(0., 2.);
464	leg->AddEntry(hTemp, names[istep-2], "p");
465	}
466	leg->Draw();
467	}