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