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