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