[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 "AliESDtrack.h"
#include "AliLog.h"
#include "AliTPCpidESD.h"
#include "AliTOFpidESD.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)
  , fTPCpid(NULL)
  , fTOFpid(NULL)
  , fCutTPCclusterRatio(0.)
  , fMinMomentumTOF(0.5)
  , fHcutStatistics(NULL)
  , fHcutCorrelation(NULL)
{
  //
  // Default constructor
  //
  
  fTPCpid = new AliTPCpidESD;
  fTOFpid = new AliTOFpidESD;
  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
    fCutTPCnSigma[ispec] = -1;
    fCutTOFnSigma[ispec] = -1;
  }
  //memset(fCutTPCnSigma, -1, sizeof(Float_t) * AliPID::kSPECIES);
  //memset(fCutTOFnSigma, -1, sizeof(Float_t) * AliPID::kSPECIES);

  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)
  , fTPCpid(NULL)
  , fTOFpid(NULL)
  , fCutTPCclusterRatio(ref.fCutTPCclusterRatio)
  , fMinMomentumTOF(ref.fMinMomentumTOF)
  , fHcutStatistics(NULL)
  , fHcutCorrelation(NULL)
{
  //
  // Copy constructor
  //
  fTPCpid = new AliTPCpidESD(*ref.fTPCpid);
  fTOFpid = new AliTOFpidESD(*ref.fTOFpid);
  memcpy(fCutTPCnSigma, ref.fCutTPCnSigma, sizeof(Float_t) * AliPID::kSPECIES);
  memcpy(fCutTOFnSigma, ref.fCutTOFnSigma, sizeof(Float_t) * AliPID::kSPECIES);
  
  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 fTPCpid;
  delete fTOFpid;

  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 *o){
  //
  // Select Track
  // 
  if(TString(o->IsA()->GetName()).CompareTo("AliESDtrack")){
    Char_t errormessage[256];
    sprintf(errormessage, "Provided object not an AliESDtrack: Type %s", o->IsA()->GetName());
    AliError(errormessage);
    return kFALSE;
  }
  return AcceptTrack(const_cast<const AliESDtrack *>(dynamic_cast<AliESDtrack *>(o)));
}

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

  target.fTPCpid = new AliTPCpidESD;
  target.fTOFpid = new AliTOFpidESD;

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