[u/mrichter/AliRoot.git] / STEER / STEERBase / AliEMCALPIDResponse.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.                  *
 **************************************************************************/
//////////////////////////////////////////////////////////////////////////
//                                                                      //
// AliEMCALPIDResponse                                                  //
//                                                                      //
// EMCAL class to perfom PID                                            //
// This is a prototype and still under development                      //
//                                                                      //
// ---------------------------------------------------------------------//
// GetNumberOfSigmas():                                                 //
//                                                                      //
// Electrons:  Number of Sigmas for E/p value                           //
//             Parametrization of LHC11a (after recalibration)          //
//                                                                      //
// NON electrons:                                                       //
//             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
//             --> return +/- 99                                        //
//             Otherwise                                                //
//             --> return nsigma (parametrization of LHC10e)            //
//                                                                      //
// NO Parametrization (outside pT range): --> return -999               //
//                                                                      //
// ---------------------------------------------------------------------//
// ComputeEMCALProbability():                                           //
//                                                                      //
// Electrons:  Probability from Gaussian distribution                   //
//                                                                      //
// NON electrons:                                                       //
//             Below or above E/p thresholds ( E/p < 0.5 || E/p > 1.5)  //
//             --> probability to find particles below or above thr.    //
//             Otherwise                                                //
//             -->  Probability from Gaussian distribution              //
//                  (proper normalization to each other?)               //
//                                                                      //
// NO Parametrization (outside pT range): --> return kFALSE             //
//////////////////////////////////////////////////////////////////////////

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

#include "AliEMCALPIDResponse.h"       //class header

#include "AliLog.h"   

ClassImp(AliEMCALPIDResponse)

//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
AliEMCALPIDResponse::AliEMCALPIDResponse():
  TObject(),
  fNorm(NULL),
  fCurrCentrality(-1.),
  fkPIDParams(NULL)
{
  //
  //  The default constructor
  //


  fNorm = new TF1("fNorm","gaus",-20,20); 
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
  TObject(other),
  fNorm(other.fNorm),
  fCurrCentrality(other.fCurrCentrality),
  fkPIDParams(other.fkPIDParams)
{
  //
  //  The copy constructor
  //

}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
{
  //
  //  The assignment operator
  //

  if(this == &other) return *this;
  
  // Make copy
  TObject::operator=(other);
  fNorm = other.fNorm;
  fCurrCentrality = other.fCurrCentrality;
  fkPIDParams = other.fkPIDParams;

 
  return *this;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
AliEMCALPIDResponse::~AliEMCALPIDResponse() {

  delete fNorm;

}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n,  Int_t charge) const  {
  //
  // Calculates the expected sigma of the PID signal as the function of 
  // the information stored in the track, for the specified particle type 
  //  
  //
  
  Double_t norm = 1.;

  // Check the charge
  if( charge != -1 && charge != 1){
    return norm;
  }

  // Get the parameters for this particle type and pt
  const TVectorD *params = GetParams(n, pt, charge);

  // IF not in momentum range, NULL is returned --> return default value
  if(!params) return norm;

  Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
  Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
  Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
  Double_t eopMax   = (*params)[5];   // max E/p value for parametrization
  Double_t probLow  = (*params)[6];   // probability to be below eopMin
  Double_t probHigh = (*params)[7];   // probability to be above eopMax

  // Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
  fNorm->SetParameters(1./TMath::Sqrt(2*TMath::Pi()*sigma*sigma),mean,sigma);
  norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);

  return norm;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliPID::EParticleType n,  Int_t charge) const {
      
  Double_t nsigma = -999.;

  // Check the charge
  if( charge != -1 && charge != 1){
    return nsigma;
  }

  // Get the parameters for this particle type and pt
  const TVectorD *params = GetParams(n, pt, charge);

  // IF not in momentum range, NULL is returned --> return default value
  if(!params) return nsigma;

  Double_t mean     = (*params)[2];   // mean value of Gausiian parametrization
  Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
  Double_t eopMin   = (*params)[4];   // min E/p value for parametrization
  Double_t eopMax   = (*params)[5];   // max E/p value for parametrization

  // if electron
  if(n == AliPID::kElectron){
    if(sigma != 0) nsigma = (eop - mean) / sigma;
  }

  // if NON electron
  else{
    if ( eop < eopMin )
      nsigma = -99;    // not parametrized (smaller than eopMin)
    else if ( eop > eopMax )
      nsigma = 99.;     // not parametrized (bigger than eopMax)
    else{
      if(sigma != 0) nsigma = (eop - mean) / sigma; 
    }
  }

  return nsigma;

}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
  //
  //
  Double_t fRange  = 5.0;   // hardcoded (???)
  Double_t nsigma  = 0.0;
  

  // Check the charge
  if( charge != -1 && charge != 1){
    return kFALSE;
  }

 
  // default value (will be returned, if pt below threshold)
  for (Int_t species = 0; species < nSpecies; species++) {
    pEMCAL[species] = 1./nSpecies;
  }

  // set E/p range
  if(eop < 0.05) eop = 0.05;
  if(eop > 2.00) eop = 2.00;
  
  for (Int_t species = 0; species < nSpecies; species++) {
    
    AliPID::EParticleType type = AliPID::EParticleType(species);

    // Get the parameters for this particle type and pt
    const TVectorD *params = GetParams(species, pt, charge);
    
    // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
    if(!params) return kFALSE;

    Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
    Double_t probLow  = (*params)[6];   // probability to be below eopMin
    Double_t probHigh = (*params)[7];   // probability to be above eopMax

    // get nsigma value for each particle type at this E/p value
    nsigma = GetNumberOfSigmas(pt,eop,type,charge);

    // electrons (standard Gaussian calculation of probabilities)
    if(type == AliPID::kElectron){
      if (TMath::Abs(nsigma) > fRange) {
	pEMCAL[species]=TMath::Exp(-0.5*fRange*fRange)/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
      }
      else{
	pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
      }
    }
    //NON electrons
    else{
      // E/p < eopMin  -->  return probability below E/p = eopMin
      if ( nsigma == -99){
	pEMCAL[species] = probLow;
      }
      // E/p > eopMax  -->  return probability above E/p = eopMax
      else if ( nsigma == 99){
	pEMCAL[species] = probHigh;
      }
      // in parametrized region --> calculate probability for corresponding Gauss curve
      else{
	pEMCAL[species]=TMath::Exp(-0.5*(nsigma)*(nsigma))/TMath::Sqrt(2*TMath::Pi()*sigma*sigma);
	
	// normalize to total probability == 1
	pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
      }
    }
  }

  return kTRUE;

}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
  //
  // returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
  //
  // 0 = momMin
  // 1 = momMax
  // 2 = mean of Gaus
  // 3 = sigma of Gaus
  // 4 = eopLow   
  // 5 = eopHig   
  // 6 = probLow  (not used for electrons)
  // 7 = probHigh (not used for electrons)
  //
  // for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
  // so first the correct centrality bin has to be found

  // **** Centrality bins (hard coded for the moment)
  const Int_t nCent = 7;
  Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};

  if(nParticle > AliPID::kSPECIES || nParticle <0) return NULL;
  if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons

  TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
  if(!particlePar) return NULL;

  const TVectorD *parameters = NULL;
  Double_t momMin = 0.;
  Double_t momMax = 0.;

  // is the centrality dependent parametrization used
  TString arrayName = particlePar->GetName();

  // centrality dependent parametrization
  if(arrayName.Contains("Centrality")){
    
    for(Int_t iCent = 0; iCent < nCent; iCent++){

      if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
	
	TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
	if(!centPar) return NULL;
	
	TIter centIter(centPar);
	parameters = NULL;
	momMin = 0.;
	momMax = 0.;
	
	while((parameters = static_cast<const TVectorD *>(centIter()))){
	  
	  momMin = (*parameters)[0];
	  momMax = (*parameters)[1];
	  
	  if( fPt > momMin && fPt < momMax ) return parameters;
	 
	} 
      }
    }
  }

  // NO centrality dependent parametrization
  else{

    TIter parIter(particlePar);
    while((parameters = static_cast<const TVectorD *>(parIter()))){
      
      momMin = (*parameters)[0];
      momMax = (*parameters)[1];
      
      if( fPt > momMin && fPt < momMax ) return parameters;
      
    }  
  }
  AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));

  return parameters;
}
Commit	Line	Data
b2138b40	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	// //
	17	// AliEMCALPIDResponse //
	18	// //
	19	// EMCAL class to perfom PID //
	20	// This is a prototype and still under development //
	21	// //
	22	// ---------------------------------------------------------------------//
	23	// GetNumberOfSigmas(): //
	24	// //
	25	// Electrons: Number of Sigmas for E/p value //
	26	// Parametrization of LHC11a (after recalibration) //
	27	// //
	28	// NON electrons: //
	29	// Below or above E/p thresholds ( E/p < 0.5 \|\| E/p > 1.5) //
	30	// --> return +/- 99 //
	31	// Otherwise //
	32	// --> return nsigma (parametrization of LHC10e) //
	33	// //
	34	// NO Parametrization (outside pT range): --> return -999 //
	35	// //
	36	// ---------------------------------------------------------------------//
	37	// ComputeEMCALProbability(): //
	38	// //
	39	// Electrons: Probability from Gaussian distribution //
	40	// //
	41	// NON electrons: //
	42	// Below or above E/p thresholds ( E/p < 0.5 \|\| E/p > 1.5) //
	43	// --> probability to find particles below or above thr. //
	44	// Otherwise //
	45	// --> Probability from Gaussian distribution //
	46	// (proper normalization to each other?) //
	47	// //
00a38d07	48	// NO Parametrization (outside pT range): --> return kFALSE //
b2138b40	49	//////////////////////////////////////////////////////////////////////////
	50
	51	#include <TF1.h>
	52	#include <TMath.h>
	53
	54	#include "AliEMCALPIDResponse.h" //class header
	55
	56	#include "AliLog.h"
	57
	58	ClassImp(AliEMCALPIDResponse)
	59
	60	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	61	AliEMCALPIDResponse::AliEMCALPIDResponse():
	62	TObject(),
	63	fNorm(NULL),
87da0205	64	fCurrCentrality(-1.),
b2138b40	65	fkPIDParams(NULL)
	66	{
	67	//
	68	// The default constructor
	69	//
	70
	71
	72	fNorm = new TF1("fNorm","gaus",-20,20);
	73	}
	74	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	75	AliEMCALPIDResponse::AliEMCALPIDResponse(const AliEMCALPIDResponse &other):
	76	TObject(other),
	77	fNorm(other.fNorm),
87da0205	78	fCurrCentrality(other.fCurrCentrality),
b2138b40	79	fkPIDParams(other.fkPIDParams)
	80	{
	81	//
	82	// The copy constructor
	83	//
	84
	85	}
	86	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	87	AliEMCALPIDResponse & AliEMCALPIDResponse::operator=( const AliEMCALPIDResponse& other)
	88	{
	89	//
	90	// The assignment operator
	91	//
	92
	93	if(this == &other) return *this;
	94
	95	// Make copy
	96	TObject::operator=(other);
	97	fNorm = other.fNorm;
87da0205	98	fCurrCentrality = other.fCurrCentrality;
b2138b40	99	fkPIDParams = other.fkPIDParams;
	100
	101
	102	return *this;
	103	}
	104	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	105	AliEMCALPIDResponse::~AliEMCALPIDResponse() {
	106
	107	delete fNorm;
	108
	109	}
	110	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	111	Double_t AliEMCALPIDResponse::GetExpectedNorm( Float_t pt, AliPID::EParticleType n, Int_t charge) const {
	112	//
	113	// Calculates the expected sigma of the PID signal as the function of
	114	// the information stored in the track, for the specified particle type
	115	//
	116	//
	117
	118	Double_t norm = 1.;
	119
	120	// Check the charge
	121	if( charge != -1 && charge != 1){
	122	return norm;
	123	}
	124
	125	// Get the parameters for this particle type and pt
1f631618	126	const TVectorD *params = GetParams(n, pt, charge);
b2138b40	127
	128	// IF not in momentum range, NULL is returned --> return default value
	129	if(!params) return norm;
	130
	131	Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
	132	Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
	133	Double_t eopMin = (*params)[4]; // min E/p value for parametrization
	134	Double_t eopMax = (*params)[5]; // max E/p value for parametrization
	135	Double_t probLow = (*params)[6]; // probability to be below eopMin
	136	Double_t probHigh = (*params)[7]; // probability to be above eopMax
	137
	138	// Get the normalization factor ( Probability in the parametrized area / Integral of parametrized Gauss function in this area )
	139	fNorm->SetParameters(1./TMath::Sqrt(2TMath::Pi()sigma*sigma),mean,sigma);
	140	norm = 1./fNorm->Integral(eopMin,eopMax)*(1-probLow-probHigh);
	141
	142	return norm;
	143	}
	144	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	145	Double_t AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt, Float_t eop, AliPID::EParticleType n, Int_t charge) const {
	146
	147	Double_t nsigma = -999.;
	148
	149	// Check the charge
	150	if( charge != -1 && charge != 1){
	151	return nsigma;
	152	}
	153
	154	// Get the parameters for this particle type and pt
1f631618	155	const TVectorD *params = GetParams(n, pt, charge);
b2138b40	156
	157	// IF not in momentum range, NULL is returned --> return default value
	158	if(!params) return nsigma;
	159
	160	Double_t mean = (*params)[2]; // mean value of Gausiian parametrization
	161	Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
	162	Double_t eopMin = (*params)[4]; // min E/p value for parametrization
	163	Double_t eopMax = (*params)[5]; // max E/p value for parametrization
	164
	165	// if electron
	166	if(n == AliPID::kElectron){
	167	if(sigma != 0) nsigma = (eop - mean) / sigma;
	168	}
	169
	170	// if NON electron
	171	else{
	172	if ( eop < eopMin )
	173	nsigma = -99; // not parametrized (smaller than eopMin)
	174	else if ( eop > eopMax )
	175	nsigma = 99.; // not parametrized (bigger than eopMax)
	176	else{
	177	if(sigma != 0) nsigma = (eop - mean) / sigma;
	178	}
	179	}
	180
	181	return nsigma;
	182
	183	}
	184	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00a38d07	185	Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
b2138b40	186	//
	187	//
	188	Double_t fRange = 5.0; // hardcoded (???)
	189	Double_t nsigma = 0.0;
b2138b40	190
	191
	192	// Check the charge
	193	if( charge != -1 && charge != 1){
00a38d07	194	return kFALSE;
b2138b40	195	}
	196
	197
	198	// default value (will be returned, if pt below threshold)
00a38d07	199	for (Int_t species = 0; species < nSpecies; species++) {
00a38d07	200	pEMCAL[species] = 1./nSpecies;
b2138b40	201	}
	202
	203	// set E/p range
	204	if(eop < 0.05) eop = 0.05;
	205	if(eop > 2.00) eop = 2.00;
	206
00a38d07	207	for (Int_t species = 0; species < nSpecies; species++) {
b2138b40	208
	209	AliPID::EParticleType type = AliPID::EParticleType(species);
	210
	211	// Get the parameters for this particle type and pt
1f631618	212	const TVectorD *params = GetParams(species, pt, charge);
b2138b40	213
00a38d07	214	// IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
00a38d07	215	if(!params) return kFALSE;
b2138b40	216
	217	Double_t sigma = (*params)[3]; // sigma value of Gausiian parametrization
	218	Double_t probLow = (*params)[6]; // probability to be below eopMin
	219	Double_t probHigh = (*params)[7]; // probability to be above eopMax
	220
	221	// get nsigma value for each particle type at this E/p value
	222	nsigma = GetNumberOfSigmas(pt,eop,type,charge);
	223
	224	// electrons (standard Gaussian calculation of probabilities)
	225	if(type == AliPID::kElectron){
	226	if (TMath::Abs(nsigma) > fRange) {
	227	pEMCAL[species]=TMath::Exp(-0.5fRangefRange)/TMath::Sqrt(2TMath::Pi()sigma*sigma);
	228	}
	229	else{
	230	pEMCAL[species]=TMath::Exp(-0.5(nsigma)(nsigma))/TMath::Sqrt(2TMath::Pi()sigma*sigma);
	231	}
	232	}
	233	//NON electrons
	234	else{
	235	// E/p < eopMin --> return probability below E/p = eopMin
	236	if ( nsigma == -99){
	237	pEMCAL[species] = probLow;
	238	}
	239	// E/p > eopMax --> return probability above E/p = eopMax
	240	else if ( nsigma == 99){
	241	pEMCAL[species] = probHigh;
	242	}
	243	// in parametrized region --> calculate probability for corresponding Gauss curve
	244	else{
	245	pEMCAL[species]=TMath::Exp(-0.5(nsigma)(nsigma))/TMath::Sqrt(2TMath::Pi()sigma*sigma);
	246
	247	// normalize to total probability == 1
	248	pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
	249	}
	250	}
b2138b40	251	}
b2138b40	252
00a38d07	253	return kTRUE;
b2138b40	254
	255	}
	256	//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1f631618	257	const TVectorD* AliEMCALPIDResponse::GetParams(Int_t nParticle, Float_t fPt, Int_t charge) const {
b2138b40	258	//
	259	// returns the PID parameters (mean, sigma, probabilities for Hadrons) for a certain particle and pt
	260	//
	261	// 0 = momMin
	262	// 1 = momMax
	263	// 2 = mean of Gaus
	264	// 3 = sigma of Gaus
	265	// 4 = eopLow
	266	// 5 = eopHig
	267	// 6 = probLow (not used for electrons)
	268	// 7 = probHigh (not used for electrons)
	269	//
87da0205	270	// for PbPb the parametrization is done centrality dependent (marked by TString "Centrality")
	271	// so first the correct centrality bin has to be found
	272
	273	// **** Centrality bins (hard coded for the moment)
	274	const Int_t nCent = 7;
	275	Int_t centBins[nCent+1] = {0,10,20,30,40,50,70,90};
b2138b40	276
b2138b40	277	if(nParticle > AliPID::kSPECIES \|\| nParticle <0) return NULL;
1f631618	278	if(nParticle == AliPID::kProton && charge == -1) nParticle = AliPID::kSPECIES; // special case for antiprotons
b2138b40	279
	280	TObjArray * particlePar = dynamic_cast<TObjArray *>(fkPIDParams->At(nParticle));
	281	if(!particlePar) return NULL;
1f631618	282
b2138b40	283	const TVectorD *parameters = NULL;
	284	Double_t momMin = 0.;
	285	Double_t momMax = 0.;
	286
87da0205	287	// is the centrality dependent parametrization used
	288	TString arrayName = particlePar->GetName();
	289
	290	// centrality dependent parametrization
	291	if(arrayName.Contains("Centrality")){
	292
	293	for(Int_t iCent = 0; iCent < nCent; iCent++){
b2138b40	294
87da0205	295	if( fCurrCentrality > centBins[iCent] && fCurrCentrality < centBins[iCent+1] ){
	296
	297	TObjArray * centPar = dynamic_cast<TObjArray *>(particlePar->At(iCent));
	298	if(!centPar) return NULL;
	299
	300	TIter centIter(centPar);
	301	parameters = NULL;
	302	momMin = 0.;
	303	momMax = 0.;
	304
	305	while((parameters = static_cast<const TVectorD *>(centIter()))){
	306
	307	momMin = (*parameters)[0];
	308	momMax = (*parameters)[1];
	309
	310	if( fPt > momMin && fPt < momMax ) return parameters;
	311
	312	}
	313	}
	314	}
	315	}
b2138b40	316
87da0205	317	// NO centrality dependent parametrization
87da0205	318	else{
b2138b40	319
87da0205	320	TIter parIter(particlePar);
	321	while((parameters = static_cast<const TVectorD *>(parIter()))){
	322
	323	momMin = (*parameters)[0];
	324	momMax = (*parameters)[1];
	325
	326	if( fPt > momMin && fPt < momMax ) return parameters;
	327
	328	}
	329	}
b2138b40	330	AliDebug(2, Form("NO params for particle %d and momentum %f \n", nParticle, fPt));
	331
	332	return parameters;
	333	}