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