[u/mrichter/AliRoot.git] / ITS / AliITSPident.cxx

//////////////////////////////////////////////////////////////////////////
// USER Class for PID in the ITS                                          //
//The PID is based on the likelihood of all the four ITS' layers,         //
//without using the truncated mean for the dE/dx. The response            //
//functions for each layer are convoluted Landau-Gaussian functions.      // 
// Origin: Elena Bruna bruna@to.infn.it,, Massimo Masera masera@to.infn.it//
//////////////////////////////////////////////////////////////////////////
#include "AliITSPident.h"

ClassImp(AliITSPident)
  //_______________________________________________________________________
AliITSPident::AliITSPident():
fMom(0),
fdEdx(0),
fPBayesp(0),
fPBayesk(0),
fPBayespi(0),
fPPriorip(0),
fPPriorik(0),
fPPrioripi(0),
fPPriorie(0),
fInvPt(0)
{
  // default constructor
  for (Int_t i=0;i<4;i++){
    fCondFunProLay[i]=0;
    fCondFunKLay[i]=0;
    fCondFunPiLay[i]=0;
  }
}
//_______________________________________________________________________
AliITSPident::~AliITSPident(){
  // destructor
}
//______________________________________________________________________
AliITSPident::AliITSPident(const AliITSPident &ob) :TObject(ob),
fMom(ob.fMom),
fdEdx(ob.fdEdx),
fPBayesp(ob.fPBayesp),
fPBayesk(ob.fPBayesk),
fPBayespi(ob.fPBayespi),
fPPriorip(ob.fPPriorip),
fPPriorik(ob.fPPriorik),
fPPrioripi(ob.fPPrioripi),
fPPriorie(ob.fPPriorie),
fInvPt(ob.fInvPt)  
{
  // Copy constructor
}

//______________________________________________________________________
AliITSPident& AliITSPident::operator=(const AliITSPident&  ob){
  // Assignment operator
  this->~AliITSPident();
  new(this) AliITSPident(ob);
  return *this;
}


//_______________________________________________________________________
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):
fMom(mom),
fdEdx(dEdx),
fPBayesp(0),
fPBayesk(0),
fPBayespi(0),
fPPriorip(priorip),
fPPriorik(priorik),
fPPrioripi(prioripi),
fPPriorie(priorie),
fInvPt(invPt){

  //test

  for(Int_t la=0;la<4;la++){//loop on layers
    Double_t parp[3];Double_t park[3];Double_t parpi[3];
    sp->GetParFitLayer(la,fMom,parp,park,parpi);

    Double_t range[6];
    range[0]=0.3*parp[1];
    range[1]=2.*parp[1];

    range[2]=0.3*park[1];
    range[3]=2.*park[1];
    
    range[4]=0.3*parpi[1];
    range[5]=2.*parpi[1];
      CookFunItsLay(la,0,parp,Qlay[la],fMom,range[0],range[1],"fPro");
      CookFunItsLay(la,1,park,Qlay[la],fMom,range[2],range[3],"fKao");
      CookFunItsLay(la,2,parpi,Qlay[la],fMom,range[4],range[5],"fPi");
    
  }

  Float_t prior[4];Double_t condFun[4][3];

  prior[0]=fPPriorip;
  prior[1]=fPPriorik;
  prior[2]=fPPrioripi;
  prior[3]=fPPriorie;
  for(Int_t la=0;la<4;la++){
    condFun[la][0]= fCondFunProLay[la];
    condFun[la][1]= fCondFunKLay[la]; 
    condFun[la][2]= fCondFunPiLay[la];
    
  }

  fPBayesp=CookCombinedBayes(condFun,prior,0);
  fPBayesk=CookCombinedBayes(condFun,prior,1); 
  fPBayespi=CookCombinedBayes(condFun,prior,2); 
 }


//_______________________________________________________________________
AliITSPident::AliITSPident(AliITStrackV2 *trackITS,AliITSSteerPid *sp,Float_t *Qlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie):
fMom(0),
fdEdx(0),
fPBayesp(0),
fPBayesk(0),
fPBayespi(0),
fPPriorip(priorip),
fPPriorik(priorik),
fPPrioripi(prioripi),
fPPriorie(priorie),
fInvPt(0)
{
  //
  Double_t xr;
  Double_t par[5];
  trackITS->GetExternalParameters(xr,par);
  if (par[4]!=0) {
    Float_t lamb=TMath::ATan(par[3]);
    fMom=1/(TMath::Abs(par[4])*TMath::Cos(lamb));
  }
 
  fInvPt=par[4];
 
  for(Int_t la=0;la<4;la++){//loop on layers
    Double_t parp[3];Double_t park[3];Double_t parpi[3];
    sp->GetParFitLayer(la,fMom,parp,park,parpi);
    Double_t range[8];
    range[0]=0.3*parp[1];
    range[1]=2.*parp[1];

    range[2]=0.3*park[1];
    range[3]=2.*park[1];
    
    range[4]=0.3*parpi[1];
    range[5]=2.*parpi[1];
      CookFunItsLay(la,0,parp,Qlay[la],fMom,range[0],range[1],"fPro");
      CookFunItsLay(la,1,park,Qlay[la],fMom,range[2],range[3],"fKao");
      CookFunItsLay(la,2,parpi,Qlay[la],fMom,range[4],range[5],"fPi");
  }

  Float_t prior[4];Double_t condFun[4][3];

  prior[0]=fPPriorip;
  prior[1]=fPPriorik;
  prior[2]=fPPrioripi;
  prior[3]=fPPriorie;
  for(Int_t la=0;la<4;la++){
    condFun[la][0]= fCondFunProLay[la];
    condFun[la][1]= fCondFunKLay[la]; 
    condFun[la][2]= fCondFunPiLay[la];
    
  }

  fPBayesp=CookCombinedBayes(condFun,prior,0);
  fPBayesk=CookCombinedBayes(condFun,prior,1); 
  fPBayespi=CookCombinedBayes(condFun,prior,2); 
  fdEdx=trackITS->GetdEdx();

}
//_______________________________________________________________________
void AliITSPident::PrintParameters() const{
 //print parameters
  cout<<"___________________________\n";
  cout<<"Track Local Momentum = "<<"  "<<fMom<<endl;
  cout<<"Track dEdx = "<<"  "<<fdEdx<<endl;
  cout<<"Inv Pt = "<<fInvPt<<endl;
}

//_______________________________________________________________________
Double_t AliITSPident::Langaufun(Double_t *x, Double_t *par) {

  //Fit parameters:
  //par[0]=Width (scale) parameter of Landau density
  //par[1]=Most Probable (MP, location) parameter of Landau density
  //par[2]=Total area (integral -inf to inf, normalization constant)
  //par[3]=Width (sigma) of convoluted Gaussian function
  //
  //In the Landau distribution (represented by the CERNLIB approximation), 
  //the maximum is located at x=-0.22278298 with the location parameter=0.
  //This shift is corrected within this function, so that the actual
  //maximum is identical to the MP parameter.

  // Numeric constants
  Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
  Double_t mpshift  = -0.22278298;       // Landau maximum location

  // Control constants
  Double_t np = 100.0;      // number of convolution steps
  Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

  // Variables
  Double_t xx;
  Double_t mpc;
  Double_t fland;
  Double_t sum = 0.0;
  Double_t xlow,xupp;
  Double_t step;
  Double_t i;


  // MP shift correction
  mpc = par[1] - mpshift * par[0]; 

  // Range of convolution integral
  xlow = x[0] - sc * par[3];
  xupp = x[0] + sc * par[3];

  step = (xupp-xlow) / np;

  // Convolution integral of Landau and Gaussian by sum
  for(i=1.0; i<=np/2; i++) {
    xx = xlow + (i-.5) * step;
    fland = TMath::Landau(xx,mpc,par[0]) / par[0];
    sum += fland * TMath::Gaus(x[0],xx,par[3]);

    xx = xupp - (i-.5) * step;
    fland = TMath::Landau(xx,mpc,par[0]) / par[0];
    sum += fland * TMath::Gaus(x[0],xx,par[3]);
  }

  return (par[2] * step * sum * invsq2pi / par[3]);
}

//_______________________________________________________________________
Double_t AliITSPident::Langaufun2(Double_t *x, Double_t *par){
  //normalized langaufun
  return 1/par[4]*Langaufun(x,par);
}
//_______________________________________________________________________
Double_t AliITSPident::Langaufunnorm(Double_t *x, Double_t *par){
   //Fit parameters:
  //par[0]=Width (scale) parameter of Landau density
  //par[1]=Most Probable (MP, location) parameter of Landau density
  
  //par[2]=Width (sigma) of convoluted Gaussian function
  //
  //In the Landau distribution (represented by the CERNLIB approximation), 
  //the maximum is located at x=-0.22278298 with the location parameter=0.
  //This shift is corrected within this function, so that the actual
  //maximum is identical to the MP parameter.

  // Numeric constants
  Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
  Double_t mpshift  = -0.22278298;       // Landau maximum location

  // Control constants
  Double_t np = 100.0;      // number of convolution steps
  Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

  // Variables
  Double_t xx;
  Double_t mpc;
  Double_t fland;
  Double_t sum = 0.0;
  Double_t xlow,xupp;
  Double_t step;
  Double_t i;


  // MP shift correction
  mpc = par[1] - mpshift * par[0]; 

  // Range of convolution integral
  xlow = x[0] - sc * par[2];
  xupp = x[0] + sc * par[2];

  step = (xupp-xlow) / np;

  // Convolution integral of Landau and Gaussian by sum
  for(i=1.0; i<=np/2; i++) {
    xx = xlow + (i-.5) * step;
    fland = TMath::Landau(xx,mpc,par[0]) / par[0];
    sum += fland * TMath::Gaus(x[0],xx,par[2]);

    xx = xupp - (i-.5) * step;
    fland = TMath::Landau(xx,mpc,par[0]) / par[0];
    sum += fland * TMath::Gaus(x[0],xx,par[2]);
  }

  return (step * sum * invsq2pi / par[2]);
}
//_______________________________________________________________________
Double_t AliITSPident::Gaus2(Double_t *x, Double_t *par){
  //normalized gaussian function
  return 1/(sqrt(2*TMath::Pi())*par[1])*TMath::Gaus(x[0],par[0],par[1]);
}
//_______________________________________________________________________
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){
  //it gives the response functions
 TF1 funLay(comment,Langaufunnorm,rangei,rangef,3);
 funLay.SetParameters(par);
 Double_t condFun=funLay.Eval(dedx);
  if(opt==0){
    fCondFunProLay[lay]=condFun;
    if(mom<0.4 && dedx<100)fCondFunProLay[lay]=0.00001;
    if(mom<0.4 &&dedx<50)fCondFunProLay[lay]=0.0000001;
  }
  if(opt==1){
    fCondFunKLay[lay]=condFun; 
    if(mom<0.25 && dedx<100)fCondFunKLay[lay]=0.00001;
    if(mom<0.4 &&dedx<30)fCondFunKLay[lay]=0.0000001;   
  }
  if(opt==2){
    fCondFunPiLay[lay]=condFun;
    if(mom<0.6 &&dedx<20)fCondFunPiLay[lay]=0.001;
  }

}
//_______________________________________________________________________
Float_t AliITSPident::CookCombinedBayes(Double_t condfun[][3],Float_t *prior,Int_t part)const {
  //Bayesian combined PID in the ITS
  Int_t test=0; 
  Float_t bayes;
  Float_t pprior[4]={0,0,0,0};
  for(Int_t j=0;j<4;j++)pprior[j]=prior[j];
  pprior[2]+=pprior[3];//prior for electrons summed to priors for pions
  for(Int_t i=0;i<4;i++){//layer
    if (condfun[i][0]>0 || condfun[i][1]>0 ||condfun[i][2]>0) test++; 
  }
  if(test==4){
    if ((pprior[0]!=0 || pprior[1]!=0 ||pprior[2]!=0)&&CookSum(condfun,pprior)!=0){

      bayes=pprior[part]*CookProd(condfun,part)*1/CookSum(condfun,pprior);
  
    }
    else bayes=-100;
  }
  else bayes=-100;
  return bayes;
}
//_______________________________________________________________________
Float_t AliITSPident::CookProd(Double_t condfun[][3],Int_t part)const{
  //
  Float_t p=1;
  for(Int_t lay=0;lay<4;lay++){
  
    p=p*condfun[lay][part];
  }

  return p;
}
//_______________________________________________________________________
Float_t AliITSPident::CookSum(Double_t condfun[][3],Float_t *prior)const{
  //
  Float_t sum=0;
  Float_t pprior[4]={0,0,0,0};
  for(Int_t j=0;j<4;j++)pprior[j]=prior[j];
  pprior[2]+=pprior[3];//prior for electrons summed to priors for pions
  for(Int_t i=0;i<3;i++){//sum over the particles--electrons excluded
    sum+=pprior[i]*CookProd(condfun,i);
  }
  return sum;

}
Commit	Line	Data
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	//////////////////////////////////////////////////////////////////////////
	8	#include "AliITSPident.h"
	9
	10	ClassImp(AliITSPident)
	11	//_______________________________________________________________________
94631b2f	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	{
e62c1aea	24	// default constructor
e62c1aea	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	//______________________________________________________________________
94631b2f	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	{
e62c1aea	48	// Copy constructor
e62c1aea	49	}
	50
	51	//______________________________________________________________________
94631b2f	52	AliITSPident& AliITSPident::operator=(const AliITSPident& ob){
e62c1aea	53	// Assignment operator
94631b2f	54	this->~AliITSPident();
94631b2f	55	new(this) AliITSPident(ob);
e62c1aea	56	return *this;
	57	}
	58
	59
	60	//_______________________________________________________________________
94631b2f	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){
e62c1aea	72
	73	//test
	74
e62c1aea	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);
149c4c59	78
e62c1aea	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	}
e62c1aea	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
e62c1aea	107	fPBayesp=CookCombinedBayes(condFun,prior,0);
e62c1aea	108	fPBayesk=CookCombinedBayes(condFun,prior,1);
e62c1aea	109	fPBayespi=CookCombinedBayes(condFun,prior,2);
94631b2f	110	}
e62c1aea	111
	112
	113	//_______________________________________________________________________
94631b2f	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	{
e62c1aea	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	}
e62c1aea	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
e62c1aea	167	fPBayesp=CookCombinedBayes(condFun,prior,0);
e62c1aea	168	fPBayesk=CookCombinedBayes(condFun,prior,1);
e62c1aea	169	fPBayespi=CookCombinedBayes(condFun,prior,2);
e62c1aea	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(2TMath::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	}