[u/mrichter/AliRoot.git] / ANALYSIS / AliESDpidCuts.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.                  *
 **************************************************************************/
#include <TCanvas.h>
#include <TClass.h>
#include <TCollection.h>
#include <TDirectory.h>
#include <TH1F.h>
#include <TH1I.h>
#include <TH2I.h>
#include <TMath.h>
#include <TIterator.h>
#include <TString.h>
#include <TList.h>

#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliESDtrack.h"
#include "AliESDEvent.h"
#include "AliLog.h"
#include "AliPIDResponse.h"

#include "AliESDpidCuts.h"

ClassImp(AliESDpidCuts)

const Int_t AliESDpidCuts::kNcuts = 3;

//_____________________________________________________________________
AliESDpidCuts::AliESDpidCuts(const Char_t *name, const Char_t *title):
    AliAnalysisCuts(name, title)
  , fPIDresponse(NULL)
  , fTPCsigmaCutRequired(0)
  , fTOFsigmaCutRequired(0)
  , fCutTPCclusterRatio(0.)
  , fMinMomentumTOF(0.5)
  , fHcutStatistics(NULL)
  , fHcutCorrelation(NULL)
{
  //
  // Default constructor
  //
  
  memset(fCutTPCnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);
  memset(fCutTOFnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);

  memset(fHclusterRatio, 0, sizeof(TH1F *) * 2);
  memset(fHnSigmaTPC, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);
  memset(fHnSigmaTOF, 0, sizeof(TH1F *) * AliPID::kSPECIES * 2);
}

//_____________________________________________________________________
AliESDpidCuts::AliESDpidCuts(const AliESDpidCuts &ref):
    AliAnalysisCuts(ref)
  , fPIDresponse(NULL)
  , fTPCsigmaCutRequired(ref.fTPCsigmaCutRequired)
  , fTOFsigmaCutRequired(ref.fTOFsigmaCutRequired)
  , fCutTPCclusterRatio(ref.fCutTPCclusterRatio)
  , fMinMomentumTOF(ref.fMinMomentumTOF)
  , fHcutStatistics(NULL)
  , fHcutCorrelation(NULL)
{
  //
  // Copy constructor
  //
  fPIDresponse = ref.fPIDresponse;
  memcpy(fCutTPCnSigma, ref.fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
  memcpy(fCutTOFnSigma, ref.fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
  
  if(ref.fHcutStatistics) fHcutStatistics = dynamic_cast<TH1I *>(ref.fHcutStatistics->Clone());
  if(ref.fHcutCorrelation) fHcutCorrelation = dynamic_cast<TH2I *>(ref.fHcutCorrelation->Clone());
  for(Int_t imode = 0; imode < 2; imode++){
    if(ref.fHclusterRatio[imode]) fHclusterRatio[imode] = dynamic_cast<TH1F *>(ref.fHclusterRatio[imode]->Clone());
    for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
      if(fHnSigmaTPC[ispec][imode]) fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
      if(fHnSigmaTOF[ispec][imode]) fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
    }
  }
}

//_____________________________________________________________________
AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){
  //
  // Assignment operator
  //
  if(this != &ref)
    ref.Copy(*this);
  return *this;
}

//_____________________________________________________________________
AliESDpidCuts::~AliESDpidCuts(){
  //
  // Destructor
  //

  delete fHcutCorrelation;
  for(Int_t imode = 0; imode < 2; imode++){
    delete fHclusterRatio[imode];
    for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
      delete fHnSigmaTPC[ispec][imode];
      delete fHnSigmaTOF[ispec][imode];
    }
  }
}

//_____________________________________________________________________
void AliESDpidCuts::Init(){
  //
  // Init function, get PID response from the Analysis Manager
  //
  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
  if(mgr){
    AliInputEventHandler *handler = dynamic_cast<AliInputEventHandler*>(mgr->GetInputEventHandler());
    if(handler){
      fPIDresponse = handler->GetPIDResponse();
    }
  }
}

//_____________________________________________________________________
Bool_t AliESDpidCuts::IsSelected(TObject *obj){
  //
  // Select Track
  // 
  AliESDtrack * trk = dynamic_cast<AliESDtrack*>(obj);
  if(!trk){
    AliError("Provided object is not AliESDtrack!");
    return kFALSE;
  }
  const AliESDEvent* evt = trk->GetESDEvent();
  if(!evt){
    AliError("No AliESDEvent!");
    return kFALSE;
  }
  return AcceptTrack(trk, evt);
}

//_____________________________________________________________________
void AliESDpidCuts::Copy(TObject &c) const {
  //
  // Copy function
  //
  AliESDpidCuts &target = dynamic_cast<AliESDpidCuts &>(c);

  target.fPIDresponse = fPIDresponse;

  target.fCutTPCclusterRatio = fCutTPCclusterRatio;
  target.fMinMomentumTOF = fMinMomentumTOF;

  target.fTPCsigmaCutRequired = fTPCsigmaCutRequired;
  target.fTOFsigmaCutRequired = fTOFsigmaCutRequired;
  
  if(fHcutStatistics) target.fHcutStatistics = dynamic_cast<TH1I *>(fHcutStatistics->Clone());
  if(fHcutCorrelation) target.fHcutCorrelation = dynamic_cast<TH2I *>(fHcutCorrelation->Clone());
  for(Int_t imode = 0; imode < 2; imode++){
    if(fHclusterRatio[imode]) target.fHclusterRatio[imode] = dynamic_cast<TH1F *>(fHclusterRatio[imode]->Clone());
    for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
      if(fHnSigmaTPC[ispec][imode]) target.fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
      if(fHnSigmaTOF[ispec][imode]) target.fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTOF[ispec][imode]->Clone());
    }
  }
 
  memcpy(target.fCutTPCnSigma, fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
  memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
 
  AliESDpidCuts::Copy(c);
}

//_____________________________________________________________________
Long64_t AliESDpidCuts::Merge(TCollection *coll){
  //
  // Merge Cut objects
  //
  if(!coll) return 0;
  if(coll->IsEmpty())   return 1;
  if(!HasHistograms())  return 0;
  
  TIterator *iter = coll->MakeIterator();
  TObject *o = NULL; 
  AliESDpidCuts *ref = NULL;
  Int_t counter = 0;
  while((o = iter->Next())){
    ref = dynamic_cast<AliESDpidCuts *>(o);
    if(!ref) continue;
    if(!ref->HasHistograms()) continue;

    fHcutStatistics->Add(ref->fHcutStatistics);
    fHcutCorrelation->Add(ref->fHcutCorrelation);
    for(Int_t imode = 0; imode < 2; imode++){
      fHclusterRatio[imode]->Add(ref->fHclusterRatio[imode]);
      for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
        fHnSigmaTPC[ispec][imode]->Add(ref->fHnSigmaTPC[ispec][imode]);
        fHnSigmaTOF[ispec][imode]->Add(ref->fHnSigmaTOF[ispec][imode]);
      }
    }
    ++counter;
  }
  return ++counter;
}

//_____________________________________________________________________
void AliESDpidCuts::DefineHistograms(Color_t color){
  //
  // Swich on QA and create the histograms
  //
  SetBit(kHasHistograms, kTRUE);
  fHcutStatistics = new TH1I("fHcutStatistics", "Cut Statistics", kNcuts, 0, kNcuts);
  fHcutStatistics->SetLineColor(color);
  fHcutCorrelation = new TH2I("fHcutCorrelation", "Cut Correlation", kNcuts, 0, kNcuts, kNcuts, 0, kNcuts);
  TString cutname[kNcuts] = {"TPCclusterRatio", "TPC sigma", "TOF sigma"};
  for(Int_t icut = 0; icut < kNcuts; icut++){
    fHcutStatistics->GetXaxis()->SetBinLabel(fHcutStatistics->GetXaxis()->GetFirst() + icut, cutname[icut].Data());
    fHcutCorrelation->GetXaxis()->SetBinLabel(fHcutCorrelation->GetXaxis()->GetFirst() + icut, cutname[icut].Data());
    fHcutCorrelation->GetYaxis()->SetBinLabel(fHcutCorrelation->GetYaxis()->GetFirst() + icut, cutname[icut].Data());
  }
  Char_t hname[256], htitle[256];
  for(Int_t imode = 0; imode < 2; imode++){
    snprintf(hname, 256, "fHclusterRatio%s", imode ? "After" : "Before");
    snprintf(htitle, 256, "TPC cluster Ratio %s cuts;Ratio;Entries", imode ? "after" : "before");
    fHclusterRatio[imode] = new TH1F(hname, htitle, 20, 0., 1.);
    for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
      snprintf(hname, 256, "fHnSigma%sTPC%s", AliPID::ParticleName(ispec), imode ? "after" : "before");
      snprintf(htitle, 256, "TPC sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");
      fHnSigmaTPC[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);
      snprintf(hname, 256, "fHnSigma%sTOF%s", AliPID::ParticleName(ispec), imode ? "after" : "before");
      snprintf(htitle, 256, "TOF sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");
      fHnSigmaTOF[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);
    }
  }
}

//_____________________________________________________________________
Bool_t AliESDpidCuts::AcceptTrack(const AliESDtrack *track, const AliESDEvent *event){
  //
  // Check whether the tracks survived the cuts
  //
  if(!fPIDresponse){
    Init();
    if (!fPIDresponse)
      AliError("PID Response not available");
    return 0;
  }
  enum{
    kCutClusterRatioTPC,
    kCutNsigmaTPC,
    kCutNsigmaTOF
  };
  Long64_t cutRequired=0, cutFullfiled = 0;
  if(fTOFsigmaCutRequired && event == 0)  {
    AliError("No event pointer. Need event pointer for T0 for TOF cut");
    return (0);
  }
  Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
  if(fCutTPCclusterRatio > 0.){
    SETBIT(cutRequired, kCutClusterRatioTPC);
    if(clusterRatio >= fCutTPCclusterRatio) 
      SETBIT(cutFullfiled, kCutClusterRatioTPC);
  }
  // check TPC nSigma cut
  Float_t nsigmaTPC[AliPID::kSPECIES], nsigma;   // need all sigmas for QA plotting
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    nsigmaTPC[ispec] = nsigma = fPIDresponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
    if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;
    SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required
    if(nsigma >= fCutTPCnSigma[2*ispec] && nsigma <= fCutTPCnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC);    // Fullfiled for at least one species
  }
  // check TOF nSigma cut
  Float_t nsigmaTOF[AliPID::kSPECIES];    // see above
  Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set
  Double_t times[AliPID::kSPECIES];
  track->GetIntegratedTimes(times);
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    
    if(hasTOFpid && event) nsigmaTOF[ispec] = nsigma = fPIDresponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);
    if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;
    SETBIT(cutRequired, kCutNsigmaTOF);
    if(track->GetOuterParam()->P() >= fMinMomentumTOF){
      if(hasTOFpid && nsigma <= fCutTOFnSigma[2*ispec] && nsigma >= fCutTOFnSigma[2*ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);
    }
  }

  // Fill Histograms before cuts
  if(HasHistograms()){
    fHclusterRatio[0]->Fill(clusterRatio);
    for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
      fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);
      if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);
    }
  }
  if(cutRequired != cutFullfiled){
    // Fill cut statistics
    if(HasHistograms()){
      for(Int_t icut = 0; icut < kNcuts; icut++){
	if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){
	  // cut not fullfiled
	  fHcutStatistics->Fill(icut);
	  for(Int_t jcut = 0; jcut <= icut; jcut++)
	    if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);
	}
      }
    }
    return kFALSE;    // At least one cut is not fullfiled
  }

  // Fill Histograms after cuts
  if(HasHistograms()){
    fHclusterRatio[1]->Fill(clusterRatio);
    for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
      fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);
      if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);
    }
  }

  return kTRUE;
}

//_____________________________________________________________________
void AliESDpidCuts::SaveHistograms(const Char_t * location){
  //
  // Save the histograms to a file
  //
  if(!HasHistograms()){
    AliError("Histograms not on - Exiting");
    return;
  }
  if(!location) location = GetName();
  gDirectory->mkdir(location);
  gDirectory->cd(location);
  fHcutStatistics->Write();
  fHcutCorrelation->Write();

  gDirectory->mkdir("before_cuts");
  gDirectory->mkdir("after_cuts");

  gDirectory->cd("before_cuts");
  fHclusterRatio[0]->Write();
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    fHnSigmaTPC[ispec][0]->Write();
    fHnSigmaTOF[ispec][0]->Write();
  }

  gDirectory->cd("../after_cuts");
  fHclusterRatio[1]->Write();
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    fHnSigmaTPC[ispec][1]->Write();
    fHnSigmaTOF[ispec][1]->Write();
  }

  gDirectory->cd("..");
}

//_____________________________________________________________________
void AliESDpidCuts::DrawHistograms(){
  //
  // Draw the Histograms
  //
  TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);
  stat->cd();
  fHcutStatistics->SetStats(kFALSE);
  fHcutStatistics->Draw();
  stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));

  TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);
  correl->cd();
  fHcutCorrelation->SetStats(kFALSE);
  fHcutCorrelation->Draw("colz");
  correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));

  TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);
  cRatio->cd();
  fHclusterRatio[0]->SetLineColor(kRed);
  fHclusterRatio[0]->SetStats(kFALSE);
  fHclusterRatio[0]->Draw();
  fHclusterRatio[1]->SetLineColor(kBlue);
  fHclusterRatio[1]->SetStats(kFALSE);
  fHclusterRatio[1]->Draw("same");
  cRatio->SaveAs(Form("%s_%s.gif",  GetName(), cRatio->GetName()));

  TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);
  cNsigTPC->Divide(3,2);
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    cNsigTPC->cd(ispec + 1);
    fHnSigmaTPC[ispec][0]->SetLineColor(kRed);
    fHnSigmaTPC[ispec][0]->SetStats(kFALSE);
    fHnSigmaTPC[ispec][0]->Draw();
    fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);
    fHnSigmaTPC[ispec][1]->SetStats(kFALSE);
    fHnSigmaTPC[ispec][1]->Draw("same");
  }
  cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));

  TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);
  cNsigTOF->Divide(3,2);
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    cNsigTOF->cd(ispec + 1);
    fHnSigmaTOF[ispec][0]->SetLineColor(kRed);
    fHnSigmaTOF[ispec][0]->SetStats(kFALSE);
    fHnSigmaTOF[ispec][0]->Draw();
    fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);
    fHnSigmaTOF[ispec][1]->SetStats(kFALSE);
    fHnSigmaTOF[ispec][1]->Draw("same");
  }
  cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));
}
Commit	Line	Data
a65a7e70	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	#include <TCanvas.h>
	16	#include <TClass.h>
	17	#include <TCollection.h>
	18	#include <TDirectory.h>
	19	#include <TH1F.h>
	20	#include <TH1I.h>
	21	#include <TH2I.h>
	22	#include <TMath.h>
	23	#include <TIterator.h>
	24	#include <TString.h>
	25	#include <TList.h>
	26
	27	#include "AliAnalysisManager.h"
	28	#include "AliInputEventHandler.h"
	29	#include "AliESDtrack.h"
	30	#include "AliESDEvent.h"
	31	#include "AliLog.h"
	32	#include "AliPIDResponse.h"
	33
	34	#include "AliESDpidCuts.h"
	35
	36	ClassImp(AliESDpidCuts)
	37
	38	const Int_t AliESDpidCuts::kNcuts = 3;
	39
	40	//_____________________________________________________________________
	41	AliESDpidCuts::AliESDpidCuts(const Char_t name, const Char_t title):
	42	AliAnalysisCuts(name, title)
	43	, fPIDresponse(NULL)
	44	, fTPCsigmaCutRequired(0)
	45	, fTOFsigmaCutRequired(0)
	46	, fCutTPCclusterRatio(0.)
	47	, fMinMomentumTOF(0.5)
	48	, fHcutStatistics(NULL)
	49	, fHcutCorrelation(NULL)
	50	{
	51	//
	52	// Default constructor
	53	//
	54
	55	memset(fCutTPCnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);
	56	memset(fCutTOFnSigma, 0, sizeof(Float_t)* AliPID::kSPECIES * 2);
	57
	58	memset(fHclusterRatio, 0, sizeof(TH1F ) 2);
	59	memset(fHnSigmaTPC, 0, sizeof(TH1F ) AliPID::kSPECIES * 2);
	60	memset(fHnSigmaTOF, 0, sizeof(TH1F ) AliPID::kSPECIES * 2);
	61	}
	62
	63	//_____________________________________________________________________
	64	AliESDpidCuts::AliESDpidCuts(const AliESDpidCuts &ref):
65	AliAnalysisCuts(ref)
66	, fPIDresponse(NULL)
67	, fTPCsigmaCutRequired(ref.fTPCsigmaCutRequired)
68	, fTOFsigmaCutRequired(ref.fTOFsigmaCutRequired)
69	, fCutTPCclusterRatio(ref.fCutTPCclusterRatio)
70	, fMinMomentumTOF(ref.fMinMomentumTOF)
71	, fHcutStatistics(NULL)
72	, fHcutCorrelation(NULL)
73	{
74	//
75	// Copy constructor
76	//
77	fPIDresponse = ref.fPIDresponse;
78	memcpy(fCutTPCnSigma, ref.fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
79	memcpy(fCutTOFnSigma, ref.fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
80
81	if(ref.fHcutStatistics) fHcutStatistics = dynamic_cast<TH1I *>(ref.fHcutStatistics->Clone());
82	if(ref.fHcutCorrelation) fHcutCorrelation = dynamic_cast<TH2I *>(ref.fHcutCorrelation->Clone());
83	for(Int_t imode = 0; imode < 2; imode++){
84	if(ref.fHclusterRatio[imode]) fHclusterRatio[imode] = dynamic_cast<TH1F *>(ref.fHclusterRatio[imode]->Clone());
85	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
86	if(fHnSigmaTPC[ispec][imode]) fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
87	if(fHnSigmaTOF[ispec][imode]) fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
88	}
89	}
90	}
91
92	//_____________________________________________________________________
93	AliESDpidCuts &AliESDpidCuts::operator=(const AliESDpidCuts &ref){
94	//
95	// Assignment operator
96	//
97	if(this != &ref)
98	ref.Copy(*this);
99	return *this;
100	}
101
102	//_____________________________________________________________________
103	AliESDpidCuts::~AliESDpidCuts(){
104	//
105	// Destructor
106	//
107
108	delete fHcutCorrelation;
109	for(Int_t imode = 0; imode < 2; imode++){
110	delete fHclusterRatio[imode];
111	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
112	delete fHnSigmaTPC[ispec][imode];
113	delete fHnSigmaTOF[ispec][imode];
114	}
115	}
116	}
117
118	//_____________________________________________________________________
119	void AliESDpidCuts::Init(){
120	//
121	// Init function, get PID response from the Analysis Manager
122	//
123	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
124	if(mgr){
125	AliInputEventHandler handler = dynamic_cast<AliInputEventHandler>(mgr->GetInputEventHandler());
126	if(handler){
127	fPIDresponse = handler->GetPIDResponse();
128	}
129	}
130	}
131
132	//_____________________________________________________________________
133	Bool_t AliESDpidCuts::IsSelected(TObject *obj){
134	//
135	// Select Track
136	//
137	AliESDtrack * trk = dynamic_cast<AliESDtrack*>(obj);
138	if(!trk){
139	AliError("Provided object is not AliESDtrack!");
140	return kFALSE;
141	}
142	const AliESDEvent* evt = trk->GetESDEvent();
143	if(!evt){
144	AliError("No AliESDEvent!");
145	return kFALSE;
146	}
147	return AcceptTrack(trk, evt);
148	}
149
150	//_____________________________________________________________________
151	void AliESDpidCuts::Copy(TObject &c) const {
152	//
153	// Copy function
154	//
155	AliESDpidCuts &target = dynamic_cast<AliESDpidCuts &>(c);
156
157	target.fPIDresponse = fPIDresponse;
158
159	target.fCutTPCclusterRatio = fCutTPCclusterRatio;
160	target.fMinMomentumTOF = fMinMomentumTOF;
161
162	target.fTPCsigmaCutRequired = fTPCsigmaCutRequired;
163	target.fTOFsigmaCutRequired = fTOFsigmaCutRequired;
164
165	if(fHcutStatistics) target.fHcutStatistics = dynamic_cast<TH1I *>(fHcutStatistics->Clone());
166	if(fHcutCorrelation) target.fHcutCorrelation = dynamic_cast<TH2I *>(fHcutCorrelation->Clone());
167	for(Int_t imode = 0; imode < 2; imode++){
168	if(fHclusterRatio[imode]) target.fHclusterRatio[imode] = dynamic_cast<TH1F *>(fHclusterRatio[imode]->Clone());
169	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
170	if(fHnSigmaTPC[ispec][imode]) target.fHnSigmaTPC[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTPC[ispec][imode]->Clone());
171	if(fHnSigmaTOF[ispec][imode]) target.fHnSigmaTOF[ispec][imode] = dynamic_cast<TH1F *>(fHnSigmaTOF[ispec][imode]->Clone());
172	}
173	}
174
175	memcpy(target.fCutTPCnSigma, fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
176	memcpy(target.fCutTOFnSigma, fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES * 2);
177
178	AliESDpidCuts::Copy(c);
179	}
180
181	//_____________________________________________________________________
182	Long64_t AliESDpidCuts::Merge(TCollection *coll){
183	//
184	// Merge Cut objects
185	//
186	if(!coll) return 0;
187	if(coll->IsEmpty()) return 1;
188	if(!HasHistograms()) return 0;
189
190	TIterator *iter = coll->MakeIterator();
191	TObject *o = NULL;
192	AliESDpidCuts *ref = NULL;
193	Int_t counter = 0;
194	while((o = iter->Next())){
195	ref = dynamic_cast<AliESDpidCuts *>(o);
196	if(!ref) continue;
197	if(!ref->HasHistograms()) continue;
198
199	fHcutStatistics->Add(ref->fHcutStatistics);
200	fHcutCorrelation->Add(ref->fHcutCorrelation);
201	for(Int_t imode = 0; imode < 2; imode++){
202	fHclusterRatio[imode]->Add(ref->fHclusterRatio[imode]);
203	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
204	fHnSigmaTPC[ispec][imode]->Add(ref->fHnSigmaTPC[ispec][imode]);
205	fHnSigmaTOF[ispec][imode]->Add(ref->fHnSigmaTOF[ispec][imode]);
206	}
207	}
208	++counter;
209	}
210	return ++counter;
211	}
212
213	//_____________________________________________________________________
214	void AliESDpidCuts::DefineHistograms(Color_t color){
215	//
216	// Swich on QA and create the histograms
217	//
218	SetBit(kHasHistograms, kTRUE);
219	fHcutStatistics = new TH1I("fHcutStatistics", "Cut Statistics", kNcuts, 0, kNcuts);
220	fHcutStatistics->SetLineColor(color);
221	fHcutCorrelation = new TH2I("fHcutCorrelation", "Cut Correlation", kNcuts, 0, kNcuts, kNcuts, 0, kNcuts);
222	TString cutname[kNcuts] = {"TPCclusterRatio", "TPC sigma", "TOF sigma"};
223	for(Int_t icut = 0; icut < kNcuts; icut++){
224	fHcutStatistics->GetXaxis()->SetBinLabel(fHcutStatistics->GetXaxis()->GetFirst() + icut, cutname[icut].Data());
225	fHcutCorrelation->GetXaxis()->SetBinLabel(fHcutCorrelation->GetXaxis()->GetFirst() + icut, cutname[icut].Data());
226	fHcutCorrelation->GetYaxis()->SetBinLabel(fHcutCorrelation->GetYaxis()->GetFirst() + icut, cutname[icut].Data());
227	}
228	Char_t hname[256], htitle[256];
229	for(Int_t imode = 0; imode < 2; imode++){
230	snprintf(hname, 256, "fHclusterRatio%s", imode ? "After" : "Before");
231	snprintf(htitle, 256, "TPC cluster Ratio %s cuts;Ratio;Entries", imode ? "after" : "before");
232	fHclusterRatio[imode] = new TH1F(hname, htitle, 20, 0., 1.);
233	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
234	snprintf(hname, 256, "fHnSigma%sTPC%s", AliPID::ParticleName(ispec), imode ? "after" : "before");
235	snprintf(htitle, 256, "TPC sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");
236	fHnSigmaTPC[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);
237	snprintf(hname, 256, "fHnSigma%sTOF%s", AliPID::ParticleName(ispec), imode ? "after" : "before");
238	snprintf(htitle, 256, "TOF sigma for %s %s cuts;sigma;Entries", AliPID::ParticleName(ispec), imode ? "after" : "before");
239	fHnSigmaTOF[ispec][imode] = new TH1F(hname, htitle, 200, -10., 10.);
240	}
241	}
242	}
243
244	//_____________________________________________________________________
245	Bool_t AliESDpidCuts::AcceptTrack(const AliESDtrack track, const AliESDEvent event){
246	//
247	// Check whether the tracks survived the cuts
248	//
249	if(!fPIDresponse){
250	Init();
251	if (!fPIDresponse)
252	AliError("PID Response not available");
253	return 0;
254	}
255	enum{
256	kCutClusterRatioTPC,
257	kCutNsigmaTPC,
258	kCutNsigmaTOF
259	};
260	Long64_t cutRequired=0, cutFullfiled = 0;
261	if(fTOFsigmaCutRequired && event == 0) {
262	AliError("No event pointer. Need event pointer for T0 for TOF cut");
263	return (0);
264	}
265	Double_t clusterRatio = track->GetTPCNclsF() ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t>(track->GetTPCNclsF()) : 1.;
266	if(fCutTPCclusterRatio > 0.){
267	SETBIT(cutRequired, kCutClusterRatioTPC);
268	if(clusterRatio >= fCutTPCclusterRatio)
269	SETBIT(cutFullfiled, kCutClusterRatioTPC);
270	}
271	// check TPC nSigma cut
272	Float_t nsigmaTPC[AliPID::kSPECIES], nsigma; // need all sigmas for QA plotting
273	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
274	nsigmaTPC[ispec] = nsigma = fPIDresponse->NumberOfSigmasTPC(track,(AliPID::EParticleType)ispec);
275	if(!(fTPCsigmaCutRequired & 1 << ispec)) continue;
276	SETBIT(cutRequired, kCutNsigmaTPC); // We found at least one species where the n-Sigma Cut is required
277	if(nsigma >= fCutTPCnSigma[2ispec] && nsigma <= fCutTPCnSigma[2ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTPC); // Fullfiled for at least one species
278	}
279	// check TOF nSigma cut
280	Float_t nsigmaTOF[AliPID::kSPECIES]; // see above
281	Bool_t hasTOFpid = track->GetStatus() & AliESDtrack::kTOFpid; // only apply TOF n-sigma cut when PID Status Bit is set
282	Double_t times[AliPID::kSPECIES];
283	track->GetIntegratedTimes(times);
284	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
285
286	if(hasTOFpid && event) nsigmaTOF[ispec] = nsigma = fPIDresponse->NumberOfSigmasTOF(track,(AliPID::EParticleType)ispec);
287	if(!(fTOFsigmaCutRequired && 1 << ispec)) continue;
288	SETBIT(cutRequired, kCutNsigmaTOF);
289	if(track->GetOuterParam()->P() >= fMinMomentumTOF){
290	if(hasTOFpid && nsigma <= fCutTOFnSigma[2ispec] && nsigma >= fCutTOFnSigma[2ispec+1]) SETBIT(cutFullfiled, kCutNsigmaTOF);
291	}
292	}
293
294	// Fill Histograms before cuts
295	if(HasHistograms()){
296	fHclusterRatio[0]->Fill(clusterRatio);
297	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
298	fHnSigmaTPC[ispec][0]->Fill(nsigmaTPC[ispec]);
299	if(hasTOFpid) fHnSigmaTOF[ispec][0]->Fill(nsigmaTOF[ispec]);
300	}
301	}
302	if(cutRequired != cutFullfiled){
303	// Fill cut statistics
304	if(HasHistograms()){
305	for(Int_t icut = 0; icut < kNcuts; icut++){
306	if(TESTBIT(cutRequired, icut) && !TESTBIT(cutFullfiled, icut)){
307	// cut not fullfiled
308	fHcutStatistics->Fill(icut);
309	for(Int_t jcut = 0; jcut <= icut; jcut++)
310	if(TESTBIT(cutRequired, jcut) && !TESTBIT(cutFullfiled, jcut)) fHcutCorrelation->Fill(jcut, icut);
311	}
312	}
313	}
314	return kFALSE; // At least one cut is not fullfiled
315	}
316
317	// Fill Histograms after cuts
318	if(HasHistograms()){
319	fHclusterRatio[1]->Fill(clusterRatio);
320	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
321	fHnSigmaTPC[ispec][1]->Fill(nsigmaTPC[ispec]);
322	if(hasTOFpid) fHnSigmaTOF[ispec][1]->Fill(nsigmaTOF[ispec]);
323	}
324	}
325
326	return kTRUE;
327	}
328
329	//_____________________________________________________________________
330	void AliESDpidCuts::SaveHistograms(const Char_t * location){
331	//
332	// Save the histograms to a file
333	//
334	if(!HasHistograms()){
335	AliError("Histograms not on - Exiting");
336	return;
337	}
338	if(!location) location = GetName();
339	gDirectory->mkdir(location);
340	gDirectory->cd(location);
341	fHcutStatistics->Write();
342	fHcutCorrelation->Write();
343
344	gDirectory->mkdir("before_cuts");
345	gDirectory->mkdir("after_cuts");
346
347	gDirectory->cd("before_cuts");
348	fHclusterRatio[0]->Write();
349	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
350	fHnSigmaTPC[ispec][0]->Write();
351	fHnSigmaTOF[ispec][0]->Write();
352	}
353
354	gDirectory->cd("../after_cuts");
355	fHclusterRatio[1]->Write();
356	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
357	fHnSigmaTPC[ispec][1]->Write();
358	fHnSigmaTOF[ispec][1]->Write();
359	}
360
361	gDirectory->cd("..");
362	}
363
364	//_____________________________________________________________________
365	void AliESDpidCuts::DrawHistograms(){
366	//
367	// Draw the Histograms
368	//
369	TCanvas *stat = new TCanvas("cutStat", "Cut Statistics", 640, 480);
370	stat->cd();
371	fHcutStatistics->SetStats(kFALSE);
372	fHcutStatistics->Draw();
373	stat->SaveAs(Form("%s_%s.gif", GetName(), stat->GetName()));
374
375	TCanvas *correl = new TCanvas("cutCorrelation", "Cut Correlation", 640, 480);
376	correl->cd();
377	fHcutCorrelation->SetStats(kFALSE);
378	fHcutCorrelation->Draw("colz");
379	correl->SaveAs(Form("%s_%s.gif", GetName(), correl->GetName()));
380
381	TCanvas *cRatio = new TCanvas("ClusterRatioTPC", "TPC cluster Ratio", 640, 480);
382	cRatio->cd();
383	fHclusterRatio[0]->SetLineColor(kRed);
384	fHclusterRatio[0]->SetStats(kFALSE);
385	fHclusterRatio[0]->Draw();
386	fHclusterRatio[1]->SetLineColor(kBlue);
387	fHclusterRatio[1]->SetStats(kFALSE);
388	fHclusterRatio[1]->Draw("same");
389	cRatio->SaveAs(Form("%s_%s.gif", GetName(), cRatio->GetName()));
390
391	TCanvas *cNsigTPC = new TCanvas("NsigmaTPC", "TPC n-sigma", 640, 480);
392	cNsigTPC->Divide(3,2);
393	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
394	cNsigTPC->cd(ispec + 1);
395	fHnSigmaTPC[ispec][0]->SetLineColor(kRed);
396	fHnSigmaTPC[ispec][0]->SetStats(kFALSE);
397	fHnSigmaTPC[ispec][0]->Draw();
398	fHnSigmaTPC[ispec][1]->SetLineColor(kBlue);
399	fHnSigmaTPC[ispec][1]->SetStats(kFALSE);
400	fHnSigmaTPC[ispec][1]->Draw("same");
401	}
402	cNsigTPC->SaveAs(Form("%s_%s.gif", GetName(), cNsigTPC->GetName()));
403
404	TCanvas *cNsigTOF = new TCanvas("NsigmaTOF", "TOF n-sigma", 640, 480);
405	cNsigTOF->Divide(3,2);
406	for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
407	cNsigTOF->cd(ispec + 1);
408	fHnSigmaTOF[ispec][0]->SetLineColor(kRed);
409	fHnSigmaTOF[ispec][0]->SetStats(kFALSE);
410	fHnSigmaTOF[ispec][0]->Draw();
411	fHnSigmaTOF[ispec][1]->SetLineColor(kBlue);
412	fHnSigmaTOF[ispec][1]->SetStats(kFALSE);
413	fHnSigmaTOF[ispec][1]->Draw("same");
414	}
415	cNsigTOF->SaveAs(Form("%s_%s.gif", GetName(), cNsigTOF->GetName()));
416	}
417