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 "AliITSPident.h"
   9
  10 ClassImp(AliITSPident)
  11   //_______________________________________________________________________
  12 AliITSPident::AliITSPident():
  13 fMom(0),
  14 fdEdx(0),
  15 fPBayesp(0),
  16 fPBayesk(0),
  17 fPBayespi(0),
  18 fPPriorip(0),
  19 fPPriorik(0),
  20 fPPrioripi(0),
  21 fPPriorie(0),
  22 fInvPt(0)
  23 {
  24   // default constructor
  25   for (Int_t i=0;i<4;i++){
  26     fCondFunProLay[i]=0;
  27     fCondFunKLay[i]=0;
  28     fCondFunPiLay[i]=0;
  29   }
  30 }
  31 //_______________________________________________________________________
  32 AliITSPident::~AliITSPident(){
  33   // destructor
  34 }
  35 //______________________________________________________________________
  36 AliITSPident::AliITSPident(const AliITSPident &ob) :TObject(ob),
  37 fMom(ob.fMom),
  38 fdEdx(ob.fdEdx),
  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),
  46 fInvPt(ob.fInvPt)
  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,Double_t invPt,Double_t dEdx,AliITSSteerPid *sp,Float_t *Qlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie):
  62 fMom(mom),
  63 fdEdx(dEdx),
  64 fPBayesp(0),
  65 fPBayesk(0),
  66 fPBayespi(0),
  67 fPPriorip(priorip),
  68 fPPriorik(priorik),
  69 fPPrioripi(prioripi),
  70 fPPriorie(priorie),
  71 fInvPt(invPt){
  72
  73   //test
  74
  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);
  78
  79     Double_t range[6];
  80     range[0]=0.3*parp[1];
  81     range[1]=2.*parp[1];
  82
  83     range[2]=0.3*park[1];
  84     range[3]=2.*park[1];
  85
  86     range[4]=0.3*parpi[1];
  87     range[5]=2.*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");
  91
  92   }
  93
  94   Float_t prior[4];Double_t condFun[4][3];
  95
  96   prior[0]=fPPriorip;
  97   prior[1]=fPPriorik;
  98   prior[2]=fPPrioripi;
  99   prior[3]=fPPriorie;
 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];
 104
 105   }
 106
 107   fPBayesp=CookCombinedBayes(condFun,prior,0);
 108   fPBayesk=CookCombinedBayes(condFun,prior,1);
 109   fPBayespi=CookCombinedBayes(condFun,prior,2);
 110  }
 111
 112
 113 //_______________________________________________________________________
 114 AliITSPident::AliITSPident(AliITStrackV2 *trackITS,AliITSSteerPid *sp,Float_t *Qlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie):
 115 fMom(0),
 116 fdEdx(0),
 117 fPBayesp(0),
 118 fPBayesk(0),
 119 fPBayespi(0),
 120 fPPriorip(priorip),
 121 fPPriorik(priorik),
 122 fPPrioripi(prioripi),
 123 fPPriorie(priorie),
 124 fInvPt(0)
 125 {
 126   //
 127   Double_t xr;
 128   Double_t par[5];
 129   trackITS->GetExternalParameters(xr,par);
 130   if (par[4]!=0) {
 131     Float_t lamb=TMath::ATan(par[3]);
 132     fMom=1/(TMath::Abs(par[4])*TMath::Cos(lamb));
 133   }
 134
 135   fInvPt=par[4];
 136
 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);
 140     Double_t range[8];
 141     range[0]=0.3*parp[1];
 142     range[1]=2.*parp[1];
 143
 144     range[2]=0.3*park[1];
 145     range[3]=2.*park[1];
 146
 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");
 152   }
 153
 154   Float_t prior[4];Double_t condFun[4][3];
 155
 156   prior[0]=fPPriorip;
 157   prior[1]=fPPriorik;
 158   prior[2]=fPPrioripi;
 159   prior[3]=fPPriorie;
 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];
 164
 165   }
 166
 167   fPBayesp=CookCombinedBayes(condFun,prior,0);
 168   fPBayesk=CookCombinedBayes(condFun,prior,1);
 169   fPBayespi=CookCombinedBayes(condFun,prior,2);
 170   fdEdx=trackITS->GetdEdx();
 171
 172 }
 173 //_______________________________________________________________________
 174 void AliITSPident::PrintParameters() const{
 175  //print parameters
 176   cout<<"___________________________\n";
 177   cout<<"Track Local Momentum = "<<"  "<<fMom<<endl;
 178   cout<<"Track dEdx = "<<"  "<<fdEdx<<endl;
 179   cout<<"Inv Pt = "<<fInvPt<<endl;
 180 }
 181
 182 //_______________________________________________________________________
 183 Double_t AliITSPident::Langaufun(Double_t *x, Double_t *par) {
 184
 185   //Fit parameters:
 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
 190   //
 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.
 195
 196   // Numeric constants
 197   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 198   Double_t mpshift  = -0.22278298;       // Landau maximum location
 199
 200   // Control constants
 201   Double_t np = 100.0;      // number of convolution steps
 202   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 203
 204   // Variables
 205   Double_t xx;
 206   Double_t mpc;
 207   Double_t fland;
 208   Double_t sum = 0.0;
 209   Double_t xlow,xupp;
 210   Double_t step;
 211   Double_t i;
 212
 213
 214   // MP shift correction
 215   mpc = par[1] - mpshift * par[0];
 216
 217   // Range of convolution integral
 218   xlow = x[0] - sc * par[3];
 219   xupp = x[0] + sc * par[3];
 220
 221   step = (xupp-xlow) / np;
 222
 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]);
 228
 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]);
 232   }
 233
 234   return (par[2] * step * sum * invsq2pi / par[3]);
 235 }
 236
 237 //_______________________________________________________________________
 238 Double_t AliITSPident::Langaufun2(Double_t *x, Double_t *par){
 239   //normalized langaufun
 240   return 1/par[4]*Langaufun(x,par);
 241 }
 242 //_______________________________________________________________________
 243 Double_t AliITSPident::Langaufunnorm(Double_t *x, Double_t *par){
 244    //Fit parameters:
 245   //par[0]=Width (scale) parameter of Landau density
 246   //par[1]=Most Probable (MP, location) parameter of Landau density
 247
 248   //par[2]=Width (sigma) of convoluted Gaussian function
 249   //
 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.
 254
 255   // Numeric constants
 256   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
 257   Double_t mpshift  = -0.22278298;       // Landau maximum location
 258
 259   // Control constants
 260   Double_t np = 100.0;      // number of convolution steps
 261   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
 262
 263   // Variables
 264   Double_t xx;
 265   Double_t mpc;
 266   Double_t fland;
 267   Double_t sum = 0.0;
 268   Double_t xlow,xupp;
 269   Double_t step;
 270   Double_t i;
 271
 272
 273   // MP shift correction
 274   mpc = par[1] - mpshift * par[0];
 275
 276   // Range of convolution integral
 277   xlow = x[0] - sc * par[2];
 278   xupp = x[0] + sc * par[2];
 279
 280   step = (xupp-xlow) / np;
 281
 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]);
 287
 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]);
 291   }
 292
 293   return (step * sum * invsq2pi / par[2]);
 294 }
 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]);
 299 }
 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);
 306   if(opt==0){
 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;
 310   }
 311   if(opt==1){
 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;
 315   }
 316   if(opt==2){
 317     fCondFunPiLay[lay]=condFun;
 318     if(mom<0.6 &&dedx<20)fCondFunPiLay[lay]=0.001;
 319   }
 320
 321 }
 322 //_______________________________________________________________________
 323 Float_t AliITSPident::CookCombinedBayes(Double_t condfun[][3],Float_t *prior,Int_t part)const {
 324   //Bayesian combined PID in the ITS
 325   Int_t test=0;
 326   Float_t bayes;
 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++;
 332   }
 333   if(test==4){
 334     if ((pprior[0]!=0 || pprior[1]!=0 ||pprior[2]!=0)&&CookSum(condfun,pprior)!=0){
 335
 336       bayes=pprior[part]*CookProd(condfun,part)*1/CookSum(condfun,pprior);
 337
 338     }
 339     else bayes=-100;
 340   }
 341   else bayes=-100;
 342   return bayes;
 343 }
 344 //_______________________________________________________________________
 345 Float_t AliITSPident::CookProd(Double_t condfun[][3],Int_t part)const{
 346   //
 347   Float_t p=1;
 348   for(Int_t lay=0;lay<4;lay++){
 349
 350     p=p*condfun[lay][part];
 351   }
 352
 353   return p;
 354 }
 355 //_______________________________________________________________________
 356 Float_t AliITSPident::CookSum(Double_t condfun[][3],Float_t *prior)const{
 357   //
 358   Float_t sum=0;
 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);
 364   }
 365   return sum;
 366
 367 }