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