[u/mrichter/AliRoot.git] / ANALYSIS / AliESDpidCuts.cxx

#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 "AliESDtrack.h"
#include "AliESDEvent.h"
#include "AliLog.h"
#include "AliESDpid.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)
  , fESDpid(NULL)
  , fTPCsigmaCutRequired(0)
  , fTOFsigmaCutRequired(0)
  , fCutTPCclusterRatio(0.)
  , fMinMomentumTOF(0.5)
  , fHcutStatistics(NULL)
  , fHcutCorrelation(NULL)
{
  //
  // Default constructor
  //
  
  fESDpid = new AliESDpid;
  memset(fCutTPCnSigma, 0, sizeof(Float_t) * 2);
  memset(fCutTOFnSigma, 0, sizeof(Float_t) * 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)
  , fESDpid(NULL)
  , fTPCsigmaCutRequired(ref.fTPCsigmaCutRequired)
  , fTOFsigmaCutRequired(ref.fTOFsigmaCutRequired)
  , fCutTPCclusterRatio(ref.fCutTPCclusterRatio)
  , fMinMomentumTOF(ref.fMinMomentumTOF)
  , fHcutStatistics(NULL)
  , fHcutCorrelation(NULL)
{
  //
  // Copy constructor
  //
  fESDpid = new AliESDpid(*ref.fESDpid);
  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 fESDpid;

  delete fHcutStatistics;
  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];
    }
  }
}

//_____________________________________________________________________
Bool_t AliESDpidCuts::IsSelected(TObject *obj){
  //
  // Select Track
  // 
  AliESDtrack * trk = dynamic_cast<AliESDtrack*>(obj);
  if(!trk){
    AliError("Provided object is not AliESDtrack!");
    return kFALSE;
  }
  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.fESDpid = new AliESDpid(*fESDpid);

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