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 | |
51a0ce25 |
38 | #include "AliTRDPIDReference.h" |
ffb1ee30 |
39 | #include "AliTRDPIDResponse.h" |
40 | |
41 | ClassImp(AliTRDPIDResponse) |
42 | |
ffb1ee30 |
43 | //____________________________________________________________ |
e0de37e9 |
44 | AliTRDPIDResponse::AliTRDPIDResponse(): |
45 | TObject() |
51a0ce25 |
46 | ,fkPIDReference(NULL) |
e0de37e9 |
47 | ,fGainNormalisationFactor(1.) |
48 | ,fPIDmethod(kLQ1D) |
ffb1ee30 |
49 | { |
e0de37e9 |
50 | // |
51 | // Default constructor |
52 | // |
ffb1ee30 |
53 | } |
54 | |
55 | //____________________________________________________________ |
56 | AliTRDPIDResponse::AliTRDPIDResponse(const AliTRDPIDResponse &ref): |
e0de37e9 |
57 | TObject(ref) |
51a0ce25 |
58 | ,fkPIDReference(ref.fkPIDReference) |
e0de37e9 |
59 | ,fGainNormalisationFactor(ref.fGainNormalisationFactor) |
60 | ,fPIDmethod(ref.fPIDmethod) |
ffb1ee30 |
61 | { |
e0de37e9 |
62 | // |
63 | // Copy constructor |
e0de37e9 |
64 | // |
ffb1ee30 |
65 | } |
66 | |
67 | //____________________________________________________________ |
68 | AliTRDPIDResponse &AliTRDPIDResponse::operator=(const AliTRDPIDResponse &ref){ |
e0de37e9 |
69 | // |
70 | // Assignment operator |
e0de37e9 |
71 | // |
72 | if(this == &ref) return *this; |
73 | |
74 | // Make copy |
75 | TObject::operator=(ref); |
51a0ce25 |
76 | fGainNormalisationFactor = ref.fGainNormalisationFactor; |
77 | fkPIDReference = ref.fkPIDReference; |
e0de37e9 |
78 | fPIDmethod = ref.fPIDmethod; |
79 | |
80 | return *this; |
ffb1ee30 |
81 | } |
82 | |
83 | //____________________________________________________________ |
84 | AliTRDPIDResponse::~AliTRDPIDResponse(){ |
e0de37e9 |
85 | // |
86 | // Destructor |
87 | // |
51a0ce25 |
88 | if(IsOwner()) delete fkPIDReference; |
ffb1ee30 |
89 | } |
90 | |
91 | //____________________________________________________________ |
51a0ce25 |
92 | Bool_t AliTRDPIDResponse::Load(const Char_t * filename, const Char_t *refName){ |
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(); |
51a0ce25 |
108 | fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(in->Get(refName)->Clone()); |
109 | in->Close(); delete in; |
e0de37e9 |
110 | owd->cd(); |
e0de37e9 |
111 | SetBit(kIsOwner, kTRUE); |
51a0ce25 |
112 | AliDebug(2, Form("Successfully loaded References for %d Momentum bins", fkPIDReference->GetNumberOfMomentumBins())); |
e0de37e9 |
113 | return kTRUE; |
ffb1ee30 |
114 | } |
115 | |
116 | //____________________________________________________________ |
51a0ce25 |
117 | 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 |
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 |
134 | // true if calculation success |
135 | // |
ffb1ee30 |
136 | |
51a0ce25 |
137 | if(!fkPIDReference){ |
e0de37e9 |
138 | AliWarning("Missing reference data. PID calculation not possible."); |
139 | return kFALSE; |
140 | } |
ffb1ee30 |
141 | |
e0de37e9 |
142 | for(Int_t is(AliPID::kSPECIES); is--;) prob[is]=.2; |
143 | Double_t prLayer[AliPID::kSPECIES]; |
51a0ce25 |
144 | Double_t dE[10], s(0.); |
e0de37e9 |
145 | for(Int_t il(kNlayer); il--;){ |
146 | memset(prLayer, 0, AliPID::kSPECIES*sizeof(Double_t)); |
51a0ce25 |
147 | if(!CookdEdx(n, &dedx[il*n], &dE[0])) continue; |
ffb1ee30 |
148 | |
e0de37e9 |
149 | s=0.; |
150 | for(Int_t is(AliPID::kSPECIES); is--;){ |
51a0ce25 |
151 | if((dE[0] > 0.) && (p[il] > 0.)) prLayer[is] = GetProbabilitySingleLayer(is, p[il], dE[0]); |
e0de37e9 |
152 | AliDebug(3, Form("Probability for Species %d in Layer %d: %f", is, il, prLayer[is])); |
153 | s+=prLayer[is]; |
154 | } |
155 | if(s<1.e-30){ |
156 | AliDebug(2, Form("Null all species prob layer[%d].", il)); |
157 | continue; |
158 | } |
159 | for(Int_t is(AliPID::kSPECIES); is--;){ |
160 | if(kNorm) prLayer[is] /= s; |
161 | prob[is] *= prLayer[is]; |
162 | } |
163 | } |
164 | if(!kNorm) return kTRUE; |
ffb1ee30 |
165 | |
e0de37e9 |
166 | s=0.; |
167 | for(Int_t is(AliPID::kSPECIES); is--;) s+=prob[is]; |
168 | if(s<1.e-30){ |
169 | AliDebug(2, "Null total prob."); |
170 | return kFALSE; |
171 | } |
172 | for(Int_t is(AliPID::kSPECIES); is--;) prob[is]/=s; |
173 | return kTRUE; |
ffb1ee30 |
174 | } |
175 | |
ffb1ee30 |
176 | //____________________________________________________________ |
b439f460 |
177 | Double_t AliTRDPIDResponse::GetProbabilitySingleLayer(Int_t species, Double_t plocal, Double_t dEdx) const { |
e0de37e9 |
178 | // |
179 | // Get the non-normalized probability for a certain particle species as coming |
180 | // from the reference histogram |
181 | // Interpolation between momentum bins |
182 | // |
183 | AliDebug(1, Form("Make Probability for Species %d with Momentum %f", species, plocal)); |
51a0ce25 |
184 | Float_t pLower, pUpper; |
185 | Double_t probLayer = 0.; |
186 | TH1 *refUpper = dynamic_cast<TH1 *>(fkPIDReference->GetUpperReference((AliPID::EParticleType)species, plocal, pUpper)), |
187 | *refLower = dynamic_cast<TH1 *>(fkPIDReference->GetLowerReference((AliPID::EParticleType)species, plocal, pLower)); |
e0de37e9 |
188 | // Do Interpolation exept for underflow and overflow |
51a0ce25 |
189 | if(refLower && refUpper){ |
190 | Double_t probLower = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx)); |
191 | Double_t probUpper = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx)); |
e0de37e9 |
192 | |
51a0ce25 |
193 | probLayer = probLower + (probUpper - probLower)/(pUpper-pLower) * (plocal - pLower); |
194 | } else if(refLower){ |
195 | // underflow |
196 | probLayer = refLower->GetBinContent(refLower->GetXaxis()->FindBin(dEdx)); |
197 | } else if(refUpper){ |
198 | // overflow |
199 | probLayer = refUpper->GetBinContent(refUpper->GetXaxis()->FindBin(dEdx)); |
200 | } else { |
201 | AliError("No references available"); |
e0de37e9 |
202 | } |
51a0ce25 |
203 | AliDebug(1, Form("Probability %f", probLayer)); |
204 | return probLayer; |
ffb1ee30 |
205 | } |
206 | |
207 | //____________________________________________________________ |
208 | void AliTRDPIDResponse::SetOwner(){ |
e0de37e9 |
209 | // |
210 | // Make Deep Copy of the Reference Histograms |
211 | // |
51a0ce25 |
212 | if(!fkPIDReference || IsOwner()) return; |
213 | const AliTRDPIDReference *tmp = fkPIDReference; |
214 | fkPIDReference = dynamic_cast<const AliTRDPIDReference *>(tmp->Clone()); |
e0de37e9 |
215 | SetBit(kIsOwner, kTRUE); |
ffb1ee30 |
216 | } |
217 | |
218 | //____________________________________________________________ |
51a0ce25 |
219 | Bool_t AliTRDPIDResponse::CookdEdx(Int_t nSlice, const Double_t * const in, Double_t *out) const |
ffb1ee30 |
220 | { |
51a0ce25 |
221 | // |
222 | // Recalculate dE/dx |
223 | // |
e0de37e9 |
224 | switch(fPIDmethod){ |
225 | case kNN: // NN |
226 | break; |
227 | case kLQ2D: // 2D LQ |
228 | break; |
229 | case kLQ1D: // 1D LQ |
230 | out[0]= 0.; |
231 | for(Int_t islice = 0; islice < nSlice; islice++) out[0] += in[islice] * fGainNormalisationFactor; |
232 | if(out[0] < 1e-6) return kFALSE; |
233 | break; |
234 | default: |
235 | return kFALSE; |
236 | } |
237 | return kTRUE; |
ffb1ee30 |
238 | } |
239 | |