coverity fixed
[u/mrichter/AliRoot.git] / STEER / STEERBase / 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//
ffb1ee30 25#include <TAxis.h>
ea235c90 26#include <TClass.h>
27#include <TDirectory.h>
ffb1ee30 28#include <TFile.h>
29#include <TH1.h>
ce5d6908 30#include <TKey.h>
ea235c90 31#include <TMath.h>
ffb1ee30 32#include <TObjArray.h>
ce5d6908 33#include <TROOT.h>
ea235c90 34#include <TString.h>
ffb1ee30 35#include <TSystem.h>
ea235c90 36#include <TVectorT.h>
ffb1ee30 37
38#include "AliLog.h"
39
51a0ce25 40#include "AliTRDPIDReference.h"
ffb1ee30 41#include "AliTRDPIDResponse.h"
42
43ClassImp(AliTRDPIDResponse)
44
ffb1ee30 45//____________________________________________________________
e0de37e9 46AliTRDPIDResponse::AliTRDPIDResponse():
47 TObject()
51a0ce25 48 ,fkPIDReference(NULL)
ea235c90 49 ,fkPIDParams(NULL)
e0de37e9 50 ,fGainNormalisationFactor(1.)
51 ,fPIDmethod(kLQ1D)
ffb1ee30 52{
e0de37e9 53 //
54 // Default constructor
55 //
ffb1ee30 56}
57
58//____________________________________________________________
59AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref):
e0de37e9 60 TObject(ref)
51a0ce25 61 ,fkPIDReference(ref.fkPIDReference)
ea235c90 62 ,fkPIDParams(ref.fkPIDParams)
e0de37e9 63 ,fGainNormalisationFactor(ref.fGainNormalisationFactor)
64 ,fPIDmethod(ref.fPIDmethod)
ffb1ee30 65{
e0de37e9 66 //
67 // Copy constructor
e0de37e9 68 //
ffb1ee30 69}
70
71//____________________________________________________________
72AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){
e0de37e9 73 //
74 // Assignment operator
e0de37e9 75 //
76 if(this == &ref) return *this;
77
78 // Make copy
79 TObject::operator=(ref);
51a0ce25 80 fGainNormalisationFactor = ref.fGainNormalisationFactor;
81 fkPIDReference = ref.fkPIDReference;
ea235c90 82 fkPIDParams = ref.fkPIDParams;
e0de37e9 83 fPIDmethod = ref.fPIDmethod;
84
85 return *this;
ffb1ee30 86}
87
88//____________________________________________________________
89AliTRDPIDResponse::~AliTRDPIDResponse(){
e0de37e9 90 //
91 // Destructor
92 //
51a0ce25 93 if(IsOwner()) delete fkPIDReference;
ffb1ee30 94}
95
96//____________________________________________________________
51a0ce25 97Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){
e0de37e9 98 //
99 // Load References into the toolkit
100 //
101 AliDebug(1, "Loading reference histos from root file");
2ad33025 102 TDirectory *owd = gDirectory;// store old working directory
e0de37e9 103
104 if(!filename)
105 filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT"));
106 TFile *in = TFile::Open(filename);
107 if(!in){
108 AliError("Ref file not available.");
109 return kFALSE;
110 }
111
112 gROOT->cd();
51a0ce25 113 fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone());
114 in->Close(); delete in;
e0de37e9 115 owd->cd();
e0de37e9 116 SetBit(kIsOwner, kTRUE);
51a0ce25 117 AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins()));
e0de37e9 118 return kTRUE;
ffb1ee30 119}
120
121//____________________________________________________________
51a0ce25 122Bool_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const
ffb1ee30 123{
e0de37e9 124 //
125 // Calculate TRD likelihood values for the track based on dedx and
126 // momentum values. The likelihoods are calculated by query the
127 // reference data depending on the PID method selected
128 //
129 // Input parameter :
130 // n - number of dedx slices/chamber
131 // dedx - array of dedx slices organized layer wise
132 // p - array of momentum measurements organized layer wise
133 //
134 // Return parameters
135 // prob - probabilities allocated by TRD for particle specis
136 // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob.
137 //
138 // Return value
139 // true if calculation success
140 //
ffb1ee30 141
51a0ce25 142 if(!fkPIDReference){
e0de37e9 143 AliWarning("Missing reference data. PID calculation not possible.");
144 return kFALSE;
145 }
ffb1ee30 146
e0de37e9 147 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2;
148 Double_t prLayer[AliPID::kSPECIES];
51a0ce25 149 Double_t dE[10], s(0.);
e0de37e9 150 for(Int_t il(kNlayer); il--;){
151 memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t));
51a0ce25 152 if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue;
ffb1ee30 153
e0de37e9 154 s=0.;
155 for(Int_t is(AliPID::kSPECIES); is--;){
51a0ce25 156 if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]);
e0de37e9 157 AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is]));
158 s+=prLayer[is];
159 }
160 if(s<1.e-30){
161 AliDebug(2, Form("Null all species prob layer[%d].", il));
162 continue;
163 }
164 for(Int_t is(AliPID::kSPECIES); is--;){
165 if(kNorm) prLayer[is] /= s;
166 prob[is] *= prLayer[is];
167 }
168 }
169 if(!kNorm) return kTRUE;
ffb1ee30 170
e0de37e9 171 s=0.;
172 for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is];
173 if(s<1.e-30){
174 AliDebug(2, "Null total prob.");
175 return kFALSE;
176 }
177 for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s;
178 return kTRUE;
ffb1ee30 179}
180
ffb1ee30 181//____________________________________________________________
b439f460 182Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const {
e0de37e9 183 //
184 // Get the non-normalized probability for a certain particle species as coming
185 // from the reference histogram
186 // Interpolation between momentum bins
187 //
188 AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal));
51a0ce25 189 Float_t pLower, pUpper;
190 Double_t probLayer = 0.;
191 TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)),
192 *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower));
e0de37e9 193 // Do Interpolation exept for underflow and overflow
51a0ce25 194 if(refLower && refUpper){
195 Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
196 Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
e0de37e9 197
51a0ce25 198 probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower);
199 } else if(refLower){
200 // underflow
201 probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx));
202 } else if(refUpper){
203 // overflow
204 probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx));
205 } else {
206 AliError("No references available");
e0de37e9 207 }
51a0ce25 208 AliDebug(1, Form("Probability %f", probLayer));
209 return probLayer;
ffb1ee30 210}
211
212//____________________________________________________________
213void AliTRDPIDResponse::SetOwner(){
e0de37e9 214 //
215 // Make Deep Copy of the Reference Histograms
216 //
51a0ce25 217 if(!fkPIDReference || IsOwner()) return;
218 const AliTRDPIDReference *tmp = fkPIDReference;
219 fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone());
e0de37e9 220 SetBit(kIsOwner, kTRUE);
ffb1ee30 221}
222
223//____________________________________________________________
51a0ce25 224Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const
ffb1ee30 225{
51a0ce25 226 //
227 // Recalculate dE/dx
228 //
e0de37e9 229 switch(fPIDmethod){
230 case kNN: // NN
231 break;
232 case kLQ2D: // 2D LQ
233 break;
234 case kLQ1D: // 1D LQ
235 out[0]= 0.;
99e9d5ec 236 for(Int_t islice = 0; islice < nSlice; islice++)
237 if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information
e0de37e9 238 if(out[0] < 1e-6) return kFALSE;
239 break;
240 default:
241 return kFALSE;
242 }
243 return kTRUE;
ffb1ee30 244}
245
ea235c90 246//____________________________________________________________
247Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const {
248 //
249 // Check whether particle is identified as electron assuming a certain electron efficiency level
250 // Only electron and pion hypothesis is taken into account
251 //
252 // Inputs:
253 // Number of tracklets
254 // Likelihood values
255 // Momentum
256 // Electron efficiency level
257 //
258 // If the function fails when the params are not accessible, the function returns true
259 //
260 if(!fkPIDParams){
261 AliError("No PID Param object available");
262 return kTRUE;
263 }
264 Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]);
265 const TVectorD *params = GetParams(nTracklets, level);
266 if(!params){
267 AliError("No Params found for the given configuration");
268 return kTRUE;
269 }
270 Double_t threshold = 1. - (*params)[0] - (*params)[1] * p - (*params)[2] * TMath::Exp(-(*params)[3] * p);
271 if(probEle > TMath::Max(TMath::Min(threshold, 0.99), 0.2)) return kTRUE; // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
272 return kFALSE;
273}
274
275//____________________________________________________________
276const TVectorD* AliTRDPIDResponse::GetParams(Int_t ntracklets, Double_t level) const {
277 //
278 // returns the threshold for a given number of tracklets and a given efficiency level
279 //tby definition the lower of step is given.
280 //
53d016dc 281 if(ntracklets > 6 || ntracklets <=0) return NULL;
ea235c90 282 TObjArray * entry = dynamic_cast<TObjArray *>(fkPIDParams->At(ntracklets - 1));
283 if(!entry) return NULL;
284
285 TObjArray*cut = NULL;
286 TVectorF *effLevel = NULL; const TVectorD *parameters = NULL;
287 Float_t currentLower = 0.;
288 TIter cutIter(entry);
289 while((cut = dynamic_cast<TObjArray *>(cutIter()))){
290 effLevel = static_cast<TVectorF *>(cut->At(0));
53d016dc 291 if((*effLevel)[0] > currentLower && (*effLevel)[0] <= level){
ea235c90 292 // New Lower entry found
293 parameters = static_cast<const TVectorD *>(cut->At(1));
53d016dc 294 currentLower = (*effLevel)[0];
ea235c90 295 }
296 }
99e9d5ec 297 AliDebug(2, Form("Taking params for %d tracklets and %f electron Efficiency\n", ntracklets, currentLower));
ea235c90 298
299 return parameters;
300}