STEER/STEERBase/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
  37 #include "AliLog.h"
  38
  39 #include "AliTRDPIDResponseObject.h"
  40 //#include "AliTRDPIDParams.h"
  41 //#include "AliTRDPIDReference.h"
  42 #include "AliTRDPIDResponse.h"
  43 #include "AliTRDTKDInterpolator.h"
  44
  45 ClassImp(AliTRDPIDResponse)
  46
  47 //____________________________________________________________
  48 AliTRDPIDResponse::AliTRDPIDResponse():
  49   TObject()
  50   ,fkPIDResponseObject(NULL)
  51   ,fGainNormalisationFactor(1.)
  52   ,fPIDmethod(kLQ1D)
  53 {
  54   //
  55   // Default constructor
  56   //
  57 }
  58
  59 //____________________________________________________________
  60 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
  61   TObject(ref)
  62   ,fkPIDResponseObject(NULL)
  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   fkPIDResponseObject = ref.fkPIDResponseObject;
  82   fPIDmethod = ref.fPIDmethod;
  83
  84   return *this;
  85 }
  86
  87 //____________________________________________________________
  88 AliTRDPIDResponse::~AliTRDPIDResponse(){
  89   //
  90   // Destructor
  91   //
  92     if(IsOwner()) delete fkPIDResponseObject;
  93 }
  94
  95 //____________________________________________________________
  96 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
  97   //
  98   // Load References into the toolkit
  99   //
 100   AliDebug(1, "Loading reference histos from root file");
 101   TDirectory *owd = gDirectory;// store old working directory
 102
 103   if(!filename)
 104     filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
 105   TFile *in = TFile::Open(filename);
 106   if(!in){
 107     AliError("Ref file not available.");
 108     return kFALSE;
 109   }
 110
 111   gROOT->cd();
 112   fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone());
 113   in->Close(); delete in;
 114   owd->cd();
 115   SetBit(kIsOwner, kTRUE);
 116   AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins()));
 117   return kTRUE;
 118 }
 119
 120 //____________________________________________________________
 121 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
 122 {
 123   //
 124   // Calculate TRD likelihood values for the track based on dedx and
 125   // momentum values. The likelihoods are calculated by query the
 126   // reference data depending on the PID method selected
 127   //
 128   // Input parameter :
 129   //   n - number of dedx slices/chamber
 130   //   dedx - array of dedx slices organized layer wise
 131   //   p - array of momentum measurements organized layer wise
 132   //
 133   // Return parameters
 134   //   prob - probabilities allocated by TRD for particle specis
 135   //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
 136   //
 137   // Return value
 138   //   true if calculation success
 139   //
 140
 141     if(!fkPIDResponseObject){
 142         AliWarning("Missing reference data. PID calculation not possible.");
 143         return kFALSE;
 144     }
 145
 146     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
 147     Double_t prLayer[AliPID::kSPECIES];
 148     Double_t dE[10], s(0.);
 149     for(Int_t il(kNlayer); il--;){
 150         memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
 151         if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
 152
 153         s=0.;
 154         Bool_t filled=kTRUE;
 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             if(prLayer[is]<1.e-30){
 159                 AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il));
 160                 filled=kFALSE;
 161                 break;
 162             }
 163             s+=prLayer[is];
 164         }
 165         if(!filled){
 166             continue;
 167         }
 168         if(s<1.e-30){
 169             AliDebug(2, Form("Null all species prob layer[%d].", il));
 170             continue;
 171         }
 172         for(Int_t is(AliPID::kSPECIES); is--;){
 173             if(kNorm) prLayer[is] /= s;
 174             prob[is] *= prLayer[is];
 175         }
 176     }
 177     if(!kNorm) return kTRUE;
 178
 179     s=0.;
 180     for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
 181     if(s<1.e-30){
 182         AliDebug(2, "Null total prob.");
 183         return kFALSE;
 184     }
 185     for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
 186     return kTRUE;
 187 }
 188
 189 //____________________________________________________________
 190 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx) const {
 191   //
 192   // Get the non-normalized probability for a certain particle species as coming
 193   // from the reference histogram
 194   // Interpolation between momentum bins
 195   //
 196   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
 197
 198   Double_t probLayer = 0.;
 199   Float_t pLower, pUpper;
 200
 201   switch(fPIDmethod){
 202   case kNN: // NN
 203       break;
 204   case kLQ2D: // 2D LQ
 205       {
 206           Double_t error;
 207           Double_t point[kNslicesLQ2D];
 208           for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];}
 209
 210           AliTRDTKDInterpolator *refLower = dynamic_cast<AliTRDTKDInterpolator*>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D));
 211
 212           if(refLower){
 213               refLower->Eval(point,probLayer,error);
 214           }
 215           else {
 216               AliError("No references available");
 217           }
 218       }
 219       break;
 220   case kLQ1D: // 1D LQ
 221       {
 222           TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)),
 223               *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D));
 224           // Do Interpolation exept for underflow and overflow
 225           if(refLower && refUpper){
 226               Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
 227               Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
 228
 229               probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
 230           } else if(refLower){
 231               // underflow
 232               probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0]));
 233           } else if(refUpper){
 234               // overflow
 235               probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0]));
 236           } else {
 237               AliError("No references available");
 238           }
 239           AliDebug(1, Form("Probability %f", probLayer));
 240       }
 241       break;
 242   default:
 243       break;
 244   }
 245   return probLayer;
 246 }
 247
 248 //____________________________________________________________
 249 void AliTRDPIDResponse::SetOwner(){
 250   //
 251   // Make Deep Copy of the Reference Histograms
 252   //
 253     if(!fkPIDResponseObject || IsOwner()) return;
 254     const AliTRDPIDResponseObject *tmp = fkPIDResponseObject;
 255     fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone());
 256     SetBit(kIsOwner, kTRUE);
 257 }
 258
 259 //____________________________________________________________
 260 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
 261 {
 262         //
 263         // Recalculate dE/dx
 264         //
 265   switch(fPIDmethod){
 266   case kNN: // NN
 267       break;
 268   case kLQ2D: // 2D LQ
 269       out[0]=0;
 270       out[1]=0;
 271       for(Int_t islice = 0; islice < nSlice; islice++){
 272           if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice];
 273           else out[1]+= in[islice];
 274       }
 275       break;
 276   case kLQ1D: // 1D LQ
 277       out[0]= 0.;
 278       for(Int_t islice = 0; islice < nSlice; islice++)
 279           if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor;   // Protect against negative values for slices having no dE/dx information
 280       if(out[0] < 1e-6) return kFALSE;
 281       break;
 282
 283   default:
 284       return kFALSE;
 285   }
 286   return kTRUE;
 287 }
 288
 289 //____________________________________________________________
 290 Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const {
 291     //
 292     // Check whether particle is identified as electron assuming a certain electron efficiency level
 293     // Only electron and pion hypothesis is taken into account
 294     //
 295     // Inputs:
 296     //         Number of tracklets
 297     //         Likelihood values
 298     //         Momentum
 299     //         Electron efficiency level
 300     //
 301     // If the function fails when the params are not accessible, the function returns true
 302     //
 303   if(!fkPIDResponseObject){
 304     AliError("No PID Param object available");
 305     return kTRUE;
 306   }
 307   Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
 308   Double_t params[4];
 309   if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,fPIDmethod)){
 310     AliError("No Params found for the given configuration");
 311     return kTRUE;
 312   }
 313   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
 314   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
 315   return kFALSE;
 316 }
 317