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 "AliTRDPIDResponse.h"
  39
  40 ClassImp(AliTRDPIDResponse)
  41
  42 const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
  43
  44 //____________________________________________________________
  45 AliTRDPIDResponse::AliTRDPIDResponse(const Char_t * filename):
  46  TObject()
  47  ,fReferences(NULL)
  48  ,fGainNormalisationFactor(1.)
  49  ,fPIDmethod(2)
  50 {
  51  //
  52  // Default constructor
  53  //
  54   Int_t size  = (AliPID::kSPECIES)*(kNPBins);
  55   //memset(fMapRefHists, 0, sizeof(Int_t) * size);
  56   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++)
  57     for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++)
  58       fMapRefHists[ispec][ipbin] = -1;
  59   Load(filename);
  60 }
  61
  62 //____________________________________________________________
  63 AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
  64  TObject(ref)
  65  ,fReferences(ref.fReferences)
  66  ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
  67  ,fPIDmethod(ref.fPIDmethod)
  68 {
  69  //
  70  // Copy constructor
  71  // Flat copy of the reference histos
  72  // For creating a deep copy call SetOwner
  73  //
  74   Int_t size  = (AliPID::kSPECIES)*(kNPBins);
  75  memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
  76 }
  77
  78 //____________________________________________________________
  79 AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
  80  //
  81  // Assignment operator
  82  // Performs a flat copy of the reference histos
  83  //
  84  if(this == &ref) return *this;
  85
  86  // Make copy
  87  TObject::operator=(ref);
  88  if(IsOwner()){
  89    if(fReferences){
  90      fReferences->Delete();
  91      delete fReferences;
  92    }
  93    SetBit(kIsOwner, kFALSE);
  94  } else if(fReferences) delete fReferences;
  95  fReferences = ref.fReferences;
  96  Int_t size  = (AliPID::kSPECIES)*(kNPBins);
  97  memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
  98  fPIDmethod = ref.fPIDmethod;
  99
 100  return *this;
 101 }
 102
 103 //____________________________________________________________
 104 AliTRDPIDResponse::~AliTRDPIDResponse(){
 105  //
 106  // Destructor
 107  //
 108  if(IsOwner()){
 109    // Destroy histos
 110    if(fReferences) fReferences->Delete();
 111  }
 112  if(fReferences) delete fReferences;
 113 }
 114
 115 //____________________________________________________________
 116 Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
 117  //
 118  // Load References into the toolkit
 119  //
 120
 121   TDirectory *owd = gDirectory;// store old working directory
 122
 123  if(!filename)
 124    filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
 125  TFile *in = TFile::Open(filename);
 126  if(!in){
 127    AliError("Ref file not available.");
 128    return kFALSE;
 129  }
 130
 131  gROOT->cd();
 132  fReferences = new TObjArray();
 133  fReferences->SetOwner();
 134  TIter iter(in->GetListOfKeys());
 135  TKey *histkey = NULL;
 136  TObject *tmp = NULL;
 137  Int_t arrayPos = 0, pbin = 0;
 138  AliPID::EParticleType species;
 139  TString histname;
 140  while((histkey = dynamic_cast<TKey *>(iter()))){
 141    tmp = histkey->ReadObj();
 142    histname = tmp->GetName();
 143    if(histname.BeginsWith("fHQel")){ // Electron histogram
 144      histname.ReplaceAll("fHQel_p","");
 145      species = AliPID::kElectron;
 146    } else { // Pion histogram
 147      histname.ReplaceAll("fHQpi_p","");
 148      species = AliPID::kPion;
 149    }
 150    pbin = histname.Atoi() - 1;
 151    fMapRefHists[species][pbin] = arrayPos;
 152    fReferences->AddAt(new TH1F(*dynamic_cast<TH1F *>(tmp)), arrayPos);
 153    arrayPos++;
 154  }
 155  in->Close();
 156  owd->cd();
 157  AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
 158  SetBit(kIsOwner, kTRUE);
 159  delete in;
 160  return kTRUE;
 161 }
 162
 163 //____________________________________________________________
 164 Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
 165 {
 166 //
 167 // Calculate TRD likelihood values for the track based on dedx and
 168 // momentum values. The likelihoods are calculated by query the
 169 // reference data depending on the PID method selected
 170 //
 171 // Input parameter :
 172 //   n - number of dedx slices/chamber
 173 //   dedx - array of dedx slices organized layer wise
 174 //   p - array of momentum measurements organized layer wise
 175 //
 176 // Return parameters
 177 //   prob - probabilities allocated by TRD for particle specis
 178 //   kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
 179 //
 180 // Return value
 181 //   true if calculation success
 182 //
 183
 184  if(!fReferences){
 185    AliWarning("Missing reference data. PID calculation not possible.");
 186    return kFALSE;
 187  }
 188
 189  for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
 190  Double_t prLayer[AliPID::kSPECIES];
 191  memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
 192  Double_t DE[10], s(0.);
 193  for(Int_t il(kNlayer); il--;){
 194    if(!CookdEdx(&dedx[il*n], &DE[0])) continue;
 195
 196    s=0.;
 197    for(Int_t is(AliPID::kSPECIES); is--;){
 198      prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
 199      AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
 200      s+=prLayer[is];
 201    }
 202    if(s<1.e-30){
 203      AliDebug(2, Form("Null all species prob layer[%d].", il));
 204      continue;
 205    }
 206    for(Int_t is(AliPID::kSPECIES); is--;){
 207      if(kNorm) prLayer[is] /= s;
 208      prob[is] *= prLayer[is];
 209    }
 210  }
 211  if(!kNorm) return kTRUE;
 212
 213  s=0.;
 214  for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
 215  if(s<1.e-30){
 216    AliDebug(2, "Null total prob.");
 217    return kFALSE;
 218  }
 219  for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
 220  return kTRUE;
 221 }
 222
 223
 224 //____________________________________________________________
 225 Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
 226  //
 227  // Get the non-normalized probability for a certain particle species as coming
 228  // from the reference histogram
 229  // Interpolation between momentum bins
 230  //
 231  AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
 232  Int_t pbin = GetLowerMomentumBin(plocal);
 233  Double_t pLayer = 0.;
 234  // Do Interpolation exept for underflow and overflow
 235  if(pbin >= 0 && pbin < kNPBins-1){
 236    TH1 *refHistUpper = NULL, *refHistLower = NULL;
 237    if(fMapRefHists[species][pbin] > 0)
 238     refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
 239    if(fMapRefHists[species][pbin+1] > 0)
 240     refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
 241
 242    if (refHistLower && refHistUpper ) {
 243      Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
 244      Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
 245
 246      pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * plocal;
 247    }
 248  }
 249  else{
 250    TH1 *refHist = NULL;
 251    if(pbin < 0){
 252      // underflow
 253      if(fMapRefHists[species][0] >= 0){
 254        refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0]));
 255      }
 256    } else {
 257      // overflow
 258      if(fMapRefHists[species][kNPBins-1] > 0)
 259        refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1]));
 260    }
 261    if (refHist)
 262      pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
 263  }
 264  AliDebug(1, Form("Probability %f", pLayer));
 265  return pLayer;
 266 }
 267
 268 //____________________________________________________________
 269 Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
 270  //
 271  // Get the momentum bin for a given momentum value
 272  //
 273  Int_t bin = -1;
 274  if(p > fgkPBins[kNPBins-1]) return kNPBins-1;
 275  for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){
 276    if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){
 277      bin = ibin;
 278      break;
 279    }
 280  }
 281  return bin;
 282 }
 283
 284 //____________________________________________________________
 285 void AliTRDPIDResponse::SetOwner(){
 286  //
 287  // Make Deep Copy of the Reference Histograms
 288  //
 289  if(!fReferences || IsOwner()) return;
 290  TObjArray *ctmp = new TObjArray();
 291  for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++)
 292    ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien);
 293  delete fReferences;
 294  fReferences = ctmp;
 295  SetBit(kIsOwner, kTRUE);
 296 }
 297
 298 //____________________________________________________________
 299 Bool_t AliTRDPIDResponse::CookdEdx(Double_t *in, Double_t *out) const
 300 {
 301  switch(fPIDmethod){
 302  case 0: // NN
 303    break;
 304  case 1: // 2D LQ
 305    break;
 306  case 2: // 1D LQ
 307    out[0]= 0.;
 308    for(Int_t islice = 0; islice < 8; islice++) out[0] += in[islice] * fGainNormalisationFactor;
 309    if(out[0] < 1e-6) return kFALSE;
 310    break;
 311  default:
 312    return kFALSE;
 313  }
 314  return kTRUE;
 315 }
 316