new utility class for TRD PID (Alex W)

diff --git a/TRD/AliTRDpidUtil.cxx b/TRD/AliTRDpidUtil.cxx
new file mode 100644 (file)
index 0000000..9eae493
--- /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 (file)
index 0000000..34a6fa4
--- /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 <wilka@uni-muenster.de>
+//
+///////////////////////////////////////////////////////
+
+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