[u/mrichter/AliRoot.git] / STEER / STEERBase / AliITSPidParams.cxx

/**************************************************************************
 * Copyright(c) 2007-2009, 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: */

///////////////////////////////////////////////////////////////////
//                                                               //
// Implementation of the class to store parameters of ITS        //
// response funcions for dE/dx based PID                         //
// Origin: F.Prino, Torino, prino@to.infn.it                     //
//                                                               //
///////////////////////////////////////////////////////////////////

#include <TFormula.h>
#include <TNamed.h>
#include <TMath.h>
#include "AliITSPidParams.h"

ClassImp(AliITSPidParams)

//______________________________________________________________________
AliITSPidParams::AliITSPidParams(Bool_t isMC):
TNamed("default",""),
  fSDDElecLandauWidth(0),
  fSDDElecGaussWidth(0),
  fSSDElecLandauWidth(0),
  fSSDElecGaussWidth(0),
  fSDDPionLandauWidth(0),
  fSDDPionGaussWidth(0),
  fSSDPionLandauWidth(0),
  fSSDPionGaussWidth(0),
  fSDDKaonLandauWidth(0),
  fSDDKaonGaussWidth(0),
  fSSDKaonLandauWidth(0),
  fSSDKaonGaussWidth(0),
  fSDDProtLandauWidth(0),
  fSDDProtGaussWidth(0),
  fSSDProtLandauWidth(0),
  fSSDProtGaussWidth(0)
{
  // default constructor
  if (isMC) InitMC();
  else InitData();
}
//______________________________________________________________________
AliITSPidParams::AliITSPidParams(Char_t * name, Bool_t isMC):
  TNamed(name,""),
  fSDDElecLandauWidth(0),
  fSDDElecGaussWidth(0),
  fSSDElecLandauWidth(0),
  fSSDElecGaussWidth(0),
  fSDDPionLandauWidth(0),
  fSDDPionGaussWidth(0),
  fSSDPionLandauWidth(0),
  fSSDPionGaussWidth(0),
  fSDDKaonLandauWidth(0),
  fSDDKaonGaussWidth(0),
  fSSDKaonLandauWidth(0),
  fSSDKaonGaussWidth(0),
  fSDDProtLandauWidth(0),
  fSDDProtGaussWidth(0),
  fSSDProtLandauWidth(0), 
  fSSDProtGaussWidth(0)
{
  // standard constructor
  if (isMC) InitMC();
  else InitData();
}
//______________________________________________________________________
AliITSPidParams::~AliITSPidParams(){
  //
  if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
  if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
  
  if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
  if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;

  if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
  if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;

  if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
  if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
  
  if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
  if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
  
  if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
  if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;

  if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
  if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;

  if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
  if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
}
//______________________________________________________________________
Double_t AliITSPidParams::BetheBloch(Double_t mom, Double_t mass, const Double_t *p) const{
 //
  // This is the empirical parameterization of the Bethe-Bloch formula.
  // It is normalized to 1 at the minimum.
  //
  // bg - beta*gamma
  // 
  Double_t bg = mom/mass;
  Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
  Double_t aa = TMath::Power(beta,p[3]);
  Double_t bb = TMath::Power(1./bg,p[4]);
  bb=TMath::Log(p[2]+bb);
  return (p[0]-aa-bb)*p[1]/aa;
}
//______________________________________________________________________
Double_t AliITSPidParams::ExtrapolateWidth(Double_t mom, Double_t x1, Double_t y1,Double_t x2, Double_t y2) const{
  //
  // This is a linear extrapolation of Landau width and Gaussian width 
  // for low momentum.
  
 Double_t slope = (y2-y1)/(x2-x1);
 return slope*mom+(y1-slope*x1);
}//______________________________________________________________________
void AliITSPidParams::InitMC(){
  // initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
  // parameter values from LHC10d1 
  // MPV BetheBloch parameters;
  
  //sdd MC electrons parameters
  fSDDElecMPVBetheParams[0] = -0.0931934;
  fSDDElecMPVBetheParams[1] = 77.8422;
  fSDDElecMPVBetheParams[2] = -0.889085;
  fSDDElecMPVBetheParams[3] = -154.455;
  fSDDElecMPVBetheParams[4] = -0.000256076;
  
  //ssd MC electrons parameters
  fSSDElecMPVBetheParams[0] = -0.0989358;
  fSSDElecMPVBetheParams[1] = 77.8271;
  fSSDElecMPVBetheParams[2] = -0.900887;
  fSSDElecMPVBetheParams[3] =-1241.14;
  fSSDElecMPVBetheParams[4] = -0.0014204;
  
   // electrons 
  if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
  fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
  fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);

  if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
  fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
  fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);

  if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
  fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
  fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);

  if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
  fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
  fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
  
  //sdd MC hadrons parameters
  fSDDHadronMPVBetheParams[0] = 1.13531;
  fSDDHadronMPVBetheParams[1] = -156.651;
  fSDDHadronMPVBetheParams[2] = 1.87562;
  fSDDHadronMPVBetheParams[3] = 0.45819;
  fSDDHadronMPVBetheParams[4] = 2.26506;
  
  //ssd MC hadrons parameters
  fSSDHadronMPVBetheParams[0] = -0.451908;
  fSSDHadronMPVBetheParams[1] = -55.4368;
  fSSDHadronMPVBetheParams[2] = 0.984636;
  fSSDHadronMPVBetheParams[3] = 0.97078;
  fSSDHadronMPVBetheParams[4] = 2.50883;

  // pions 
  if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
  fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
  fSDDPionLandauWidth->SetParameters(0.08026, 5.87922);

  if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
  fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSDDPionGaussWidth->SetParameters(-0.090495, 7.705286);

  if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
  fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
  fSSDPionLandauWidth->SetParameters(0.083882, 5.823419);

  if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
  fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSSDPionGaussWidth->SetParameters(-0.105218, 5.650956);

  // kaons
  if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
  fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
  fSDDKaonLandauWidth->SetParameters(1.430692, 5.581389);

  if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
  fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSDDKaonGaussWidth->SetParameters(-1.6537, 8.071832);

  if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
  fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
  fSSDKaonLandauWidth->SetParameters(1.368824, 5.639291);

  if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
  fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSSDKaonGaussWidth->SetParameters(-1.901858, 6.0932281);

  // protons
  if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
  fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
  fSDDProtLandauWidth->SetParameters(6.529418, 5.049098);

  if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
  fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSDDProtGaussWidth->SetParameters(-9.360599, 11.318026);

  if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
  fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
  fSSDProtLandauWidth->SetParameters(6.419493, 4.925070);

  if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
  fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSSDProtGaussWidth->SetParameters(-9.599064, 9.358656);
}
//______________________________________________________________________
void AliITSPidParams::InitData(){
  // initialize TFormulas to Real Data values (=p-p ALICE Experiment)
  // parameter values from LHC10b 
  // MPV BetheBloch parameters;
 
 //sdd data electrons parameters
  fSDDElecMPVBetheParams[0] = -0.130204;
  fSDDElecMPVBetheParams[1] = 75.1267;
  fSDDElecMPVBetheParams[2] = -0.889337;
  fSDDElecMPVBetheParams[3] = 388.372;
  fSDDElecMPVBetheParams[4] = 0.00134649;
  
  //ssd data electrons parameters
  fSSDElecMPVBetheParams[0] = -0.162773;
  fSSDElecMPVBetheParams[1] = 72.9393;
  fSSDElecMPVBetheParams[2] = -0.896944;
  fSSDElecMPVBetheParams[3] = 3233.02;
  fSSDElecMPVBetheParams[4] = 0.00146896;
  
   // electrons 
  if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
  fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
  fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);

  if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
  fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
  fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);

  if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
  fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
  fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);

  if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
  fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
  fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
  
  //sdd data hadrons parameters
  fSDDHadronMPVBetheParams[0] = -18.1867;
  fSDDHadronMPVBetheParams[1] = -3.45806;
  fSDDHadronMPVBetheParams[2] =  2.23635;
  fSDDHadronMPVBetheParams[3] =  2.08328;
  fSDDHadronMPVBetheParams[4] = -1.92331;
 
  //ssd data hadrons parameters
  fSSDHadronMPVBetheParams[0] = -12.6459;
  fSSDHadronMPVBetheParams[1] = -4.84598;
  fSSDHadronMPVBetheParams[2] =  1.50253;
  fSSDHadronMPVBetheParams[3] =  2.08328;
  fSSDHadronMPVBetheParams[4] = -1.3719;
  
  // pions
  if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
  fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
  fSDDPionLandauWidth->SetParameters(0.07694, 5.468728);

  if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
  fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSDDPionGaussWidth->SetParameters(-0.098209, 10.409441);

  if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
  fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
  fSSDPionLandauWidth->SetParameters(0.071602, 5.365442);

  if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
  fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(x*x)*TMath::Log(x)+[1]");
  fSSDPionGaussWidth->SetParameters(-0.098045, 7.617583);

  // kaons
  if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
  fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
  fSDDKaonLandauWidth->SetParameters(0.998191, 5.461668);

  if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
  fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]");
  fSDDKaonGaussWidth->SetParameters(1.629308, 9.851873);

  if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
  fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
  fSSDKaonLandauWidth->SetParameters(0.773113, 5.618683);

  if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
  fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]");
  fSSDKaonGaussWidth->SetParameters(1.510713, 6.862774);

  // protons
  if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
  fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
  fSDDProtLandauWidth->SetParameters(3.561429, 5.372105);

  if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
  fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]");
  fSDDProtGaussWidth->SetParameters(5.395926, 10.044613);

  if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
  fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
  fSSDProtLandauWidth->SetParameters(2.647428, 5.678460);

  if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
  fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
  fSSDProtGaussWidth->SetParameters(4.91025, 6.779763);
}
//_______________________________________________________________________
Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
  // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
  if(TMath::Abs(pdgCode)==11) return GetLandauGausNorm(dedx,AliPID::kElectron,mom,lay);
  else if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
  else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
  else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
  else return 0.;
}
//_______________________________________________________________________
Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
  // Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
  
  Double_t par[4];
  Bool_t isSet=kFALSE;
  if(partType==AliPID::kElectron){
    if(lay==3 || lay==4){
      par[0]=GetSDDElecLandauWidth(mom);
      par[1]=GetSDDElecMPV(mom);
      par[2]=GetSDDElecGaussWidth(mom);
      isSet=kTRUE;
    }
    else if(lay==5 || lay==6){
      par[0]=GetSSDElecLandauWidth(mom);
      par[1]=GetSSDElecMPV(mom);
      par[2]=GetSSDElecGaussWidth(mom);
      isSet=kTRUE;
    }
  }else if(partType==AliPID::kPion){
    if(lay==3 || lay==4){
      par[0]=GetSDDPionLandauWidth(mom);
      par[1]=GetSDDPionMPV(mom);
      par[2]=GetSDDPionGaussWidth(mom);
      isSet=kTRUE;
    }
    else if(lay==5 || lay==6){
      par[0]=GetSSDPionLandauWidth(mom);
      par[1]=GetSSDPionMPV(mom);
      par[2]=GetSSDPionGaussWidth(mom);
      isSet=kTRUE;
    }
  }else if(partType==AliPID::kKaon){
    if(lay==3 || lay==4){
      par[0]=GetSDDKaonLandauWidth(mom);
      par[1]=GetSDDKaonMPV(mom);
      par[2]=GetSDDKaonGaussWidth(mom);
      isSet=kTRUE;
    }
    else if(lay==5 || lay==6){
      par[0]=GetSSDKaonLandauWidth(mom);
      par[1]=GetSSDKaonMPV(mom);
      par[2]=GetSSDKaonGaussWidth(mom);
      isSet=kTRUE;
    }
  }else if(partType==AliPID::kProton){
    if(lay==3 || lay==4){
      par[0]=GetSDDProtLandauWidth(mom);
      par[1]=GetSDDProtMPV(mom);
      par[2]=GetSDDProtGaussWidth(mom);
      isSet=kTRUE;
    }
    else if(lay==5 || lay==6){
      par[0]=GetSSDProtLandauWidth(mom);
      par[1]=GetSSDProtMPV(mom);
      par[2]=GetSSDProtGaussWidth(mom);
      isSet=kTRUE;
    }
  }
  if(!isSet) return 0.;
  // Numeric constants
  const Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
  const Double_t mpshift  = -0.22278298;       // Landau maximum location
  // Control constants
  const Double_t np = 100.0;      // number of convolution steps
  const 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 = dedx - sc * par[2];
  xupp = dedx + 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(dedx,xx,par[2]);
    
    xx = xupp - (i-.5) * step;
    fland = TMath::Landau(xx,mpc,par[0]) / par[0];
    sum += fland * TMath::Gaus(dedx,xx,par[2]);
  }
  
  return (step * sum * invsq2pi / par[2]);
}
Commit	Line	Data
b536a002	1	/**************************************************************************
	2	* Copyright(c) 2007-2009, 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
	18	///////////////////////////////////////////////////////////////////
	19	// //
	20	// Implementation of the class to store parameters of ITS //
	21	// response funcions for dE/dx based PID //
	22	// Origin: F.Prino, Torino, prino@to.infn.it //
	23	// //
	24	///////////////////////////////////////////////////////////////////
	25
	26	#include <TFormula.h>
	27	#include <TNamed.h>
	28	#include <TMath.h>
	29	#include "AliITSPidParams.h"
b536a002	30
	31	ClassImp(AliITSPidParams)
	32
	33	//______________________________________________________________________
b52bfc67	34	AliITSPidParams::AliITSPidParams(Bool_t isMC):
b52bfc67	35	TNamed("default",""),
2ca1f4ee	36	fSDDElecLandauWidth(0),
	37	fSDDElecGaussWidth(0),
	38	fSSDElecLandauWidth(0),
	39	fSSDElecGaussWidth(0),
b536a002	40	fSDDPionLandauWidth(0),
b536a002	41	fSDDPionGaussWidth(0),
b536a002	42	fSSDPionLandauWidth(0),
b536a002	43	fSSDPionGaussWidth(0),
b536a002	44	fSDDKaonLandauWidth(0),
b536a002	45	fSDDKaonGaussWidth(0),
b536a002	46	fSSDKaonLandauWidth(0),
b536a002	47	fSSDKaonGaussWidth(0),
b536a002	48	fSDDProtLandauWidth(0),
b536a002	49	fSDDProtGaussWidth(0),
b536a002	50	fSSDProtLandauWidth(0),
	51	fSSDProtGaussWidth(0)
	52	{
	53	// default constructor
b52bfc67	54	if (isMC) InitMC();
b52bfc67	55	else InitData();
b536a002	56	}
b536a002	57	//______________________________________________________________________
b52bfc67	58	AliITSPidParams::AliITSPidParams(Char_t * name, Bool_t isMC):
b536a002	59	TNamed(name,""),
2ca1f4ee	60	fSDDElecLandauWidth(0),
	61	fSDDElecGaussWidth(0),
	62	fSSDElecLandauWidth(0),
	63	fSSDElecGaussWidth(0),
63a0c396	64	fSDDPionLandauWidth(0),
63a0c396	65	fSDDPionGaussWidth(0),
b536a002	66	fSSDPionLandauWidth(0),
b536a002	67	fSSDPionGaussWidth(0),
b536a002	68	fSDDKaonLandauWidth(0),
b536a002	69	fSDDKaonGaussWidth(0),
b536a002	70	fSSDKaonLandauWidth(0),
b536a002	71	fSSDKaonGaussWidth(0),
b536a002	72	fSDDProtLandauWidth(0),
b536a002	73	fSDDProtGaussWidth(0),
b52bfc67	74	fSSDProtLandauWidth(0),
b536a002	75	fSSDProtGaussWidth(0)
	76	{
	77	// standard constructor
b52bfc67	78	if (isMC) InitMC();
b52bfc67	79	else InitData();
b536a002	80	}
	81	//______________________________________________________________________
	82	AliITSPidParams::~AliITSPidParams(){
2ca1f4ee	83	//
	84	if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
	85	if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
	86
	87	if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
	88	if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
	89
b536a002	90	if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
	91	if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
	92
b536a002	93	if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
b536a002	94	if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
b52bfc67	95
b536a002	96	if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
b536a002	97	if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
b52bfc67	98
b536a002	99	if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
	100	if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
	101
b536a002	102	if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
	103	if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
	104
b536a002	105	if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
	106	if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
	107	}
b536a002	108	//______________________________________________________________________
2ca1f4ee	109	Double_t AliITSPidParams::BetheBloch(Double_t mom, Double_t mass, const Double_t *p) const{
	110	//
	111	// This is the empirical parameterization of the Bethe-Bloch formula.
	112	// It is normalized to 1 at the minimum.
	113	//
	114	// bg - beta*gamma
	115	//
	116	Double_t bg = mom/mass;
	117	Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
	118	Double_t aa = TMath::Power(beta,p[3]);
	119	Double_t bb = TMath::Power(1./bg,p[4]);
	120	bb=TMath::Log(p[2]+bb);
	121	return (p[0]-aa-bb)*p[1]/aa;
	122	}
	123	//______________________________________________________________________
	124	Double_t AliITSPidParams::ExtrapolateWidth(Double_t mom, Double_t x1, Double_t y1,Double_t x2, Double_t y2) const{
	125	//
	126	// This is a linear extrapolation of Landau width and Gaussian width
	127	// for low momentum.
	128
	129	Double_t slope = (y2-y1)/(x2-x1);
	130	return slopemom+(y1-slopex1);
	131	}//______________________________________________________________________
8abeb05b	132	void AliITSPidParams::InitMC(){
	133	// initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
	134	// parameter values from LHC10d1
2ca1f4ee	135	// MPV BetheBloch parameters;
	136
	137	//sdd MC electrons parameters
	138	fSDDElecMPVBetheParams[0] = -0.0931934;
	139	fSDDElecMPVBetheParams[1] = 77.8422;
	140	fSDDElecMPVBetheParams[2] = -0.889085;
	141	fSDDElecMPVBetheParams[3] = -154.455;
	142	fSDDElecMPVBetheParams[4] = -0.000256076;
	143
	144	//ssd MC electrons parameters
	145	fSSDElecMPVBetheParams[0] = -0.0989358;
	146	fSSDElecMPVBetheParams[1] = 77.8271;
	147	fSSDElecMPVBetheParams[2] = -0.900887;
	148	fSSDElecMPVBetheParams[3] =-1241.14;
	149	fSSDElecMPVBetheParams[4] = -0.0014204;
	150
	151	// electrons
	152	if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
	153	fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
	154	fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);
	155
	156	if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
	157	fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
	158	fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);
	159
	160	if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
	161	fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
	162	fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);
	163
	164	if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
	165	fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
	166	fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
	167
	168	//sdd MC hadrons parameters
	169	fSDDHadronMPVBetheParams[0] = 1.13531;
	170	fSDDHadronMPVBetheParams[1] = -156.651;
	171	fSDDHadronMPVBetheParams[2] = 1.87562;
	172	fSDDHadronMPVBetheParams[3] = 0.45819;
	173	fSDDHadronMPVBetheParams[4] = 2.26506;
	174
	175	//ssd MC hadrons parameters
	176	fSSDHadronMPVBetheParams[0] = -0.451908;
	177	fSSDHadronMPVBetheParams[1] = -55.4368;
	178	fSSDHadronMPVBetheParams[2] = 0.984636;
	179	fSSDHadronMPVBetheParams[3] = 0.97078;
	180	fSSDHadronMPVBetheParams[4] = 2.50883;
	181
	182	// pions
b536a002	183	if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
8abeb05b	184	fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	185	fSDDPionLandauWidth->SetParameters(0.08026, 5.87922);
b536a002	186
b536a002	187	if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
8abeb05b	188	fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	189	fSDDPionGaussWidth->SetParameters(-0.090495, 7.705286);
b536a002	190
b536a002	191	if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
8abeb05b	192	fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	193	fSSDPionLandauWidth->SetParameters(0.083882, 5.823419);
b536a002	194
b536a002	195	if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
8abeb05b	196	fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	197	fSSDPionGaussWidth->SetParameters(-0.105218, 5.650956);
b536a002	198
b536a002	199	// kaons
b536a002	200	if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
8abeb05b	201	fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	202	fSDDKaonLandauWidth->SetParameters(1.430692, 5.581389);
b52bfc67	203
	204	if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
	205	fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	206	fSDDKaonGaussWidth->SetParameters(-1.6537, 8.071832);
b52bfc67	207
	208	if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
	209	fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	210	fSSDKaonLandauWidth->SetParameters(1.368824, 5.639291);
b52bfc67	211
	212	if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
	213	fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	214	fSSDKaonGaussWidth->SetParameters(-1.901858, 6.0932281);
b52bfc67	215
b52bfc67	216	// protons
b52bfc67	217	if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
b52bfc67	218	fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	219	fSDDProtLandauWidth->SetParameters(6.529418, 5.049098);
b52bfc67	220
	221	if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
	222	fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	223	fSDDProtGaussWidth->SetParameters(-9.360599, 11.318026);
b52bfc67	224
	225	if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
	226	fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	227	fSSDProtLandauWidth->SetParameters(6.419493, 4.925070);
b52bfc67	228
	229	if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
	230	fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	231	fSSDProtGaussWidth->SetParameters(-9.599064, 9.358656);
b52bfc67	232	}
	233	//______________________________________________________________________
	234	void AliITSPidParams::InitData(){
2ca1f4ee	235	// initialize TFormulas to Real Data values (=p-p ALICE Experiment)
b52bfc67	236	// parameter values from LHC10b
2ca1f4ee	237	// MPV BetheBloch parameters;
	238
	239	//sdd data electrons parameters
	240	fSDDElecMPVBetheParams[0] = -0.130204;
	241	fSDDElecMPVBetheParams[1] = 75.1267;
	242	fSDDElecMPVBetheParams[2] = -0.889337;
	243	fSDDElecMPVBetheParams[3] = 388.372;
	244	fSDDElecMPVBetheParams[4] = 0.00134649;
	245
	246	//ssd data electrons parameters
	247	fSSDElecMPVBetheParams[0] = -0.162773;
	248	fSSDElecMPVBetheParams[1] = 72.9393;
	249	fSSDElecMPVBetheParams[2] = -0.896944;
	250	fSSDElecMPVBetheParams[3] = 3233.02;
	251	fSSDElecMPVBetheParams[4] = 0.00146896;
	252
	253	// electrons
	254	if(fSDDElecLandauWidth) delete fSDDElecLandauWidth;
	255	fSDDElecLandauWidth=new TFormula("fSDDElecLandauWidth","[0]/(x*x)+[1]");
	256	fSDDElecLandauWidth->SetParameters(-0.002702, 6.224960);
	257
	258	if(fSDDElecGaussWidth) delete fSDDElecGaussWidth;
	259	fSDDElecGaussWidth=new TFormula("fSDDElecGaussWidth","[0]/(x*x)+[1]");
	260	fSDDElecGaussWidth->SetParameters(0.012402, 7.993106);
	261
	262	if(fSSDElecLandauWidth) delete fSSDElecLandauWidth;
	263	fSSDElecLandauWidth=new TFormula("fSSDElecLandauWidth","[0]/(x*x)+[1]");
	264	fSSDElecLandauWidth->SetParameters(-0.002144, 6.231089);
	265
	266	if(fSSDElecGaussWidth) delete fSSDElecGaussWidth;
	267	fSSDElecGaussWidth=new TFormula("fSSDElecGaussWidth","[0]/(x*x)+[1]");
	268	fSSDElecGaussWidth->SetParameters(0.014530, 6.217153);
	269
	270	//sdd data hadrons parameters
	271	fSDDHadronMPVBetheParams[0] = -18.1867;
	272	fSDDHadronMPVBetheParams[1] = -3.45806;
	273	fSDDHadronMPVBetheParams[2] = 2.23635;
	274	fSDDHadronMPVBetheParams[3] = 2.08328;
	275	fSDDHadronMPVBetheParams[4] = -1.92331;
	276
	277	//ssd data hadrons parameters
	278	fSSDHadronMPVBetheParams[0] = -12.6459;
	279	fSSDHadronMPVBetheParams[1] = -4.84598;
	280	fSSDHadronMPVBetheParams[2] = 1.50253;
	281	fSSDHadronMPVBetheParams[3] = 2.08328;
	282	fSSDHadronMPVBetheParams[4] = -1.3719;
b52bfc67	283
b52bfc67	284	// pions
b52bfc67	285	if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
b52bfc67	286	fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	287	fSDDPionLandauWidth->SetParameters(0.07694, 5.468728);
b52bfc67	288
	289	if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
	290	fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	291	fSDDPionGaussWidth->SetParameters(-0.098209, 10.409441);
b52bfc67	292
	293	if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
	294	fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	295	fSSDPionLandauWidth->SetParameters(0.071602, 5.365442);
b52bfc67	296
	297	if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
	298	fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
2ca1f4ee	299	fSSDPionGaussWidth->SetParameters(-0.098045, 7.617583);
b52bfc67	300
b52bfc67	301	// kaons
b52bfc67	302	if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
b52bfc67	303	fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	304	fSDDKaonLandauWidth->SetParameters(0.998191, 5.461668);
b536a002	305
b536a002	306	if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
8abeb05b	307	fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]");
2ca1f4ee	308	fSDDKaonGaussWidth->SetParameters(1.629308, 9.851873);
b536a002	309
b536a002	310	if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
8abeb05b	311	fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	312	fSSDKaonLandauWidth->SetParameters(0.773113, 5.618683);
b536a002	313
b536a002	314	if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
8abeb05b	315	fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]");
2ca1f4ee	316	fSSDKaonGaussWidth->SetParameters(1.510713, 6.862774);
b536a002	317
b536a002	318	// protons
b536a002	319	if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
8abeb05b	320	fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	321	fSDDProtLandauWidth->SetParameters(3.561429, 5.372105);
b536a002	322
b536a002	323	if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
8abeb05b	324	fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]");
2ca1f4ee	325	fSDDProtGaussWidth->SetParameters(5.395926, 10.044613);
8abeb05b	326
b536a002	327	if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
8abeb05b	328	fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
2ca1f4ee	329	fSSDProtLandauWidth->SetParameters(2.647428, 5.678460);
b536a002	330
b536a002	331	if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
8abeb05b	332	fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
2ca1f4ee	333	fSSDProtGaussWidth->SetParameters(4.91025, 6.779763);
b536a002	334	}
	335	//_______________________________________________________________________
	336	Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
	337	// Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
2ca1f4ee	338	if(TMath::Abs(pdgCode)==11) return GetLandauGausNorm(dedx,AliPID::kElectron,mom,lay);
2ca1f4ee	339	else if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
b536a002	340	else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
	341	else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
	342	else return 0.;
	343	}
	344	//_______________________________________________________________________
b52bfc67	345	Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
b536a002	346	// Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
b52bfc67	347
2ca1f4ee	348	Double_t par[4];
b536a002	349	Bool_t isSet=kFALSE;
2ca1f4ee	350	if(partType==AliPID::kElectron){
	351	if(lay==3 \|\| lay==4){
	352	par[0]=GetSDDElecLandauWidth(mom);
	353	par[1]=GetSDDElecMPV(mom);
	354	par[2]=GetSDDElecGaussWidth(mom);
	355	isSet=kTRUE;
	356	}
	357	else if(lay==5 \|\| lay==6){
	358	par[0]=GetSSDElecLandauWidth(mom);
	359	par[1]=GetSSDElecMPV(mom);
	360	par[2]=GetSSDElecGaussWidth(mom);
	361	isSet=kTRUE;
	362	}
	363	}else if(partType==AliPID::kPion){
b536a002	364	if(lay==3 \|\| lay==4){
	365	par[0]=GetSDDPionLandauWidth(mom);
	366	par[1]=GetSDDPionMPV(mom);
	367	par[2]=GetSDDPionGaussWidth(mom);
	368	isSet=kTRUE;
	369	}
	370	else if(lay==5 \|\| lay==6){
	371	par[0]=GetSSDPionLandauWidth(mom);
	372	par[1]=GetSSDPionMPV(mom);
	373	par[2]=GetSSDPionGaussWidth(mom);
	374	isSet=kTRUE;
	375	}
	376	}else if(partType==AliPID::kKaon){
	377	if(lay==3 \|\| lay==4){
	378	par[0]=GetSDDKaonLandauWidth(mom);
	379	par[1]=GetSDDKaonMPV(mom);
	380	par[2]=GetSDDKaonGaussWidth(mom);
	381	isSet=kTRUE;
	382	}
	383	else if(lay==5 \|\| lay==6){
	384	par[0]=GetSSDKaonLandauWidth(mom);
	385	par[1]=GetSSDKaonMPV(mom);
	386	par[2]=GetSSDKaonGaussWidth(mom);
	387	isSet=kTRUE;
	388	}
	389	}else if(partType==AliPID::kProton){
	390	if(lay==3 \|\| lay==4){
	391	par[0]=GetSDDProtLandauWidth(mom);
	392	par[1]=GetSDDProtMPV(mom);
	393	par[2]=GetSDDProtGaussWidth(mom);
	394	isSet=kTRUE;
	395	}
	396	else if(lay==5 \|\| lay==6){
	397	par[0]=GetSSDProtLandauWidth(mom);
	398	par[1]=GetSSDProtMPV(mom);
	399	par[2]=GetSSDProtGaussWidth(mom);
	400	isSet=kTRUE;
	401	}
	402	}
	403	if(!isSet) return 0.;
	404	// Numeric constants
63a0c396	405	const Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
63a0c396	406	const Double_t mpshift = -0.22278298; // Landau maximum location
b536a002	407	// Control constants
63a0c396	408	const Double_t np = 100.0; // number of convolution steps
63a0c396	409	const Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
b536a002	410	// Variables
	411	Double_t xx;
	412	Double_t mpc;
	413	Double_t fland;
	414	Double_t sum = 0.0;
	415	Double_t xlow,xupp;
2ca1f4ee	416	Double_t step;
b536a002	417	Double_t i;
	418
	419	// MP shift correction
	420	mpc = par[1] - mpshift * par[0];
	421	// Range of convolution integral
	422	xlow = dedx - sc * par[2];
	423	xupp = dedx + sc * par[2];
63a0c396	424	step = (xupp-xlow) / np;
b536a002	425
	426	// Convolution integral of Landau and Gaussian by sum
	427	for(i=1.0; i<=np/2; i++) {
	428	xx = xlow + (i-.5) * step;
	429
	430	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
	431	sum += fland * TMath::Gaus(dedx,xx,par[2]);
	432
	433	xx = xupp - (i-.5) * step;
	434	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
	435	sum += fland * TMath::Gaus(dedx,xx,par[2]);
	436	}
	437
	438	return (step * sum * invsq2pi / par[2]);
	439	}
	440