]>
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 | } |
accb90b5 | 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 | } | |
fa8f5c53 | 96 | return *this; |
accb90b5 | 97 | } |
98 | ||
9189f38b | 99 | //______________________________________________________________________ |
92329a66 | 100 | void AliITSdEdxSamples::SetdESamples(Int_t nSamples, Double_t* samples){ |
9189f38b | 101 | // Set the samples |
102 | ||
92329a66 | 103 | if(nSamples>kMaxSamples){ |
104 | AliWarning(Form("Too many dE samples,only first %d will be used",kMaxSamples)); | |
105 | fNSamples=kMaxSamples; | |
9189f38b | 106 | }else{ |
107 | fNSamples=nSamples; | |
108 | } | |
92329a66 | 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.; | |
9189f38b | 111 | return; |
112 | } | |
113 | //______________________________________________________________________ | |
92329a66 | 114 | void AliITSdEdxSamples::SetdxSamples(Int_t nSamples, Double_t* samples){ |
9189f38b | 115 | // Set the samples |
92329a66 | 116 | |
117 | if(nSamples>kMaxSamples){ | |
13242232 | 118 | AliWarning(Form("Too many dx samples,only first %d will be used",kMaxSamples)); |
92329a66 | 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); | |
9189f38b | 133 | for(Int_t i=0; i<fNSamples; i++) fPAtSample[i]=mom[i]; |
92329a66 | 134 | for(Int_t i=fNSamples; i<kMaxSamples; i++) fPAtSample[i]=0.; |
9189f38b | 135 | return; |
136 | } | |
137 | //______________________________________________________________________ | |
0d66557a | 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 | //______________________________________________________________________ | |
9189f38b | 154 | Double_t AliITSdEdxSamples::GetTruncatedMean(Double_t frac, Double_t mindedx) const { |
155 | // compute truncated mean | |
156 | ||
157 | Int_t nc=0; | |
92329a66 | 158 | Double_t dedx[kMaxSamples]; |
9189f38b | 159 | for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values |
92329a66 | 160 | Double_t dedxsamp=GetdEdxSample(il); |
0d66557a | 161 | if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){ |
92329a66 | 162 | dedx[nc]= dedxsamp; |
9189f38b | 163 | nc++; |
92329a66 | 164 | } |
9189f38b | 165 | } |
166 | if(nc<1) return 0.; | |
0d66557a | 167 | |
9189f38b | 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; | |
92329a66 | 181 | Double_t weight[kMaxSamples]; |
182 | for(Int_t iw=0; iw<kMaxSamples; iw++) weight[iw]=0.; | |
9189f38b | 183 | Int_t lastUsed=(Int_t)(frac*nc+0.00001); |
184 | if(lastUsed==0) lastUsed=1; | |
92329a66 | 185 | if(lastUsed>kMaxSamples) lastUsed=kMaxSamples; |
9189f38b | 186 | for(Int_t iw=0; iw<lastUsed; iw++) weight[iw]=1.; |
92329a66 | 187 | if((frac*nc-lastUsed)>0.4 && lastUsed<kMaxSamples) weight[lastUsed]=0.5; |
9189f38b | 188 | for (Int_t i=0; i<nc; i++) { |
92329a66 | 189 | // AliDebug(5,Form("dE/dx %f weight %f",dedx[i],weight[i])); |
9189f38b | 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) | |
9189f38b | 200 | Int_t nc=0; |
92329a66 | 201 | Double_t dedx[kMaxSamples]; |
9189f38b | 202 | for (Int_t il=0; il<fNSamples; il++) { // count good (>0) dE/dx values |
92329a66 | 203 | Double_t dedxsamp=GetdEdxSample(il); |
0d66557a | 204 | if(HasPointOnLayer(il) && UseLayerForPid(il) && dedxsamp>mindedx){ |
92329a66 | 205 | dedx[nc]= dedxsamp; |
206 | nc++; | |
9189f38b | 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 | //______________________________________________________________________ | |
b34fed4e | 221 | void AliITSdEdxSamples::GetConditionalProbabilities(const AliITSPidParams* pars, Double_t condprob[AliPID::kSPECIES], Double_t mindedx) const { |
9189f38b | 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++){ | |
0d66557a | 227 | if(!HasPointOnLayer(iS)) continue; |
228 | if(!UseLayerForPid(iS)) continue; | |
9189f38b | 229 | Int_t iLayer=iS+3; // to match with present ITS |
230 | if(iLayer>6) iLayer=6; // all extra points are treated as SSD | |
92329a66 | 231 | Float_t dedx = GetdEdxSample(iS); |
232 | if(dedx<mindedx) continue; | |
9189f38b | 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 | } | |
0d66557a | 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 | } |