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 | // |
25 | #include <TClass.h> |
26 | #include <TAxis.h> |
27 | #include <TFile.h> |
28 | #include <TH1.h> |
ce5d6908 |
29 | #include <TKey.h> |
ffb1ee30 |
30 | #include <TObjArray.h> |
31 | #include <TString.h> |
ce5d6908 |
32 | #include <TROOT.h> |
ffb1ee30 |
33 | #include <TSystem.h> |
2ad33025 |
34 | #include <TDirectory.h> |
ffb1ee30 |
35 | |
36 | #include "AliLog.h" |
37 | |
38 | #include "AliTRDPIDResponse.h" |
39 | |
40 | ClassImp(AliTRDPIDResponse) |
41 | |
42 | const Double_t AliTRDPIDResponse::fgkPBins[kNPBins] = {1, 2, 3, 4, 5, 6}; |
43 | |
44 | //____________________________________________________________ |
45 | AliTRDPIDResponse::AliTRDPIDResponse(const Char_t * filename): |
46 | TObject() |
47 | ,fReferences(NULL) |
ce5d6908 |
48 | ,fGainNormalisationFactor(1.) |
ffb1ee30 |
49 | ,fPIDmethod(2) |
50 | { |
51 | // |
52 | // Default constructor |
53 | // |
ce5d6908 |
54 | for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++) |
55 | for(Int_t ipbin = 0; ipbin < kNPBins; ipbin++) |
56 | fMapRefHists[ispec][ipbin] = -1; |
57 | Load(filename); |
ffb1ee30 |
58 | } |
59 | |
60 | //____________________________________________________________ |
61 | AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref): |
62 | TObject(ref) |
63 | ,fReferences(ref.fReferences) |
ce5d6908 |
64 | ,fGainNormalisationFactor(ref.fGainNormalisationFactor) |
ffb1ee30 |
65 | ,fPIDmethod(ref.fPIDmethod) |
66 | { |
67 | // |
68 | // Copy constructor |
69 | // Flat copy of the reference histos |
70 | // For creating a deep copy call SetOwner |
71 | // |
ce5d6908 |
72 | Int_t size = (AliPID::kSPECIES)*(kNPBins); |
338c9979 |
73 | memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size); |
ffb1ee30 |
74 | } |
75 | |
76 | //____________________________________________________________ |
77 | AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){ |
78 | // |
79 | // Assignment operator |
80 | // Performs a flat copy of the reference histos |
81 | // |
82 | if(this == &ref) return *this; |
83 | |
84 | // Make copy |
85 | TObject::operator=(ref); |
86 | if(IsOwner()){ |
87 | if(fReferences){ |
88 | fReferences->Delete(); |
89 | delete fReferences; |
90 | } |
91 | SetBit(kIsOwner, kFALSE); |
f59c1465 |
92 | } else if(fReferences) delete fReferences; |
ffb1ee30 |
93 | fReferences = ref.fReferences; |
ce5d6908 |
94 | Int_t size = (AliPID::kSPECIES)*(kNPBins); |
338c9979 |
95 | memcpy(fMapRefHists, ref.fMapRefHists, sizeof(Double_t) * size); |
ffb1ee30 |
96 | fPIDmethod = ref.fPIDmethod; |
97 | |
98 | return *this; |
99 | } |
100 | |
101 | //____________________________________________________________ |
102 | AliTRDPIDResponse::~AliTRDPIDResponse(){ |
103 | // |
104 | // Destructor |
105 | // |
f59c1465 |
106 | if(IsOwner()){ |
107 | // Destroy histos |
108 | if(fReferences) fReferences->Delete(); |
ffb1ee30 |
109 | } |
f59c1465 |
110 | if(fReferences) delete fReferences; |
ffb1ee30 |
111 | } |
112 | |
113 | //____________________________________________________________ |
114 | Bool_t AliTRDPIDResponse::Load(const Char_t * filename){ |
115 | // |
116 | // Load References into the toolkit |
117 | // |
2ad33025 |
118 | |
119 | TDirectory *owd = gDirectory;// store old working directory |
120 | |
ffb1ee30 |
121 | if(!filename) |
122 | filename = Form("%s/STEER/LQ1dRef_v1.root",gSystem->ExpandPathName("$ALICE_ROOT")); |
123 | TFile *in = TFile::Open(filename); |
124 | if(!in){ |
125 | AliError("Ref file not available."); |
126 | return kFALSE; |
127 | } |
128 | |
ce5d6908 |
129 | gROOT->cd(); |
ffb1ee30 |
130 | fReferences = new TObjArray(); |
131 | fReferences->SetOwner(); |
132 | TIter iter(in->GetListOfKeys()); |
ce5d6908 |
133 | TKey *histkey = NULL; |
ffb1ee30 |
134 | TObject *tmp = NULL; |
135 | Int_t arrayPos = 0, pbin = 0; |
136 | AliPID::EParticleType species; |
137 | TString histname; |
ce5d6908 |
138 | while((histkey = dynamic_cast<TKey *>(iter()))){ |
139 | tmp = histkey->ReadObj(); |
ffb1ee30 |
140 | histname = tmp->GetName(); |
141 | if(histname.BeginsWith("fHQel")){ // Electron histogram |
142 | histname.ReplaceAll("fHQel_p",""); |
143 | species = AliPID::kElectron; |
144 | } else { // Pion histogram |
145 | histname.ReplaceAll("fHQpi_p",""); |
146 | species = AliPID::kPion; |
147 | } |
ce5d6908 |
148 | pbin = histname.Atoi() - 1; |
1bbef4ba |
149 | AliDebug(1, Form("Species %d, Histname %s, Pbin %d, Position in container %d", species, histname.Data(), pbin, arrayPos)); |
ffb1ee30 |
150 | fMapRefHists[species][pbin] = arrayPos; |
ce5d6908 |
151 | fReferences->AddAt(new TH1F(*dynamic_cast<TH1F *>(tmp)), arrayPos); |
ffb1ee30 |
152 | arrayPos++; |
153 | } |
154 | in->Close(); |
2ad33025 |
155 | owd->cd(); |
ffb1ee30 |
156 | AliDebug(2, Form("Successfully loaded %d Reference Histograms", arrayPos)); |
f59c1465 |
157 | SetBit(kIsOwner, kTRUE); |
ffb1ee30 |
158 | delete in; |
159 | return kTRUE; |
160 | } |
161 | |
162 | //____________________________________________________________ |
b439f460 |
163 | Bool_t AliTRDPIDResponse::GetResponse(Int_t n, Double_t *dedx, Float_t *p, Double_t prob[AliPID::kSPECIES], Bool_t kNorm) const |
ffb1ee30 |
164 | { |
165 | // |
166 | // Calculate TRD likelihood values for the track based on dedx and |
167 | // momentum values. The likelihoods are calculated by query the |
168 | // reference data depending on the PID method selected |
169 | // |
170 | // Input parameter : |
171 | // n - number of dedx slices/chamber |
172 | // dedx - array of dedx slices organized layer wise |
173 | // p - array of momentum measurements organized layer wise |
174 | // |
175 | // Return parameters |
176 | // prob - probabilities allocated by TRD for particle specis |
177 | // kNorm - switch to normalize probabilities to 1. By default true. If false return not normalized prob. |
178 | // |
179 | // Return value |
180 | // true if calculation success |
181 | // |
182 | |
183 | if(!fReferences){ |
184 | AliWarning("Missing reference data. PID calculation not possible."); |
185 | return kFALSE; |
186 | } |
187 | |
188 | for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2; |
189 | Double_t prLayer[AliPID::kSPECIES]; |
ffb1ee30 |
190 | Double_t DE[10], s(0.); |
191 | for(Int_t il(kNlayer); il--;){ |
1bbef4ba |
192 | memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t)); |
ffb1ee30 |
193 | if(!CookdEdx(&dedx[il*n], &DE[0])) continue; |
194 | |
195 | s=0.; |
196 | for(Int_t is(AliPID::kSPECIES); is--;){ |
1bbef4ba |
197 | if((DE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], DE[0]); |
ce5d6908 |
198 | AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is])); |
ffb1ee30 |
199 | s+=prLayer[is]; |
200 | } |
201 | if(s<1.e-30){ |
202 | AliDebug(2, Form("Null all species prob layer[%d].", il)); |
203 | continue; |
204 | } |
205 | for(Int_t is(AliPID::kSPECIES); is--;){ |
206 | if(kNorm) prLayer[is] /= s; |
207 | prob[is] *= prLayer[is]; |
208 | } |
209 | } |
210 | if(!kNorm) return kTRUE; |
211 | |
212 | s=0.; |
213 | for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is]; |
214 | if(s<1.e-30){ |
215 | AliDebug(2, "Null total prob."); |
216 | return kFALSE; |
217 | } |
218 | for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s; |
219 | return kTRUE; |
220 | } |
221 | |
222 | |
223 | //____________________________________________________________ |
b439f460 |
224 | Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const { |
ffb1ee30 |
225 | // |
226 | // Get the non-normalized probability for a certain particle species as coming |
227 | // from the reference histogram |
228 | // Interpolation between momentum bins |
229 | // |
ce5d6908 |
230 | AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal)); |
1bbef4ba |
231 | Int_t pbin = GetLowerMomentumBin(plocal); |
232 | AliDebug(1, Form("Bin %d", pbin)); |
ffb1ee30 |
233 | Double_t pLayer = 0.; |
234 | // Do Interpolation exept for underflow and overflow |
235 | if(pbin >= 0 && pbin < kNPBins-1){ |
236 | TH1 *refHistUpper = NULL, *refHistLower = NULL; |
1bbef4ba |
237 | if(fMapRefHists[species][pbin] >= 0) |
ce5d6908 |
238 | refHistLower = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin])); |
1bbef4ba |
239 | if(fMapRefHists[species][pbin+1] >= 0) |
ce5d6908 |
240 | refHistUpper = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][pbin+1])); |
1bbef4ba |
241 | AliDebug(1, Form("Reference Histos (Upper/Lower): [%p|%p]", refHistUpper, refHistLower)); |
ffb1ee30 |
242 | |
338c9979 |
243 | if (refHistLower && refHistUpper ) { |
244 | Double_t pLower = refHistLower->GetBinContent(refHistLower->GetXaxis()->FindBin(dEdx)); |
245 | Double_t pUpper = refHistUpper->GetBinContent(refHistUpper->GetXaxis()->FindBin(dEdx)); |
ffb1ee30 |
246 | |
9590970b |
247 | pLayer = pLower + (pUpper - pLower)/(fgkPBins[pbin+1]-fgkPBins[pbin]) * (plocal - fgkPBins[pbin]); |
338c9979 |
248 | } |
ffb1ee30 |
249 | } |
250 | else{ |
251 | TH1 *refHist = NULL; |
252 | if(pbin < 0){ |
253 | // underflow |
ce5d6908 |
254 | if(fMapRefHists[species][0] >= 0){ |
255 | refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][0])); |
256 | } |
ffb1ee30 |
257 | } else { |
258 | // overflow |
ce5d6908 |
259 | if(fMapRefHists[species][kNPBins-1] > 0) |
260 | refHist = dynamic_cast<TH1 *>(fReferences->UncheckedAt(fMapRefHists[species][kNPBins-1])); |
ffb1ee30 |
261 | } |
338c9979 |
262 | if (refHist) |
263 | pLayer = refHist->GetBinContent(refHist->GetXaxis()->FindBin(dEdx)); |
ffb1ee30 |
264 | } |
ce5d6908 |
265 | AliDebug(1, Form("Probability %f", pLayer)); |
ffb1ee30 |
266 | return pLayer; |
267 | } |
268 | |
269 | //____________________________________________________________ |
b439f460 |
270 | Int_t AliTRDPIDResponse::GetLowerMomentumBin(Double_t p) const { |
ffb1ee30 |
271 | // |
272 | // Get the momentum bin for a given momentum value |
273 | // |
274 | Int_t bin = -1; |
275 | if(p > fgkPBins[kNPBins-1]) return kNPBins-1; |
276 | for(Int_t ibin = 0; ibin < kNPBins - 1; ibin++){ |
277 | if(p >= fgkPBins[ibin] && p < fgkPBins[ibin+1]){ |
278 | bin = ibin; |
279 | break; |
280 | } |
281 | } |
282 | return bin; |
283 | } |
284 | |
285 | //____________________________________________________________ |
286 | void AliTRDPIDResponse::SetOwner(){ |
287 | // |
288 | // Make Deep Copy of the Reference Histograms |
289 | // |
290 | if(!fReferences || IsOwner()) return; |
291 | TObjArray *ctmp = new TObjArray(); |
292 | for(Int_t ien = 0; ien < fReferences->GetEntriesFast(); ien++) |
293 | ctmp->AddAt(fReferences->UncheckedAt(ien)->Clone(), ien); |
f59c1465 |
294 | delete fReferences; |
ffb1ee30 |
295 | fReferences = ctmp; |
296 | SetBit(kIsOwner, kTRUE); |
297 | } |
298 | |
299 | //____________________________________________________________ |
b439f460 |
300 | Bool_t AliTRDPIDResponse::CookdEdx(Double_t *in, Double_t *out) const |
ffb1ee30 |
301 | { |
302 | switch(fPIDmethod){ |
303 | case 0: // NN |
304 | break; |
305 | case 1: // 2D LQ |
306 | break; |
307 | case 2: // 1D LQ |
308 | out[0]= 0.; |
ce5d6908 |
309 | for(Int_t islice = 0; islice < 8; islice++) out[0] += in[islice] * fGainNormalisationFactor; |
ffb1ee30 |
310 | if(out[0] < 1e-6) return kFALSE; |
311 | break; |
312 | default: |
313 | return kFALSE; |
314 | } |
315 | return kTRUE; |
316 | } |
317 | |