Fixing coverity warnings (Matthieu)
[u/mrichter/AliRoot.git] / STEER / AliTRDPIDResponse.cxx
CommitLineData
ffb1ee30 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>
ce5d6908 29#include <TKey.h>
ffb1ee30 30#include <TObjArray.h>
31#include <TString.h>
ce5d6908 32#include <TROOT.h>
ffb1ee30 33#include <TSystem.h>
2ad33025 34#include <TDirectory.h>
ffb1ee30 35
36#include "AliLog.h"
37
38#include "AliTRDPIDResponse.h"
39
40ClassImp(AliTRDPIDResponse)
41
42const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6};
43
44//____________________________________________________________
45AliTRDPIDResponse::AliTRDPIDResponse(const Char_t * filename):
46 TObject()
47 ,fReferences(NULL)
ce5d6908 48 ,fGainNormalisationFactor(1.)
ffb1ee30 49 ,fPIDmethod(2)
50{
51 //
52 // Default constructor
53 //
ce5d6908 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);
ffb1ee30 58}
59
60//____________________________________________________________
61AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
62 TObject(ref)
63 ,fReferences(ref.fReferences)
ce5d6908 64 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
ffb1ee30 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 //
ce5d6908 72 Int_t size = (AliPID::kSPECIES)*(kNPBins);
338c9979 73 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
ffb1ee30 74}
75
76//____________________________________________________________
77AliTRDPIDResponse &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);
f59c1465 92 } else if(fReferences) delete fReferences;
ffb1ee30 93 fReferences = ref.fReferences;
ce5d6908 94 Int_t size = (AliPID::kSPECIES)*(kNPBins);
338c9979 95 memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size);
ffb1ee30 96 fPIDmethod = ref.fPIDmethod;
97
98 return *this;
99}
100
101//____________________________________________________________
102AliTRDPIDResponse::~AliTRDPIDResponse(){
103 //
104 // Destructor
105 //
f59c1465 106 if(IsOwner()){
107 // Destroy histos
108 if(fReferences) fReferences->Delete();
ffb1ee30 109 }
f59c1465 110 if(fReferences) delete fReferences;
ffb1ee30 111}
112
113//____________________________________________________________
114Bool_t AliTRDPIDResponse::Load(const Char_t * filename){
115 //
116 // Load References into the toolkit
117 //
2ad33025 118
119 TDirectory *owd = gDirectory;// store old working directory
120
ffb1ee30 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
ce5d6908 129 gROOT->cd();
ffb1ee30 130 fReferences = new TObjArray();
131 fReferences->SetOwner();
132 TIter iter(in->GetListOfKeys());
ce5d6908 133 TKey *histkey = NULL;
ffb1ee30 134 TObject *tmp = NULL;
135 Int_t arrayPos = 0, pbin = 0;
136 AliPID::EParticleType species;
137 TString histname;
ce5d6908 138 while((histkey = dynamic_cast<TKey *>(iter()))){
139 tmp = histkey->ReadObj();
ffb1ee30 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 }
ce5d6908 148 pbin = histname.Atoi() - 1;
1bbef4ba 149 AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos));
ffb1ee30 150 fMapRefHists[species][pbin] = arrayPos;
ce5d6908 151 fReferences->AddAt(new TH1F(*dynamic_cast<TH1F *>(tmp)), arrayPos);
ffb1ee30 152 arrayPos++;
153 }
154 in->Close();
2ad33025 155 owd->cd();
ffb1ee30 156 AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos));
f59c1465 157 SetBit(kIsOwner, kTRUE);
ffb1ee30 158 delete in;
159 return kTRUE;
160}
161
162//____________________________________________________________
b439f460 163Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
ffb1ee30 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];
ffb1ee30 190 Double_t DE[10], s(0.);
191 for(Int_t il(kNlayer); il--;){
1bbef4ba 192 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
ffb1ee30 193 if(!CookdEdx(&dedx[il*n], &DE[0])) continue;
194
195 s=0.;
196 for(Int_t is(AliPID::kSPECIES); is--;){
1bbef4ba 197 if((DE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]);
ce5d6908 198 AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
ffb1ee30 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//____________________________________________________________
b439f460 224Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
ffb1ee30 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 //
ce5d6908 230 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
1bbef4ba 231 Int_t pbin = GetLowerMomentumBin(plocal);
232 AliDebug(1, Form("Bin %d", pbin));
ffb1ee30 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;
1bbef4ba 237 if(fMapRefHists[species][pbin] >= 0)
ce5d6908 238 refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin]));
1bbef4ba 239 if(fMapRefHists[species][pbin+1] >= 0)
ce5d6908 240 refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1]));
1bbef4ba 241 AliDebug(1, Form("Reference Histos (Upper/Lower): [%p|%p]", refHistUpper, refHistLower));
ffb1ee30 242
338c9979 243 if (refHistLower && refHistUpper ) {
244 Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx));
245 Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx));
ffb1ee30 246
9590970b 247 pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * (plocal - fgkPBins[pbin]);
338c9979 248 }
ffb1ee30 249 }
250 else{
251 TH1 *refHist = NULL;
252 if(pbin < 0){
253 // underflow
ce5d6908 254 if(fMapRefHists[species][0] >= 0){
255 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0]));
256 }
ffb1ee30 257 } else {
258 // overflow
ce5d6908 259 if(fMapRefHists[species][kNPBins-1] > 0)
260 refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1]));
ffb1ee30 261 }
338c9979 262 if (refHist)
263 pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx));
ffb1ee30 264 }
ce5d6908 265 AliDebug(1, Form("Probability %f", pLayer));
ffb1ee30 266 return pLayer;
267}
268
269//____________________________________________________________
b439f460 270Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const {
ffb1ee30 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//____________________________________________________________
286void 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);
f59c1465 294 delete fReferences;
ffb1ee30 295 fReferences = ctmp;
296 SetBit(kIsOwner, kTRUE);
297}
298
299//____________________________________________________________
b439f460 300Bool_t AliTRDPIDResponse::CookdEdx(Double_t *in, Double_t *out) const
ffb1ee30 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.;
ce5d6908 309 for(Int_t islice = 0; islice < 8; islice++) out[0] += in[islice] * fGainNormalisationFactor;
ffb1ee30 310 if(out[0] < 1e-6) return kFALSE;
311 break;
312 default:
313 return kFALSE;
314 }
315 return kTRUE;
316}
317