From 0e60774c6eb3f98709e5258ec7f4eab1f7588408 Mon Sep 17 00:00:00 2001 From: abercuci Date: Mon, 29 Sep 2008 15:00:15 +0000 Subject: [PATCH] new utility class for TRD PID (Alex W) --- TRD/AliTRDpidUtil.cxx | 136 ++++++++++++++++++++++++++++++++++++++++++ TRD/AliTRDpidUtil.h | 51 ++++++++++++++++ 2 files changed, 187 insertions(+) create mode 100644 TRD/AliTRDpidUtil.cxx create mode 100644 TRD/AliTRDpidUtil.h diff --git a/TRD/AliTRDpidUtil.cxx b/TRD/AliTRDpidUtil.cxx new file mode 100644 index 00000000000..9eae4931833 --- /dev/null +++ b/TRD/AliTRDpidUtil.cxx @@ -0,0 +1,136 @@ +#include "TObject.h" +#include "TMath.h" +#include "TF1.h" +#include "TH1F.h" +#include "TGraph.h" +#include "TGraphErrors.h" + +#include "AliLog.h" +#include "Cal/AliTRDCalPID.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 + // +} + + + +//________________________________________________________________________ +void AliTRDpidUtil::CalculatePionEffi(TH1F* histo1, TH1F* histo2) +// Double_t AliTRDpidUtil::GetError() +{ + // + // Calculates the pion efficiency + // + + 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; + } + 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)); + + + // calculate electron efficiency of each bin + for (abinE = (histo1 -> GetNbinsX())-2; 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())-2; 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 = new TGraph(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))); + + // 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; + } + fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1 -> Derivative(fEleEffi))*(f1 -> Derivative(fEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi))); +} + + +//__________________________________________________________________________ +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; +} + diff --git a/TRD/AliTRDpidUtil.h b/TRD/AliTRDpidUtil.h new file mode 100644 index 00000000000..34a6fa46999 --- /dev/null +++ b/TRD/AliTRDpidUtil.h @@ -0,0 +1,51 @@ +#ifndef ALITRDPIDUTIL_H +#define ALITRDPIDUTIL_H + +////////////////////////////////////////////////////// +// +// Class to calculate PID performance of the TRD +// +// Author : Alex Wilk +// +/////////////////////////////////////////////////////// + +class TH1F; +class AliTRDpidUtil : public TObject { +public: + + enum { + kBins = 10001 + }; + + AliTRDpidUtil(); + virtual ~AliTRDpidUtil(){;} + + void CalculatePionEffi(TH1F* histo1, TH1F* histo2); + + static Float_t ElectronEfficiency() { return fEleEffi;}; + + Double_t GetCalcElectronEfficiency() {return fCalcEleEffi;}; + Double_t GetPionEfficiency() {return fPionEffi;}; + Double_t GetError() {return fError;}; + Double_t GetThreshold() {return fThreshold;}; + + Int_t GetMomentumBin(Double_t p); + + static void SetElectronEfficiency(Float_t eleeffi) {fEleEffi = eleeffi;}; + +private: + AliTRDpidUtil(const AliTRDpidUtil&); // not implemented + AliTRDpidUtil& operator=(const AliTRDpidUtil&); // not implemented + + static Float_t fEleEffi; // electron efficiency + + Double_t fCalcEleEffi; // electron efficiency after calculation + Double_t fPionEffi; // pion efficiency + Double_t fError; // pion efficiency error + Double_t fThreshold; // threshold for calculated electron efficiency + + ClassDef(AliTRDpidUtil, 1) // TRD PID efficiency calculator + +}; + +#endif -- 2.31.1