]>
Commit | Line | Data |
---|---|---|
1 | /************************************************************************** | |
2 | * Copyright(c) 2009-2012, 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 | /* $Id$ */ | |
16 | ||
17 | ||
18 | ||
19 | /////////////////////////////////////////////////////////////////// | |
20 | // // | |
21 | // Class to store information for PID with ITS // | |
22 | // and truncated mean computation methods // | |
23 | // Origin: F.Prino, Torino, prino@to.infn.it // | |
24 | // // | |
25 | /////////////////////////////////////////////////////////////////// | |
26 | ||
27 | #include "AliITSPidParams.h" | |
28 | #include "AliITSdEdxSamples.h" | |
29 | #include "AliLog.h" | |
30 | #include <TMath.h> | |
31 | ||
32 | ClassImp(AliITSdEdxSamples) | |
33 | ||
34 | //______________________________________________________________________ | |
35 | AliITSdEdxSamples::AliITSdEdxSamples():TObject(), | |
36 | fNSamples(0), | |
37 | fClusterMap(0), | |
38 | fP(0.), | |
39 | fParticleSpecie(0), | |
40 | fLayersForPid(0xFFFF) | |
41 | { | |
42 | // Default constructor | |
43 | for(Int_t i=0; i<kMaxSamples; i++){ | |
44 | fdESamples[i]=0.; | |
45 | fdxSamples[i]=0.; | |
46 | fPAtSample[i]=0.; | |
47 | } | |
48 | } | |
49 | ||
50 | //______________________________________________________________________ | |
51 | AliITSdEdxSamples::AliITSdEdxSamples(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t mom, Int_t specie) : | |
52 | TObject(), | |
53 | fNSamples(nSamples), | |
54 | fClusterMap(0), | |
55 | fP(mom), | |
56 | fParticleSpecie(specie), | |
57 | fLayersForPid(0xFFFF) | |
58 | { | |
59 | // Standard constructor | |
60 | SetdESamples(nSamples,esamples); | |
61 | SetdxSamples(nSamples,xsamples); | |
62 | SetClusterMapFromdE(); | |
63 | } | |
64 | ||
65 | //______________________________________________________________________ | |
66 | AliITSdEdxSamples::AliITSdEdxSamples(const AliITSdEdxSamples& source) : | |
67 | TObject(), | |
68 | fNSamples(source.fNSamples), | |
69 | fClusterMap(source.fClusterMap), | |
70 | fP(source.fP), | |
71 | fParticleSpecie(source.fParticleSpecie), | |
72 | fLayersForPid(source.fLayersForPid) | |
73 | { | |
74 | // Copy constructor | |
75 | for(Int_t i=0; i<kMaxSamples; i++){ | |
76 | fdESamples[i]=source.GetdESample(i); | |
77 | fdxSamples[i]=source.GetdxSample(i); | |
78 | fPAtSample[i]=source.GetMomentumAtSample(i); | |
79 | } | |
80 | } | |
81 | //_____________________________________________________________________________ | |
82 | AliITSdEdxSamples& AliITSdEdxSamples::operator=(const AliITSdEdxSamples &source){ | |
83 | // Assignment operator | |
84 | if(this==&source) return *this; | |
85 | ((TObject *)this)->operator=(source); | |
86 | fNSamples = source.fNSamples; | |
87 | fClusterMap = source.fClusterMap; | |
88 | fP = source.fP; | |
89 | fParticleSpecie = source.fParticleSpecie; | |
90 | fLayersForPid = source.fLayersForPid; | |
91 | for(Int_t i=0; i<kMaxSamples; i++){ | |
92 | fdESamples[i]=source.GetdESample(i); | |
93 | fdxSamples[i]=source.GetdxSample(i); | |
94 | fPAtSample[i]=source.GetMomentumAtSample(i); | |
95 | } | |
96 | return *this; | |
97 | } | |
98 | ||
99 | //______________________________________________________________________ | |
100 | void AliITSdEdxSamples::SetdESamples(Int_t nSamples, Double_t* samples){ | |
101 | // Set the samples | |
102 | ||
103 | if(nSamples>kMaxSamples){ | |
104 | AliWarning(Form("Too many dE samples,only first %d will be used",kMaxSamples)); | |
105 | fNSamples=kMaxSamples; | |
106 | }else{ | |
107 | fNSamples=nSamples; | |
108 | } | |
109 | for(Int_t i=0; i<fNSamples; i++) fdESamples[i]=samples[i]; | |
110 | for(Int_t i=fNSamples; i<kMaxSamples; i++) fdESamples[i]=0.; | |
111 | return; | |
112 | } | |
113 | //______________________________________________________________________ | |
114 | void AliITSdEdxSamples::SetdxSamples(Int_t nSamples, Double_t* samples){ | |
115 | // Set the samples | |
116 | ||
117 | if(nSamples>kMaxSamples){ | |
118 | AliWarning(Form("Too many dx samples,only first %d will be used",kMaxSamples)); | |
119 | fNSamples=kMaxSamples; | |
120 | }else{ | |
121 | fNSamples=nSamples; | |
122 | } | |
123 | for(Int_t i=0; i<fNSamples; i++) fdxSamples[i]=samples[i]; | |
124 | for(Int_t i=fNSamples; i<kMaxSamples; i++) fdxSamples[i]=0.; | |
125 | return; | |
126 | } | |
127 | ||
128 | //______________________________________________________________________ | |
129 | void AliITSdEdxSamples::SetSamplesAndMomenta(Int_t nSamples, Double_t* esamples, Double_t* xsamples, Double_t* mom){ | |
130 | // Set the samples | |
131 | SetdESamples(nSamples,esamples); | |
132 | SetdxSamples(nSamples,xsamples); | |
133 | for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i]; | |
134 | for(Int_t i=fNSamples; i<kMaxSamples; i++) fPAtSample[i]=0.; | |
135 | return; | |
136 | } | |
137 | //______________________________________________________________________ | |
138 | void AliITSdEdxSamples::SetLayerSample(Int_t iLayer, Bool_t haspoint, Double_t dE, Double_t dx, Double_t p){ | |
139 | // set info from single layer | |
140 | if(haspoint){ | |
141 | SetPointOnLayer(iLayer); | |
142 | fdESamples[iLayer]=dE; | |
143 | fdxSamples[iLayer]=dx; | |
144 | fPAtSample[iLayer]=p; | |
145 | }else{ | |
146 | if(HasPointOnLayer(iLayer)) fClusterMap-=(1<<iLayer); | |
147 | fdESamples[iLayer]=0.; | |
148 | fdxSamples[iLayer]=0.; | |
149 | fPAtSample[iLayer]=0.; | |
150 | ||
151 | } | |
152 | } | |
153 | //______________________________________________________________________ | |
154 | Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const { | |
155 | // compute truncated mean | |
156 | ||
157 | Int_t nc=0; | |
158 | Double_t dedx[kMaxSamples]; | |
159 | for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values | |
160 | Double_t dedxsamp=GetdEdxSample(il); | |
161 | if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){ | |
162 | dedx[nc]= dedxsamp; | |
163 | nc++; | |
164 | } | |
165 | } | |
166 | if(nc<1) return 0.; | |
167 | ||
168 | Int_t swap; // sort in ascending order | |
169 | do { | |
170 | swap=0; | |
171 | for (Int_t i=0; i<nc-1; i++) { | |
172 | if (dedx[i]<=dedx[i+1]) continue; | |
173 | Double_t tmp=dedx[i]; | |
174 | dedx[i]=dedx[i+1]; | |
175 | dedx[i+1]=tmp; | |
176 | swap++; | |
177 | } | |
178 | } while (swap); | |
179 | ||
180 | Double_t sumamp=0,sumweight=0; | |
181 | Double_t weight[kMaxSamples]; | |
182 | for(Int_t iw=0; iw<kMaxSamples; iw++) weight[iw]=0.; | |
183 | Int_t lastUsed=(Int_t)(frac*nc+0.00001); | |
184 | if(lastUsed==0) lastUsed=1; | |
185 | if(lastUsed>kMaxSamples) lastUsed=kMaxSamples; | |
186 | for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.; | |
187 | if((frac*nc-lastUsed)>0.4 && lastUsed<kMaxSamples) weight[lastUsed]=0.5; | |
188 | for (Int_t i=0; i<nc; i++) { | |
189 | // AliDebug(5,Form("dE/dx %f weight %f",dedx[i],weight[i])); | |
190 | sumamp+= dedx[i]*weight[i]; | |
191 | sumweight+=weight[i]; | |
192 | } | |
193 | if(sumweight>0.) return sumamp/sumweight; | |
194 | else return 0.; | |
195 | ||
196 | } | |
197 | //______________________________________________________________________ | |
198 | Double_t AliITSdEdxSamples::GetWeightedMean(Double_t mindedx) const { | |
199 | // compute generalized mean with k=-2 (used by CMS) | |
200 | Int_t nc=0; | |
201 | Double_t dedx[kMaxSamples]; | |
202 | for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values | |
203 | Double_t dedxsamp=GetdEdxSample(il); | |
204 | if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){ | |
205 | dedx[nc]= dedxsamp; | |
206 | nc++; | |
207 | } | |
208 | } | |
209 | if(nc<1) return 0.; | |
210 | ||
211 | Double_t weiSum = 0.; | |
212 | for (Int_t i=0; i<nc; i++) { | |
213 | weiSum+=TMath::Power(dedx[i],-2); | |
214 | } | |
215 | Double_t wMean=0.; | |
216 | if(weiSum>0.) wMean= TMath::Power(weiSum/nc,-0.5); | |
217 | return wMean; | |
218 | ||
219 | } | |
220 | //______________________________________________________________________ | |
221 | void AliITSdEdxSamples::GetConditionalProbabilities(const AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const { | |
222 | // compute conditional probablilities | |
223 | const Int_t nPart = 3; | |
224 | Double_t itsProb[nPart] = {1,1,1}; // p, K, pi | |
225 | ||
226 | for(Int_t iS=0; iS<fNSamples; iS++){ | |
227 | if(!HasPointOnLayer(iS)) continue; | |
228 | if(!UseLayerForPid(iS)) continue; | |
229 | Int_t iLayer=iS+3; // to match with present ITS | |
230 | if(iLayer>6) iLayer=6; // all extra points are treated as SSD | |
231 | Float_t dedx = GetdEdxSample(iS); | |
232 | if(dedx<mindedx) continue; | |
233 | Float_t layProb = pars->GetLandauGausNorm(dedx,AliPID::kProton,fP,iLayer); | |
234 | itsProb[0] *= layProb; | |
235 | ||
236 | layProb = pars->GetLandauGausNorm(dedx,AliPID::kKaon,fP,iLayer); | |
237 | if (fP < 0.16) layProb=0.00001; | |
238 | itsProb[1] *= layProb; | |
239 | ||
240 | layProb = pars->GetLandauGausNorm(dedx,AliPID::kPion,fP,iLayer); | |
241 | itsProb[2] *= layProb; | |
242 | } | |
243 | ||
244 | // Normalise probabilities | |
245 | Double_t sumProb = 0; | |
246 | for (Int_t iPart = 0; iPart < nPart; iPart++) { | |
247 | sumProb += itsProb[iPart]; | |
248 | } | |
249 | sumProb += 2*itsProb[2]; // muon and electron cannot be distinguished from pions | |
250 | ||
251 | for (Int_t iPart = 0; iPart < nPart; iPart++) { | |
252 | itsProb[iPart]/=sumProb; | |
253 | } | |
254 | ||
255 | condprob[AliPID::kElectron] = itsProb[2]; | |
256 | condprob[AliPID::kMuon] = itsProb[2]; | |
257 | condprob[AliPID::kPion] = itsProb[2]; | |
258 | condprob[AliPID::kKaon] = itsProb[1]; | |
259 | condprob[AliPID::kProton] = itsProb[0]; | |
260 | return; | |
261 | ||
262 | } | |
263 | //______________________________________________________________________ | |
264 | void AliITSdEdxSamples::PrintAll() const{ | |
265 | // print all the infos | |
266 | printf("Particle %d momentum %f GeV/c, number of points %d\n", | |
267 | GetParticleSpecieMC(), | |
268 | fP, | |
269 | GetNumberOfEffectiveSamples()); | |
270 | for(Int_t iLay=0; iLay<fNSamples; iLay++){ | |
271 | printf(" Layer %d Point %d dE %f keV dx %f cm mom %f GeV/c\n",iLay, | |
272 | HasPointOnLayer(iLay), | |
273 | GetdESample(iLay), | |
274 | GetdxSample(iLay), | |
275 | GetMomentumAtSample(iLay)); | |
276 | } | |
277 | ||
278 | printf("Layers used for PID:\n"); | |
279 | printf("Layer "); | |
280 | for(Int_t iLay=0; iLay<fNSamples; iLay++){ | |
281 | printf("%d ",iLay); | |
282 | } | |
283 | printf("\n"); | |
284 | ||
285 | printf("Use "); | |
286 | for(Int_t iLay=0; iLay<fNSamples; iLay++){ | |
287 | printf("%d ",UseLayerForPid(iLay)); | |
288 | } | |
289 | printf("\n"); | |
290 | printf("Truncated mean = %f\n",GetTruncatedMean()); | |
291 | } | |
292 | //______________________________________________________________________ | |
293 | void AliITSdEdxSamples::PrintClusterMap() const{ | |
294 | // print the cluster map | |
295 | ||
296 | printf("Layer "); | |
297 | for(Int_t iLay=0; iLay<fNSamples; iLay++){ | |
298 | printf("%d ",iLay); | |
299 | } | |
300 | printf("\n"); | |
301 | ||
302 | printf("Point "); | |
303 | for(Int_t iLay=0; iLay<fNSamples; iLay++){ | |
304 | printf("%d ",HasPointOnLayer(iLay)); | |
305 | } | |
306 | printf("\n"); | |
307 | } |