[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(){
  // default constructor
  fMom=0;
  fdEdx=0;
  fInvPt=0;
  fPBayesp=0;
  fPBayesk=0;
  fPBayespi=0;
  fPPriorip=0;
  fPPriorik=0;
  fPPrioripi=0;
  fPPriorie=0;
  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) {
  // Copy constructor
  // Copies are not allowed. The method is protected to avoid misuse.
  Error("AliITSPident","Copy constructor not allowed\n");
}

//______________________________________________________________________
AliITSPident& AliITSPident::operator=(const AliITSPident& /* ob */){
  // Assignment operator
  // Assignment is not allowed. The method is protected to avoid misuse.
  Error("= operator","Assignment operator not allowed\n");
  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){

  //test

  fMom=mom;
  fInvPt=invPt;
  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");
    
  }
  fPPriorip=priorip;
  fPPriorik=priorik;
  fPPrioripi=prioripi;
  fPPriorie=priorie;


  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=dEdx;
}


//_______________________________________________________________________
AliITSPident::AliITSPident(AliITStrackV2 *trackITS,AliITSSteerPid *sp,Float_t *Qlay,Float_t priorip,Float_t priorik,Float_t prioripi,Float_t priorie){
  //
  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");
  }
  fPPriorip=priorip;
  fPPriorik=priorik;
  fPPrioripi=prioripi;
  fPPriorie=priorie;


  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	//_______________________________________________________________________
	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);
149c4c59	60
e62c1aea	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
e62c1aea	94	fPBayesp=CookCombinedBayes(condFun,prior,0);
e62c1aea	95	fPBayesk=CookCombinedBayes(condFun,prior,1);
e62c1aea	96	fPBayespi=CookCombinedBayes(condFun,prior,2);
e62c1aea	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
e62c1aea	149	fPBayesp=CookCombinedBayes(condFun,prior,0);
e62c1aea	150	fPBayesk=CookCombinedBayes(condFun,prior,1);
e62c1aea	151	fPBayespi=CookCombinedBayes(condFun,prior,2);
e62c1aea	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(2TMath::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	}