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 <TClass.h>
  26 #include <TAxis.h>
  27 #include <TFile.h>
  28 #include <TH1.h>
  29 #include <TKey.h>
  30 #include <TObjArray.h>
  31 #include <TString.h>
  32 #include <TROOT.h>
  33 #include <TSystem.h>
  34 #include <TDirectory.h>
  35
  36 #include "AliLog.h"
  37
  38 #include "AliTRDPIDReference.h"
  39 #include "AliTRDPIDResponse.h"
  40
  41 ClassImp(AliTRDPIDResponse)
  42
  43 //____________________________________________________________
  44 AliTRDPIDResponse::AliTRDPIDResponse():
  45   TObject()
  46   ,fkPIDReference(NULL)
  47   ,fGainNormalisationFactor(1.)
  48   ,fPIDmethod(kLQ1D)
  49 {
  50   //
  51   // Default constructor
  52   //
  53 }
  54
  55 //____________________________________________________________
  56 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
  57   TObject(ref)
  58   ,fkPIDReference(ref.fkPIDReference)
  59   ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
  60   ,fPIDmethod(ref.fPIDmethod)
  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   fkPIDReference = ref.fkPIDReference;
  78   fPIDmethod = ref.fPIDmethod;
  79
  80   return *this;
  81 }
  82
  83 //____________________________________________________________
  84 AliTRDPIDResponse::~AliTRDPIDResponse(){
  85   //
  86   // Destructor
  87   //
  88   if(IsOwner()) delete fkPIDReference;
  89 }
  90
  91 //____________________________________________________________
  92 Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
  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   fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
 109   in->Close(); delete in;
 110   owd->cd();
 111   SetBit(kIsOwner, kTRUE);
 112   AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins()));
 113   return kTRUE;
 114 }
 115
 116 //____________________________________________________________
 117 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
 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   //   true if calculation success
 135   //
 136
 137   if(!fkPIDReference){
 138     AliWarning("Missing reference data. PID calculation not possible.");
 139     return kFALSE;
 140   }
 141
 142   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
 143   Double_t prLayer[AliPID::kSPECIES];
 144   Double_t dE[10], s(0.);
 145   for(Int_t il(kNlayer); il--;){
 146     memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
 147     if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
 148
 149     s=0.;
 150     for(Int_t is(AliPID::kSPECIES); is--;){
 151       if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]);
 152       AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
 153       s+=prLayer[is];
 154     }
 155     if(s<1.e-30){
 156       AliDebug(2, Form("Null all species prob layer[%d].", il));
 157       continue;
 158     }
 159     for(Int_t is(AliPID::kSPECIES); is--;){
 160       if(kNorm) prLayer[is] /= s;
 161       prob[is] *= prLayer[is];
 162     }
 163   }
 164   if(!kNorm) return kTRUE;
 165
 166   s=0.;
 167   for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
 168   if(s<1.e-30){
 169     AliDebug(2, "Null total prob.");
 170     return kFALSE;
 171   }
 172   for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
 173   return kTRUE;
 174 }
 175
 176 //____________________________________________________________
 177 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
 178   //
 179   // Get the non-normalized probability for a certain particle species as coming
 180   // from the reference histogram
 181   // Interpolation between momentum bins
 182   //
 183   AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
 184   Float_t pLower, pUpper;
 185   Double_t probLayer = 0.;
 186   TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)),
 187       *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower));
 188   // Do Interpolation exept for underflow and overflow
 189   if(refLower && refUpper){
 190     Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
 191     Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
 192
 193     probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
 194   } else if(refLower){
 195     // underflow
 196     probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
 197   } else if(refUpper){
 198     // overflow
 199     probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
 200   } else {
 201     AliError("No references available");
 202   }
 203   AliDebug(1, Form("Probability %f", probLayer));
 204   return probLayer;
 205 }
 206
 207 //____________________________________________________________
 208 void AliTRDPIDResponse::SetOwner(){
 209   //
 210   // Make Deep Copy of the Reference Histograms
 211   //
 212   if(!fkPIDReference || IsOwner()) return;
 213   const AliTRDPIDReference *tmp = fkPIDReference;
 214   fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
 215   SetBit(kIsOwner, kTRUE);
 216 }
 217
 218 //____________________________________________________________
 219 Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
 220 {
 221         //
 222         // Recalculate dE/dx
 223         //
 224   switch(fPIDmethod){
 225   case kNN: // NN
 226     break;
 227   case kLQ2D: // 2D LQ
 228     break;
 229   case kLQ1D: // 1D LQ
 230     out[0]= 0.;
 231     for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor;
 232     if(out[0] < 1e-6) return kFALSE;
 233     break;
 234   default:
 235     return kFALSE;
 236   }
 237   return kTRUE;
 238 }
 239