changes from fzhou
[u/mrichter/AliRoot.git] / TRD / AliTRDpidUtil.cxx
index fa2668f..85d5eaa 100644 (file)
@@ -1,3 +1,34 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: AliTRDdigitizer.cxx 44182 2010-10-10 16:23:39Z cblume $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// Helper class for TRD PID efficiency calculation.                       //
+// Calculation of the hadron efficiency dependent on momentum and of      //
+// the errors implemented in function CalculatePionEff. The pion          //
+// efficiency is based on a predefined electron efficiency.               //
+// The default is 90%. To change the, one has to call the function        //
+// SetElectronEfficiency.                                                 //
+// Other Helper functions decide based on 90% electron efficiency         //
+// whether a certain track is accepted as Electron Track.                 //
+// The reference data is stored in the TRD OCDB.                          //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
 #include "TObject.h"
 #include "TObjArray.h"
 #include "TMath.h"
@@ -17,8 +48,7 @@
 
 ClassImp(AliTRDpidUtil)
 
-
-Float_t AliTRDpidUtil::fEleEffi = 0.9;
+Float_t AliTRDpidUtil::fgEleEffi = 0.9;
 
 //________________________________________________________________________
 AliTRDpidUtil::AliTRDpidUtil() 
@@ -33,8 +63,6 @@ AliTRDpidUtil::AliTRDpidUtil()
   //
 }
 
-
-
 //________________________________________________________________________
 Bool_t  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
 // Double_t  AliTRDpidUtil::GetError()
@@ -45,65 +73,62 @@ Bool_t  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
 
   histo1 -> SetLineColor(kRed);
   histo2 -> SetLineColor(kBlue); 
-  if(!histo1 -> GetEntries() || !histo2 -> GetEntries()){
-    AliWarning("Histo with no entries !");
+  if(!histo1->GetEntries() || !histo2 -> GetEntries()){
+    AliError("Probability histos empty !");
+    return kFALSE;
+  }
+  if(histo1->GetNbinsX() != histo2->GetNbinsX()){
+    AliError(Form("Electron probability discretized differently from pions [%d %d] !", histo1->GetNbinsX(), histo2->GetNbinsX()));
     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;
+  // array of the integrated sum in each bin
+  Double_t sumElecE[kBins+2], sumPionsE[kBins+2];  
+  memset(sumElecE, 0, (kBins+2)*sizeof(Double_t));
+  memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t));
+
+  Int_t nbinE(histo1->GetNbinsX()),
+        abinE(nbinE),
+        bbinE(nbinE),
+        cbinE(-1);
+  // calculate electron efficiency for each bin
+  // and also integral distribution
+  for(Bool_t bFirst(kTRUE); abinE--;){
+    sumElecE[abinE] = sumElecE[abinE+1] + histo1->GetBinContent(abinE+1);
+    if((sumElecE[abinE] >= fgEleEffi) && bFirst){
       cbinE = abinE;
-      fCalcEleEffi = (SumElecsE[cbinE]); 
+      fCalcEleEffi = sumElecE[cbinE];
+      bFirst = kFALSE;
     }
   }
-  
-  fThreshold = histo1 -> GetBinCenter(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]);
-    }
+  // and also integral distribution
+  for (;bbinE--;){
+    sumPionsE[bbinE] = sumPionsE[bbinE+1] + histo2->GetBinContent(bbinE+1);
+    if(bbinE == cbinE) fPionEffi = sumPionsE[cbinE];
   }
   
 
   // pion efficiency vs electron efficiency
-  TGraph gEffis(kBins, SumElecsE, SumPionsE);
+  TGraph gEffis(nbinE, sumElecE, 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);
+  TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05);
+  gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.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!");
+    AliError(" 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)));
+  fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(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)));
+  AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
+  AliDebug(2, Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi)));
   return kTRUE;
 }
 
@@ -154,15 +179,18 @@ Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method)
   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);
+  AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
+  if (!cdbThresholds) return kFALSE;
+  TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
+  if (!histos) return kFALSE;
+  TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
+  if (!thresholdHist) return kFALSE;
+  Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
   
   // Do Decision
-  Double_t pid_probs[5];
-  track->GetTRDpid(pid_probs);
-  if(pid_probs[AliPID::kElectron] >= threshold) return kTRUE;
+  Double_t pidProbs[5];
+  track->GetTRDpid(pidProbs);
+  if(pidProbs[AliPID::kElectron] >= threshold) return kTRUE;
   return kFALSE; 
 }
 
@@ -176,16 +204,19 @@ Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMeth
   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);
+  AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
+  if (!cdbThresholds) return kFALSE;
+  TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
+  if (!histos) return kFALSE;
+  TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
+  if (!thresholdHist) return kFALSE;
+  return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
 }
 
 //________________________________________________________________________
 Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
   //
-  // Private Helper function to get the paticle species (ALICE notation) 
+  // Private Helper function to get the paticle species (ALICE notation)
   // from the Pdg code
   //
   Int_t species;
@@ -211,3 +242,14 @@ Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
   return species;
 }
 
+//________________________________________________________________________
+Int_t AliTRDpidUtil::Mass2Pid(Float_t m){
+  //
+  // Private Helper function to get the paticle species (ALICE notation)
+  // from the Pdg mass
+  //
+
+  for(Int_t is(0); is<AliPID::kSPECIES; is++) if(TMath::Abs(m-AliPID::ParticleMass(is))<1.e-4) return is;
+  return -1;
+}
+