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

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

/////////////////////////////////////////////////////////
//Class for PID in the ITS                             //
//                                                     //
//                                                     //
/////////////////////////////////////////////////////////

#include <Riostream.h>
#include <TF1.h>
#include <TMath.h>

#include "AliITSPidParItem.h"

ClassImp(AliITSPidParItem)
//____________________________________________________________________
AliITSPidParItem::AliITSPidParItem():
fPCenter(0),
fPWidth(0)
{
  // default constructor
  for(Int_t i=0;i<39;i++){
  fBuff[i]=0;
  }
}//____________________________________________________________________
AliITSPidParItem::AliITSPidParItem(Float_t center,Float_t width,Double_t *buff):
fPCenter(center),
fPWidth(width)
{
  // standard constructor
  for (Int_t i=0;i<39;i++) fBuff[i]=buff[i];

}

//____________________________________________________________________
void AliITSPidParItem::GetParameters(Double_t *buff) const{
  //get all the parameters
  for (Int_t i=0;i<39;i++) buff[i]=fBuff[i];

}
//____________________________________________________________________
void AliITSPidParItem::GetProtonPar(Double_t *buffp) const{
  //get the protons' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
  buffp[0]=fBuff[0];
  buffp[1]=fBuff[1];
  buffp[2]=fBuff[2];
  buffp[3]=fBuff[3];
  buffp[4]=fBuff[8];
  buffp[5]=fBuff[9];
  buffp[6]=fBuff[10];
}
//____________________________________________________________________
void AliITSPidParItem::GetKaonPar(Double_t *buffk) const{
  //get the kaons' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
  buffk[0]=fBuff[13];
  buffk[1]=fBuff[14];
  buffk[2]=fBuff[15];
  buffk[3]=fBuff[16];
  buffk[4]=fBuff[21];
  buffk[5]=fBuff[22];
  buffk[6]=fBuff[23];
}
//____________________________________________________________________
void AliITSPidParItem::GetPionPar(Double_t *buffpi) const{
  //get the pions' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
  buffpi[0]=fBuff[26];
  buffpi[1]=fBuff[27];
  buffpi[2]=fBuff[28];
  buffpi[3]=fBuff[29];
  buffpi[4]=fBuff[34];
  buffpi[5]=fBuff[35];
  buffpi[6]=fBuff[36];
}
//____________________________________________________________________
void AliITSPidParItem::GetPar0(Double_t *buff0) const{
  //Width Landau for protons, kaons, pions.
  buff0[0]=fBuff[0];
  buff0[1]=fBuff[13];
  buff0[2]=fBuff[26];
}
//____________________________________________________________________
void AliITSPidParItem::GetPar1(Double_t *buff1) const{
  //Most Probable for protons, kaons, pions.
  buff1[0]=fBuff[1];
  buff1[1]=fBuff[14];
  buff1[2]=fBuff[27];
}
//____________________________________________________________________
void AliITSPidParItem::GetPar2(Double_t *buff2) const{
  //Area for protons, kaons, pions.
  buff2[0]=fBuff[2];
  buff2[1]=fBuff[15];
  buff2[2]=fBuff[28];
}
//____________________________________________________________________
void AliITSPidParItem::GetPar3(Double_t *buff3) const{
  //Width Gaussian for protons, kaons, pions.
  buff3[0]=fBuff[3];
  buff3[1]=fBuff[16];
  buff3[2]=fBuff[29];
}
//____________________________________________________________________
void AliITSPidParItem::GetChisquare(Double_t *buffchi) const{
  //Chi2 of the fit for protons, kaons, pions.
  buffchi[0]=fBuff[8];
  buffchi[1]=fBuff[21];
  buffchi[2]=fBuff[34];
}
//____________________________________________________________________
void AliITSPidParItem::GetNDF(Double_t *buffndf) const{
  //NDF of the fit for protons, kaons, pions.
  buffndf[0]=fBuff[9];
  buffndf[1]=fBuff[22];
  buffndf[2]=fBuff[35];
}
//____________________________________________________________________
void AliITSPidParItem::GetProParFun(Double_t *pfun) const{
  //some Protons parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
  pfun[0]=fBuff[0];
  pfun[1]=fBuff[1];
  pfun[2]=fBuff[2];
  pfun[3]=fBuff[3];
  pfun[4]=fBuff[10];
}
//____________________________________________________________________
void AliITSPidParItem::GetKaoParFun(Double_t *kfun) const{
  //some Kaons parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
  kfun[0]=fBuff[13];
  kfun[1]=fBuff[14];
  kfun[2]=fBuff[15];
  kfun[3]=fBuff[16];
  kfun[4]=fBuff[23];
}
//____________________________________________________________________
void AliITSPidParItem::GetPiParFun(Double_t *pifun) const{
  //some Pions parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
  pifun[0]=fBuff[26];
  pifun[1]=fBuff[27];
  pifun[2]=fBuff[28];
  pifun[3]=fBuff[29];
  pifun[4]=fBuff[36];
}
//____________________________________________________________________
void AliITSPidParItem::GetRangeLim(Double_t *range) const{
  //Range limits for the response functions
  range[0]=fBuff[11];//proton low
  range[1]=fBuff[12];//proton high
  range[2]=fBuff[24];//kaon low
  range[3]=fBuff[25];//kaon high
  range[4]=fBuff[37];//pion low
  range[5]=fBuff[38];//pion high
}
//____________________________________________________________________
void AliITSPidParItem::GetProtonParErr(Double_t *bufferp)const{
  //errors on the protons' parameters
  bufferp[0]=fBuff[4];
  bufferp[1]=fBuff[5];
  bufferp[2]=fBuff[6];
  bufferp[3]=fBuff[7];
}
//____________________________________________________________________
void AliITSPidParItem::GetKaonParErr(Double_t *bufferk)const{
  //errors on the kaons' parameters
  bufferk[0]=fBuff[17];
  bufferk[1]=fBuff[18];
  bufferk[2]=fBuff[19];
  bufferk[3]=fBuff[20];
}
//____________________________________________________________________
void AliITSPidParItem::GetPionParErr(Double_t *bufferpi)const{
  //errors on the pions' parameters
  bufferpi[0]=fBuff[30];
  bufferpi[1]=fBuff[31];
  bufferpi[2]=fBuff[32];
  bufferpi[3]=fBuff[33];
}
//____________________________________________________________________
void AliITSPidParItem::PrintParameters() const {
  // Prints the data members of this class
cout<<"==============================***************"<<endl;
cout<<"Momentum (GeV/c) - center of the bin - "<<fPCenter<<endl;
cout<<" Width of momentum bin (GeV/c) "<<fPWidth<<endl;
for (Int_t i=0;i<39;i++) cout<<"Parameter"<<i<<" = "<<fBuff[i]<<endl;

}

//_______________________________________________________________________
TF1* AliITSPidParItem::CookFunIts(TString namefun,Double_t *par,Double_t rangei,Double_t rangef,TString comment){
  //
  TF1 *fun;
    if (par[4]!=0) {
      fun=new TF1(comment,Langaufun2,rangei,rangef,5);
      fun->SetParameters(par);
     
    }
    else {fun=new TF1(namefun,"0");}
    return fun;
}

//_________________________________________________________________________
Double_t AliITSPidParItem::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 AliITSPidParItem::Langaufun2(Double_t *x, Double_t *par){
  //
  return 1/par[4]*Langaufun(x,par);
}
Commit	Line	Data
090026bf	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
e62c1aea	18	/////////////////////////////////////////////////////////
	19	//Class for PID in the ITS //
	20	// //
	21	// //
	22	/////////////////////////////////////////////////////////
	23
	24	#include <Riostream.h>
	25	#include <TF1.h>
090026bf	26	#include <TMath.h>
090026bf	27
e62c1aea	28	#include "AliITSPidParItem.h"
	29
	30	ClassImp(AliITSPidParItem)
	31	//____________________________________________________________________
94631b2f	32	AliITSPidParItem::AliITSPidParItem():
	33	fPCenter(0),
	34	fPWidth(0)
	35	{
e62c1aea	36	// default constructor
e62c1aea	37	for(Int_t i=0;i<39;i++){
	38	fBuff[i]=0;
	39	}
	40	}//____________________________________________________________________
94631b2f	41	AliITSPidParItem::AliITSPidParItem(Float_t center,Float_t width,Double_t *buff):
	42	fPCenter(center),
	43	fPWidth(width)
	44	{
e62c1aea	45	// standard constructor
e62c1aea	46	for (Int_t i=0;i<39;i++) fBuff[i]=buff[i];
	47
	48	}
	49
	50	//____________________________________________________________________
	51	void AliITSPidParItem::GetParameters(Double_t *buff) const{
	52	//get all the parameters
	53	for (Int_t i=0;i<39;i++) buff[i]=fBuff[i];
	54
	55	}
	56	//____________________________________________________________________
	57	void AliITSPidParItem::GetProtonPar(Double_t *buffp) const{
	58	//get the protons' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
	59	buffp[0]=fBuff[0];
	60	buffp[1]=fBuff[1];
	61	buffp[2]=fBuff[2];
	62	buffp[3]=fBuff[3];
	63	buffp[4]=fBuff[8];
	64	buffp[5]=fBuff[9];
	65	buffp[6]=fBuff[10];
	66	}
	67	//____________________________________________________________________
	68	void AliITSPidParItem::GetKaonPar(Double_t *buffk) const{
	69	//get the kaons' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
	70	buffk[0]=fBuff[13];
	71	buffk[1]=fBuff[14];
	72	buffk[2]=fBuff[15];
	73	buffk[3]=fBuff[16];
	74	buffk[4]=fBuff[21];
	75	buffk[5]=fBuff[22];
	76	buffk[6]=fBuff[23];
	77	}
	78	//____________________________________________________________________
	79	void AliITSPidParItem::GetPionPar(Double_t *buffpi) const{
	80	//get the pions' parameters (Width Landau, Most Probable, Area, Width Gaussian, Chi2 fit, NDF fit, Integral fit)
	81	buffpi[0]=fBuff[26];
	82	buffpi[1]=fBuff[27];
	83	buffpi[2]=fBuff[28];
	84	buffpi[3]=fBuff[29];
	85	buffpi[4]=fBuff[34];
	86	buffpi[5]=fBuff[35];
	87	buffpi[6]=fBuff[36];
	88	}
	89	//____________________________________________________________________
	90	void AliITSPidParItem::GetPar0(Double_t *buff0) const{
	91	//Width Landau for protons, kaons, pions.
	92	buff0[0]=fBuff[0];
	93	buff0[1]=fBuff[13];
	94	buff0[2]=fBuff[26];
	95	}
	96	//____________________________________________________________________
	97	void AliITSPidParItem::GetPar1(Double_t *buff1) const{
	98	//Most Probable for protons, kaons, pions.
	99	buff1[0]=fBuff[1];
	100	buff1[1]=fBuff[14];
	101	buff1[2]=fBuff[27];
	102	}
	103	//____________________________________________________________________
	104	void AliITSPidParItem::GetPar2(Double_t *buff2) const{
	105	//Area for protons, kaons, pions.
	106	buff2[0]=fBuff[2];
	107	buff2[1]=fBuff[15];
	108	buff2[2]=fBuff[28];
	109	}
110	//____________________________________________________________________
111	void AliITSPidParItem::GetPar3(Double_t *buff3) const{
112	//Width Gaussian for protons, kaons, pions.
113	buff3[0]=fBuff[3];
114	buff3[1]=fBuff[16];
115	buff3[2]=fBuff[29];
116	}
117	//____________________________________________________________________
118	void AliITSPidParItem::GetChisquare(Double_t *buffchi) const{
119	//Chi2 of the fit for protons, kaons, pions.
120	buffchi[0]=fBuff[8];
121	buffchi[1]=fBuff[21];
122	buffchi[2]=fBuff[34];
123	}
124	//____________________________________________________________________
125	void AliITSPidParItem::GetNDF(Double_t *buffndf) const{
126	//NDF of the fit for protons, kaons, pions.
127	buffndf[0]=fBuff[9];
128	buffndf[1]=fBuff[22];
129	buffndf[2]=fBuff[35];
130	}
131	//____________________________________________________________________
132	void AliITSPidParItem::GetProParFun(Double_t *pfun) const{
133	//some Protons parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
134	pfun[0]=fBuff[0];
135	pfun[1]=fBuff[1];
136	pfun[2]=fBuff[2];
137	pfun[3]=fBuff[3];
138	pfun[4]=fBuff[10];
139	}
140	//____________________________________________________________________
141	void AliITSPidParItem::GetKaoParFun(Double_t *kfun) const{
142	//some Kaons parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
143	kfun[0]=fBuff[13];
144	kfun[1]=fBuff[14];
145	kfun[2]=fBuff[15];
146	kfun[3]=fBuff[16];
147	kfun[4]=fBuff[23];
148	}
149	//____________________________________________________________________
150	void AliITSPidParItem::GetPiParFun(Double_t *pifun) const{
151	//some Pions parameters: Width Landau, Most Probable, Area, Width Gaussian, Integral fit
152	pifun[0]=fBuff[26];
153	pifun[1]=fBuff[27];
154	pifun[2]=fBuff[28];
155	pifun[3]=fBuff[29];
156	pifun[4]=fBuff[36];
157	}
158	//____________________________________________________________________
159	void AliITSPidParItem::GetRangeLim(Double_t *range) const{
160	//Range limits for the response functions
161	range[0]=fBuff[11];//proton low
162	range[1]=fBuff[12];//proton high
163	range[2]=fBuff[24];//kaon low
164	range[3]=fBuff[25];//kaon high
165	range[4]=fBuff[37];//pion low
166	range[5]=fBuff[38];//pion high
167	}
168	//____________________________________________________________________
169	void AliITSPidParItem::GetProtonParErr(Double_t *bufferp)const{
170	//errors on the protons' parameters
171	bufferp[0]=fBuff[4];
172	bufferp[1]=fBuff[5];
173	bufferp[2]=fBuff[6];
174	bufferp[3]=fBuff[7];
175	}
176	//____________________________________________________________________
177	void AliITSPidParItem::GetKaonParErr(Double_t *bufferk)const{
178	//errors on the kaons' parameters
179	bufferk[0]=fBuff[17];
180	bufferk[1]=fBuff[18];
181	bufferk[2]=fBuff[19];
182	bufferk[3]=fBuff[20];
183	}
184	//____________________________________________________________________
185	void AliITSPidParItem::GetPionParErr(Double_t *bufferpi)const{
186	//errors on the pions' parameters
187	bufferpi[0]=fBuff[30];
188	bufferpi[1]=fBuff[31];
189	bufferpi[2]=fBuff[32];
190	bufferpi[3]=fBuff[33];
191	}
192	//____________________________________________________________________
193	void AliITSPidParItem::PrintParameters() const {
194	// Prints the data members of this class
195	cout<<"==============================***************"<<endl;
196	cout<<"Momentum (GeV/c) - center of the bin - "<<fPCenter<<endl;
197	cout<<" Width of momentum bin (GeV/c) "<<fPWidth<<endl;
198	for (Int_t i=0;i<39;i++) cout<<"Parameter"<<i<<" = "<<fBuff[i]<<endl;
199
200	}
201
202	//_______________________________________________________________________
203	TF1* AliITSPidParItem::CookFunIts(TString namefun,Double_t *par,Double_t rangei,Double_t rangef,TString comment){
204	//
205	TF1 *fun;
206	if (par[4]!=0) {
207	fun=new TF1(comment,Langaufun2,rangei,rangef,5);
208	fun->SetParameters(par);
209
210	}
211	else {fun=new TF1(namefun,"0");}
212	return fun;
213	}
214
215	//_________________________________________________________________________
216	Double_t AliITSPidParItem::Langaufun(Double_t x, Double_t par) {
217
218	//Fit parameters:
219	//par[0]=Width (scale) parameter of Landau density
220	//par[1]=Most Probable (MP, location) parameter of Landau density
221	//par[2]=Total area (integral -inf to inf, normalization constant)
222	//par[3]=Width (sigma) of convoluted Gaussian function
223	//
224	//In the Landau distribution (represented by the CERNLIB approximation),
225	//the maximum is located at x=-0.22278298 with the location parameter=0.
226	//This shift is corrected within this function, so that the actual
227	//maximum is identical to the MP parameter.
228
229	// Numeric constants
230	Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
231	Double_t mpshift = -0.22278298; // Landau maximum location
232
233	// Control constants
234	Double_t np = 100.0; // number of convolution steps
235	Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
236
237	// Variables
238	Double_t xx;
239	Double_t mpc;
240	Double_t fland;
241	Double_t sum = 0.0;
242	Double_t xlow,xupp;
243	Double_t step;
244	Double_t i;
245
246
247	// MP shift correction
248	mpc = par[1] - mpshift * par[0];
249
250	// Range of convolution integral
251	xlow = x[0] - sc * par[3];
252	xupp = x[0] + sc * par[3];
253
254	step = (xupp-xlow) / np;
255
256	// Convolution integral of Landau and Gaussian by sum
257	for(i=1.0; i<=np/2; i++) {
258	xx = xlow + (i-.5) * step;
259	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
260	sum += fland * TMath::Gaus(x[0],xx,par[3]);
261
262	xx = xupp - (i-.5) * step;
263	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
264	sum += fland * TMath::Gaus(x[0],xx,par[3]);
265	}
266
267	return (par[2] * step * sum * invsq2pi / par[3]);
268	}
269	//_______________________________________________________________________
270	Double_t AliITSPidParItem::Langaufun2(Double_t x, Double_t par){
271	//
272	return 1/par[4]*Langaufun(x,par);
273	}