[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"
#include "AliPID.h"

ClassImp(AliITSPidParams)

//______________________________________________________________________
AliITSPidParams::AliITSPidParams(Bool_t isMC):
TNamed("default",""),
  fSDDPionMPV(0),
  fSDDPionLandauWidth(0),
  fSDDPionGaussWidth(0),
  fSSDPionMPV(0),
  fSSDPionLandauWidth(0),
  fSSDPionGaussWidth(0),
  fSDDKaonMPV(0),
  fSDDKaonLandauWidth(0),
  fSDDKaonGaussWidth(0),
  fSSDKaonMPV(0),
  fSSDKaonLandauWidth(0),
  fSSDKaonGaussWidth(0),
  fSDDProtMPV(0),
  fSDDProtLandauWidth(0),
  fSDDProtGaussWidth(0),
  fSSDProtMPV(0),
  fSSDProtLandauWidth(0),
  fSSDProtGaussWidth(0)
{
  // default constructor
  if (isMC) InitMC();
  else InitData();
}
//______________________________________________________________________
AliITSPidParams::AliITSPidParams(Char_t * name, Bool_t isMC):
  TNamed(name,""),
  fSDDPionMPV(0),
  fSDDPionLandauWidth(0),
  fSDDPionGaussWidth(0),
  fSSDPionMPV(0),
  fSSDPionLandauWidth(0),
  fSSDPionGaussWidth(0),
  fSDDKaonMPV(0),
  fSDDKaonLandauWidth(0),
  fSDDKaonGaussWidth(0),
  fSSDKaonMPV(0),
  fSSDKaonLandauWidth(0),
  fSSDKaonGaussWidth(0),
  fSDDProtMPV(0),
  fSDDProtLandauWidth(0),
  fSDDProtGaussWidth(0),
  fSSDProtMPV(0),
  fSSDProtLandauWidth(0), 
  fSSDProtGaussWidth(0)
{
  // standard constructor
  if (isMC) InitMC();
  else InitData();
}
//______________________________________________________________________
AliITSPidParams::~AliITSPidParams(){
  // 
  if(fSDDPionMPV) delete fSDDPionMPV;
  if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
  if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;

  if(fSSDPionMPV) delete fSSDPionMPV;
  if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
  if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
  
  if(fSDDKaonMPV) delete fSDDKaonMPV;
  if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
  if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
  
  if(fSSDKaonMPV) delete fSSDKaonMPV;
  if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
  if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;

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

  if(fSSDProtMPV) delete fSSDProtMPV;
  if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
  if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
}

//______________________________________________________________________
void AliITSPidParams::InitMC(){
  // initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
  // parameter values from LHC10d1 

  // pions
  if(fSDDPionMPV) delete fSDDPionMPV;
  fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
  fSDDPionMPV->SetParameters(-0.892291, 0.003630, 1.866484, 78.378179);

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

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

  if(fSSDPionMPV) delete fSSDPionMPV;
  fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
  fSSDPionMPV->SetParameters(-0.896507, 0.003173, 2.017155, 80.682567);

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

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

  // kaons
  if(fSDDKaonMPV) delete fSDDKaonMPV;
  fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSDDKaonMPV->SetParameters(17.581590, -0.120134, 72.550701);

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

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

  if(fSSDKaonMPV) delete fSSDKaonMPV;
  fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSSDKaonMPV->SetParameters(16.238778, 0.039318, 75.863719);

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

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

  // protons
  if(fSDDProtMPV) delete fSDDProtMPV;
  fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSDDProtMPV->SetParameters(64.482762, -1.667823, 71.850731);

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

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

  if(fSSDProtMPV) delete fSSDProtMPV;
  fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSSDProtMPV->SetParameters(63.817375, -1.221779, 73.233644);

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

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

}
//______________________________________________________________________
void AliITSPidParams::InitData(){
  // initialize TFormulas to Real Data values (=p-p simulations PYTHIA+GEANT)
  // parameter values from LHC10b 
  
  // pions
  if(fSDDPionMPV) delete fSDDPionMPV;
  fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
  fSDDPionMPV->SetParameters(-0.907609, 0.006521, 3.340347, 81.297942);

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

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

  if(fSSDPionMPV) delete fSSDPionMPV;
  fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(x*x)*TMath::Log(x)+[1]/(x*x*x*x)*TMath::Log(x)+[2]*TMath::Log(x)+[3]");
  fSSDPionMPV->SetParameters(-0.920046, 0.006061, 3.428578, 81.401816);

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

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

  // kaons
  if(fSDDKaonMPV) delete fSDDKaonMPV;
  fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSDDKaonMPV->SetParameters(15.429146, 0.178251, 74.640293);

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

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

  if(fSSDKaonMPV) delete fSSDKaonMPV;
  fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSSDKaonMPV->SetParameters(15.170715, 0.181379, 74.951884);

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

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

  // protons
  if(fSDDProtMPV) delete fSDDProtMPV;
  fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSDDProtMPV->SetParameters(61.452534, 0.372908, 71.668352);

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

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

  if(fSSDProtMPV) delete fSSDProtMPV;
  fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(x*x)+[1]/(x*x*x*x)+[2]");
  fSSDProtMPV->SetParameters(60.246538, 0.000323, 71.992031);

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

  if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
  fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
  fSSDProtGaussWidth->SetParameters(5.050541, 6.290964);

}
//_______________________________________________________________________
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)==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[3];
  Bool_t isSet=kFALSE;
  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
  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 = 0.;
  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];
  if(np!=0) 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"
	30	#include "AliPID.h"
	31
	32	ClassImp(AliITSPidParams)
	33
	34	//______________________________________________________________________
b52bfc67	35	AliITSPidParams::AliITSPidParams(Bool_t isMC):
b52bfc67	36	TNamed("default",""),
b536a002	37	fSDDPionMPV(0),
	38	fSDDPionLandauWidth(0),
	39	fSDDPionGaussWidth(0),
	40	fSSDPionMPV(0),
	41	fSSDPionLandauWidth(0),
	42	fSSDPionGaussWidth(0),
	43	fSDDKaonMPV(0),
	44	fSDDKaonLandauWidth(0),
	45	fSDDKaonGaussWidth(0),
	46	fSSDKaonMPV(0),
	47	fSSDKaonLandauWidth(0),
	48	fSSDKaonGaussWidth(0),
	49	fSDDProtMPV(0),
	50	fSDDProtLandauWidth(0),
	51	fSDDProtGaussWidth(0),
	52	fSSDProtMPV(0),
	53	fSSDProtLandauWidth(0),
	54	fSSDProtGaussWidth(0)
	55	{
	56	// default constructor
b52bfc67	57	if (isMC) InitMC();
b52bfc67	58	else InitData();
b536a002	59	}
b536a002	60	//______________________________________________________________________
b52bfc67	61	AliITSPidParams::AliITSPidParams(Char_t * name, Bool_t isMC):
b536a002	62	TNamed(name,""),
	63	fSDDPionMPV(0),
	64	fSDDPionLandauWidth(0),
	65	fSDDPionGaussWidth(0),
	66	fSSDPionMPV(0),
	67	fSSDPionLandauWidth(0),
	68	fSSDPionGaussWidth(0),
	69	fSDDKaonMPV(0),
	70	fSDDKaonLandauWidth(0),
	71	fSDDKaonGaussWidth(0),
	72	fSSDKaonMPV(0),
	73	fSSDKaonLandauWidth(0),
	74	fSSDKaonGaussWidth(0),
	75	fSDDProtMPV(0),
	76	fSDDProtLandauWidth(0),
	77	fSDDProtGaussWidth(0),
	78	fSSDProtMPV(0),
b52bfc67	79	fSSDProtLandauWidth(0),
b536a002	80	fSSDProtGaussWidth(0)
	81	{
	82	// standard constructor
b52bfc67	83	if (isMC) InitMC();
b52bfc67	84	else InitData();
b536a002	85	}
	86	//______________________________________________________________________
	87	AliITSPidParams::~AliITSPidParams(){
	88	//
	89	if(fSDDPionMPV) delete fSDDPionMPV;
	90	if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
	91	if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
	92
	93	if(fSSDPionMPV) delete fSSDPionMPV;
	94	if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
	95	if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
b52bfc67	96
b536a002	97	if(fSDDKaonMPV) delete fSDDKaonMPV;
	98	if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
	99	if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
b52bfc67	100
b536a002	101	if(fSSDKaonMPV) delete fSSDKaonMPV;
	102	if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
	103	if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
	104
	105	if(fSDDProtMPV) delete fSDDProtMPV;
	106	if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
	107	if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
	108
	109	if(fSSDProtMPV) delete fSSDProtMPV;
	110	if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
	111	if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
	112	}
	113
	114	//______________________________________________________________________
8abeb05b	115	void AliITSPidParams::InitMC(){
	116	// initialize TFormulas to Monte Carlo values (=p-p simulations PYTHIA+GEANT)
	117	// parameter values from LHC10d1
b536a002	118
	119	// pions
	120	if(fSDDPionMPV) delete fSDDPionMPV;
8abeb05b	121	fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(xx)TMath::Log(x)+[1]/(xxxx)TMath::Log(x)+[2]*TMath::Log(x)+[3]");
b52bfc67	122	fSDDPionMPV->SetParameters(-0.892291, 0.003630, 1.866484, 78.378179);
b536a002	123
b536a002	124	if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
8abeb05b	125	fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
b52bfc67	126	fSDDPionLandauWidth->SetParameters(0.080999, 5.917715);
b536a002	127
b536a002	128	if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
8abeb05b	129	fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
b52bfc67	130	fSDDPionGaussWidth->SetParameters(-0.092822, 7.839729);
b536a002	131
b536a002	132	if(fSSDPionMPV) delete fSSDPionMPV;
8abeb05b	133	fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(xx)TMath::Log(x)+[1]/(xxxx)TMath::Log(x)+[2]*TMath::Log(x)+[3]");
b52bfc67	134	fSSDPionMPV->SetParameters(-0.896507, 0.003173, 2.017155, 80.682567);
b536a002	135
b536a002	136	if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
8abeb05b	137	fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
b52bfc67	138	fSSDPionLandauWidth->SetParameters(0.087182, 5.843610);
b536a002	139
b536a002	140	if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
8abeb05b	141	fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
b52bfc67	142	fSSDPionGaussWidth->SetParameters(-0.110444, 5.837737);
b536a002	143
	144	// kaons
	145	if(fSDDKaonMPV) delete fSDDKaonMPV;
b52bfc67	146	fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(xx)+[1]/(xxxx)+[2]");
b52bfc67	147	fSDDKaonMPV->SetParameters(17.581590, -0.120134, 72.550701);
b536a002	148
b536a002	149	if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
8abeb05b	150	fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
b52bfc67	151	fSDDKaonLandauWidth->SetParameters(1.271756, 5.778888);
	152
	153	if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
	154	fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
	155	fSDDKaonGaussWidth->SetParameters(-1.650298, 8.322084);
	156
	157	if(fSSDKaonMPV) delete fSSDKaonMPV;
	158	fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(xx)+[1]/(xxxx)+[2]");
	159	fSSDKaonMPV->SetParameters(16.238778, 0.039318, 75.863719);
	160
	161	if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
	162	fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
	163	fSSDKaonLandauWidth->SetParameters(1.179541, 5.961353);
	164
	165	if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
	166	fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
	167	fSSDKaonGaussWidth->SetParameters(-2.019126, 6.155977);
	168
	169	// protons
	170	if(fSDDProtMPV) delete fSDDProtMPV;
	171	fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(xx)+[1]/(xxxx)+[2]");
	172	fSDDProtMPV->SetParameters(64.482762, -1.667823, 71.850731);
	173
	174	if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
	175	fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
	176	fSDDProtLandauWidth->SetParameters(6.948997, 3.928018);
	177
	178	if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
	179	fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
	180	fSDDProtGaussWidth->SetParameters(-6.522760, 12.673959);
	181
	182	if(fSSDProtMPV) delete fSSDProtMPV;
	183	fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(xx)+[1]/(xxxx)+[2]");
	184	fSSDProtMPV->SetParameters(63.817375, -1.221779, 73.233644);
	185
	186	if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
	187	fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
	188	fSSDProtLandauWidth->SetParameters(7.286942, 3.581451);
	189
	190	if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
	191	fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
	192	fSSDProtGaussWidth->SetParameters(-8.327867, 9.723422);
	193
	194	}
	195	//______________________________________________________________________
	196	void AliITSPidParams::InitData(){
	197	// initialize TFormulas to Real Data values (=p-p simulations PYTHIA+GEANT)
	198	// parameter values from LHC10b
	199
	200	// pions
	201	if(fSDDPionMPV) delete fSDDPionMPV;
	202	fSDDPionMPV=new TFormula("fSDDPionMPV","[0]/(xx)TMath::Log(x)+[1]/(xxxx)TMath::Log(x)+[2]*TMath::Log(x)+[3]");
	203	fSDDPionMPV->SetParameters(-0.907609, 0.006521, 3.340347, 81.297942);
	204
	205	if(fSDDPionLandauWidth) delete fSDDPionLandauWidth;
	206	fSDDPionLandauWidth=new TFormula("fSDDPionLandauWidth","[0]/(x*x)+[1]");
	207	fSDDPionLandauWidth->SetParameters(0.077272, 5.478557);
	208
	209	if(fSDDPionGaussWidth) delete fSDDPionGaussWidth;
	210	fSDDPionGaussWidth=new TFormula("fSDDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
	211	fSDDPionGaussWidth->SetParameters(-0.098529, 10.265711);
	212
	213	if(fSSDPionMPV) delete fSSDPionMPV;
	214	fSSDPionMPV=new TFormula("fSSDPionMPV","[0]/(xx)TMath::Log(x)+[1]/(xxxx)TMath::Log(x)+[2]*TMath::Log(x)+[3]");
215	fSSDPionMPV->SetParameters(-0.920046, 0.006061, 3.428578, 81.401816);
216
217	if(fSSDPionLandauWidth) delete fSSDPionLandauWidth;
218	fSSDPionLandauWidth=new TFormula("fSSDPionLandauWidth","[0]/(x*x)+[1]");
219	fSSDPionLandauWidth->SetParameters(0.071243, 5.388830);
220
221	if(fSSDPionGaussWidth) delete fSSDPionGaussWidth;
222	fSSDPionGaussWidth=new TFormula("fSSDPionGaussWidth","[0]/(xx)TMath::Log(x)+[1]");
223	fSSDPionGaussWidth->SetParameters(-0.099189, 7.412309);
224
225	// kaons
226	if(fSDDKaonMPV) delete fSDDKaonMPV;
227	fSDDKaonMPV=new TFormula("fSDDKaonMPV","[0]/(xx)+[1]/(xxxx)+[2]");
228	fSDDKaonMPV->SetParameters(15.429146, 0.178251, 74.640293);
229
230	if(fSDDKaonLandauWidth) delete fSDDKaonLandauWidth;
231	fSDDKaonLandauWidth=new TFormula("fSDDKaonLandauWidth","[0]/(x*x)+[1]");
232	fSDDKaonLandauWidth->SetParameters(0.975202, 5.699311);
b536a002	233
b536a002	234	if(fSDDKaonGaussWidth) delete fSDDKaonGaussWidth;
8abeb05b	235	fSDDKaonGaussWidth=new TFormula("fSDDKaonGaussWidth","[0]/(x*x)+[1]");
b52bfc67	236	fSDDKaonGaussWidth->SetParameters(1.660840, 9.389192);
b536a002	237
b536a002	238	if(fSSDKaonMPV) delete fSSDKaonMPV;
b52bfc67	239	fSSDKaonMPV=new TFormula("fSSDKaonMPV","[0]/(xx)+[1]/(xxxx)+[2]");
b52bfc67	240	fSSDKaonMPV->SetParameters(15.170715, 0.181379, 74.951884);
b536a002	241
b536a002	242	if(fSSDKaonLandauWidth) delete fSSDKaonLandauWidth;
8abeb05b	243	fSSDKaonLandauWidth=new TFormula("fSSDKaonLandauWidth","[0]/(x*x)+[1]");
b52bfc67	244	fSSDKaonLandauWidth->SetParameters(0.756466, 5.818274);
b536a002	245
b536a002	246	if(fSSDKaonGaussWidth) delete fSSDKaonGaussWidth;
8abeb05b	247	fSSDKaonGaussWidth=new TFormula("fSSDKaonGaussWidth","[0]/(x*x)+[1]");
b52bfc67	248	fSSDKaonGaussWidth->SetParameters(1.546693, 6.389872);
b536a002	249
	250	// protons
	251	if(fSDDProtMPV) delete fSDDProtMPV;
b52bfc67	252	fSDDProtMPV=new TFormula("fSDDProtMPV","[0]/(xx)+[1]/(xxxx)+[2]");
b52bfc67	253	fSDDProtMPV->SetParameters(61.452534, 0.372908, 71.668352);
b536a002	254
b536a002	255	if(fSDDProtLandauWidth) delete fSDDProtLandauWidth;
8abeb05b	256	fSDDProtLandauWidth=new TFormula("fSDDProtLandauWidth","[0]/(x*x)+[1]");
b52bfc67	257	fSDDProtLandauWidth->SetParameters(3.667023, 5.430721);
b536a002	258
b536a002	259	if(fSDDProtGaussWidth) delete fSDDProtGaussWidth;
8abeb05b	260	fSDDProtGaussWidth=new TFormula("fSDDProtGaussWidth","[0]/(x*x)+[1]");
b52bfc67	261	fSDDProtGaussWidth->SetParameters(5.503814, 9.657439);
b536a002	262
b536a002	263	if(fSSDProtMPV) delete fSSDProtMPV;
b52bfc67	264	fSSDProtMPV=new TFormula("fSSDProtMPV","[0]/(xx)+[1]/(xxxx)+[2]");
b52bfc67	265	fSSDProtMPV->SetParameters(60.246538, 0.000323, 71.992031);
8abeb05b	266
b536a002	267	if(fSSDProtLandauWidth) delete fSSDProtLandauWidth;
8abeb05b	268	fSSDProtLandauWidth=new TFormula("fSSDProtLandauWidth","[0]/(x*x)+[1]");
b52bfc67	269	fSSDProtLandauWidth->SetParameters(2.568323, 5.939774);
b536a002	270
b536a002	271	if(fSSDProtGaussWidth) delete fSSDProtGaussWidth;
8abeb05b	272	fSSDProtGaussWidth=new TFormula("fSSDProtGaussWidth","[0]/(x*x)+[1]");
b52bfc67	273	fSSDProtGaussWidth->SetParameters(5.050541, 6.290964);
b536a002	274
	275	}
	276	//_______________________________________________________________________
	277	Double_t AliITSPidParams::GetLandauGausNormPdgCode(Double_t dedx, Int_t pdgCode, Double_t mom, Int_t lay) const {
	278	// Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
	279	if(TMath::Abs(pdgCode)==211) return GetLandauGausNorm(dedx,AliPID::kPion,mom,lay);
	280	else if(TMath::Abs(pdgCode)==321) return GetLandauGausNorm(dedx,AliPID::kKaon,mom,lay);
	281	else if(TMath::Abs(pdgCode)==2212) return GetLandauGausNorm(dedx,AliPID::kProton,mom,lay);
	282	else return 0.;
	283	}
	284	//_______________________________________________________________________
b52bfc67	285	Double_t AliITSPidParams::GetLandauGausNorm(Double_t dedx, Int_t partType, Double_t mom, Int_t lay) const{
b536a002	286	// Computes Landau Gauss convolution for given particle specie and given momentum in a given ITS layer
b52bfc67	287
b536a002	288	Double_t par[3];
	289	Bool_t isSet=kFALSE;
	290	if(partType==AliPID::kPion){
	291	if(lay==3 \|\| lay==4){
	292	par[0]=GetSDDPionLandauWidth(mom);
	293	par[1]=GetSDDPionMPV(mom);
	294	par[2]=GetSDDPionGaussWidth(mom);
	295	isSet=kTRUE;
	296	}
	297	else if(lay==5 \|\| lay==6){
	298	par[0]=GetSSDPionLandauWidth(mom);
	299	par[1]=GetSSDPionMPV(mom);
	300	par[2]=GetSSDPionGaussWidth(mom);
	301	isSet=kTRUE;
	302	}
	303	}else if(partType==AliPID::kKaon){
	304	if(lay==3 \|\| lay==4){
	305	par[0]=GetSDDKaonLandauWidth(mom);
	306	par[1]=GetSDDKaonMPV(mom);
	307	par[2]=GetSDDKaonGaussWidth(mom);
	308	isSet=kTRUE;
	309	}
	310	else if(lay==5 \|\| lay==6){
	311	par[0]=GetSSDKaonLandauWidth(mom);
	312	par[1]=GetSSDKaonMPV(mom);
	313	par[2]=GetSSDKaonGaussWidth(mom);
	314	isSet=kTRUE;
	315	}
	316	}else if(partType==AliPID::kProton){
	317	if(lay==3 \|\| lay==4){
	318	par[0]=GetSDDProtLandauWidth(mom);
	319	par[1]=GetSDDProtMPV(mom);
	320	par[2]=GetSDDProtGaussWidth(mom);
	321	isSet=kTRUE;
	322	}
	323	else if(lay==5 \|\| lay==6){
	324	par[0]=GetSSDProtLandauWidth(mom);
	325	par[1]=GetSSDProtMPV(mom);
	326	par[2]=GetSSDProtGaussWidth(mom);
	327	isSet=kTRUE;
	328	}
	329	}
	330	if(!isSet) return 0.;
	331	// Numeric constants
	332	Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
	333	Double_t mpshift = -0.22278298; // Landau maximum location
	334	// Control constants
	335	Double_t np = 100.0; // number of convolution steps
	336	Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
	337	// Variables
	338	Double_t xx;
	339	Double_t mpc;
	340	Double_t fland;
	341	Double_t sum = 0.0;
	342	Double_t xlow,xupp;
14f1402f	343	Double_t step = 0.;
b536a002	344	Double_t i;
	345
	346	// MP shift correction
	347	mpc = par[1] - mpshift * par[0];
	348	// Range of convolution integral
	349	xlow = dedx - sc * par[2];
	350	xupp = dedx + sc * par[2];
	351	if(np!=0) step = (xupp-xlow) / np;
	352
	353	// Convolution integral of Landau and Gaussian by sum
	354	for(i=1.0; i<=np/2; i++) {
	355	xx = xlow + (i-.5) * step;
	356
	357	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
	358	sum += fland * TMath::Gaus(dedx,xx,par[2]);
	359
	360	xx = xupp - (i-.5) * step;
	361	fland = TMath::Landau(xx,mpc,par[0]) / par[0];
	362	sum += fland * TMath::Gaus(dedx,xx,par[2]);
	363	}
	364
	365	return (step * sum * invsq2pi / par[2]);
	366	}
	367