#ifndef PndMvdPidLikelihoodfCC
#define PndMvdPidLikelihoodfCC

#include "MvdPid/PndMvdPidLikelihood.h"
#include "MvdData/PndMvdPidCand.h"
//public

void PndMvdPidLikelihood::CalcLikelihood(Double_t momentum, Double_t energyloss, std::map<Int_t, Double_t>& likelihood) {
  //computes pid likelihoods for pi, k, p, mu and stores them in the likelihood map 
  //with the map key used as pdg id. This is basically done by integrating the
  //respective energy loss distribution within the chosen selector limits
  
  Double_t weightP;
  Double_t weightK;
  Double_t weightPi;

  //Proton selector
  if (energyloss>=LowerProtonBoundary(momentum)) {
    weightP=1-pRemainder;
    weightK=1-LandauIntegral((LowerProtonBoundary(momentum)-LowerBoundary(momentum, kMass)-kShift)/kScale);
    weightPi=1-LandauIntegral((LowerProtonBoundary(momentum)-LowerBoundary(momentum, piMass)-piShift)/piScale);

  //Kaon selector
  } else if (energyloss>=LowerKaonBoundary(momentum)) {
    weightP=pRemainder;
    weightK=LandauIntegral((LowerProtonBoundary(momentum)-LowerBoundary(momentum, kMass)-kShift)/kScale)
           -kRemainder;
    weightPi=LandauIntegral((LowerProtonBoundary(momentum)-LowerBoundary(momentum, piMass)-piShift)/piScale)
            -LandauIntegral((LowerKaonBoundary(momentum)-LowerBoundary(momentum, piMass)-piShift)/piScale);
  
  //Pion selector
  } else if (energyloss>=LowerPionBoundary(momentum)) {
    weightP=0;
    weightK=kRemainder;
    weightPi=LandauIntegral((LowerKaonBoundary(momentum)-LowerBoundary(momentum, piMass)-piShift)/piScale)
            -piRemainder;

  //Electron selector
  } else {
    weightP=0;
    weightK=0;
    weightPi=piRemainder;
  }
  //Muons have approximately the same distribution as pions.
  Double_t sum=weightP+weightK+2*weightPi;
  weightP/=sum;
  weightK/=sum;
  weightPi/=sum;

  likelihood[211]=weightPi;
  likelihood[2212]=weightP;
  likelihood[321]=weightK;
  likelihood[13]=weightPi;
}

void PndMvdPidLikelihood::CalcLikelihood(PndMvdPidCand* cand) {
  cand->fmeanMomentum=0;
  Double_t dE=0;
  Double_t dx=0;
  for (Int_t k=0;k<cand->fdE.size(); k++) {
    dE+=cand->fdE[k];
    dx+=cand->fdx[k];
    cand->fmeanMomentum+=cand->fmomentum[k];
  }
  cand->fmeanEnergyloss=dE/dx;
  cand->fmeanMomentum/=cand->fdE.size();

  PndMvdPidLikelihood::CalcLikelihood(cand->fmeanMomentum, cand->fmeanEnergyloss, cand->flikelihood);
};


//private

Double_t PndMvdPidLikelihood::LowerProtonBoundary(Double_t p) {
  return LowerBoundary(p+0.02, pMass)-5e-4;
};

Double_t PndMvdPidLikelihood::LowerKaonBoundary(Double_t p) {
  return LowerBoundary(p+0.02, kMass)-3e-4;
};

Double_t PndMvdPidLikelihood::LowerPionBoundary(Double_t p) {
  return LowerBoundary(p+0.01, piMass)-3e-4;
};


Double_t PndMvdPidLikelihood::LowerBoundary(Double_t p, Double_t m) {
//Calculate the lower boundary of the energy loss distribution.
  
  Double_t sqrfBeta=1/(1+pow(m/p,2));
  return 4.9312e-05 * (log(2*eMass*c*c/eb*sqrfBeta/(1-sqrfBeta))-sqrfBeta)/sqrfBeta;
};

Double_t PndMvdPidLikelihood::LandauIntegral(Double_t x) {
//compute integral approximation of a Landau-p.d.f from -INF to x
//with expected value 0 and deviation 1. Hard coded numbers have
//been fitted within their respective intervals.

  if (x<-0.8)
  //gauss
    return 0.199785*exp(-pow(x+0.149198,2)*0.769779);
  else if (x<0.5)
    //linear
    return 0.177214*(x+1.61437);
  else
    //power function
    return 1-19.0054*pow(x+5.860003,-1.84611);
}

Double_t PndMvdPidLikelihood::EvalPolynom(Double_t x, Double_t* c) {
//Evaluate Polynom at x -- HONER STYLE!!!

  Double_t v=c[6];
   for (Int_t k=6;k>0;) 
    v=v*x+c[--k];
  return v;
}

//[GeV/c**2]
float PndMvdPidLikelihood::piMass=0.1396;
float PndMvdPidLikelihood::kMass=0.4937;
float PndMvdPidLikelihood::pMass=0.9383;
float PndMvdPidLikelihood::eMass=0.511e-3;
//[GeV]
float PndMvdPidLikelihood::eb=0.14e-6;
//[m/s]
float PndMvdPidLikelihood::c=2.99792458e8;

//Landau distribution parameters; these parameters have been 
//obtained by fitting the Landau-p.d.f on the energy loss shifted
//by the Bethe-Bloch energy loss at momentum >= 0.4 GeV
Double_t PndMvdPidLikelihood::kShift=3.38598e-4;
Double_t PndMvdPidLikelihood::kScale=1.36362e-4;

Double_t PndMvdPidLikelihood::piShift=3.09159e-4;
Double_t PndMvdPidLikelihood::piScale=1.58696e-4;

//General remainders of integrals for integrating from -INF to 
//below the Bethe-Bloch lower boundary as obtained by MC Simulation
Double_t PndMvdPidLikelihood::pRemainder=0.01;
Double_t PndMvdPidLikelihood::kRemainder=0.004;
Double_t PndMvdPidLikelihood::piRemainder=0.001;

#endif