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