TRD/AliTRDpidUtil.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliTRDdigitizer.cxx 44182 2010-10-10 16:23:39Z cblume $ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 // Helper class for TRD PID efficiency calculation.                       //
  21 // Calculation of the hadron efficiency dependent on momentum and of      //
  22 // the errors implemented in function CalculatePionEff. The pion          //
  23 // efficiency is based on a predefined electron efficiency.               //
  24 // The default is 90%. To change the, one has to call the function        //
  25 // SetElectronEfficiency.                                                 //
  26 // Other Helper functions decide based on 90% electron efficiency         //
  27 // whether a certain track is accepted as Electron Track.                 //
  28 // The reference data is stored in the TRD OCDB.                          //
  29 //                                                                        //
  30 ////////////////////////////////////////////////////////////////////////////
  31
  32 #include "TObject.h"
  33 #include "TObjArray.h"
  34 #include "TMath.h"
  35 #include "TF1.h"
  36 #include "TH1F.h"
  37 #include "TGraph.h"
  38 #include "TGraphErrors.h"
  39 #include "TPDGCode.h"
  40
  41 #include "AliLog.h"
  42 #include "Cal/AliTRDCalPID.h"
  43 #include "AliCDBManager.h"
  44 #include "AliCDBEntry.h"
  45 #include "AliESDtrack.h"
  46 #include "AliPID.h"
  47 #include "AliTRDpidUtil.h"
  48
  49 ClassImp(AliTRDpidUtil)
  50
  51 Float_t AliTRDpidUtil::fgEleEffi = 0.9;
  52
  53 //________________________________________________________________________
  54 AliTRDpidUtil::AliTRDpidUtil()
  55   : TObject()
  56   ,fCalcEleEffi(0.)
  57   ,fPionEffi(-1.)
  58   ,fError(-1.)
  59   ,fThreshold(-1.)
  60 {
  61   //
  62   // Default constructor
  63   //
  64 }
  65
  66 //________________________________________________________________________
  67 Bool_t  AliTRDpidUtil::CalculatePionEffi(TH1* histo1, TH1* histo2)
  68 // Double_t  AliTRDpidUtil::GetError()
  69 {
  70   //
  71   // Calculates the pion efficiency
  72   //
  73
  74   histo1 -> SetLineColor(kRed);
  75   histo2 -> SetLineColor(kBlue);
  76   if(!histo1->GetEntries() || !histo2 -> GetEntries()){
  77     AliError("Probability histos empty !");
  78     return kFALSE;
  79   }
  80   if(histo1->GetNbinsX() != histo2->GetNbinsX()){
  81     AliError(Form("Electron probability discretized differently from pions [%d %d] !", histo1->GetNbinsX(), histo2->GetNbinsX()));
  82     return kFALSE;
  83   }
  84   histo1 -> Scale(1./histo1->GetEntries());
  85   histo2 -> Scale(1./histo2->GetEntries());
  86
  87   // array of the integrated sum in each bin
  88   Double_t sumElecE[kBins+2], sumPionsE[kBins+2];
  89   memset(sumElecE, 0, (kBins+2)*sizeof(Double_t));
  90   memset(sumPionsE, 0, (kBins+2)*sizeof(Double_t));
  91
  92   Int_t nbinE(histo1->GetNbinsX()),
  93         abinE(nbinE),
  94         bbinE(nbinE),
  95         cbinE(-1);
  96   // calculate electron efficiency for each bin
  97   // and also integral distribution
  98   for(Bool_t bFirst(kTRUE); abinE--;){
  99     sumElecE[abinE] = sumElecE[abinE+1] + histo1->GetBinContent(abinE+1);
 100     if((sumElecE[abinE] >= fgEleEffi) && bFirst){
 101       cbinE = abinE;
 102       fCalcEleEffi = sumElecE[cbinE];
 103       bFirst = kFALSE;
 104     }
 105   }
 106   fThreshold = histo1->GetBinCenter(cbinE);
 107
 108   // calculate pion efficiency of each bin
 109   // and also integral distribution
 110   for (;bbinE--;){
 111     sumPionsE[bbinE] = sumPionsE[bbinE+1] + histo2->GetBinContent(bbinE+1);
 112     if(bbinE == cbinE) fPionEffi = sumPionsE[cbinE];
 113   }
 114
 115
 116   // pion efficiency vs electron efficiency
 117   TGraph gEffis(nbinE, sumElecE, sumPionsE);
 118
 119   // use fit function to get derivate of the TGraph for error calculation
 120   TF1 f1("f1","[0]*x*x+[1]*x+[2]", fgEleEffi-.05, fgEleEffi+.05);
 121   gEffis.Fit(&f1, "Q", "", fgEleEffi-.05, fgEleEffi+.05);
 122
 123   // return the error of the pion efficiency
 124   if(((1.-fPionEffi) < 0) || ((1.-fCalcEleEffi) < 0)){
 125     AliError(" EleEffi or PioEffi > 1. Error can not be calculated. Please increase statistics or check your simulation!");
 126     return kFALSE;
 127   }
 128   fError = sqrt(((1/histo2 -> GetEntries())*fPionEffi*(1-fPionEffi))+((f1.Derivative(fgEleEffi))*(f1.Derivative(fgEleEffi))*(1/histo1 -> GetEntries())*fCalcEleEffi*(1-fCalcEleEffi)));
 129
 130   AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", fCalcEleEffi, fPionEffi, fError, fThreshold));
 131   AliDebug(2, Form("Derivative at %4.2f : %f\n", fgEleEffi, f1.Derivative(fgEleEffi)));
 132   return kTRUE;
 133 }
 134
 135
 136 //__________________________________________________________________________
 137 Int_t AliTRDpidUtil::GetMomentumBin(Double_t p)
 138 {
 139   //
 140   // returns the momentum bin for a given momentum
 141   //
 142
 143   Int_t pBin1 = -1;                                    // check bin1
 144   Int_t pBin2 = -1;                                    // check bin2
 145
 146   if(p < 0) return -1;                                 // return -1 if momentum < 0
 147   if(p < AliTRDCalPID::GetMomentum(0)) return 0;                                      // smallest momentum bin
 148   if(p >= AliTRDCalPID::GetMomentum(AliTRDCalPID::kNMom-1)) return AliTRDCalPID::kNMom-1; // largest momentum bin
 149
 150
 151   // calculate momentum bin for non extremal momenta
 152   for(Int_t iMomBin = 1; iMomBin < AliTRDCalPID::kNMom; iMomBin++){
 153     if(p < AliTRDCalPID::GetMomentum(iMomBin)){
 154       pBin1 = iMomBin - 1;
 155       pBin2 = iMomBin;
 156     }
 157     else
 158       continue;
 159
 160     if(p - AliTRDCalPID::GetMomentum(pBin1) >= AliTRDCalPID::GetMomentum(pBin2) - p){
 161        return pBin2;
 162     }
 163     else{
 164       return pBin1;
 165     }
 166   }
 167
 168   return -1;
 169 }
 170
 171
 172 //__________________________________________________________________________
 173 Bool_t AliTRDpidUtil::IsElectron(const AliESDtrack *track, ETRDPIDMethod method){
 174   //
 175   // Do PID decision for the TRD based on 90% Electron efficiency threshold
 176   //
 177   // Author: Markus Fasel (M.Fasel@gsi.de)
 178   //
 179   if(method == kESD) method = kNN;
 180   TString histname[2] = {"fHistThreshLQ", "fHistThreshNN"};
 181   AliCDBManager *cdb = AliCDBManager::Instance();
 182   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
 183   if (!cdbThresholds) return kFALSE;
 184   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
 185   if (!histos) return kFALSE;
 186   TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
 187   if (!thresholdHist) return kFALSE;
 188   Double_t threshold = thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
 189
 190   // Do Decision
 191   Double_t pidProbs[5];
 192   track->GetTRDpid(pidProbs);
 193   if(pidProbs[AliPID::kElectron] >= threshold) return kTRUE;
 194   return kFALSE;
 195 }
 196
 197 //__________________________________________________________________________
 198 Double_t AliTRDpidUtil::GetSystematicError(const AliESDtrack *track, ETRDPIDMethod method){
 199   //
 200   // Returns the pion efficiency at 90% electron efficiency from the OCDB
 201   //
 202   // Author: Markus Fasel (M.Fasel@gsi.de)
 203   //
 204   if(method == kESD) method = kNN;
 205   TString histname[2] = {"fHistPionEffLQ", "fHistPionEffNN"};
 206   AliCDBManager *cdb = AliCDBManager::Instance();
 207   AliCDBEntry *cdbThresholds = cdb->Get("TRD/Calib/PIDThresholds");
 208   if (!cdbThresholds) return kFALSE;
 209   TObjArray *histos = dynamic_cast<TObjArray *>(cdbThresholds->GetObject());
 210   if (!histos) return kFALSE;
 211   TH1 *thresholdHist = dynamic_cast<TH1F *>(histos->FindObject(histname[method].Data()));
 212   if (!thresholdHist) return kFALSE;
 213   return thresholdHist->GetBinContent(GetMomentumBin(track->P()) + 1);
 214 }
 215
 216 //________________________________________________________________________
 217 Int_t AliTRDpidUtil::Pdg2Pid(Int_t pdg){
 218   //
 219   // Private Helper function to get the paticle species (ALICE notation)
 220   // from the Pdg code
 221   //
 222   Int_t species;
 223   switch(TMath::Abs(pdg)){
 224   case kElectron:
 225     species = AliPID::kElectron;
 226     break;
 227   case kMuonMinus:
 228     species = AliPID::kMuon;
 229     break;
 230   case kPiPlus:
 231     species = AliPID::kPion;
 232     break;
 233   case kKPlus:
 234     species = AliPID::kKaon;
 235     break;
 236   case kProton:
 237     species = AliPID::kProton;
 238     break;
 239   default:
 240     species = -1;
 241   }
 242   return species;
 243 }
 244
 245 //________________________________________________________________________
 246 Int_t AliTRDpidUtil::Mass2Pid(Float_t m){
 247   //
 248   // Private Helper function to get the paticle species (ALICE notation)
 249   // from the Pdg mass
 250   //
 251
 252   for(Int_t is(0); is<AliPID::kSPECIES; is++) if(TMath::Abs(m-AliPID::ParticleMass(is))<1.e-4) return is;
 253   return -1;
 254 }
 255