[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 "AliITStrackV2.h"
#include "AliITSPident.h"
#include "AliITSSteerPid.h"
#include <Riostream.h>

ClassImp(AliITSPident)
  //_______________________________________________________________________
AliITSPident::AliITSPident():
fMom(0),
fdEdx(0),
fPBayesp(0),
fPBayesk(0),
fPBayespi(0),
fPPriorip(0),
fPPriorik(0),
fPPrioripi(0),
fPPriorie(0)
{
  // default constructor
  for (Int_t i=0;i<8;i++){
    fCondFunProLay[i]=0;
    fCondFunKLay[i]=0;
    fCondFunPiLay[i]=0;
  }
  for (Int_t i=0;i<4;i++)fNcls[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)
{
  // 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 dEdx,AliITSSteerPid *sp,Float_t *Qlay,Float_t *nlay,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)
{
  //
  for (Int_t i=0;i<8;i++){
    fCondFunProLay[i]=-1;
    fCondFunKLay[i]=-1;
    fCondFunPiLay[i]=-1;
  }
  for(Int_t la=0;la<4;la++){//loop on layers
    Double_t parp[3];Double_t park[3];Double_t parpi[3];
    fNcls[la]=0;
    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];
    Int_t layer=la+2;
    for(Int_t ii=0;ii<8;ii++){
      if(nlay[ii]==layer){
	fNcls[la]++;
	if(Qlay[ii]>0){
	  sp->GetParFitLayer(la,fMom,parp,park,parpi);
	  CookFunItsLay(ii,0,parp,Qlay[ii],fMom,range[0],range[1],"fPro");
	  CookFunItsLay(ii,1,park,Qlay[ii],fMom,range[2],range[3],"fKao");
	  CookFunItsLay(ii,2,parpi,Qlay[ii],fMom,range[4],range[5],"fPi");
	}  
      }
    }
    
  }
 
  Float_t prior[4];Double_t condFun[8][3];

  prior[0]=fPPriorip;
  prior[1]=fPPriorik;
  prior[2]=fPPrioripi;
  prior[3]=fPPriorie;
  for(Int_t la=0;la<8;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 *nlay,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)
{
  //
  for (Int_t i=0;i<8;i++){
    fCondFunProLay[i]=-1;
    fCondFunKLay[i]=-1;
    fCondFunPiLay[i]=-1;
  }
  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));
  }
 
 
  for(Int_t la=0;la<4;la++){//loop on layers
    Double_t parp[3];Double_t park[3];Double_t parpi[3];
    fNcls[la]=0;
    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];

    Int_t layer=la+2;
    for(Int_t ii=0;ii<8;ii++){
      if(nlay[ii]==layer){
	fNcls[la]++;
	if(Qlay[ii]>0){
	  sp->GetParFitLayer(la,fMom,parp,park,parpi);
	  CookFunItsLay(ii,0,parp,Qlay[ii],fMom,range[0],range[1],"fPro");
	  CookFunItsLay(ii,1,park,Qlay[ii],fMom,range[2],range[3],"fKao");
	  CookFunItsLay(ii,2,parpi,Qlay[ii],fMom,range[4],range[5],"fPi");
	}  
      }
    }
    
  }

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

  prior[0]=fPPriorip;
  prior[1]=fPPriorik;
  prior[2]=fPPrioripi;
  prior[3]=fPPriorie;
  for(Int_t la=0;la<8;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::GetNclsPerLayer(Int_t *ncls) const{
  //number of clusters for each layer (sdd1,sdd2,ssd1,ssd2)
 for(Int_t la=0;la<4;la++){
   ncls[la]=fNcls[la];
 }

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