[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEefficiency.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.                  *
**************************************************************************/
//
// Task for Efficiency studies
// Used for testing classes AliHFEcontainer and AliHFEfilter
// Creates Efficiency Histograms
//
// Author:
//   Markus Fasel <M.Fasel@gsi.de>
//
#include <TAxis.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TGraphErrors.h>
#include <TH1F.h>
#include <TH2.h>
#include <TIterator.h>
#include <TLegend.h>
#include <TObjArray.h>
#include <TPad.h>

#include "AliAnalysisManager.h"
#include "AliCFAcceptanceCuts.h"
#include "AliCFTrackIsPrimaryCuts.h"
#include "AliCFContainer.h"
#include "AliCFEffGrid.h"
#include "AliESDEvent.h"
#include "AliHFEcollection.h"
#include "AliHFEcontainer.h"
#include "AliHFEcutStep.h"
#include "AliHFEefficiency.h"
#include "AliHFEextraCuts.h"
#include "AliHFEtrackFilter.h"
#include "AliHFEtools.h"
#include "AliMCEvent.h"
#include "AliMCEventHandler.h"

ClassImp(AliHFEefficiency)

AliHFEefficiency::AliHFEefficiency():
  AliAnalysisTaskSE(),
  fFilter(NULL),
  fMCcut(NULL),
  fAcceptanceCuts(NULL),
  fEfficiency(NULL),
  fOutput(NULL),
  fCutTRD(kFALSE)
{
  //
  // Default constructor
  //
}

AliHFEefficiency::AliHFEefficiency(const Char_t *name):
  AliAnalysisTaskSE(name),
  fFilter(NULL),
  fMCcut(NULL),
  fAcceptanceCuts(NULL),
  fEfficiency(NULL),
  fOutput(NULL),
  fCutTRD(kFALSE)
{
  //
  // Main Constructor
  //
  DefineOutput(1, AliHFEcontainer::Class());
  DefineOutput(2, TList::Class());

  SetRunTerminate();
}

AliHFEefficiency::~AliHFEefficiency(){
  //
  // Destructor
  //
  if(fFilter) delete fFilter;
  if(fEfficiency) delete fEfficiency;
  if(fAcceptanceCuts) delete fAcceptanceCuts;
  if(fOutput) delete fOutput;
}

void AliHFEefficiency::UserCreateOutputObjects(){
  //
  // Create the output objects
  //
  fEfficiency = new AliHFEcontainer("Efficiency");
  fEfficiency->SetNumberOfVariables(4);
  // InitializeVariables
  fEfficiency->SetVariableName(0, "pt");
  fEfficiency->MakeLogarithmicBinning(0, 40, 0.1, 10);
  fEfficiency->SetVariableName(1, "Eta");
  fEfficiency->MakeLinearBinning(1, 20, -0.9, 0.9);
  fEfficiency->SetVariableName(2, "phi");
  fEfficiency->MakeLinearBinning(2 , 18, 0, 2 * TMath::Pi());
  fEfficiency->SetVariableName(3, "Charge");
  fEfficiency->MakeLinearBinning(3, 2, -1.1, 1.1);
  // New container
  fEfficiency->CreateContainer("MCFilter", "Efficiency for MC Studies", 2);
  AliCFContainer *cont = fEfficiency->GetCFContainer("MCFilter");
  cont->SetStepTitle(0, "MC Signal");
  cont->SetStepTitle(1, "MC Signal in acceptance");

  fFilter = new AliHFEtrackFilter("testfilter");
  fFilter->GenerateCutSteps();
  fFilter->InitCF(fEfficiency);
  fMCcut = fFilter->GetMCSignalCuts();
  if(fCutTRD){
    AliHFEcutStep *hfeTRD = fFilter->GetCutStep("HFETRD");
    AliHFEextraCuts *hfecuts = dynamic_cast<AliHFEextraCuts *>(hfeTRD->GetCut("HFETRDCuts"));
    if(hfecuts) hfecuts->SetMinTrackletsTRD(4);
  }
  AliHFEcutStep *hfeITS = fFilter->GetCutStep("HFEITS");
  AliHFEextraCuts *hfeitscuts = dynamic_cast<AliHFEextraCuts *>(hfeITS->GetCut("HFEPixelsCuts"));
  if(hfeitscuts)hfeitscuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);

  AliHFEcutStep *primary = fFilter->GetCutStep("Primary");
  AliCFTrackIsPrimaryCuts *primcuts = dynamic_cast<AliCFTrackIsPrimaryCuts *>(primary->GetCut("PrimaryCuts"));
  if(primcuts){ 
    primcuts->SetMaxDCAToVertexXY(3);
    primcuts->SetMaxDCAToVertexZ(5);
  }

  fAcceptanceCuts = new AliCFAcceptanceCuts("Acceptance", "MC Acceptance Cuts");
  fAcceptanceCuts->SetMinNHitITS(3);
  fAcceptanceCuts->SetMinNHitTPC(2);
  fAcceptanceCuts->SetMinNHitTRD(0);

  // create additional histos for testing purpose
  fOutput = new AliHFEcollection("histos", "QA histograms");
  fOutput->CreateTH1F("hNtracks","Number of Tracks per Event", 100, 0, 100);
  fOutput->CreateTH1F("hPt","Pt of the tracks", 40, 0.1, 10, 0);
  fOutput->CreateTH1F("kinkIndex", "Kink Index", 400, -200, 200);
  fOutput->CreateTH1F("itspixel", "ITS PIXEL", 2, 0, 2);
  fOutput->CreateTH1F("mcmother", "Mother PDG", 1000, -500, 500);
  fOutput->CreateTH2F("ptres", "Pt Resolution", 40, 0.1, 10, 200, -1.5, 1.5, 0);
}

void AliHFEefficiency::UserExec(Option_t *){
  //
  // Event Loop
  // Filter track, fill Efficiency container
  //
  AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fInputEvent);
  if(!esdevent){
    AliError("ESD Event required");
    return;
  }
  fEfficiency->NewEvent();
  fFilter->SetRecEvent(fInputEvent);
  if(fMCEvent){
    AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
   if ( ! mcH ) {
     AliError("Cannot get MC truth event handler");
     return;
    }  
    if(mcH &&(!mcH->InitOk())) return;
    if(mcH &&(!mcH->TreeK())) return;
    if(mcH &&(!mcH->TreeTR())) return;
    fFilter->SetMC(fMCEvent);
    FilterMC();
  }
  fFilter->FilterTracks(esdevent);
  TObjArray *tracks = fFilter->GetFilteredTracks();
  TIterator *iter = tracks->MakeIterator();
  fOutput->Fill("hNtracks", tracks->GetEntriesFast());
  AliESDtrack *track = NULL;
  while((track = dynamic_cast<AliESDtrack *>(iter->Next()))){
    fOutput->Fill("hPt", track->Pt());  
    fOutput->Fill("kinkIndex", track->GetKinkIndex(0));
    if(TESTBIT(track->GetITSClusterMap(), 0))
      fOutput->Fill("itspixel",1);
    else
      fOutput->Fill("itspixel",0);
    if(fMCEvent){
      AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
      if(mctrack){
        Int_t motherLabel = mctrack->Particle()->GetFirstMother();
        if(motherLabel){
          AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(motherLabel));
          if(mother)fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
        }
        fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
      }
    }
  }
  delete iter;
  fFilter->Flush();
  PostData(1, fEfficiency);
  PostData(2, fOutput->GetList());
}

void AliHFEefficiency::Terminate(Option_t *){
  //
  // Evaluate Results
  //
  fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
  if(!fEfficiency){
    AliError("No Output data available");
    return;
  }

  if(!IsRunningTerminate()) return;
  PostProcess();

  TList *l = dynamic_cast<TList *>(GetOutputData(2));
  if(l) DrawPtResolution(l);
}

void AliHFEefficiency::FilterMC(){
  //
  // Cut MC tracks
  //
  Double_t cont[4] = {0., 0., 0., 0.};
  AliMCParticle *track = NULL;
  for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
    track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
    if(!track) continue;
    if(!fMCcut->IsSelected(track)) continue;
    cont[0] = track->Pt();
    cont[1] = track->Eta();
    cont[2] = track->Phi();
    cont[3] = track->Charge()/3;
    //AliMCParticle *mother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(track->Particle()->GetFirstMother()));
    //if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) continue;
    fEfficiency->FillCFContainer("MCFilter", 0, cont);
    if(!fAcceptanceCuts->IsSelected(track)) continue;
    fEfficiency->FillCFContainer("MCFilter", 1, cont);
  }
}

void AliHFEefficiency::Load(const char* filename){
  //
  // Load results for post processing
  //
  TFile *input = TFile::Open(filename);
  AliHFEcontainer *cin = dynamic_cast<AliHFEcontainer *>(input->Get("Efficiency"));
  if(cin) fEfficiency = dynamic_cast<AliHFEcontainer *>(cin->Clone());
  input->Close();
  delete input;
}

void AliHFEefficiency::PostProcess(){
  //
  // Calculate Efficiencies
  //
  fEfficiency->Print();

  // Prepare containers
  AliCFContainer *crec  = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signal");
  AliCFContainer *cmc   = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signalMC");
  AliCFEffGrid *gridRec = new AliCFEffGrid("effrec", "Efficiency Calculation using rec vars", *crec);
  AliCFEffGrid *gridMC  = new AliCFEffGrid("effmc", "Efficiency Calculation using mc vars", *cmc);

  TCanvas *ctmp = NULL;
  for(Int_t ivar = 0; ivar < crec->GetNVar(); ivar++){
    ctmp = new TCanvas(Form("cEffSignal%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s)", crec->GetVarTitle(ivar)));
    ctmp->cd();
    DrawSignalEfficiency(gridRec, crec, ivar);
    
    ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
    ctmp->cd();
    DrawCutEfficiency(gridRec, crec, ivar);
  
    ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
    ctmp->cd();
    DrawSignalEfficiency(gridMC, cmc, ivar);

    ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
    ctmp->cd();
    DrawCutEfficiency(gridMC, cmc, ivar);
  } 

  CalculatePTsmearing();

  delete gridRec;
  delete gridMC;
  delete crec;
  delete cmc;
}

void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
  //
  // Efficiency of a cut step with respect to the signal
  //
  TH1 *hEff = NULL;
  Bool_t kFirst = kTRUE;
  TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
  for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
    eff->CalculateEfficiency(istep, 0);
    hEff = eff->Project(var);
    hEff->SetTitle(Form("Signal Efficiency (%s)", cont->GetVarTitle(var)));
    hEff->GetYaxis()->SetRangeUser(0., 1.6);
    hEff->GetYaxis()->SetTitle("Efficiency");
    hEff->SetStats(kFALSE);
    hEff->SetMarkerStyle(20+istep);
    hEff->SetMarkerColor(kBlue - 5 + istep);
    hEff->Draw(kFirst ? "ep" : "epsame");
    kFirst = kFALSE;
    leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
  }
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->Draw();
  gPad->Update();
}

void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid *eff, AliCFContainer *cont, Int_t var){
  //
  // Efficiency of a cut step with respect to the previous one
  //
  Bool_t kFirst = kTRUE;
  TH1 *hEff = NULL;
  TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
  TString effTitle;
  Int_t start = 0, length = 0;
  for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
    eff->CalculateEfficiency(istep, istep -1);
    hEff = eff->Project(var);
    effTitle = hEff->GetTitle();
    start = effTitle.Index("projected"); length = effTitle.Length() - start;
    effTitle.Remove(start - 1, length + 1);
    effTitle.ReplaceAll("Efficiency: ", "");
    hEff->SetTitle(Form("Cut Efficiency (%s)", cont->GetVarTitle(var)));
    hEff->GetYaxis()->SetRangeUser(0., 1.6);
    hEff->GetYaxis()->SetTitle("Efficiency");
    hEff->SetStats(kFALSE);
    hEff->SetMarkerStyle(20+istep);
    hEff->SetMarkerColor(kBlue - 5 + istep);
    hEff->Draw(kFirst ? "ep" : "epsame");
    kFirst = kFALSE;
    leg->AddEntry(hEff, effTitle.Data(), "p");
  }
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->Draw();
  gPad->Update();
}

void AliHFEefficiency::CalculatePTsmearing(){
  //
  // Efficiency of a cut step with respect to the same cut step filled with MC 
  // Information
  //
  AliCFContainer *cont = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "testfilter_container_signal:testfilter_container_signalMC");
  AliCFEffGrid *grid = new AliCFEffGrid("effc", "Efficiency Calculation", *cont);
  Bool_t kFirst = kTRUE;
  TH1 *hEff = NULL;
  TCanvas *c2 = new TCanvas("cSmear", "pT modification");
  c2->cd();
  TLegend *leg1 = new TLegend(0.5, 0.7, 0.89, 0.89);
  Int_t nStep = cont->GetNStep()/2;
  for(Int_t istep = 0; istep < nStep; istep++){
    grid->CalculateEfficiency(istep, istep + nStep);
    hEff = grid->Project(0);
    hEff->SetStats(kFALSE);
    hEff->SetMarkerStyle(20+istep);
    hEff->SetMarkerColor(kBlue - 5 + istep);
    hEff->Draw(kFirst ? "ep" : "epsame");
    kFirst = kFALSE;
    leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
  }
  leg1->Draw();
  gPad->Update();

  delete cont;
  delete grid;
}

void AliHFEefficiency::DrawPtResolution(const TList * const l){
  //
  // Draw pt resolution
  //
  TH2 *h2 = dynamic_cast<TH2 *>(l->FindObject("ptres"));
  if(!h2) return;
  TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
  TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());

  Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
  TH1 *htmp = NULL;
  for(Int_t imom = 1; imom <= h2->GetXaxis()->GetLast(); imom++){
    pt = h2->GetXaxis()->GetBinCenter(imom);
    pterr = h2->GetXaxis()->GetBinWidth(imom)/2.;
    htmp = h2->ProjectionY("py", imom, imom);
    mymean = htmp->GetMean();
    myrms = htmp->GetRMS();
    mymeanerr = htmp->GetMeanError();
    myrmserr = htmp->GetRMSError();
    delete htmp;
    mean->SetPoint(imom -1, pt, mymean);
    rms->SetPoint(imom -1, pt, myrms);
    mean->SetPointError(imom-1,pterr,mymeanerr);
    rms->SetPointError(imom-1,pterr,myrmserr);
  }

  TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
  c1->cd();
  mean->SetMarkerStyle(22);
  mean->SetMarkerColor(kBlue);
  mean->Draw("ape");
  rms->SetMarkerStyle(22);
  rms->SetMarkerColor(kBlack);
  rms->Draw("pesame");
  TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
  leg->AddEntry(mean, "Mean", "lp");
  leg->AddEntry(rms, "RMS", "lp");
  leg->Draw();
  gPad->Update();
}
Commit	Line	Data
70da6c5a	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	// Task for Efficiency studies
	17	// Used for testing classes AliHFEcontainer and AliHFEfilter
	18	// Creates Efficiency Histograms
	19	//
	20	// Author:
	21	// Markus Fasel <M.Fasel@gsi.de>
	22	//
	23	#include <TAxis.h>
	24	#include <TCanvas.h>
	25	#include <TFile.h>
	26	#include <TGraphErrors.h>
	27	#include <TH1F.h>
	28	#include <TH2.h>
	29	#include <TIterator.h>
	30	#include <TLegend.h>
	31	#include <TObjArray.h>
	32	#include <TPad.h>
	33
91c7e1ec	34	#include "AliAnalysisManager.h"
70da6c5a	35	#include "AliCFAcceptanceCuts.h"
91c7e1ec	36	#include "AliCFTrackIsPrimaryCuts.h"
70da6c5a	37	#include "AliCFContainer.h"
	38	#include "AliCFEffGrid.h"
	39	#include "AliESDEvent.h"
	40	#include "AliHFEcollection.h"
	41	#include "AliHFEcontainer.h"
	42	#include "AliHFEcutStep.h"
	43	#include "AliHFEefficiency.h"
	44	#include "AliHFEextraCuts.h"
	45	#include "AliHFEtrackFilter.h"
	46	#include "AliHFEtools.h"
	47	#include "AliMCEvent.h"
91c7e1ec	48	#include "AliMCEventHandler.h"
70da6c5a	49
	50	ClassImp(AliHFEefficiency)
	51
	52	AliHFEefficiency::AliHFEefficiency():
	53	AliAnalysisTaskSE(),
	54	fFilter(NULL),
	55	fMCcut(NULL),
	56	fAcceptanceCuts(NULL),
	57	fEfficiency(NULL),
	58	fOutput(NULL),
	59	fCutTRD(kFALSE)
	60	{
	61	//
	62	// Default constructor
	63	//
	64	}
	65
	66	AliHFEefficiency::AliHFEefficiency(const Char_t *name):
	67	AliAnalysisTaskSE(name),
	68	fFilter(NULL),
	69	fMCcut(NULL),
	70	fAcceptanceCuts(NULL),
	71	fEfficiency(NULL),
	72	fOutput(NULL),
	73	fCutTRD(kFALSE)
	74	{
	75	//
	76	// Main Constructor
	77	//
	78	DefineOutput(1, AliHFEcontainer::Class());
	79	DefineOutput(2, TList::Class());
faee3b18	80
faee3b18	81	SetRunTerminate();
70da6c5a	82	}
	83
	84	AliHFEefficiency::~AliHFEefficiency(){
	85	//
	86	// Destructor
	87	//
	88	if(fFilter) delete fFilter;
	89	if(fEfficiency) delete fEfficiency;
	90	if(fAcceptanceCuts) delete fAcceptanceCuts;
	91	if(fOutput) delete fOutput;
	92	}
	93
	94	void AliHFEefficiency::UserCreateOutputObjects(){
	95	//
	96	// Create the output objects
	97	//
	98	fEfficiency = new AliHFEcontainer("Efficiency");
	99	fEfficiency->SetNumberOfVariables(4);
	100	// InitializeVariables
	101	fEfficiency->SetVariableName(0, "pt");
	102	fEfficiency->MakeLogarithmicBinning(0, 40, 0.1, 10);
	103	fEfficiency->SetVariableName(1, "Eta");
	104	fEfficiency->MakeLinearBinning(1, 20, -0.9, 0.9);
	105	fEfficiency->SetVariableName(2, "phi");
	106	fEfficiency->MakeLinearBinning(2 , 18, 0, 2 * TMath::Pi());
	107	fEfficiency->SetVariableName(3, "Charge");
	108	fEfficiency->MakeLinearBinning(3, 2, -1.1, 1.1);
	109	// New container
	110	fEfficiency->CreateContainer("MCFilter", "Efficiency for MC Studies", 2);
	111	AliCFContainer *cont = fEfficiency->GetCFContainer("MCFilter");
	112	cont->SetStepTitle(0, "MC Signal");
	113	cont->SetStepTitle(1, "MC Signal in acceptance");
	114
	115	fFilter = new AliHFEtrackFilter("testfilter");
	116	fFilter->GenerateCutSteps();
	117	fFilter->InitCF(fEfficiency);
	118	fMCcut = fFilter->GetMCSignalCuts();
	119	if(fCutTRD){
	120	AliHFEcutStep *hfeTRD = fFilter->GetCutStep("HFETRD");
	121	AliHFEextraCuts hfecuts = dynamic_cast<AliHFEextraCuts >(hfeTRD->GetCut("HFETRDCuts"));
bf892a6a	122	if(hfecuts) hfecuts->SetMinTrackletsTRD(4);
70da6c5a	123	}
	124	AliHFEcutStep *hfeITS = fFilter->GetCutStep("HFEITS");
	125	AliHFEextraCuts hfeitscuts = dynamic_cast<AliHFEextraCuts >(hfeITS->GetCut("HFEPixelsCuts"));
bf892a6a	126	if(hfeitscuts)hfeitscuts->SetRequireITSpixel(AliHFEextraCuts::kFirst);
70da6c5a	127
91c7e1ec	128	AliHFEcutStep *primary = fFilter->GetCutStep("Primary");
91c7e1ec	129	AliCFTrackIsPrimaryCuts primcuts = dynamic_cast<AliCFTrackIsPrimaryCuts >(primary->GetCut("PrimaryCuts"));
bf892a6a	130	if(primcuts){
	131	primcuts->SetMaxDCAToVertexXY(3);
	132	primcuts->SetMaxDCAToVertexZ(5);
	133	}
91c7e1ec	134
70da6c5a	135	fAcceptanceCuts = new AliCFAcceptanceCuts("Acceptance", "MC Acceptance Cuts");
	136	fAcceptanceCuts->SetMinNHitITS(3);
	137	fAcceptanceCuts->SetMinNHitTPC(2);
	138	fAcceptanceCuts->SetMinNHitTRD(0);
	139
	140	// create additional histos for testing purpose
	141	fOutput = new AliHFEcollection("histos", "QA histograms");
	142	fOutput->CreateTH1F("hNtracks","Number of Tracks per Event", 100, 0, 100);
faee3b18	143	fOutput->CreateTH1F("hPt","Pt of the tracks", 40, 0.1, 10, 0);
70da6c5a	144	fOutput->CreateTH1F("kinkIndex", "Kink Index", 400, -200, 200);
	145	fOutput->CreateTH1F("itspixel", "ITS PIXEL", 2, 0, 2);
	146	fOutput->CreateTH1F("mcmother", "Mother PDG", 1000, -500, 500);
faee3b18	147	fOutput->CreateTH2F("ptres", "Pt Resolution", 40, 0.1, 10, 200, -1.5, 1.5, 0);
70da6c5a	148	}
	149
	150	void AliHFEefficiency::UserExec(Option_t *){
	151	//
	152	// Event Loop
	153	// Filter track, fill Efficiency container
	154	//
e3ae862b	155	AliESDEvent esdevent = dynamic_cast<AliESDEvent >(fInputEvent);
	156	if(!esdevent){
	157	AliError("ESD Event required");
	158	return;
	159	}
70da6c5a	160	fEfficiency->NewEvent();
	161	fFilter->SetRecEvent(fInputEvent);
	162	if(fMCEvent){
91c7e1ec	163	AliMCEventHandler mcH = dynamic_cast<AliMCEventHandler >(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
215ffe88	164	if ( ! mcH ) {
	165	AliError("Cannot get MC truth event handler");
	166	return;
	167	}
e3ae862b	168	if(mcH &&(!mcH->InitOk())) return;
	169	if(mcH &&(!mcH->TreeK())) return;
	170	if(mcH &&(!mcH->TreeTR())) return;
70da6c5a	171	fFilter->SetMC(fMCEvent);
	172	FilterMC();
	173	}
e3ae862b	174	fFilter->FilterTracks(esdevent);
70da6c5a	175	TObjArray *tracks = fFilter->GetFilteredTracks();
	176	TIterator *iter = tracks->MakeIterator();
	177	fOutput->Fill("hNtracks", tracks->GetEntriesFast());
	178	AliESDtrack *track = NULL;
	179	while((track = dynamic_cast<AliESDtrack *>(iter->Next()))){
	180	fOutput->Fill("hPt", track->Pt());
	181	fOutput->Fill("kinkIndex", track->GetKinkIndex(0));
	182	if(TESTBIT(track->GetITSClusterMap(), 0))
	183	fOutput->Fill("itspixel",1);
	184	else
	185	fOutput->Fill("itspixel",0);
faee3b18	186	if(fMCEvent){
	187	AliMCParticle mctrack = dynamic_cast<AliMCParticle >(fMCEvent->GetTrack(TMath::Abs(track->GetLabel())));
	188	if(mctrack){
	189	Int_t motherLabel = mctrack->Particle()->GetFirstMother();
	190	if(motherLabel){
	191	AliMCParticle mother = dynamic_cast<AliMCParticle >(fMCEvent->GetTrack(motherLabel));
	192	if(mother)fOutput->Fill("mcmother", mother->Particle()->GetPdgCode());
	193	}
	194	fOutput->Fill("ptres", mctrack->Pt(), (track->Pt() - mctrack->Pt())/mctrack->Pt());
91c7e1ec	195	}
91c7e1ec	196	}
70da6c5a	197	}
	198	delete iter;
	199	fFilter->Flush();
	200	PostData(1, fEfficiency);
	201	PostData(2, fOutput->GetList());
	202	}
	203
	204	void AliHFEefficiency::Terminate(Option_t *){
	205	//
	206	// Evaluate Results
	207	//
	208	fEfficiency = dynamic_cast<AliHFEcontainer *>(GetOutputData(1));
faee3b18	209	if(!fEfficiency){
	210	AliError("No Output data available");
	211	return;
	212	}
	213
	214	if(!IsRunningTerminate()) return;
70da6c5a	215	PostProcess();
	216
	217	TList l = dynamic_cast<TList >(GetOutputData(2));
bf892a6a	218	if(l) DrawPtResolution(l);
70da6c5a	219	}
	220
	221	void AliHFEefficiency::FilterMC(){
	222	//
	223	// Cut MC tracks
	224	//
	225	Double_t cont[4] = {0., 0., 0., 0.};
	226	AliMCParticle *track = NULL;
	227	for(Int_t itrack = 0; itrack < fMCEvent->GetNumberOfTracks(); itrack++){
	228	track = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(itrack));
e3ae862b	229	if(!track) continue;
70da6c5a	230	if(!fMCcut->IsSelected(track)) continue;
	231	cont[0] = track->Pt();
	232	cont[1] = track->Eta();
	233	cont[2] = track->Phi();
	234	cont[3] = track->Charge()/3;
	235	//AliMCParticle mother = dynamic_cast<AliMCParticle >(fMCEvent->GetTrack(track->Particle()->GetFirstMother()));
	236	//if(TMath::Abs(mother->Particle()->GetPdgCode()) != 443) continue;
	237	fEfficiency->FillCFContainer("MCFilter", 0, cont);
	238	if(!fAcceptanceCuts->IsSelected(track)) continue;
	239	fEfficiency->FillCFContainer("MCFilter", 1, cont);
	240	}
	241	}
	242
	243	void AliHFEefficiency::Load(const char* filename){
3a72645a	244	//
	245	// Load results for post processing
	246	//
70da6c5a	247	TFile *input = TFile::Open(filename);
70da6c5a	248	AliHFEcontainer cin = dynamic_cast<AliHFEcontainer >(input->Get("Efficiency"));
bf892a6a	249	if(cin) fEfficiency = dynamic_cast<AliHFEcontainer *>(cin->Clone());
70da6c5a	250	input->Close();
	251	delete input;
	252	}
	253
	254	void AliHFEefficiency::PostProcess(){
	255	//
	256	// Calculate Efficiencies
	257	//
	258	fEfficiency->Print();
	259
	260	// Prepare containers
	261	AliCFContainer *crec = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signal");
	262	AliCFContainer *cmc = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "MCFilter:testfilter_container_signalMC");
	263	AliCFEffGrid gridRec = new AliCFEffGrid("effrec", "Efficiency Calculation using rec vars", crec);
	264	AliCFEffGrid gridMC = new AliCFEffGrid("effmc", "Efficiency Calculation using mc vars", cmc);
	265
	266	TCanvas *ctmp = NULL;
	267	for(Int_t ivar = 0; ivar < crec->GetNVar(); ivar++){
	268	ctmp = new TCanvas(Form("cEffSignal%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s)", crec->GetVarTitle(ivar)));
	269	ctmp->cd();
	270	DrawSignalEfficiency(gridRec, crec, ivar);
	271
	272	ctmp = new TCanvas(Form("cCutEff%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s)", crec->GetVarTitle(ivar)));
	273	ctmp->cd();
	274	DrawCutEfficiency(gridRec, crec, ivar);
	275
	276	ctmp = new TCanvas(Form("cEffSignalMC%s", crec->GetVarTitle(ivar)), Form("Signal Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
	277	ctmp->cd();
	278	DrawSignalEfficiency(gridMC, cmc, ivar);
	279
	280	ctmp = new TCanvas(Form("cCutEffMC%s", crec->GetVarTitle(ivar)), Form("Cut Step Efficiency (%s, MC)", crec->GetVarTitle(ivar)));
	281	ctmp->cd();
	282	DrawCutEfficiency(gridMC, cmc, ivar);
	283	}
	284
	285	CalculatePTsmearing();
	286
	287	delete gridRec;
	288	delete gridMC;
	289	delete crec;
	290	delete cmc;
	291	}
	292
	293	void AliHFEefficiency::DrawSignalEfficiency(AliCFEffGrid eff, AliCFContainer cont, Int_t var){
	294	//
	295	// Efficiency of a cut step with respect to the signal
	296	//
	297	TH1 *hEff = NULL;
	298	Bool_t kFirst = kTRUE;
	299	TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
	300	for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
	301	eff->CalculateEfficiency(istep, 0);
	302	hEff = eff->Project(var);
	303	hEff->SetTitle(Form("Signal Efficiency (%s)", cont->GetVarTitle(var)));
	304	hEff->GetYaxis()->SetRangeUser(0., 1.6);
	305	hEff->GetYaxis()->SetTitle("Efficiency");
	306	hEff->SetStats(kFALSE);
	307	hEff->SetMarkerStyle(20+istep);
	308	hEff->SetMarkerColor(kBlue - 5 + istep);
	309	hEff->Draw(kFirst ? "ep" : "epsame");
	310	kFirst = kFALSE;
	311	leg->AddEntry(hEff, cont->GetStepTitle(istep), "p");
	312	}
	313	leg->SetBorderSize(0);
314	leg->SetFillStyle(0);
315	leg->Draw();
316	gPad->Update();
317	}
318
319	void AliHFEefficiency::DrawCutEfficiency(AliCFEffGrid eff, AliCFContainer cont, Int_t var){
320	//
321	// Efficiency of a cut step with respect to the previous one
322	//
323	Bool_t kFirst = kTRUE;
324	TH1 *hEff = NULL;
325	TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
326	TString effTitle;
327	Int_t start = 0, length = 0;
328	for(Int_t istep = 1; istep < cont->GetNStep(); istep++){
329	eff->CalculateEfficiency(istep, istep -1);
330	hEff = eff->Project(var);
331	effTitle = hEff->GetTitle();
332	start = effTitle.Index("projected"); length = effTitle.Length() - start;
333	effTitle.Remove(start - 1, length + 1);
334	effTitle.ReplaceAll("Efficiency: ", "");
335	hEff->SetTitle(Form("Cut Efficiency (%s)", cont->GetVarTitle(var)));
336	hEff->GetYaxis()->SetRangeUser(0., 1.6);
337	hEff->GetYaxis()->SetTitle("Efficiency");
338	hEff->SetStats(kFALSE);
339	hEff->SetMarkerStyle(20+istep);
340	hEff->SetMarkerColor(kBlue - 5 + istep);
341	hEff->Draw(kFirst ? "ep" : "epsame");
342	kFirst = kFALSE;
343	leg->AddEntry(hEff, effTitle.Data(), "p");
344	}
345	leg->SetBorderSize(0);
346	leg->SetFillStyle(0);
347	leg->Draw();
348	gPad->Update();
349	}
350
351	void AliHFEefficiency::CalculatePTsmearing(){
352	//
353	// Efficiency of a cut step with respect to the same cut step filled with MC
354	// Information
355	//
356	AliCFContainer *cont = fEfficiency->MakeMergedCFContainer("merged", "Merged Container", "testfilter_container_signal:testfilter_container_signalMC");
357	AliCFEffGrid grid = new AliCFEffGrid("effc", "Efficiency Calculation", cont);
358	Bool_t kFirst = kTRUE;
359	TH1 *hEff = NULL;
360	TCanvas *c2 = new TCanvas("cSmear", "pT modification");
361	c2->cd();
362	TLegend *leg1 = new TLegend(0.5, 0.7, 0.89, 0.89);
363	Int_t nStep = cont->GetNStep()/2;
364	for(Int_t istep = 0; istep < nStep; istep++){
365	grid->CalculateEfficiency(istep, istep + nStep);
366	hEff = grid->Project(0);
367	hEff->SetStats(kFALSE);
368	hEff->SetMarkerStyle(20+istep);
369	hEff->SetMarkerColor(kBlue - 5 + istep);
370	hEff->Draw(kFirst ? "ep" : "epsame");
371	kFirst = kFALSE;
372	leg1->AddEntry(hEff, cont->GetStepTitle(istep), "p");
373	}
374	leg1->Draw();
375	gPad->Update();
376
377	delete cont;
378	delete grid;
379	}
380
3a72645a	381	void AliHFEefficiency::DrawPtResolution(const TList * const l){
70da6c5a	382	//
	383	// Draw pt resolution
	384	//
	385	TH2 h2 = dynamic_cast<TH2 >(l->FindObject("ptres"));
bf892a6a	386	if(!h2) return;
70da6c5a	387	TGraphErrors *mean = new TGraphErrors(h2->GetNbinsX());
	388	TGraphErrors *rms = new TGraphErrors(h2->GetNbinsX());
	389
	390	Double_t pt = 0., pterr = 0., mymean = 0., myrms = 0., mymeanerr = 0., myrmserr = 0.;
	391	TH1 *htmp = NULL;
	392	for(Int_t imom = 1; imom <= h2->GetXaxis()->GetLast(); imom++){
	393	pt = h2->GetXaxis()->GetBinCenter(imom);
	394	pterr = h2->GetXaxis()->GetBinWidth(imom)/2.;
	395	htmp = h2->ProjectionY("py", imom, imom);
	396	mymean = htmp->GetMean();
	397	myrms = htmp->GetRMS();
	398	mymeanerr = htmp->GetMeanError();
	399	myrmserr = htmp->GetRMSError();
	400	delete htmp;
	401	mean->SetPoint(imom -1, pt, mymean);
	402	rms->SetPoint(imom -1, pt, myrms);
	403	mean->SetPointError(imom-1,pterr,mymeanerr);
	404	rms->SetPointError(imom-1,pterr,myrmserr);
	405	}
	406
	407	TCanvas *c1 = new TCanvas("cPtShift", "pT Shift");
	408	c1->cd();
	409	mean->SetMarkerStyle(22);
	410	mean->SetMarkerColor(kBlue);
	411	mean->Draw("ape");
	412	rms->SetMarkerStyle(22);
	413	rms->SetMarkerColor(kBlack);
	414	rms->Draw("pesame");
	415	TLegend *leg = new TLegend(0.5, 0.7, 0.89, 0.89);
	416	leg->AddEntry(mean, "Mean", "lp");
	417	leg->AddEntry(rms, "RMS", "lp");
	418	leg->Draw();
	419	gPad->Update();
	420	}
	421