[u/mrichter/AliRoot.git] / PWGHF / 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 "AliHFEcontainer.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
  //
  if(this == &ref) return *this;
  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::SetTaskQA(const TList *input){
  //
  // Publish the results to the post analysis
  //
  Int_t nFound = 0;
  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());
  if(hEff){
    TString hname, cname;
    hname = mode ? "sigToBack" : "sigEff";
    Color_t mycolor = kBlack;
    switch(charge){
      case 0: mycolor = kBlue; cname = "All"; break;
      case 1: mycolor = kBlack; cname = "Neg"; break;
      case 2: mycolor = kRed; cname ="Pos"; break;
      default: break;
    }
    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;
  AliCFContainer *tracks = fEfficiencyContainer->MakeMergedCFContainer("trackContCombined", "MC + Rec(reco) Track Information", "MCTrackCont:recTrackContReco");
  AliCFEffGrid *effCalc = new AliCFEffGrid("effCalc", "Efficiency Calculation Grid", *tracks);
  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::kStepMCGeneratedZOutNoPileUpCentralityFine, 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(tracks->GetNStep() - 1, 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(tracks->GetNStep() - 1, 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;
  TString hname, htitle, cname;
  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);
      hname = "electronPurity";
      htitle = "Electron Purity";
      break;
    case 1: // Signal purity
      fPIDperformance->GetAxis(4)->SetRange(3,3);   // here signal not divided into charm and beauty
      hname = "signalPurity";
      htitle = "Signal Purity";
      break;
    case 2: // Fake contamination
      fPIDperformance->GetAxis(4)->SetRange(1,1);
      hname = "fakeContamination";
      htitle = "Contamination of misidentified hadrons";
      break;
    default: break;
  }
  switch(charge){
    case 0: 
      cname = "All"; 
      mycolor = kBlue;
      break;
    case 1: 
      cname = "Neg"; 
      mycolor = kBlack;
      break;
    case 2: 
      cname = "Pos"; 
      mycolor = kRed;
      break;
  }
  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());
  if(hEff){
    hEff->SetName(hname.Data());
    hEff->SetTitle(htitle.Data());
    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(Bool_t MC, Int_t source){
  //
  // 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);

  AliCFContainer *tracks = fEfficiencyContainer->MakeMergedCFContainer("mergedTracks", "Container for MC and reconstructed Track information", "MCTrackCont:recTrackContReco");
  Int_t nStepMC = fEfficiencyContainer->GetCFContainer("MCTrackCont")->GetNStep();
  AliDebug(1, Form("Number of MC Cut Steps: %d", nStepMC));
  const Int_t markerStart = 24;
  const Int_t colorStart = 1;
  TH1 *hTemp = NULL;

  if(MC){
    if(source > -1 && source < 4){
      AliInfo(Form("Setting source to %d", source));
      for(Int_t istep = 0; istep < tracks->GetNStep(); istep++) tracks->GetAxis(4, istep)->SetRange(source + 1, source + 1);
    }
  }
  AliCFEffGrid effcalc("cutEfficiency", "Cut step efficiency calculation", *tracks);
  Bool_t first = kTRUE;
  // Calculate efficiency for MC Steps
  Int_t histcounter = 0;
  if(MC){
    //
    // Draw plots related to MC values
    //
    effcalc.CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
    hTemp = effcalc.Project(0);
    hTemp->SetName(Form("hEff%d", histcounter));
    hTemp->SetTitle("Cut Step Efficiency");
    hTemp->SetMarkerColor(colorStart + histcounter);
    hTemp->SetMarkerStyle(markerStart + histcounter);
    hTemp->GetXaxis()->SetTitle("p / GeV/c");
    hTemp->GetYaxis()->SetTitle("Efficiency");
    hTemp->GetYaxis()->SetRangeUser(0., 2.);
    hTemp->SetStats(kFALSE);
    hTemp->Draw("ep");
    leg->AddEntry(hTemp, tracks->GetStepTitle(AliHFEcuts::kStepMCInAcceptance), "p");
    histcounter++;
    effcalc.CalculateEfficiency(nStepMC, AliHFEcuts::kStepMCGenerated);
    hTemp = effcalc.Project(0);
    hTemp->SetName("hEff2");
    hTemp->SetTitle("Cut Step Efficiency");
    hTemp->SetMarkerColor(colorStart + histcounter);
    hTemp->SetMarkerStyle(markerStart + histcounter);
    hTemp->GetXaxis()->SetTitle("p / GeV/c");
    hTemp->GetYaxis()->SetTitle("Efficiency");
    hTemp->GetYaxis()->SetRangeUser(0., 2.);
    hTemp->SetStats(kFALSE);
    hTemp->Draw("epsame");
    leg->AddEntry(hTemp, tracks->GetStepTitle(nStepMC), "p");
    first = kFALSE;
    histcounter++;
  }
  for(Int_t istep = nStepMC+1; istep < tracks->GetNStep(); istep++){
    effcalc.CalculateEfficiency(istep, istep - 1); 
    hTemp = effcalc.Project(0);
    hTemp->SetName(Form("hEff%d", istep));
    hTemp->SetTitle("Cut Step Efficiency");
    hTemp->SetMarkerColor(colorStart + histcounter);
    hTemp->SetMarkerStyle(markerStart + histcounter);
    hTemp->SetStats(kFALSE);
    hTemp->Draw(first ? "ep" : "epsame");
    hTemp->GetXaxis()->SetTitle("P / GeV/c");
    hTemp->GetYaxis()->SetTitle("Efficiency");
    hTemp->GetYaxis()->SetRangeUser(0., 2.);
    leg->AddEntry(hTemp, tracks->GetStepTitle(istep), "p");
    first = kFALSE;
    histcounter++;
  }
  leg->Draw();
  delete tracks;
}
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	//
3a72645a	25
d2af20c5	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
3a72645a	35	#include "AliHFEcontainer.h"
d2af20c5	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	//
aa1417ee	74	if(this == &ref) return *this;
d2af20c5	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	//____________________________________________________________
3a72645a	93	Int_t AliHFEpostAnalysis::SetTaskQA(const TList *input){
d2af20c5	94	//
	95	// Publish the results to the post analysis
	96	//
	97	Int_t nFound = 0;
d2af20c5	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
70da6c5a	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
d2af20c5	152	// Prepare Canvas
70da6c5a	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());
bf892a6a	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);
70da6c5a	205
bf892a6a	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	}
d2af20c5	215
70da6c5a	216	delete hNom; delete hDenom;
70da6c5a	217	return hEff;
d2af20c5	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;
3a72645a	233	AliCFContainer *tracks = fEfficiencyContainer->MakeMergedCFContainer("trackContCombined", "MC + Rec(reco) Track Information", "MCTrackCont:recTrackContReco");
3a72645a	234	AliCFEffGrid effCalc = new AliCFEffGrid("effCalc", "Efficiency Calculation Grid", tracks);
d2af20c5	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");
c2690925	248	effCalc->CalculateEfficiency(AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, AliHFEcuts::kStepMCGenerated);
d2af20c5	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");
3a72645a	261	effCalc->CalculateEfficiency(tracks->GetNStep() - 1, AliHFEcuts::kStepMCGenerated);
d2af20c5	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");
3a72645a	274	effCalc->CalculateEfficiency(tracks->GetNStep() - 1, AliHFEcuts::kStepMCInAcceptance);
d2af20c5	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)
d2af20c5	296	//
	297
	298	if(!fPIDperformance) return;
	299	// Make projection
70da6c5a	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);
d2af20c5	304
70da6c5a	305	signalPurity[0] = CreateHistoPIDperformance(1, 0);
	306	signalPurity[1] = CreateHistoPIDperformance(1, 1);
	307	signalPurity[2] = CreateHistoPIDperformance(1, 2);
d2af20c5	308
70da6c5a	309	fakeContamination[0] = CreateHistoPIDperformance(2, 0);
	310	fakeContamination[1] = CreateHistoPIDperformance(2, 1);
	311	fakeContamination[2] = CreateHistoPIDperformance(2, 2);
d2af20c5	312
d2af20c5	313	// Draw output
d2af20c5	314	TCanvas *cRatios = new TCanvas("cRatios", "Ratio Plots", 800, 600);
70da6c5a	315	const char *chargename[3] = {"All Charges", "Negative Charge", "Positive Charge"};
d2af20c5	316	cRatios->Divide(2,2);
70da6c5a	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	}
d2af20c5	331	}
d2af20c5	332
70da6c5a	333	//____________________________________________________________
	334	TH1 *AliHFEpostAnalysis::CreateHistoPIDperformance(Int_t mode, Int_t charge){
	335	//
	336	// Make Histograms for PID performance plots
	337	//
faee3b18	338	fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
	339	fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
	340
70da6c5a	341	TH1 hNom = NULL, hDenom = NULL;
bf892a6a	342	TString hname, htitle, cname;
70da6c5a	343	Color_t mycolor = kBlack;
faee3b18	344	if(charge) fPIDperformance->GetAxis(3)->SetRange(charge, charge);
70da6c5a	345	// Normalisation by all candidates - no restriction in axis 4 - only for mode == 1
faee3b18	346	if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(2,3);
70da6c5a	347	hDenom = fPIDperformance->Projection(0);
	348	hDenom->Sumw2();
	349	hDenom->SetName("hDenom");
faee3b18	350	if(mode == 1) fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins() + 1);
70da6c5a	351	// Nominator need a restriction in the 4th axis
	352	switch(mode){
	353	case 0: // Electron purity
	354	fPIDperformance->GetAxis(4)->SetRange(2,3);
bf892a6a	355	hname = "electronPurity";
bf892a6a	356	htitle = "Electron Purity";
70da6c5a	357	break;
	358	case 1: // Signal purity
	359	fPIDperformance->GetAxis(4)->SetRange(3,3); // here signal not divided into charm and beauty
bf892a6a	360	hname = "signalPurity";
bf892a6a	361	htitle = "Signal Purity";
70da6c5a	362	break;
	363	case 2: // Fake contamination
	364	fPIDperformance->GetAxis(4)->SetRange(1,1);
bf892a6a	365	hname = "fakeContamination";
bf892a6a	366	htitle = "Contamination of misidentified hadrons";
70da6c5a	367	break;
	368	default: break;
	369	}
	370	switch(charge){
	371	case 0:
bf892a6a	372	cname = "All";
70da6c5a	373	mycolor = kBlue;
	374	break;
	375	case 1:
bf892a6a	376	cname = "Neg";
70da6c5a	377	mycolor = kBlack;
	378	break;
	379	case 2:
bf892a6a	380	cname = "Pos";
70da6c5a	381	mycolor = kRed;
	382	break;
	383	}
bf892a6a	384	hname += cname;
70da6c5a	385	hNom = fPIDperformance->Projection(0);
	386	hNom->Sumw2();
	387	hNom->SetName("hNom");
	388	// Reset axis
faee3b18	389	fPIDperformance->GetAxis(4)->SetRange(0, fPIDperformance->GetAxis(4)->GetNbins()+1);
faee3b18	390	if(charge) fPIDperformance->GetAxis(3)->SetRange(0, fPIDperformance->GetAxis(3)->GetNbins() + 1);
70da6c5a	391
	392	// Create Efficiency histogram
	393	TH1 hEff = dynamic_cast<TH1D >(hNom->Clone());
bf892a6a	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	}
70da6c5a	409	return hEff;
	410	}
	411
	412	//____________________________________________________________
67fe7bd0	413	void AliHFEpostAnalysis::DrawCutEfficiency(Bool_t MC, Int_t source){
70da6c5a	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);
3a72645a	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;
70da6c5a	431	TH1 *hTemp = NULL;
67fe7bd0	432
	433	if(MC){
	434	if(source > -1 && source < 4){
13242232	435	AliInfo(Form("Setting source to %d", source));
3a72645a	436	for(Int_t istep = 0; istep < tracks->GetNStep(); istep++) tracks->GetAxis(4, istep)->SetRange(source + 1, source + 1);
67fe7bd0	437	}
67fe7bd0	438	}
3a72645a	439	AliCFEffGrid effcalc("cutEfficiency", "Cut step efficiency calculation", *tracks);
67fe7bd0	440	Bool_t first = kTRUE;
70da6c5a	441	// Calculate efficiency for MC Steps
3a72645a	442	Int_t histcounter = 0;
67fe7bd0	443	if(MC){
	444	//
	445	// Draw plots related to MC values
	446	//
	447	effcalc.CalculateEfficiency(AliHFEcuts::kStepMCInAcceptance, AliHFEcuts::kStepMCGenerated);
	448	hTemp = effcalc.Project(0);
3a72645a	449	hTemp->SetName(Form("hEff%d", histcounter));
67fe7bd0	450	hTemp->SetTitle("Cut Step Efficiency");
3a72645a	451	hTemp->SetMarkerColor(colorStart + histcounter);
	452	hTemp->SetMarkerStyle(markerStart + histcounter);
	453	hTemp->GetXaxis()->SetTitle("p / GeV/c");
67fe7bd0	454	hTemp->GetYaxis()->SetTitle("Efficiency");
	455	hTemp->GetYaxis()->SetRangeUser(0., 2.);
	456	hTemp->SetStats(kFALSE);
	457	hTemp->Draw("ep");
3a72645a	458	leg->AddEntry(hTemp, tracks->GetStepTitle(AliHFEcuts::kStepMCInAcceptance), "p");
	459	histcounter++;
	460	effcalc.CalculateEfficiency(nStepMC, AliHFEcuts::kStepMCGenerated);
67fe7bd0	461	hTemp = effcalc.Project(0);
	462	hTemp->SetName("hEff2");
	463	hTemp->SetTitle("Cut Step Efficiency");
3a72645a	464	hTemp->SetMarkerColor(colorStart + histcounter);
	465	hTemp->SetMarkerStyle(markerStart + histcounter);
	466	hTemp->GetXaxis()->SetTitle("p / GeV/c");
67fe7bd0	467	hTemp->GetYaxis()->SetTitle("Efficiency");
	468	hTemp->GetYaxis()->SetRangeUser(0., 2.);
	469	hTemp->SetStats(kFALSE);
	470	hTemp->Draw("epsame");
3a72645a	471	leg->AddEntry(hTemp, tracks->GetStepTitle(nStepMC), "p");
67fe7bd0	472	first = kFALSE;
3a72645a	473	histcounter++;
67fe7bd0	474	}
3a72645a	475	for(Int_t istep = nStepMC+1; istep < tracks->GetNStep(); istep++){
3a72645a	476	effcalc.CalculateEfficiency(istep, istep - 1);
70da6c5a	477	hTemp = effcalc.Project(0);
70da6c5a	478	hTemp->SetName(Form("hEff%d", istep));
67fe7bd0	479	hTemp->SetTitle("Cut Step Efficiency");
3a72645a	480	hTemp->SetMarkerColor(colorStart + histcounter);
3a72645a	481	hTemp->SetMarkerStyle(markerStart + histcounter);
70da6c5a	482	hTemp->SetStats(kFALSE);
67fe7bd0	483	hTemp->Draw(first ? "ep" : "epsame");
70da6c5a	484	hTemp->GetXaxis()->SetTitle("P / GeV/c");
	485	hTemp->GetYaxis()->SetTitle("Efficiency");
	486	hTemp->GetYaxis()->SetRangeUser(0., 2.);
3a72645a	487	leg->AddEntry(hTemp, tracks->GetStepTitle(istep), "p");
67fe7bd0	488	first = kFALSE;
3a72645a	489	histcounter++;
70da6c5a	490	}
70da6c5a	491	leg->Draw();
3a72645a	492	delete tracks;
70da6c5a	493	}