#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)
//________________________________________________________________________
-void AliTRDpidUtil::CalculatePionEffi(TH1F* histo1, TH1F* histo2)
+Bool_t AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
// Double_t AliTRDpidUtil::GetError()
{
//
histo1 -> SetLineColor(kRed);
histo2 -> SetLineColor(kBlue);
- AliInfo(Form("Histo1[%d] Histo2[%d]", (Int_t)histo1 -> GetEntries(), (Int_t)histo2 -> GetEntries()));
- if(!(histo1 -> GetEntries() > 0 && histo2 -> GetEntries() > 0)){
- AliWarning("Histo has no Entries !");
- return;
+ 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], SumPionsE[kBins]; // array of the integrated sum in each bin
- memset(SumElecsE, 0, kBins*sizeof(Double_t));
- memset(SumPionsE, 0, kBins*sizeof(Double_t));
+ 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())-2; abinE >= 0; abinE--){
- aBinSumE = 0;
+ 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;
fThreshold = histo1 -> GetBinCenter(cbinE);
// calculate pion efficiency of each bin
- for (bbinE = (histo2 -> GetNbinsX())-2; bbinE >= abinE; bbinE--){
+ for (bbinE = (histo2 -> GetNbinsX()); bbinE >= abinE; bbinE--){
bBinSumE = 0;
bBinSumE = histo2 -> GetBinContent(bbinE);
// pion efficiency vs electron efficiency
- TGraph *gEffis = new TGraph(kBins, SumElecsE, SumPionsE);
+ TGraph gEffis(kBins, SumElecsE, SumPionsE);
// use fit function to get derivate of the TGraph for error calculation
- TF1 *f1 = new TF1("f1","[0]*x*x+[1]*x+[2]", fEleEffi-.05, fEleEffi+.05);
- gEffis -> Fit("f1","Q","",fEleEffi-.05, fEleEffi+.05);
- AliInfo(Form("Derivative at %4.2f : %f", fEleEffi, f1 -> Derivative(fEleEffi)));
-
+ 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;
+ return kFALSE;
}
- fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1 -> Derivative(fEleEffi))*(f1 -> Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
+ 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;
}
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;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff