ITS/AliITSPident.cxx

   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 "AliITStrackV2.h"
   9 #include "AliITSPident.h"
  10 #include "AliITSSteerPid.h"
  11 #include <Riostream.h>
  12
  13 ClassImp(AliITSPident)
  14   //_______________________________________________________________________
  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),
  24 fPPriorie(0)
  25 {
  26   // default constructor
  27   for (Int_t i=0;i<8;i++){
  28     fCondFunProLay[i]=0;
  29     fCondFunKLay[i]=0;
  30     fCondFunPiLay[i]=0;
  31   }
  32   for (Int_t i=0;i<4;i++)fNcls[i]=0;
  33 }
  34 //_______________________________________________________________________
  35 AliITSPident::~AliITSPident(){
  36   // destructor
  37 }
  38 //______________________________________________________________________
  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),
  48 fPPriorie(ob.fPPriorie)
  49 {
  50   // Copy constructor
  51 }
  52
  53 //______________________________________________________________________
  54 AliITSPident& AliITSPident::operator=(const AliITSPident&  ob){
  55   // Assignment operator
  56   this->~AliITSPident();
  57   new(this) AliITSPident(ob);
  58   return *this;
  59 }
  60
  61
  62 //_______________________________________________________________________
  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):
  64 fMom(mom),
  65 fdEdx(dEdx),
  66 fPBayesp(0),
  67 fPBayesk(0),
  68 fPBayespi(0),
  69 fPPriorip(priorip),
  70 fPPriorik(priorik),
  71 fPPrioripi(prioripi),
  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   }
  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];
  82     fNcls[la]=0;
  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];
  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     }
 104
 105   }
 106
 107   Float_t prior[4];Double_t condFun[8][3];
 108
 109   prior[0]=fPPriorip;
 110   prior[1]=fPPriorik;
 111   prior[2]=fPPrioripi;
 112   prior[3]=fPPriorie;
 113   for(Int_t la=0;la<8;la++){
 114     condFun[la][0]= fCondFunProLay[la];
 115     condFun[la][1]= fCondFunKLay[la];
 116     condFun[la][2]= fCondFunPiLay[la];
 117
 118   }
 119
 120   fPBayesp=CookCombinedBayes(condFun,prior,0);
 121   fPBayesk=CookCombinedBayes(condFun,prior,1);
 122   fPBayespi=CookCombinedBayes(condFun,prior,2);
 123 }
 124
 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):
 127 fMom(0),
 128 fdEdx(0),
 129 fPBayesp(0),
 130 fPBayesk(0),
 131 fPBayespi(0),
 132 fPPriorip(priorip),
 133 fPPriorik(priorik),
 134 fPPrioripi(prioripi),
 135 fPPriorie(priorie)
 136 {
 137   //
 138   for (Int_t i=0;i<8;i++){
 139     fCondFunProLay[i]=-1;
 140     fCondFunKLay[i]=-1;
 141     fCondFunPiLay[i]=-1;
 142   }
 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
 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];
 154     fNcls[la]=0;
 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];
 161
 162     range[4]=0.3*parpi[1];
 163     range[5]=2.*parpi[1];
 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
 178   }
 179
 180   Float_t prior[4];Double_t condFun[8][3];
 181
 182   prior[0]=fPPriorip;
 183   prior[1]=fPPriorik;
 184   prior[2]=fPPrioripi;
 185   prior[3]=fPPriorie;
 186   for(Int_t la=0;la<8;la++){
 187     condFun[la][0]= fCondFunProLay[la];
 188     condFun[la][1]= fCondFunKLay[la];
 189     condFun[la][2]= fCondFunPiLay[la];
 190
 191   }
 192
 193   fPBayesp=CookCombinedBayes(condFun,prior,0);
 194   fPBayesk=CookCombinedBayes(condFun,prior,1);
 195   fPBayespi=CookCombinedBayes(condFun,prior,2);
 196   fdEdx=trackITS->GetdEdx();
 197
 198 }
 199 //_______________________________________________________________________
 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 //_______________________________________________________________________
 236 void AliITSPident::PrintParameters() const{
 237  //print parameters
 238   cout<<"___________________________\n";
 239   cout<<"Track Local Momentum = "<<"  "<<fMom<<endl;
 240   cout<<"Track dEdx = "<<"  "<<fdEdx<<endl;
 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
 391   for(Int_t i=0;i<8;i++){//layer
 392     if (condfun[i][0]>0 || condfun[i][1]>0 ||condfun[i][2]>0) test++;
 393   }
 394   if(test>0){
 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;
 409   for(Int_t lay=0;lay<8;lay++){
 410     if(condfun[lay][part]>=0)p=p*condfun[lay][part];
 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 }