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