[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx

#include "TObject.h"
#include "TObjArray.h"
#include "TMath.h"
#include "TF1.h"
#include "TH1F.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TPDGCode.h"

#include "AliLog.h"
#include "Cal/AliTRDCalPID.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliESDtrack.h"
#include "AliPID.h"
#include "AliTRDpidUtil.h"

ClassImp(AliTRDpidUtil)


Float_t AliTRDpidUtil::fEleEffi = 0.9;

//________________________________________________________________________
AliTRDpidUtil::AliTRDpidUtil() 
  : TObject()
  ,fCalcEleEffi(0.)
  ,fPionEffi(-1.)
  ,fError(-1.)
  ,fThreshold(-1.)
{
  //
  // Default constructor
  //
}


//________________________________________________________________________
Bool_t  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
// Double_t  AliTRDpidUtil::GetError()
{
  //
  // Calculates the pion efficiency
  //

  histo1 -> SetLineColor(kRed);
  histo2 -> SetLineColor(kBlue); 
  if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
    AliWarning("Histo with no entries !");
    return kFALSE;
  }
  histo1 -> Scale(1./histo1->GetEntries());
  histo2 -> Scale(1./histo2->GetEntries());

  Int_t abinE, bbinE, cbinE = -1;                    
  Double_t aBinSumE, bBinSumE;                  // content of a single bin
  Bool_t bFirst = 1;                            // checks if threshold is crossed for the first time
  Double_t SumElecsE[kBins+2], SumPionsE[kBins+2];  // array of the integrated sum in each bin
  memset(SumElecsE, 0, (kBins+2)*sizeof(Double_t));
  memset(SumPionsE, 0, (kBins+2)*sizeof(Double_t));


  // calculate electron efficiency of each bin
  for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){  
   aBinSumE = 0;
    aBinSumE = histo1 -> GetBinContent(abinE);
    
    SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;

    if((SumElecsE[abinE] >= fEleEffi) && (bFirst == 1)){
      bFirst = 0;
      cbinE = abinE;
      fCalcEleEffi = (SumElecsE[cbinE]); 
    }
  }
  
  fThreshold = histo1 -> GetBinCenter(cbinE);

  // calculate pion efficiency of each bin
  for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){	
    bBinSumE = 0;
    bBinSumE = histo2 -> GetBinContent(bbinE);

    SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
    if(bbinE == cbinE){
      fPionEffi = (SumPionsE[cbinE]);
    }
  }
  

  // pion efficiency vs electron efficiency
  TGraph gEffis(kBins, SumElecsE, SumPionsE);

  // use fit function to get derivate of the TGraph for error calculation
  TF1 f1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
  gEffis.Fit(&f1, "Q", "", fEleEffi-.05, fEleEffi+.05);
  
  // return the error of the pion efficiency
  if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
    AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
    return kFALSE;
  }
  fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fEleEffi))*(f1.Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));

//   AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
//   AliInfo(Form("Derivative at %4.2f : %f\n", fEleEffi, f1.Derivative(fEleEffi)));
  return kTRUE;
}


//__________________________________________________________________________
Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
{
  //
  // returns the momentum bin for a given momentum
  //

  Int_t pBin1 = -1;                                    // check bin1
  Int_t pBin2 = -1;                                    // check bin2

  if(p < 0) return -1;                                 // return -1 if momentum < 0
  if(p < AliTRDCalPID::GetMomentum(0)) return 0;                                      // smallest momentum bin
  if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin


  // calculate momentum bin for non extremal momenta
  for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
    if(p < AliTRDCalPID::GetMomentum(iMomBin)){
      pBin1 = iMomBin - 1;
      pBin2 = iMomBin;
    }
    else
      continue;

    if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
       return pBin2;
    }
    else{
      return pBin1;
    }
  }

  return -1;
}


//__________________________________________________________________________
Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
  //
  // Do PID decision for the TRD based on 90% Electron efficiency threshold
  //
  // Author: Markus Fasel (M.Fasel@gsi.de)
  //
  if(method == kESD) method = kNN;
  TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"};
  AliCDBManager *cdb = AliCDBManager::Instance(); 
  AliCDBEntry *cdb_thresholds = cdb->Get("TRD/Calib/PIDThresholds");
  TObjArray *histos = dynamic_cast<TObjArray *>(cdb_thresholds->GetObject());
  TH1 * threshold_hist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
  Double_t threshold = threshold_hist->GetBinContent(GetMomentumBin(track->P()) + 1);
  
  // Do Decision
  Double_t pid_probs[5];
  track->GetTRDpid(pid_probs);
  if(pid_probs[AliPID::kElectron] >= threshold) return kTRUE;
  return kFALSE; 
}

//__________________________________________________________________________
Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
  //
  // Returns the pion efficiency at 90% electron efficiency from the OCDB
  //
  // Author: Markus Fasel (M.Fasel@gsi.de)
  //
  if(method == kESD) method = kNN;
  TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"};
  AliCDBManager *cdb = AliCDBManager::Instance(); 
  AliCDBEntry *cdb_thresholds = cdb->Get("TRD/Calib/PIDThresholds");
  TObjArray *histos = dynamic_cast<TObjArray *>(cdb_thresholds->GetObject());
  TH1 * threshold_hist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
  return threshold_hist->GetBinContent(GetMomentumBin(track->P()) + 1);
}

//________________________________________________________________________
Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
  //
  // Private Helper function to get the paticle species (ALICE notation) 
  // from the Pdg code
  //
  Int_t species;
  switch(TMath::Abs(pdg)){
  case kElectron:
    species = AliPID::kElectron;
    break;
  case kMuonMinus:
    species = AliPID::kMuon;
    break;
  case kPiPlus:
    species = AliPID::kPion;
    break;
  case kKPlus:
    species = AliPID::kKaon;
    break;
  case kProton:
    species = AliPID::kProton;
    break;
  default:
    species = -1;
  }
  return species;
}
Commit	Line	Data
0e60774c	1	#include "TObject.h"
0affc9e5	2	#include "TObjArray.h"
0e60774c	3	#include "TMath.h"
	4	#include "TF1.h"
	5	#include "TH1F.h"
	6	#include "TGraph.h"
	7	#include "TGraphErrors.h"
0d83b3a5	8	#include "TPDGCode.h"
0e60774c	9
	10	#include "AliLog.h"
	11	#include "Cal/AliTRDCalPID.h"
0affc9e5	12	#include "AliCDBManager.h"
	13	#include "AliCDBEntry.h"
	14	#include "AliESDtrack.h"
	15	#include "AliPID.h"
0e60774c	16	#include "AliTRDpidUtil.h"
	17
	18	ClassImp(AliTRDpidUtil)
	19
	20
	21	Float_t AliTRDpidUtil::fEleEffi = 0.9;
	22
	23	//________________________________________________________________________
	24	AliTRDpidUtil::AliTRDpidUtil()
	25	: TObject()
	26	,fCalcEleEffi(0.)
	27	,fPionEffi(-1.)
	28	,fError(-1.)
	29	,fThreshold(-1.)
	30	{
	31	//
	32	// Default constructor
	33	//
	34	}
	35
	36
	37
	38	//________________________________________________________________________
c46a7947	39	Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
0e60774c	40	// Double_t AliTRDpidUtil::GetError()
	41	{
	42	//
	43	// Calculates the pion efficiency
	44	//
	45
	46	histo1 -> SetLineColor(kRed);
	47	histo2 -> SetLineColor(kBlue);
c46a7947	48	if(!histo1 -> GetEntries() \|\| !histo2 -> GetEntries()){
	49	AliWarning("Histo with no entries !");
	50	return kFALSE;
0e60774c	51	}
	52	histo1 -> Scale(1./histo1->GetEntries());
	53	histo2 -> Scale(1./histo2->GetEntries());
	54
	55	Int_t abinE, bbinE, cbinE = -1;
	56	Double_t aBinSumE, bBinSumE; // content of a single bin
	57	Bool_t bFirst = 1; // checks if threshold is crossed for the first time
422a2dc0	58	Double_t SumElecsE[kBins+2], SumPionsE[kBins+2]; // array of the integrated sum in each bin
	59	memset(SumElecsE, 0, (kBins+2)*sizeof(Double_t));
	60	memset(SumPionsE, 0, (kBins+2)*sizeof(Double_t));
0e60774c	61
	62
	63	// calculate electron efficiency of each bin
422a2dc0	64	for (abinE = (histo1 -> GetNbinsX()); abinE >= 0; abinE--){
422a2dc0	65	aBinSumE = 0;
0e60774c	66	aBinSumE = histo1 -> GetBinContent(abinE);
	67
	68	SumElecsE[abinE] = SumElecsE[abinE+1] + aBinSumE;
422a2dc0	69
0e60774c	70	if((SumElecsE[abinE] >= fEleEffi) && (bFirst == 1)){
	71	bFirst = 0;
	72	cbinE = abinE;
	73	fCalcEleEffi = (SumElecsE[cbinE]);
	74	}
	75	}
	76
	77	fThreshold = histo1 -> GetBinCenter(cbinE);
	78
	79	// calculate pion efficiency of each bin
422a2dc0	80	for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
0e60774c	81	bBinSumE = 0;
	82	bBinSumE = histo2 -> GetBinContent(bbinE);
	83
	84	SumPionsE[bbinE] = SumPionsE[bbinE+1] + bBinSumE;
	85	if(bbinE == cbinE){
	86	fPionEffi = (SumPionsE[cbinE]);
	87	}
	88	}
	89
	90
	91	// pion efficiency vs electron efficiency
c46a7947	92	TGraph gEffis(kBins, SumElecsE, SumPionsE);
0e60774c	93
0e60774c	94	// use fit function to get derivate of the TGraph for error calculation
c46a7947	95	TF1 f1("f1","[0]xx+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
c46a7947	96	gEffis.Fit(&f1, "Q", "", fEleEffi-.05, fEleEffi+.05);
422a2dc0	97
0e60774c	98	// return the error of the pion efficiency
	99	if(((1.-fPionEffi) < 0) \|\| ((1.-fCalcEleEffi) < 0)){
	100	AliWarning(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
c46a7947	101	return kFALSE;
0e60774c	102	}
c46a7947	103	fError = sqrt(((1/histo2 -> GetEntries())fPionEffi(1-fPionEffi))+((f1.Derivative(fEleEffi))(f1.Derivative(fEleEffi))(1/histo1 -> GetEntries())fCalcEleEffi(1-fCalcEleEffi)));
422a2dc0	104
c46a7947	105	// AliInfo(Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
	106	// AliInfo(Form("Derivative at %4.2f : %f\n", fEleEffi, f1.Derivative(fEleEffi)));
	107	return kTRUE;
0e60774c	108	}
	109
	110
	111	//__________________________________________________________________________
	112	Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
	113	{
	114	//
	115	// returns the momentum bin for a given momentum
	116	//
	117
	118	Int_t pBin1 = -1; // check bin1
	119	Int_t pBin2 = -1; // check bin2
	120
	121	if(p < 0) return -1; // return -1 if momentum < 0
	122	if(p < AliTRDCalPID::GetMomentum(0)) return 0; // smallest momentum bin
	123	if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
	124
	125
	126	// calculate momentum bin for non extremal momenta
	127	for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
	128	if(p < AliTRDCalPID::GetMomentum(iMomBin)){
	129	pBin1 = iMomBin - 1;
	130	pBin2 = iMomBin;
	131	}
	132	else
	133	continue;
	134
	135	if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
	136	return pBin2;
	137	}
	138	else{
	139	return pBin1;
	140	}
	141	}
	142
	143	return -1;
	144	}
	145
0affc9e5	146
0affc9e5	147	//__________________________________________________________________________
0d83b3a5	148	Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
0affc9e5	149	//
	150	// Do PID decision for the TRD based on 90% Electron efficiency threshold
	151	//
	152	// Author: Markus Fasel (M.Fasel@gsi.de)
	153	//
	154	if(method == kESD) method = kNN;
	155	TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"};
	156	AliCDBManager *cdb = AliCDBManager::Instance();
	157	AliCDBEntry *cdb_thresholds = cdb->Get("TRD/Calib/PIDThresholds");
	158	TObjArray histos = dynamic_cast<TObjArray >(cdb_thresholds->GetObject());
	159	TH1 * threshold_hist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
	160	Double_t threshold = threshold_hist->GetBinContent(GetMomentumBin(track->P()) + 1);
	161
	162	// Do Decision
	163	Double_t pid_probs[5];
	164	track->GetTRDpid(pid_probs);
	165	if(pid_probs[AliPID::kElectron] >= threshold) return kTRUE;
	166	return kFALSE;
	167	}
	168
	169	//__________________________________________________________________________
0d83b3a5	170	Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
0affc9e5	171	//
	172	// Returns the pion efficiency at 90% electron efficiency from the OCDB
	173	//
	174	// Author: Markus Fasel (M.Fasel@gsi.de)
	175	//
	176	if(method == kESD) method = kNN;
	177	TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"};
	178	AliCDBManager *cdb = AliCDBManager::Instance();
	179	AliCDBEntry *cdb_thresholds = cdb->Get("TRD/Calib/PIDThresholds");
	180	TObjArray histos = dynamic_cast<TObjArray >(cdb_thresholds->GetObject());
	181	TH1 * threshold_hist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
	182	return threshold_hist->GetBinContent(GetMomentumBin(track->P()) + 1);
	183	}
0d83b3a5	184
	185	//________________________________________________________________________
	186	Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
	187	//
	188	// Private Helper function to get the paticle species (ALICE notation)
	189	// from the Pdg code
	190	//
	191	Int_t species;
	192	switch(TMath::Abs(pdg)){
	193	case kElectron:
	194	species = AliPID::kElectron;
	195	break;
	196	case kMuonMinus:
	197	species = AliPID::kMuon;
	198	break;
	199	case kPiPlus:
	200	species = AliPID::kPion;
	201	break;
	202	case kKPlus:
	203	species = AliPID::kKaon;
	204	break;
	205	case kProton:
	206	species = AliPID::kProton;
	207	break;
	208	default:
	209	species = -1;
	210	}
	211	return species;
	212	}
	213