]>
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 | ////////////////////////////////////////////////////////////////////////// | |
8 | #include "AliITSPident.h" | |
9 | ||
10 | ClassImp(AliITSPident) | |
11 | //_______________________________________________________________________ | |
12 | AliITSPident::AliITSPident(){ | |
13 | // default constructor | |
14 | fMom=0; | |
15 | fdEdx=0; | |
16 | fInvPt=0; | |
17 | fPBayesp=0; | |
18 | fPBayesk=0; | |
19 | fPBayespi=0; | |
20 | fPPriorip=0; | |
21 | fPPriorik=0; | |
22 | fPPrioripi=0; | |
23 | fPPriorie=0; | |
24 | for (Int_t i=0;i<4;i++){ | |
25 | fCondFunProLay[i]=0; | |
26 | fCondFunKLay[i]=0; | |
27 | fCondFunPiLay[i]=0; | |
28 | } | |
29 | } | |
30 | //_______________________________________________________________________ | |
31 | AliITSPident::~AliITSPident(){ | |
32 | // destructor | |
33 | } | |
34 | //______________________________________________________________________ | |
35 | AliITSPident::AliITSPident(const AliITSPident &ob) :TObject(ob) { | |
36 | // Copy constructor | |
37 | // Copies are not allowed. The method is protected to avoid misuse. | |
38 | Error("AliITSPident","Copy constructor not allowed\n"); | |
39 | } | |
40 | ||
41 | //______________________________________________________________________ | |
42 | AliITSPident& AliITSPident::operator=(const AliITSPident& /* ob */){ | |
43 | // Assignment operator | |
44 | // Assignment is not allowed. The method is protected to avoid misuse. | |
45 | Error("= operator","Assignment operator not allowed\n"); | |
46 | return *this; | |
47 | } | |
48 | ||
49 | ||
50 | //_______________________________________________________________________ | |
51 | AliITSPident::AliITSPident(Double_t mom,Double_t invPt,Double_t dEdx,AliITSSteerPid *sp,Float_t *Qlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie){ | |
52 | ||
53 | //test | |
54 | ||
55 | fMom=mom; | |
56 | fInvPt=invPt; | |
57 | for(Int_t la=0;la<4;la++){//loop on layers | |
58 | Double_t parp[3];Double_t park[3];Double_t parpi[3]; | |
59 | sp->GetParFitLayer(la,fMom,parp,park,parpi); | |
149c4c59 | 60 | |
e62c1aea | 61 | Double_t range[6]; |
62 | range[0]=0.3*parp[1]; | |
63 | range[1]=2.*parp[1]; | |
64 | ||
65 | range[2]=0.3*park[1]; | |
66 | range[3]=2.*park[1]; | |
67 | ||
68 | range[4]=0.3*parpi[1]; | |
69 | range[5]=2.*parpi[1]; | |
70 | CookFunItsLay(la,0,parp,Qlay[la],fMom,range[0],range[1],"fPro"); | |
71 | CookFunItsLay(la,1,park,Qlay[la],fMom,range[2],range[3],"fKao"); | |
72 | CookFunItsLay(la,2,parpi,Qlay[la],fMom,range[4],range[5],"fPi"); | |
73 | ||
74 | } | |
75 | fPPriorip=priorip; | |
76 | fPPriorik=priorik; | |
77 | fPPrioripi=prioripi; | |
78 | fPPriorie=priorie; | |
79 | ||
80 | ||
81 | Float_t prior[4];Double_t condFun[4][3]; | |
82 | ||
83 | prior[0]=fPPriorip; | |
84 | prior[1]=fPPriorik; | |
85 | prior[2]=fPPrioripi; | |
86 | prior[3]=fPPriorie; | |
87 | for(Int_t la=0;la<4;la++){ | |
88 | condFun[la][0]= fCondFunProLay[la]; | |
89 | condFun[la][1]= fCondFunKLay[la]; | |
90 | condFun[la][2]= fCondFunPiLay[la]; | |
91 | ||
92 | } | |
93 | ||
e62c1aea | 94 | fPBayesp=CookCombinedBayes(condFun,prior,0); |
e62c1aea | 95 | fPBayesk=CookCombinedBayes(condFun,prior,1); |
e62c1aea | 96 | fPBayespi=CookCombinedBayes(condFun,prior,2); |
e62c1aea | 97 | fdEdx=dEdx; |
98 | } | |
99 | ||
100 | ||
101 | //_______________________________________________________________________ | |
102 | AliITSPident::AliITSPident(AliITStrackV2 *trackITS,AliITSSteerPid *sp,Float_t *Qlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie){ | |
103 | // | |
104 | Double_t xr; | |
105 | Double_t par[5]; | |
106 | trackITS->GetExternalParameters(xr,par); | |
107 | if (par[4]!=0) { | |
108 | Float_t lamb=TMath::ATan(par[3]); | |
109 | fMom=1/(TMath::Abs(par[4])*TMath::Cos(lamb)); | |
110 | } | |
111 | ||
112 | fInvPt=par[4]; | |
113 | ||
114 | for(Int_t la=0;la<4;la++){//loop on layers | |
115 | Double_t parp[3];Double_t park[3];Double_t parpi[3]; | |
116 | sp->GetParFitLayer(la,fMom,parp,park,parpi); | |
117 | Double_t range[8]; | |
118 | range[0]=0.3*parp[1]; | |
119 | range[1]=2.*parp[1]; | |
120 | ||
121 | range[2]=0.3*park[1]; | |
122 | range[3]=2.*park[1]; | |
123 | ||
124 | range[4]=0.3*parpi[1]; | |
125 | range[5]=2.*parpi[1]; | |
126 | CookFunItsLay(la,0,parp,Qlay[la],fMom,range[0],range[1],"fPro"); | |
127 | CookFunItsLay(la,1,park,Qlay[la],fMom,range[2],range[3],"fKao"); | |
128 | CookFunItsLay(la,2,parpi,Qlay[la],fMom,range[4],range[5],"fPi"); | |
129 | } | |
130 | fPPriorip=priorip; | |
131 | fPPriorik=priorik; | |
132 | fPPrioripi=prioripi; | |
133 | fPPriorie=priorie; | |
134 | ||
135 | ||
136 | Float_t prior[4];Double_t condFun[4][3]; | |
137 | ||
138 | prior[0]=fPPriorip; | |
139 | prior[1]=fPPriorik; | |
140 | prior[2]=fPPrioripi; | |
141 | prior[3]=fPPriorie; | |
142 | for(Int_t la=0;la<4;la++){ | |
143 | condFun[la][0]= fCondFunProLay[la]; | |
144 | condFun[la][1]= fCondFunKLay[la]; | |
145 | condFun[la][2]= fCondFunPiLay[la]; | |
146 | ||
147 | } | |
148 | ||
e62c1aea | 149 | fPBayesp=CookCombinedBayes(condFun,prior,0); |
e62c1aea | 150 | fPBayesk=CookCombinedBayes(condFun,prior,1); |
e62c1aea | 151 | fPBayespi=CookCombinedBayes(condFun,prior,2); |
e62c1aea | 152 | fdEdx=trackITS->GetdEdx(); |
153 | ||
154 | } | |
155 | //_______________________________________________________________________ | |
156 | void AliITSPident::PrintParameters() const{ | |
157 | //print parameters | |
158 | cout<<"___________________________\n"; | |
159 | cout<<"Track Local Momentum = "<<" "<<fMom<<endl; | |
160 | cout<<"Track dEdx = "<<" "<<fdEdx<<endl; | |
161 | cout<<"Inv Pt = "<<fInvPt<<endl; | |
162 | } | |
163 | ||
164 | //_______________________________________________________________________ | |
165 | Double_t AliITSPident::Langaufun(Double_t *x, Double_t *par) { | |
166 | ||
167 | //Fit parameters: | |
168 | //par[0]=Width (scale) parameter of Landau density | |
169 | //par[1]=Most Probable (MP, location) parameter of Landau density | |
170 | //par[2]=Total area (integral -inf to inf, normalization constant) | |
171 | //par[3]=Width (sigma) of convoluted Gaussian function | |
172 | // | |
173 | //In the Landau distribution (represented by the CERNLIB approximation), | |
174 | //the maximum is located at x=-0.22278298 with the location parameter=0. | |
175 | //This shift is corrected within this function, so that the actual | |
176 | //maximum is identical to the MP parameter. | |
177 | ||
178 | // Numeric constants | |
179 | Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) | |
180 | Double_t mpshift = -0.22278298; // Landau maximum location | |
181 | ||
182 | // Control constants | |
183 | Double_t np = 100.0; // number of convolution steps | |
184 | Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas | |
185 | ||
186 | // Variables | |
187 | Double_t xx; | |
188 | Double_t mpc; | |
189 | Double_t fland; | |
190 | Double_t sum = 0.0; | |
191 | Double_t xlow,xupp; | |
192 | Double_t step; | |
193 | Double_t i; | |
194 | ||
195 | ||
196 | // MP shift correction | |
197 | mpc = par[1] - mpshift * par[0]; | |
198 | ||
199 | // Range of convolution integral | |
200 | xlow = x[0] - sc * par[3]; | |
201 | xupp = x[0] + sc * par[3]; | |
202 | ||
203 | step = (xupp-xlow) / np; | |
204 | ||
205 | // Convolution integral of Landau and Gaussian by sum | |
206 | for(i=1.0; i<=np/2; i++) { | |
207 | xx = xlow + (i-.5) * step; | |
208 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
209 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
210 | ||
211 | xx = xupp - (i-.5) * step; | |
212 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
213 | sum += fland * TMath::Gaus(x[0],xx,par[3]); | |
214 | } | |
215 | ||
216 | return (par[2] * step * sum * invsq2pi / par[3]); | |
217 | } | |
218 | ||
219 | //_______________________________________________________________________ | |
220 | Double_t AliITSPident::Langaufun2(Double_t *x, Double_t *par){ | |
221 | //normalized langaufun | |
222 | return 1/par[4]*Langaufun(x,par); | |
223 | } | |
224 | //_______________________________________________________________________ | |
225 | Double_t AliITSPident::Langaufunnorm(Double_t *x, Double_t *par){ | |
226 | //Fit parameters: | |
227 | //par[0]=Width (scale) parameter of Landau density | |
228 | //par[1]=Most Probable (MP, location) parameter of Landau density | |
229 | ||
230 | //par[2]=Width (sigma) of convoluted Gaussian function | |
231 | // | |
232 | //In the Landau distribution (represented by the CERNLIB approximation), | |
233 | //the maximum is located at x=-0.22278298 with the location parameter=0. | |
234 | //This shift is corrected within this function, so that the actual | |
235 | //maximum is identical to the MP parameter. | |
236 | ||
237 | // Numeric constants | |
238 | Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2) | |
239 | Double_t mpshift = -0.22278298; // Landau maximum location | |
240 | ||
241 | // Control constants | |
242 | Double_t np = 100.0; // number of convolution steps | |
243 | Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas | |
244 | ||
245 | // Variables | |
246 | Double_t xx; | |
247 | Double_t mpc; | |
248 | Double_t fland; | |
249 | Double_t sum = 0.0; | |
250 | Double_t xlow,xupp; | |
251 | Double_t step; | |
252 | Double_t i; | |
253 | ||
254 | ||
255 | // MP shift correction | |
256 | mpc = par[1] - mpshift * par[0]; | |
257 | ||
258 | // Range of convolution integral | |
259 | xlow = x[0] - sc * par[2]; | |
260 | xupp = x[0] + sc * par[2]; | |
261 | ||
262 | step = (xupp-xlow) / np; | |
263 | ||
264 | // Convolution integral of Landau and Gaussian by sum | |
265 | for(i=1.0; i<=np/2; i++) { | |
266 | xx = xlow + (i-.5) * step; | |
267 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
268 | sum += fland * TMath::Gaus(x[0],xx,par[2]); | |
269 | ||
270 | xx = xupp - (i-.5) * step; | |
271 | fland = TMath::Landau(xx,mpc,par[0]) / par[0]; | |
272 | sum += fland * TMath::Gaus(x[0],xx,par[2]); | |
273 | } | |
274 | ||
275 | return (step * sum * invsq2pi / par[2]); | |
276 | } | |
277 | //_______________________________________________________________________ | |
278 | Double_t AliITSPident::Gaus2(Double_t *x, Double_t *par){ | |
279 | //normalized gaussian function | |
280 | return 1/(sqrt(2*TMath::Pi())*par[1])*TMath::Gaus(x[0],par[0],par[1]); | |
281 | } | |
282 | //_______________________________________________________________________ | |
283 | 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){ | |
284 | //it gives the response functions | |
285 | TF1 funLay(comment,Langaufunnorm,rangei,rangef,3); | |
286 | funLay.SetParameters(par); | |
287 | Double_t condFun=funLay.Eval(dedx); | |
288 | if(opt==0){ | |
289 | fCondFunProLay[lay]=condFun; | |
290 | if(mom<0.4 && dedx<100)fCondFunProLay[lay]=0.00001; | |
291 | if(mom<0.4 &&dedx<50)fCondFunProLay[lay]=0.0000001; | |
292 | } | |
293 | if(opt==1){ | |
294 | fCondFunKLay[lay]=condFun; | |
295 | if(mom<0.25 && dedx<100)fCondFunKLay[lay]=0.00001; | |
296 | if(mom<0.4 &&dedx<30)fCondFunKLay[lay]=0.0000001; | |
297 | } | |
298 | if(opt==2){ | |
299 | fCondFunPiLay[lay]=condFun; | |
300 | if(mom<0.6 &&dedx<20)fCondFunPiLay[lay]=0.001; | |
301 | } | |
302 | ||
303 | } | |
304 | //_______________________________________________________________________ | |
305 | Float_t AliITSPident::CookCombinedBayes(Double_t condfun[][3],Float_t *prior,Int_t part)const { | |
306 | //Bayesian combined PID in the ITS | |
307 | Int_t test=0; | |
308 | Float_t bayes; | |
309 | Float_t pprior[4]={0,0,0,0}; | |
310 | for(Int_t j=0;j<4;j++)pprior[j]=prior[j]; | |
311 | pprior[2]+=pprior[3];//prior for electrons summed to priors for pions | |
312 | for(Int_t i=0;i<4;i++){//layer | |
313 | if (condfun[i][0]>0 || condfun[i][1]>0 ||condfun[i][2]>0) test++; | |
314 | } | |
315 | if(test==4){ | |
316 | if ((pprior[0]!=0 || pprior[1]!=0 ||pprior[2]!=0)&&CookSum(condfun,pprior)!=0){ | |
317 | ||
318 | bayes=pprior[part]*CookProd(condfun,part)*1/CookSum(condfun,pprior); | |
319 | ||
320 | } | |
321 | else bayes=-100; | |
322 | } | |
323 | else bayes=-100; | |
324 | return bayes; | |
325 | } | |
326 | //_______________________________________________________________________ | |
327 | Float_t AliITSPident::CookProd(Double_t condfun[][3],Int_t part)const{ | |
328 | // | |
329 | Float_t p=1; | |
330 | for(Int_t lay=0;lay<4;lay++){ | |
331 | ||
332 | p=p*condfun[lay][part]; | |
333 | } | |
334 | ||
335 | return p; | |
336 | } | |
337 | //_______________________________________________________________________ | |
338 | Float_t AliITSPident::CookSum(Double_t condfun[][3],Float_t *prior)const{ | |
339 | // | |
340 | Float_t sum=0; | |
341 | Float_t pprior[4]={0,0,0,0}; | |
342 | for(Int_t j=0;j<4;j++)pprior[j]=prior[j]; | |
343 | pprior[2]+=pprior[3];//prior for electrons summed to priors for pions | |
344 | for(Int_t i=0;i<3;i++){//sum over the particles--electrons excluded | |
345 | sum+=pprior[i]*CookProd(condfun,i); | |
346 | } | |
347 | return sum; | |
348 | ||
349 | } |