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"
10 ClassImp(AliITSPident)
11 //_______________________________________________________________________
12 AliITSPident::AliITSPident():
24 // default constructor
25 for (Int_t i=0;i<4;i++){
31 //_______________________________________________________________________
32 AliITSPident::~AliITSPident(){
35 //______________________________________________________________________
36 AliITSPident::AliITSPident(const AliITSPident &ob) :TObject(ob),
39 fPBayesp(ob.fPBayesp),
40 fPBayesk(ob.fPBayesk),
41 fPBayespi(ob.fPBayespi),
42 fPPriorip(ob.fPPriorip),
43 fPPriorik(ob.fPPriorik),
44 fPPrioripi(ob.fPPrioripi),
45 fPPriorie(ob.fPPriorie),
51 //______________________________________________________________________
52 AliITSPident& AliITSPident::operator=(const AliITSPident& ob){
53 // Assignment operator
54 this->~AliITSPident();
55 new(this) AliITSPident(ob);
60 //_______________________________________________________________________
61 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):
75 for(Int_t la=0;la<4;la++){//loop on layers
76 Double_t parp[3];Double_t park[3];Double_t parpi[3];
77 sp->GetParFitLayer(la,fMom,parp,park,parpi);
86 range[4]=0.3*parpi[1];
88 CookFunItsLay(la,0,parp,Qlay[la],fMom,range[0],range[1],"fPro");
89 CookFunItsLay(la,1,park,Qlay[la],fMom,range[2],range[3],"fKao");
90 CookFunItsLay(la,2,parpi,Qlay[la],fMom,range[4],range[5],"fPi");
94 Float_t prior[4];Double_t condFun[4][3];
100 for(Int_t la=0;la<4;la++){
101 condFun[la][0]= fCondFunProLay[la];
102 condFun[la][1]= fCondFunKLay[la];
103 condFun[la][2]= fCondFunPiLay[la];
107 fPBayesp=CookCombinedBayes(condFun,prior,0);
108 fPBayesk=CookCombinedBayes(condFun,prior,1);
109 fPBayespi=CookCombinedBayes(condFun,prior,2);
113 //_______________________________________________________________________
114 AliITSPident::AliITSPident(AliITStrackV2 *trackITS,AliITSSteerPid *sp,Float_t *Qlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie):
122 fPPrioripi(prioripi),
129 trackITS->GetExternalParameters(xr,par);
131 Float_t lamb=TMath::ATan(par[3]);
132 fMom=1/(TMath::Abs(par[4])*TMath::Cos(lamb));
137 for(Int_t la=0;la<4;la++){//loop on layers
138 Double_t parp[3];Double_t park[3];Double_t parpi[3];
139 sp->GetParFitLayer(la,fMom,parp,park,parpi);
141 range[0]=0.3*parp[1];
144 range[2]=0.3*park[1];
147 range[4]=0.3*parpi[1];
148 range[5]=2.*parpi[1];
149 CookFunItsLay(la,0,parp,Qlay[la],fMom,range[0],range[1],"fPro");
150 CookFunItsLay(la,1,park,Qlay[la],fMom,range[2],range[3],"fKao");
151 CookFunItsLay(la,2,parpi,Qlay[la],fMom,range[4],range[5],"fPi");
154 Float_t prior[4];Double_t condFun[4][3];
160 for(Int_t la=0;la<4;la++){
161 condFun[la][0]= fCondFunProLay[la];
162 condFun[la][1]= fCondFunKLay[la];
163 condFun[la][2]= fCondFunPiLay[la];
167 fPBayesp=CookCombinedBayes(condFun,prior,0);
168 fPBayesk=CookCombinedBayes(condFun,prior,1);
169 fPBayespi=CookCombinedBayes(condFun,prior,2);
170 fdEdx=trackITS->GetdEdx();
173 //_______________________________________________________________________
174 void AliITSPident::PrintParameters() const{
176 cout<<"___________________________\n";
177 cout<<"Track Local Momentum = "<<" "<<fMom<<endl;
178 cout<<"Track dEdx = "<<" "<<fdEdx<<endl;
179 cout<<"Inv Pt = "<<fInvPt<<endl;
182 //_______________________________________________________________________
183 Double_t AliITSPident::Langaufun(Double_t *x, Double_t *par) {
186 //par[0]=Width (scale) parameter of Landau density
187 //par[1]=Most Probable (MP, location) parameter of Landau density
188 //par[2]=Total area (integral -inf to inf, normalization constant)
189 //par[3]=Width (sigma) of convoluted Gaussian function
191 //In the Landau distribution (represented by the CERNLIB approximation),
192 //the maximum is located at x=-0.22278298 with the location parameter=0.
193 //This shift is corrected within this function, so that the actual
194 //maximum is identical to the MP parameter.
197 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
198 Double_t mpshift = -0.22278298; // Landau maximum location
201 Double_t np = 100.0; // number of convolution steps
202 Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
214 // MP shift correction
215 mpc = par[1] - mpshift * par[0];
217 // Range of convolution integral
218 xlow = x[0] - sc * par[3];
219 xupp = x[0] + sc * par[3];
221 step = (xupp-xlow) / np;
223 // Convolution integral of Landau and Gaussian by sum
224 for(i=1.0; i<=np/2; i++) {
225 xx = xlow + (i-.5) * step;
226 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
227 sum += fland * TMath::Gaus(x[0],xx,par[3]);
229 xx = xupp - (i-.5) * step;
230 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
231 sum += fland * TMath::Gaus(x[0],xx,par[3]);
234 return (par[2] * step * sum * invsq2pi / par[3]);
237 //_______________________________________________________________________
238 Double_t AliITSPident::Langaufun2(Double_t *x, Double_t *par){
239 //normalized langaufun
240 return 1/par[4]*Langaufun(x,par);
242 //_______________________________________________________________________
243 Double_t AliITSPident::Langaufunnorm(Double_t *x, Double_t *par){
245 //par[0]=Width (scale) parameter of Landau density
246 //par[1]=Most Probable (MP, location) parameter of Landau density
248 //par[2]=Width (sigma) of convoluted Gaussian function
250 //In the Landau distribution (represented by the CERNLIB approximation),
251 //the maximum is located at x=-0.22278298 with the location parameter=0.
252 //This shift is corrected within this function, so that the actual
253 //maximum is identical to the MP parameter.
256 Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
257 Double_t mpshift = -0.22278298; // Landau maximum location
260 Double_t np = 100.0; // number of convolution steps
261 Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
273 // MP shift correction
274 mpc = par[1] - mpshift * par[0];
276 // Range of convolution integral
277 xlow = x[0] - sc * par[2];
278 xupp = x[0] + sc * par[2];
280 step = (xupp-xlow) / np;
282 // Convolution integral of Landau and Gaussian by sum
283 for(i=1.0; i<=np/2; i++) {
284 xx = xlow + (i-.5) * step;
285 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
286 sum += fland * TMath::Gaus(x[0],xx,par[2]);
288 xx = xupp - (i-.5) * step;
289 fland = TMath::Landau(xx,mpc,par[0]) / par[0];
290 sum += fland * TMath::Gaus(x[0],xx,par[2]);
293 return (step * sum * invsq2pi / par[2]);
295 //_______________________________________________________________________
296 Double_t AliITSPident::Gaus2(Double_t *x, Double_t *par){
297 //normalized gaussian function
298 return 1/(sqrt(2*TMath::Pi())*par[1])*TMath::Gaus(x[0],par[0],par[1]);
300 //_______________________________________________________________________
301 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){
302 //it gives the response functions
303 TF1 funLay(comment,Langaufunnorm,rangei,rangef,3);
304 funLay.SetParameters(par);
305 Double_t condFun=funLay.Eval(dedx);
307 fCondFunProLay[lay]=condFun;
308 if(mom<0.4 && dedx<100)fCondFunProLay[lay]=0.00001;
309 if(mom<0.4 &&dedx<50)fCondFunProLay[lay]=0.0000001;
312 fCondFunKLay[lay]=condFun;
313 if(mom<0.25 && dedx<100)fCondFunKLay[lay]=0.00001;
314 if(mom<0.4 &&dedx<30)fCondFunKLay[lay]=0.0000001;
317 fCondFunPiLay[lay]=condFun;
318 if(mom<0.6 &&dedx<20)fCondFunPiLay[lay]=0.001;
322 //_______________________________________________________________________
323 Float_t AliITSPident::CookCombinedBayes(Double_t condfun[][3],Float_t *prior,Int_t part)const {
324 //Bayesian combined PID in the ITS
327 Float_t pprior[4]={0,0,0,0};
328 for(Int_t j=0;j<4;j++)pprior[j]=prior[j];
329 pprior[2]+=pprior[3];//prior for electrons summed to priors for pions
330 for(Int_t i=0;i<4;i++){//layer
331 if (condfun[i][0]>0 || condfun[i][1]>0 ||condfun[i][2]>0) test++;
334 if ((pprior[0]!=0 || pprior[1]!=0 ||pprior[2]!=0)&&CookSum(condfun,pprior)!=0){
336 bayes=pprior[part]*CookProd(condfun,part)*1/CookSum(condfun,pprior);
344 //_______________________________________________________________________
345 Float_t AliITSPident::CookProd(Double_t condfun[][3],Int_t part)const{
348 for(Int_t lay=0;lay<4;lay++){
350 p=p*condfun[lay][part];
355 //_______________________________________________________________________
356 Float_t AliITSPident::CookSum(Double_t condfun[][3],Float_t *prior)const{
359 Float_t pprior[4]={0,0,0,0};
360 for(Int_t j=0;j<4;j++)pprior[j]=prior[j];
361 pprior[2]+=pprior[3];//prior for electrons summed to priors for pions
362 for(Int_t i=0;i<3;i++){//sum over the particles--electrons excluded
363 sum+=pprior[i]*CookProd(condfun,i);