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