Added two missing includes to allow macro compilation (thanks to Laurent for remarkin...
[u/mrichter/AliRoot.git] / ITS / AliITSPident.cxx
CommitLineData
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
13ClassImp(AliITSPident)
14 //_______________________________________________________________________
94631b2f 15AliITSPident::AliITSPident():
16fMom(0),
17fdEdx(0),
18fPBayesp(0),
19fPBayesk(0),
20fPBayespi(0),
21fPPriorip(0),
22fPPriorik(0),
23fPPrioripi(0),
b68da4a2 24fPPriorie(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//_______________________________________________________________________
35AliITSPident::~AliITSPident(){
36 // destructor
37}
38//______________________________________________________________________
94631b2f 39AliITSPident::AliITSPident(const AliITSPident &ob) :TObject(ob),
40fMom(ob.fMom),
41fdEdx(ob.fdEdx),
42fPBayesp(ob.fPBayesp),
43fPBayesk(ob.fPBayesk),
44fPBayespi(ob.fPBayespi),
45fPPriorip(ob.fPPriorip),
46fPPriorik(ob.fPPriorik),
47fPPrioripi(ob.fPPrioripi),
b68da4a2 48fPPriorie(ob.fPPriorie)
94631b2f 49{
e62c1aea 50 // Copy constructor
e62c1aea 51}
52
53//______________________________________________________________________
94631b2f 54AliITSPident& 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 63AliITSPident::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 64fMom(mom),
65fdEdx(dEdx),
66fPBayesp(0),
67fPBayesk(0),
68fPBayespi(0),
69fPPriorip(priorip),
70fPPriorik(priorik),
71fPPrioripi(prioripi),
b68da4a2 72fPPriorie(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//__________________________________________________________________________________________
126AliITSPident::AliITSPident(AliITStrackV2 *trackITS,AliITSSteerPid *sp,Float_t *Qlay,Float_t *nlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie):
94631b2f 127fMom(0),
128fdEdx(0),
129fPBayesp(0),
130fPBayesk(0),
131fPBayespi(0),
132fPPriorip(priorip),
133fPPriorik(priorik),
134fPPrioripi(prioripi),
b68da4a2 135fPPriorie(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 200void 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}//_______________________________________________________________________
207Double_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}//_______________________________________________________________________
216Double_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//_______________________________________________________________________
226Double_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 236void 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//_______________________________________________________________________
244Double_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//_______________________________________________________________________
299Double_t AliITSPident::Langaufun2(Double_t *x, Double_t *par){
300 //normalized langaufun
301 return 1/par[4]*Langaufun(x,par);
302}
303//_______________________________________________________________________
304Double_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//_______________________________________________________________________
357Double_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//_______________________________________________________________________
362void 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//_______________________________________________________________________
384Float_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//_______________________________________________________________________
406Float_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//_______________________________________________________________________
416Float_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}