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