]>
Commit | Line | Data |
---|---|---|
e62c1aea | 1 | ////////////////////////////////////////////////////////////////////////// |
2 | // USER Class for PID in the ITS // | |
3 | //The PID is based on the likelihood of all the four ITS' layers, // | |
4 | //without using the truncated mean for the dE/dx. The response // | |
5 | //functions for each layer are convoluted Landau-Gaussian functions. // | |
6 | // Origin: Elena Bruna bruna@to.infn.it,, Massimo Masera masera@to.infn.it// | |
7 | ////////////////////////////////////////////////////////////////////////// | |
b68da4a2 | 8 | #include "AliITStrackV2.h" |
e62c1aea | 9 | #include "AliITSPident.h" |
b68da4a2 | 10 | #include "AliITSSteerPid.h" |
11 | #include <Riostream.h> | |
e62c1aea | 12 | |
13 | ClassImp(AliITSPident) | |
14 | //_______________________________________________________________________ | |
94631b2f | 15 | AliITSPident::AliITSPident(): |
16 | fMom(0), | |
17 | fdEdx(0), | |
18 | fPBayesp(0), | |
19 | fPBayesk(0), | |
20 | fPBayespi(0), | |
21 | fPPriorip(0), | |
22 | fPPriorik(0), | |
23 | fPPrioripi(0), | |
b68da4a2 | 24 | fPPriorie(0) |
94631b2f | 25 | { |
e62c1aea | 26 | // default constructor |
b68da4a2 | 27 | for (Int_t i=0;i<8;i++){ |
e62c1aea | 28 | fCondFunProLay[i]=0; |
29 | fCondFunKLay[i]=0; | |
30 | fCondFunPiLay[i]=0; | |
31 | } | |
b68da4a2 | 32 | for (Int_t i=0;i<4;i++)fNcls[i]=0; |
e62c1aea | 33 | } |
34 | //_______________________________________________________________________ | |
35 | AliITSPident::~AliITSPident(){ | |
36 | // destructor | |
37 | } | |
38 | //______________________________________________________________________ | |
94631b2f | 39 | AliITSPident::AliITSPident(const AliITSPident &ob) :TObject(ob), |
40 | fMom(ob.fMom), | |
41 | fdEdx(ob.fdEdx), | |
42 | fPBayesp(ob.fPBayesp), | |
43 | fPBayesk(ob.fPBayesk), | |
44 | fPBayespi(ob.fPBayespi), | |
45 | fPPriorip(ob.fPPriorip), | |
46 | fPPriorik(ob.fPPriorik), | |
47 | fPPrioripi(ob.fPPrioripi), | |
b68da4a2 | 48 | fPPriorie(ob.fPPriorie) |
94631b2f | 49 | { |
e62c1aea | 50 | // Copy constructor |
e62c1aea | 51 | } |
52 | ||
53 | //______________________________________________________________________ | |
94631b2f | 54 | AliITSPident& AliITSPident::operator=(const AliITSPident& ob){ |
e62c1aea | 55 | // Assignment operator |
94631b2f | 56 | this->~AliITSPident(); |
57 | new(this) AliITSPident(ob); | |
e62c1aea | 58 | return *this; |
59 | } | |
60 | ||
61 | ||
62 | //_______________________________________________________________________ | |
b68da4a2 | 63 | AliITSPident::AliITSPident(Double_t mom,Double_t dEdx,AliITSSteerPid *sp,Float_t *Qlay,Float_t *nlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie): |
94631b2f | 64 | fMom(mom), |
65 | fdEdx(dEdx), | |
66 | fPBayesp(0), | |
67 | fPBayesk(0), | |
68 | fPBayespi(0), | |
69 | fPPriorip(priorip), | |
70 | fPPriorik(priorik), | |
71 | fPPrioripi(prioripi), | |
b68da4a2 | 72 | fPPriorie(priorie) |
73 | { | |
74 | // | |
75 | for (Int_t i=0;i<8;i++){ | |
76 | fCondFunProLay[i]=-1; | |
77 | fCondFunKLay[i]=-1; | |
78 | fCondFunPiLay[i]=-1; | |
79 | } | |
e62c1aea | 80 | for(Int_t la=0;la<4;la++){//loop on layers |
81 | Double_t parp[3];Double_t park[3];Double_t parpi[3]; | |
b68da4a2 | 82 | fNcls[la]=0; |
e62c1aea | 83 | Double_t range[6]; |
84 | range[0]=0.3*parp[1]; | |
85 | range[1]=2.*parp[1]; | |
86 | ||
87 | range[2]=0.3*park[1]; | |
88 | range[3]=2.*park[1]; | |
89 | ||
90 | range[4]=0.3*parpi[1]; | |
91 | range[5]=2.*parpi[1]; | |
b68da4a2 | 92 | Int_t layer=la+2; |
93 | for(Int_t ii=0;ii<8;ii++){ | |
94 | if(nlay[ii]==layer){ | |
95 | fNcls[la]++; | |
96 | if(Qlay[ii]>0){ | |
97 | sp->GetParFitLayer(la,fMom,parp,park,parpi); | |
98 | CookFunItsLay(ii,0,parp,Qlay[ii],fMom,range[0],range[1],"fPro"); | |
99 | CookFunItsLay(ii,1,park,Qlay[ii],fMom,range[2],range[3],"fKao"); | |
100 | CookFunItsLay(ii,2,parpi,Qlay[ii],fMom,range[4],range[5],"fPi"); | |
101 | } | |
102 | } | |
103 | } | |
e62c1aea | 104 | |
105 | } | |
b68da4a2 | 106 | |
107 | Float_t prior[4];Double_t condFun[8][3]; | |
e62c1aea | 108 | |
109 | prior[0]=fPPriorip; | |
110 | prior[1]=fPPriorik; | |
111 | prior[2]=fPPrioripi; | |
112 | prior[3]=fPPriorie; | |
b68da4a2 | 113 | for(Int_t la=0;la<8;la++){ |
e62c1aea | 114 | condFun[la][0]= fCondFunProLay[la]; |
115 | condFun[la][1]= fCondFunKLay[la]; | |
116 | condFun[la][2]= fCondFunPiLay[la]; | |
117 | ||
118 | } | |
119 | ||
e62c1aea | 120 | fPBayesp=CookCombinedBayes(condFun,prior,0); |
e62c1aea | 121 | fPBayesk=CookCombinedBayes(condFun,prior,1); |
e62c1aea | 122 | fPBayespi=CookCombinedBayes(condFun,prior,2); |
b68da4a2 | 123 | } |
e62c1aea | 124 | |
b68da4a2 | 125 | //__________________________________________________________________________________________ |
126 | AliITSPident::AliITSPident(AliITStrackV2 *trackITS,AliITSSteerPid *sp,Float_t *Qlay,Float_t *nlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie): | |
94631b2f | 127 | fMom(0), |
128 | fdEdx(0), | |
129 | fPBayesp(0), | |
130 | fPBayesk(0), | |
131 | fPBayespi(0), | |
132 | fPPriorip(priorip), | |
133 | fPPriorik(priorik), | |
134 | fPPrioripi(prioripi), | |
b68da4a2 | 135 | fPPriorie(priorie) |
94631b2f | 136 | { |
e62c1aea | 137 | // |
b68da4a2 | 138 | for (Int_t i=0;i<8;i++){ |
139 | fCondFunProLay[i]=-1; | |
140 | fCondFunKLay[i]=-1; | |
141 | fCondFunPiLay[i]=-1; | |
142 | } | |
e62c1aea | 143 | Double_t xr; |
144 | Double_t par[5]; | |
145 | trackITS->GetExternalParameters(xr,par); | |
146 | if (par[4]!=0) { | |
147 | Float_t lamb=TMath::ATan(par[3]); | |
148 | fMom=1/(TMath::Abs(par[4])*TMath::Cos(lamb)); | |
149 | } | |
150 | ||
e62c1aea | 151 | |
152 | for(Int_t la=0;la<4;la++){//loop on layers | |
153 | Double_t parp[3];Double_t park[3];Double_t parpi[3]; | |
b68da4a2 | 154 | fNcls[la]=0; |
e62c1aea | 155 | Double_t range[8]; |
156 | range[0]=0.3*parp[1]; | |
157 | range[1]=2.*parp[1]; | |
158 | ||
159 | range[2]=0.3*park[1]; | |
160 | range[3]=2.*park[1]; | |
b68da4a2 | 161 | |
e62c1aea | 162 | range[4]=0.3*parpi[1]; |
163 | range[5]=2.*parpi[1]; | |
b68da4a2 | 164 | |
165 | Int_t layer=la+2; | |
166 | for(Int_t ii=0;ii<8;ii++){ | |
167 | if(nlay[ii]==layer){ | |
168 | fNcls[la]++; | |
169 | if(Qlay[ii]>0){ | |
170 | sp->GetParFitLayer(la,fMom,parp,park,parpi); | |
171 | CookFunItsLay(ii,0,parp,Qlay[ii],fMom,range[0],range[1],"fPro"); | |
172 | CookFunItsLay(ii,1,park,Qlay[ii],fMom,range[2],range[3],"fKao"); | |
173 | CookFunItsLay(ii,2,parpi,Qlay[ii],fMom,range[4],range[5],"fPi"); | |
174 | } | |
175 | } | |
176 | } | |
177 | ||
e62c1aea | 178 | } |
e62c1aea | 179 | |
b68da4a2 | 180 | Float_t prior[4];Double_t condFun[8][3]; |
e62c1aea | 181 | |
182 | prior[0]=fPPriorip; | |
183 | prior[1]=fPPriorik; | |
184 | prior[2]=fPPrioripi; | |
185 | prior[3]=fPPriorie; | |
b68da4a2 | 186 | for(Int_t la=0;la<8;la++){ |
e62c1aea | 187 | condFun[la][0]= fCondFunProLay[la]; |
188 | condFun[la][1]= fCondFunKLay[la]; | |
189 | condFun[la][2]= fCondFunPiLay[la]; | |
190 | ||
191 | } | |
192 | ||
e62c1aea | 193 | fPBayesp=CookCombinedBayes(condFun,prior,0); |
e62c1aea | 194 | fPBayesk=CookCombinedBayes(condFun,prior,1); |
e62c1aea | 195 | fPBayespi=CookCombinedBayes(condFun,prior,2); |
e62c1aea | 196 | fdEdx=trackITS->GetdEdx(); |
197 | ||
198 | } | |
199 | //_______________________________________________________________________ | |
b68da4a2 | 200 | void AliITSPident::GetNclsPerLayer(Int_t *ncls) const{ |
201 | //number of clusters for each layer (sdd1,sdd2,ssd1,ssd2) | |
202 | for(Int_t la=0;la<4;la++){ | |
203 | ncls[la]=fNcls[la]; | |
204 | } | |
205 | ||
206 | }//_______________________________________________________________________ | |
207 | Double_t AliITSPident::GetProdCondFunPro() const { | |
208 | //Product of conditional probability functions for protons | |
209 | Double_t rv=1.; | |
210 | for(Int_t i=0;i<8;i++){ | |
211 | Double_t fun=GetCondFunPro(i); | |
212 | if(fun>=0)rv*=fun; | |
213 | } | |
214 | return rv; | |
215 | }//_______________________________________________________________________ | |
216 | Double_t AliITSPident::GetProdCondFunK() const { | |
217 | //Product of conditional probability functions for kaons | |
218 | Double_t rv=1.; | |
219 | for(Int_t i=0;i<8;i++){ | |
220 | Double_t fun=GetCondFunK(i); | |
221 | if(fun>=0)rv*=fun; | |
222 | } | |
223 | return rv; | |
224 | } | |
225 | //_______________________________________________________________________ | |
226 | Double_t AliITSPident::GetProdCondFunPi() const { | |
227 | //Product of conditional probability functions for pions | |
228 | Double_t rv=1.; | |
229 | for(Int_t i=0;i<8;i++){ | |
230 | Double_t fun=GetCondFunPi(i); | |
231 | if(fun>=0)rv*=fun; | |
232 | } | |
233 | return rv; | |
234 | } | |
235 | //_______________________________________________________________________ | |
e62c1aea | 236 | void AliITSPident::PrintParameters() const{ |
237 | //print parameters | |
238 | cout<<"___________________________\n"; | |
239 | cout<<"Track Local Momentum = "<<" "<<fMom<<endl; | |
240 | cout<<"Track dEdx = "<<" "<<fdEdx<<endl; | |
e62c1aea | 241 | } |
242 | ||
243 | //_______________________________________________________________________ | |
244 | Double_t AliITSPident::Langaufun(Double_t *x, Double_t *par) { | |
245 | ||
246 | //Fit parameters: | |
247 | //par[0]=Width (scale) parameter of Landau density | |
248 | //par[1]=Most Probable (MP, location) parameter of Landau density | |
249 | //par[2]=Total area (integral -inf to inf, normalization constant) | |
250 | //par[3]=Width (sigma) of convoluted Gaussian function | |
251 | // | |
252 | //In the Landau distribution (represented by the CERNLIB approximation), | |
253 | //the maximum is located at x=-0.22278298 with the location parameter=0. | |
254 | //This shift is corrected within this function, so that the actual | |
255 | //maximum is identical to the MP parameter. | |
256 | ||
257 | // Numeric constants | |
258 | Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) | |
259 | Double_t mpshift = -0.22278298; // Landau maximum location | |
260 | ||
261 | // Control constants | |
262 | Double_t np = 100.0; // number of convolution steps | |
263 | Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas | |
264 | ||
265 | // Variables | |
266 | Double_t xx; | |
267 | Double_t mpc; | |
268 | Double_t fland; | |
269 | Double_t sum = 0.0; | |
270 | Double_t xlow,xupp; | |
271 | Double_t step; | |
272 | Double_t i; | |
273 | ||
274 | ||
275 | // MP shift correction | |
276 | mpc = par[1] - mpshift * par[0]; | |
277 | ||
278 | // Range of convolution integral | |
279 | xlow = x[0] - sc * par[3]; | |
280 | xupp = x[0] + sc * par[3]; | |
281 | ||
282 | step = (xupp-xlow) / np; | |
283 | ||
284 | // Convolution integral of Landau and Gaussian by sum | |
285 | for(i=1.0; i<=np/2; i++) { | |
286 | xx = xlow + (i-.5) * step; | |
287 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
288 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
289 | ||
290 | xx = xupp - (i-.5) * step; | |
291 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
292 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
293 | } | |
294 | ||
295 | return (par[2] * step * sum * invsq2pi / par[3]); | |
296 | } | |
297 | ||
298 | //_______________________________________________________________________ | |
299 | Double_t AliITSPident::Langaufun2(Double_t *x, Double_t *par){ | |
300 | //normalized langaufun | |
301 | return 1/par[4]*Langaufun(x,par); | |
302 | } | |
303 | //_______________________________________________________________________ | |
304 | Double_t AliITSPident::Langaufunnorm(Double_t *x, Double_t *par){ | |
305 | //Fit parameters: | |
306 | //par[0]=Width (scale) parameter of Landau density | |
307 | //par[1]=Most Probable (MP, location) parameter of Landau density | |
308 | ||
309 | //par[2]=Width (sigma) of convoluted Gaussian function | |
310 | // | |
311 | //In the Landau distribution (represented by the CERNLIB approximation), | |
312 | //the maximum is located at x=-0.22278298 with the location parameter=0. | |
313 | //This shift is corrected within this function, so that the actual | |
314 | //maximum is identical to the MP parameter. | |
315 | ||
316 | // Numeric constants | |
317 | Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) | |
318 | Double_t mpshift = -0.22278298; // Landau maximum location | |
319 | ||
320 | // Control constants | |
321 | Double_t np = 100.0; // number of convolution steps | |
322 | Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas | |
323 | ||
324 | // Variables | |
325 | Double_t xx; | |
326 | Double_t mpc; | |
327 | Double_t fland; | |
328 | Double_t sum = 0.0; | |
329 | Double_t xlow,xupp; | |
330 | Double_t step; | |
331 | Double_t i; | |
332 | ||
333 | ||
334 | // MP shift correction | |
335 | mpc = par[1] - mpshift * par[0]; | |
336 | ||
337 | // Range of convolution integral | |
338 | xlow = x[0] - sc * par[2]; | |
339 | xupp = x[0] + sc * par[2]; | |
340 | ||
341 | step = (xupp-xlow) / np; | |
342 | ||
343 | // Convolution integral of Landau and Gaussian by sum | |
344 | for(i=1.0; i<=np/2; i++) { | |
345 | xx = xlow + (i-.5) * step; | |
346 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
347 | sum += fland * TMath::Gaus(x[0],xx,par[2]); | |
348 | ||
349 | xx = xupp - (i-.5) * step; | |
350 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
351 | sum += fland * TMath::Gaus(x[0],xx,par[2]); | |
352 | } | |
353 | ||
354 | return (step * sum * invsq2pi / par[2]); | |
355 | } | |
356 | //_______________________________________________________________________ | |
357 | Double_t AliITSPident::Gaus2(Double_t *x, Double_t *par){ | |
358 | //normalized gaussian function | |
359 | return 1/(sqrt(2*TMath::Pi())*par[1])*TMath::Gaus(x[0],par[0],par[1]); | |
360 | } | |
361 | //_______________________________________________________________________ | |
362 | void AliITSPident::CookFunItsLay(Int_t lay,Int_t opt,Double_t *par,Double_t dedx,Double_t mom,Double_t rangei,Double_t rangef,TString comment){ | |
363 | //it gives the response functions | |
364 | TF1 funLay(comment,Langaufunnorm,rangei,rangef,3); | |
365 | funLay.SetParameters(par); | |
366 | Double_t condFun=funLay.Eval(dedx); | |
367 | if(opt==0){ | |
368 | fCondFunProLay[lay]=condFun; | |
369 | if(mom<0.4 && dedx<100)fCondFunProLay[lay]=0.00001; | |
370 | if(mom<0.4 &&dedx<50)fCondFunProLay[lay]=0.0000001; | |
371 | } | |
372 | if(opt==1){ | |
373 | fCondFunKLay[lay]=condFun; | |
374 | if(mom<0.25 && dedx<100)fCondFunKLay[lay]=0.00001; | |
375 | if(mom<0.4 &&dedx<30)fCondFunKLay[lay]=0.0000001; | |
376 | } | |
377 | if(opt==2){ | |
378 | fCondFunPiLay[lay]=condFun; | |
379 | if(mom<0.6 &&dedx<20)fCondFunPiLay[lay]=0.001; | |
380 | } | |
381 | ||
382 | } | |
383 | //_______________________________________________________________________ | |
384 | Float_t AliITSPident::CookCombinedBayes(Double_t condfun[][3],Float_t *prior,Int_t part)const { | |
385 | //Bayesian combined PID in the ITS | |
386 | Int_t test=0; | |
387 | Float_t bayes; | |
388 | Float_t pprior[4]={0,0,0,0}; | |
389 | for(Int_t j=0;j<4;j++)pprior[j]=prior[j]; | |
390 | pprior[2]+=pprior[3];//prior for electrons summed to priors for pions | |
b68da4a2 | 391 | for(Int_t i=0;i<8;i++){//layer |
e62c1aea | 392 | if (condfun[i][0]>0 || condfun[i][1]>0 ||condfun[i][2]>0) test++; |
393 | } | |
b68da4a2 | 394 | if(test>0){ |
e62c1aea | 395 | if ((pprior[0]!=0 || pprior[1]!=0 ||pprior[2]!=0)&&CookSum(condfun,pprior)!=0){ |
396 | ||
397 | bayes=pprior[part]*CookProd(condfun,part)*1/CookSum(condfun,pprior); | |
398 | ||
399 | } | |
400 | else bayes=-100; | |
401 | } | |
402 | else bayes=-100; | |
403 | return bayes; | |
404 | } | |
405 | //_______________________________________________________________________ | |
406 | Float_t AliITSPident::CookProd(Double_t condfun[][3],Int_t part)const{ | |
407 | // | |
408 | Float_t p=1; | |
b68da4a2 | 409 | for(Int_t lay=0;lay<8;lay++){ |
410 | if(condfun[lay][part]>=0)p=p*condfun[lay][part]; | |
e62c1aea | 411 | } |
412 | ||
413 | return p; | |
414 | } | |
415 | //_______________________________________________________________________ | |
416 | Float_t AliITSPident::CookSum(Double_t condfun[][3],Float_t *prior)const{ | |
417 | // | |
418 | Float_t sum=0; | |
419 | Float_t pprior[4]={0,0,0,0}; | |
420 | for(Int_t j=0;j<4;j++)pprior[j]=prior[j]; | |
421 | pprior[2]+=pprior[3];//prior for electrons summed to priors for pions | |
422 | for(Int_t i=0;i<3;i++){//sum over the particles--electrons excluded | |
423 | sum+=pprior[i]*CookProd(condfun,i); | |
424 | } | |
425 | return sum; | |
426 | ||
427 | } |