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