STEER/AliTRDPIDResponse.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 // PID Response class for the TRD detector
  17 // Based on 1D Likelihood approach
  18 // Calculation of probabilities using Bayesian approach
  19 // Attention: This method is only used to separate electrons from pions
  20 //
  21 // Authors:
  22 //  Markus Fasel <M.Fasel@gsi.de>
  23 //  Anton Andronic <A.Andronic@gsi.de>
  24 //
  25 #include <TAxis.h>
  26 #include <TClass.h>
  27 #include <TDirectory.h>
  28 #include <TFile.h>
  29 #include <TH1.h>
  30 #include <TKey.h>
  31 #include <TMath.h>
  32 #include <TObjArray.h>
  33 #include <TROOT.h>
  34 #include <TString.h>
  35 #include <TSystem.h>
  36 #include <TVectorT.h>
  37
  38 #include "AliLog.h"
  39
  40 #include "AliTRDPIDReference.h"
  41 #include "AliTRDPIDResponse.h"
  42
  43 ClassImp(AliTRDPIDResponse)
  44
  45 //____________________________________________________________
  46 AliTRDPIDResponse::AliTRDPIDResponse():
  47   TObject()
  48   ,fkPIDReference(NULL)
  49   ,fkPIDParams(NULL)
  50   ,fGainNormalisationFactor(1.)
  51   ,fPIDmethod(kLQ1D)
  52 {
  53   //
  54   // Default constructor
  55   //
  56 }
  57
  58 //____________________________________________________________
  59 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
  60   TObject(ref)
  61   ,fkPIDReference(ref.fkPIDReference)
  62   ,fkPIDParams(ref.fkPIDParams)
  63   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
  64   ,fPIDmethod(ref.fPIDmethod)
  65 {
  66   //
  67   // Copy constructor
  68   //
  69 }
  70
  71 //____________________________________________________________
  72 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
  73   //
  74   // Assignment operator
  75   //
  76   if(this == &ref) return *this;
  77
  78   // Make copy
  79   TObject::operator=(ref);
  80   fGainNormalisationFactor = ref.fGainNormalisationFactor;
  81   fkPIDReference = ref.fkPIDReference;
  82   fkPIDParams = ref.fkPIDParams;
  83   fPIDmethod = ref.fPIDmethod;
  84
  85   return *this;
  86 }
  87
  88 //____________________________________________________________
  89 AliTRDPIDResponse::~AliTRDPIDResponse(){
  90   //
  91   // Destructor
  92   //
  93   if(IsOwner()) delete fkPIDReference;
  94 }
  95
  96 //____________________________________________________________
  97 Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
  98   //
  99   // Load References into the toolkit
 100   //
 101   AliDebug(1, "Loading reference histos from root file");
 102   TDirectory *owd = gDirectory;// store old working directory
 103
 104   if(!filename)
 105     filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
 106   TFile *in = TFile::Open(filename);
 107   if(!in){
 108     AliError("Ref file not available.");
 109     return kFALSE;
 110   }
 111
 112   gROOT->cd();
 113   fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
 114   in->Close(); delete in;
 115   owd->cd();
 116   SetBit(kIsOwner, kTRUE);
 117   AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins()));
 118   return kTRUE;
 119 }
 120
 121 //____________________________________________________________
 122 Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
 123 {
 124   //
 125   // Calculate TRD likelihood values for the track based on dedx and
 126   // momentum values. The likelihoods are calculated by query the
 127   // reference data depending on the PID method selected
 128   //
 129   // Input parameter :
 130   //   n - number of dedx slices/chamber
 131   //   dedx - array of dedx slices organized layer wise
 132   //   p - array of momentum measurements organized layer wise
 133   //
 134   // Return parameters
 135   //   prob - probabilities allocated by TRD for particle specis
 136   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
 137   //
 138   // Return value
 139   //   true if calculation success
 140   //
 141
 142   if(!fkPIDReference){
 143     AliWarning("Missing reference data. PID calculation not possible.");
 144     return kFALSE;
 145   }
 146
 147   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
 148   Double_t prLayer[AliPID::kSPECIES];
 149   Double_t dE[10], s(0.);
 150   for(Int_t il(kNlayer); il--;){
 151     memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
 152     if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
 153
 154     s=0.;
 155     for(Int_t is(AliPID::kSPECIES); is--;){
 156       if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]);
 157       AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
 158       s+=prLayer[is];
 159     }
 160     if(s<1.e-30){
 161       AliDebug(2, Form("Null all species prob layer[%d].", il));
 162       continue;
 163     }
 164     for(Int_t is(AliPID::kSPECIES); is--;){
 165       if(kNorm) prLayer[is] /= s;
 166       prob[is] *= prLayer[is];
 167     }
 168   }
 169   if(!kNorm) return kTRUE;
 170
 171   s=0.;
 172   for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
 173   if(s<1.e-30){
 174     AliDebug(2, "Null total prob.");
 175     return kFALSE;
 176   }
 177   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
 178   return kTRUE;
 179 }
 180
 181 //____________________________________________________________
 182 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
 183   //
 184   // Get the non-normalized probability for a certain particle species as coming
 185   // from the reference histogram
 186   // Interpolation between momentum bins
 187   //
 188   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
 189   Float_t pLower, pUpper;
 190   Double_t probLayer = 0.;
 191   TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)),
 192       *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower));
 193   // Do Interpolation exept for underflow and overflow
 194   if(refLower && refUpper){
 195     Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
 196     Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
 197
 198     probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
 199   } else if(refLower){
 200     // underflow
 201     probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
 202   } else if(refUpper){
 203     // overflow
 204     probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
 205   } else {
 206     AliError("No references available");
 207   }
 208   AliDebug(1, Form("Probability %f", probLayer));
 209   return probLayer;
 210 }
 211
 212 //____________________________________________________________
 213 void AliTRDPIDResponse::SetOwner(){
 214   //
 215   // Make Deep Copy of the Reference Histograms
 216   //
 217   if(!fkPIDReference || IsOwner()) return;
 218   const AliTRDPIDReference *tmp = fkPIDReference;
 219   fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
 220   SetBit(kIsOwner, kTRUE);
 221 }
 222
 223 //____________________________________________________________
 224 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
 225 {
 226         //
 227         // Recalculate dE/dx
 228         //
 229   switch(fPIDmethod){
 230   case kNN: // NN
 231     break;
 232   case kLQ2D: // 2D LQ
 233     break;
 234   case kLQ1D: // 1D LQ
 235     out[0]= 0.;
 236     for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor;
 237     if(out[0] < 1e-6) return kFALSE;
 238     break;
 239   default:
 240     return kFALSE;
 241   }
 242   return kTRUE;
 243 }
 244
 245 //____________________________________________________________
 246 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const {
 247   //
 248   // Check whether particle is identified as electron assuming a certain electron efficiency level
 249   // Only electron and pion hypothesis is taken into account
 250   //
 251   // Inputs:
 252   //         Number of tracklets
 253   //         Likelihood values
 254   //         Momentum
 255   //         Electron efficiency level
 256   //
 257   // If the function fails when the params are not accessible, the function returns true
 258   //
 259   if(!fkPIDParams){
 260     AliError("No PID Param object available");
 261     return kTRUE;
 262   }
 263   Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
 264   const TVectorD *params = GetParams(nTracklets, level);
 265   if(!params){
 266     AliError("No Params found for the given configuration");
 267     return kTRUE;
 268   }
 269   Double_t threshold = 1. - (*params)[0] - (*params)[1] * p - (*params)[2] * TMath::Exp(-(*params)[3] * p);
 270   if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
 271   return kFALSE;
 272 }
 273
 274 //____________________________________________________________
 275 const TVectorD* AliTRDPIDResponse::GetParams(Int_t ntracklets, Double_t level) const {
 276   //
 277   // returns the threshold for a given number of tracklets and a given efficiency level
 278   //tby definition the lower of step is given.
 279   //
 280   if(ntracklets > 6 || ntracklets) return NULL;
 281   TObjArray * entry = dynamic_cast<TObjArray *>(fkPIDParams->At(ntracklets - 1));
 282   if(!entry) return NULL;
 283
 284   TObjArray*cut = NULL;
 285   TVectorF *effLevel = NULL; const TVectorD *parameters = NULL;
 286   Float_t currentLower = 0.;
 287   TIter cutIter(entry);
 288   while((cut = dynamic_cast<TObjArray *>(cutIter()))){
 289     effLevel = static_cast<TVectorF *>(cut->At(0));
 290     if(effLevel[0] > currentLower && effLevel[0] <= level){
 291       // New Lower entry found
 292       parameters = static_cast<const TVectorD *>(cut->At(1));
 293     }
 294   }
 295
 296   return parameters;
 297 }