]>
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> |
36 | ||
37 | #include "AliLog.h" | |
38 | ||
db0e2c5f | 39 | #include "AliTRDPIDResponseObject.h" |
ffb1ee30 | 40 | #include "AliTRDPIDResponse.h" |
239fe91c | 41 | //#include "AliTRDTKDInterpolator.h" |
42 | #include "AliTRDNDFast.h" | |
ffb1ee30 | 43 | |
44 | ClassImp(AliTRDPIDResponse) | |
45 | ||
ffb1ee30 | 46 | //____________________________________________________________ |
e0de37e9 | 47 | AliTRDPIDResponse::AliTRDPIDResponse(): |
48 | TObject() | |
db0e2c5f | 49 | ,fkPIDResponseObject(NULL) |
e0de37e9 | 50 | ,fGainNormalisationFactor(1.) |
ffb1ee30 | 51 | { |
e0de37e9 | 52 | // |
53 | // Default constructor | |
54 | // | |
ffb1ee30 | 55 | } |
56 | ||
57 | //____________________________________________________________ | |
58 | AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref): | |
e0de37e9 | 59 | TObject(ref) |
db0e2c5f | 60 | ,fkPIDResponseObject(NULL) |
e0de37e9 | 61 | ,fGainNormalisationFactor(ref.fGainNormalisationFactor) |
ffb1ee30 | 62 | { |
e0de37e9 | 63 | // |
64 | // Copy constructor | |
e0de37e9 | 65 | // |
ffb1ee30 | 66 | } |
67 | ||
68 | //____________________________________________________________ | |
69 | AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){ | |
e0de37e9 | 70 | // |
71 | // Assignment operator | |
e0de37e9 | 72 | // |
73 | if(this == &ref) return *this; | |
74 | ||
75 | // Make copy | |
76 | TObject::operator=(ref); | |
51a0ce25 | 77 | fGainNormalisationFactor = ref.fGainNormalisationFactor; |
db0e2c5f | 78 | fkPIDResponseObject = ref.fkPIDResponseObject; |
e0de37e9 | 79 | |
80 | return *this; | |
ffb1ee30 | 81 | } |
82 | ||
83 | //____________________________________________________________ | |
84 | AliTRDPIDResponse::~AliTRDPIDResponse(){ | |
e0de37e9 | 85 | // |
86 | // Destructor | |
87 | // | |
db0e2c5f | 88 | if(IsOwner()) delete fkPIDResponseObject; |
ffb1ee30 | 89 | } |
90 | ||
91 | //____________________________________________________________ | |
db0e2c5f | 92 | Bool_t AliTRDPIDResponse::Load(const Char_t * filename){ |
e0de37e9 | 93 | // |
94 | // Load References into the toolkit | |
95 | // | |
96 | AliDebug(1, "Loading reference histos from root file"); | |
2ad33025 | 97 | TDirectory *owd = gDirectory;// store old working directory |
e0de37e9 | 98 | |
99 | if(!filename) | |
100 | filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT")); | |
101 | TFile *in = TFile::Open(filename); | |
102 | if(!in){ | |
103 | AliError("Ref file not available."); | |
104 | return kFALSE; | |
105 | } | |
106 | ||
107 | gROOT->cd(); | |
db0e2c5f | 108 | fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(in->Get("TRDPIDResponse")->Clone()); |
51a0ce25 | 109 | in->Close(); delete in; |
e0de37e9 | 110 | owd->cd(); |
e0de37e9 | 111 | SetBit(kIsOwner, kTRUE); |
db0e2c5f | 112 | AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDResponseObject->GetNumberOfMomentumBins())); |
e0de37e9 | 113 | return kTRUE; |
ffb1ee30 | 114 | } |
115 | ||
116 | //____________________________________________________________ | |
239fe91c | 117 | Int_t AliTRDPIDResponse::GetResponse(Int_t n, const Double_t * const dedx, const Float_t * const p, Double_t prob[AliPID::kSPECIES],ETRDPIDMethod PIDmethod,Bool_t kNorm) const |
ffb1ee30 | 118 | { |
e0de37e9 | 119 | // |
120 | // Calculate TRD likelihood values for the track based on dedx and | |
121 | // momentum values. The likelihoods are calculated by query the | |
122 | // reference data depending on the PID method selected | |
123 | // | |
124 | // Input parameter : | |
125 | // n - number of dedx slices/chamber | |
126 | // dedx - array of dedx slices organized layer wise | |
127 | // p - array of momentum measurements organized layer wise | |
128 | // | |
129 | // Return parameters | |
130 | // prob - probabilities allocated by TRD for particle specis | |
131 | // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob. | |
132 | // | |
133 | // Return value | |
239fe91c | 134 | // number of tracklets used for PID, 0 if no PID |
bd58d4b9 | 135 | // |
136 | AliDebug(3,Form(" Response for PID method: %d",PIDmethod)); | |
137 | ||
ffb1ee30 | 138 | |
db0e2c5f | 139 | if(!fkPIDResponseObject){ |
31746593 | 140 | AliDebug(3,"Missing reference data. PID calculation not possible."); |
239fe91c | 141 | return 0; |
db0e2c5f | 142 | } |
ffb1ee30 | 143 | |
db0e2c5f | 144 | for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2; |
145 | Double_t prLayer[AliPID::kSPECIES]; | |
146 | Double_t dE[10], s(0.); | |
239fe91c | 147 | Int_t ntrackletsPID=0; |
db0e2c5f | 148 | for(Int_t il(kNlayer); il--;){ |
149 | memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t)); | |
bd58d4b9 | 150 | if(!CookdEdx(n, &dedx[il*n], &dE[0],PIDmethod)) continue; |
239fe91c | 151 | s=0.; |
db0e2c5f | 152 | Bool_t filled=kTRUE; |
153 | for(Int_t is(AliPID::kSPECIES); is--;){ | |
239fe91c | 154 | //if((PIDmethod==kLQ2D)&&(!(is==0||is==2)))continue; |
bd58d4b9 | 155 | if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], &dE[0],PIDmethod); |
f2762b1c | 156 | AliDebug(3, Form("Probability for Species %d in Layer %d: %e", is, il, prLayer[is])); |
db0e2c5f | 157 | if(prLayer[is]<1.e-30){ |
158 | AliDebug(2, Form("Null for species %d species prob layer[%d].",is,il)); | |
159 | filled=kFALSE; | |
160 | break; | |
161 | } | |
162 | s+=prLayer[is]; | |
163 | } | |
164 | if(!filled){ | |
165 | continue; | |
166 | } | |
167 | if(s<1.e-30){ | |
168 | AliDebug(2, Form("Null all species prob layer[%d].", il)); | |
169 | continue; | |
170 | } | |
171 | for(Int_t is(AliPID::kSPECIES); is--;){ | |
172 | if(kNorm) prLayer[is] /= s; | |
173 | prob[is] *= prLayer[is]; | |
174 | } | |
239fe91c | 175 | ntrackletsPID++; |
e0de37e9 | 176 | } |
239fe91c | 177 | if(!kNorm) return ntrackletsPID; |
db0e2c5f | 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."); |
239fe91c | 183 | return 0; |
e0de37e9 | 184 | } |
db0e2c5f | 185 | for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s; |
239fe91c | 186 | return ntrackletsPID; |
ffb1ee30 | 187 | } |
188 | ||
ffb1ee30 | 189 | //____________________________________________________________ |
bd58d4b9 | 190 | Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t *dEdx,ETRDPIDMethod PIDmethod) 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.; |
239fe91c | 199 | |
db0e2c5f | 200 | Float_t pLower, pUpper; |
201 | ||
239fe91c | 202 | AliTRDNDFast *refUpper = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,PIDmethod)), |
203 | *refLower = dynamic_cast<AliTRDNDFast *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,PIDmethod)); | |
204 | ||
205 | // Do Interpolation exept for underflow and overflow | |
206 | if(refLower && refUpper){ | |
207 | Double_t probLower = refLower->Eval(dEdx); | |
208 | Double_t probUpper = refUpper->Eval(dEdx); | |
209 | ||
210 | probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower); | |
211 | } else if(refLower){ | |
212 | // underflow | |
213 | probLayer = refLower->Eval(dEdx); | |
214 | } else if(refUpper){ | |
215 | // overflow | |
216 | probLayer = refUpper->Eval(dEdx); | |
217 | } else { | |
218 | AliError("No references available"); | |
219 | } | |
220 | AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer)); | |
221 | ||
222 | return probLayer; | |
223 | ||
224 | /* old implementation | |
225 | ||
226 | switch(PIDmethod){ | |
227 | case kNN: // NN | |
db0e2c5f | 228 | break; |
229 | case kLQ2D: // 2D LQ | |
230 | { | |
f2762b1c | 231 | if(species==0||species==2){ // references only for electrons and pions |
2285168d | 232 | Double_t error = 0.; |
f2762b1c | 233 | Double_t point[kNslicesLQ2D]; |
234 | for(Int_t idim=0;idim<kNslicesLQ2D;idim++){point[idim]=dEdx[idim];} | |
db0e2c5f | 235 | |
239fe91c | 236 | AliTRDTKDInterpolator *refUpper = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ2D)), |
237 | *refLower = dynamic_cast<AliTRDTKDInterpolator *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ2D)); | |
238 | // Do Interpolation exept for underflow and overflow | |
239 | if(refLower && refUpper){ | |
240 | Double_t probLower=0,probUpper=0; | |
241 | refLower->Eval(point,probLower,error); | |
242 | refUpper->Eval(point,probUpper,error); | |
243 | probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower); | |
244 | } else if(refLower){ | |
245 | // underflow | |
f2762b1c | 246 | refLower->Eval(point,probLayer,error); |
239fe91c | 247 | } else if(refUpper){ |
248 | // overflow | |
249 | refUpper->Eval(point,probLayer,error); | |
250 | } else { | |
f2762b1c | 251 | AliError("No references available"); |
252 | } | |
253 | AliDebug(2,Form("Eval 2D Q0 %f Q1 %f P %e Err %e",point[0],point[1],probLayer,error)); | |
db0e2c5f | 254 | } |
255 | } | |
256 | break; | |
257 | case kLQ1D: // 1D LQ | |
258 | { | |
259 | TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDResponseObject->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper,kLQ1D)), | |
260 | *refLower = dynamic_cast<TH1 *>(fkPIDResponseObject->GetLowerReference((AliPID::EParticleType)species, plocal, pLower,kLQ1D)); | |
261 | // Do Interpolation exept for underflow and overflow | |
262 | if(refLower && refUpper){ | |
263 | Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0])); | |
264 | Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0])); | |
265 | ||
266 | probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower); | |
267 | } else if(refLower){ | |
268 | // underflow | |
269 | probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx[0])); | |
270 | } else if(refUpper){ | |
271 | // overflow | |
272 | probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx[0])); | |
273 | } else { | |
274 | AliError("No references available"); | |
275 | } | |
f2762b1c | 276 | AliDebug(1, Form("Eval 1D dEdx %f Probability %e", dEdx[0],probLayer)); |
db0e2c5f | 277 | } |
278 | break; | |
279 | default: | |
280 | break; | |
239fe91c | 281 | } |
282 | return probLayer; | |
283 | */ | |
284 | ||
ffb1ee30 | 285 | } |
286 | ||
287 | //____________________________________________________________ | |
288 | void AliTRDPIDResponse::SetOwner(){ | |
e0de37e9 | 289 | // |
290 | // Make Deep Copy of the Reference Histograms | |
291 | // | |
db0e2c5f | 292 | if(!fkPIDResponseObject || IsOwner()) return; |
293 | const AliTRDPIDResponseObject *tmp = fkPIDResponseObject; | |
294 | fkPIDResponseObject = dynamic_cast<const AliTRDPIDResponseObject *>(tmp->Clone()); | |
295 | SetBit(kIsOwner, kTRUE); | |
ffb1ee30 | 296 | } |
297 | ||
298 | //____________________________________________________________ | |
bd58d4b9 | 299 | Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out,ETRDPIDMethod PIDmethod) const |
ffb1ee30 | 300 | { |
f2762b1c | 301 | // |
302 | // Recalculate dE/dx | |
303 | // | |
bd58d4b9 | 304 | switch(PIDmethod){ |
e0de37e9 | 305 | case kNN: // NN |
db0e2c5f | 306 | break; |
307 | case kLQ2D: // 2D LQ | |
308 | out[0]=0; | |
309 | out[1]=0; | |
310 | for(Int_t islice = 0; islice < nSlice; islice++){ | |
239fe91c | 311 | if(in[islice]<=0){out[0]=0;out[1]=0;return kFALSE;} // Require that all slices are filled |
db0e2c5f | 312 | if(islice<fkPIDResponseObject->GetNSlicesQ0())out[0]+= in[islice]; |
313 | else out[1]+= in[islice]; | |
314 | } | |
239fe91c | 315 | // normalize signal to number of slices |
316 | out[0]*=1./Double_t(fkPIDResponseObject->GetNSlicesQ0()); | |
317 | out[1]*=1./Double_t(nSlice-fkPIDResponseObject->GetNSlicesQ0()); | |
f2762b1c | 318 | if(out[0] < 1e-6) return kFALSE; |
319 | AliDebug(3,Form("CookdEdx Q0 %f Q1 %f",out[0],out[1])); | |
db0e2c5f | 320 | break; |
321 | case kLQ1D: // 1D LQ | |
322 | out[0]= 0.; | |
323 | for(Int_t islice = 0; islice < nSlice; islice++) | |
324 | if(in[islice] > 0) out[0] += in[islice] * fGainNormalisationFactor; // Protect against negative values for slices having no dE/dx information | |
239fe91c | 325 | out[0]*=1./Double_t(nSlice); |
db0e2c5f | 326 | if(out[0] < 1e-6) return kFALSE; |
f2762b1c | 327 | AliDebug(3,Form("CookdEdx dEdx %f",out[0])); |
db0e2c5f | 328 | break; |
329 | ||
e0de37e9 | 330 | default: |
db0e2c5f | 331 | return kFALSE; |
e0de37e9 | 332 | } |
333 | return kTRUE; | |
ffb1ee30 | 334 | } |
335 | ||
ea235c90 | 336 | //____________________________________________________________ |
bd58d4b9 | 337 | Bool_t AliTRDPIDResponse::IdentifiedAsElectron(Int_t nTracklets, const Double_t *like, Double_t p, Double_t level,Double_t centrality,ETRDPIDMethod PIDmethod) const { |
db0e2c5f | 338 | // |
339 | // Check whether particle is identified as electron assuming a certain electron efficiency level | |
340 | // Only electron and pion hypothesis is taken into account | |
341 | // | |
342 | // Inputs: | |
343 | // Number of tracklets | |
344 | // Likelihood values | |
345 | // Momentum | |
346 | // Electron efficiency level | |
347 | // | |
348 | // If the function fails when the params are not accessible, the function returns true | |
349 | // | |
350 | if(!fkPIDResponseObject){ | |
31746593 | 351 | AliDebug(3,"No PID Param object available"); |
ea235c90 | 352 | return kTRUE; |
353 | } | |
354 | Double_t probEle = like[AliPID::kElectron]/(like[AliPID::kElectron] + like[AliPID::kPion]); | |
ce487a7f | 355 | Double_t params[4]; |
bd58d4b9 | 356 | if(!fkPIDResponseObject->GetThresholdParameters(nTracklets, level, params,centrality,PIDmethod)){ |
ea235c90 | 357 | AliError("No Params found for the given configuration"); |
358 | return kTRUE; | |
359 | } | |
ce487a7f | 360 | Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p); |
ea235c90 | 361 | 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 |
362 | return kFALSE; | |
363 | } | |
364 | ||
f2762b1c | 365 | //____________________________________________________________ |
366 | Bool_t AliTRDPIDResponse::SetPIDResponseObject(const AliTRDPIDResponseObject * obj){ | |
367 | ||
368 | fkPIDResponseObject = obj; | |
369 | if((AliLog::GetDebugLevel("",IsA()->GetName()))>0)fkPIDResponseObject->Print(""); | |
370 | return kTRUE; | |
371 | } |