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