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 | |
45 | ClassImp(AliTRDPIDResponse) |
46 | |
ffb1ee30 |
47 | //____________________________________________________________ |
e0de37e9 |
48 | AliTRDPIDResponse::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 | //____________________________________________________________ |
60 | AliTRDPIDResponse::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 | //____________________________________________________________ |
72 | AliTRDPIDResponse &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 | //____________________________________________________________ |
88 | AliTRDPIDResponse::~AliTRDPIDResponse(){ |
e0de37e9 |
89 | // |
90 | // Destructor |
91 | // |
db0e2c5f |
92 | if(IsOwner()) delete fkPIDResponseObject; |
ffb1ee30 |
93 | } |
94 | |
95 | //____________________________________________________________ |
db0e2c5f |
96 | Bool_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 |
121 | Bool_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 |
190 | Double_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 | //____________________________________________________________ |
249 | void 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 |
260 | Bool_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 | //____________________________________________________________ |
290 | Bool_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 | |