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