]>
Commit | Line | Data |
---|---|---|
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 | ||
43 | ClassImp(AliTRDPIDResponse) | |
44 | ||
ffb1ee30 | 45 | //____________________________________________________________ |
e0de37e9 | 46 | AliTRDPIDResponse::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 | //____________________________________________________________ | |
59 | AliTRDPIDResponse::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 | //____________________________________________________________ | |
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; |
81 | fkPIDReference = ref.fkPIDReference; | |
ea235c90 | 82 | fkPIDParams = ref.fkPIDParams; |
e0de37e9 | 83 | fPIDmethod = ref.fPIDmethod; |
84 | ||
85 | return *this; | |
ffb1ee30 | 86 | } |
87 | ||
88 | //____________________________________________________________ | |
89 | AliTRDPIDResponse::~AliTRDPIDResponse(){ | |
e0de37e9 | 90 | // |
91 | // Destructor | |
92 | // | |
51a0ce25 | 93 | if(IsOwner()) delete fkPIDReference; |
ffb1ee30 | 94 | } |
95 | ||
96 | //____________________________________________________________ | |
51a0ce25 | 97 | Bool_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 | 122 | 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 | 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 | 182 | Double_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 | //____________________________________________________________ | |
213 | void 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 | 224 | Bool_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.; | |
236 | for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor; | |
237 | if(out[0] < 1e-6) return kFALSE; | |
238 | break; | |
239 | default: | |
240 | return kFALSE; | |
241 | } | |
242 | return kTRUE; | |
ffb1ee30 | 243 | } |
244 | ||
ea235c90 | 245 | //____________________________________________________________ |
246 | Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level) const { | |
247 | // | |
248 | // Check whether particle is identified as electron assuming a certain electron efficiency level | |
249 | // Only electron and pion hypothesis is taken into account | |
250 | // | |
251 | // Inputs: | |
252 | // Number of tracklets | |
253 | // Likelihood values | |
254 | // Momentum | |
255 | // Electron efficiency level | |
256 | // | |
257 | // If the function fails when the params are not accessible, the function returns true | |
258 | // | |
259 | if(!fkPIDParams){ | |
260 | AliError("No PID Param object available"); | |
261 | return kTRUE; | |
262 | } | |
263 | Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]); | |
264 | const TVectorD *params = GetParams(nTracklets, level); | |
265 | if(!params){ | |
266 | AliError("No Params found for the given configuration"); | |
267 | return kTRUE; | |
268 | } | |
269 | Double_t threshold = 1. - (*params)[0] - (*params)[1] * p - (*params)[2] * TMath::Exp(-(*params)[3] * p); | |
270 | 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 | |
271 | return kFALSE; | |
272 | } | |
273 | ||
274 | //____________________________________________________________ | |
275 | const TVectorD* AliTRDPIDResponse::GetParams(Int_t ntracklets, Double_t level) const { | |
276 | // | |
277 | // returns the threshold for a given number of tracklets and a given efficiency level | |
278 | //tby definition the lower of step is given. | |
279 | // | |
280 | if(ntracklets > 6 || ntracklets) return NULL; | |
281 | TObjArray * entry = dynamic_cast<TObjArray *>(fkPIDParams->At(ntracklets - 1)); | |
282 | if(!entry) return NULL; | |
283 | ||
284 | TObjArray*cut = NULL; | |
285 | TVectorF *effLevel = NULL; const TVectorD *parameters = NULL; | |
286 | Float_t currentLower = 0.; | |
287 | TIter cutIter(entry); | |
288 | while((cut = dynamic_cast<TObjArray *>(cutIter()))){ | |
289 | effLevel = static_cast<TVectorF *>(cut->At(0)); | |
290 | if(effLevel[0] > currentLower && effLevel[0] <= level){ | |
291 | // New Lower entry found | |
292 | parameters = static_cast<const TVectorD *>(cut->At(1)); | |
293 | } | |
294 | } | |
295 | ||
296 | return parameters; | |
297 | } |