#ifndef MARS_MCalibrationFits
#define MARS_MCalibrationFits

#ifndef ROOT_TMath
#include <TMath.h>
#endif

#define GIMMEABREAK     10000000000.0

inline Double_t gfKto5(Double_t *x, Double_t *par)
{

  Double_t lambda = par[0];  
  
  Double_t sum = 0.;
  Double_t arg = 0.;

  Double_t mu0 = par[1];
  Double_t mu1 = par[2];

  if (mu1 < mu0)
    return GIMMEABREAK;

  Double_t sigma0 = par[3];
  Double_t sigma1 = par[4];

  if (sigma1 < sigma0)
    return GIMMEABREAK;


  Double_t mu2 = (2.*mu1)-mu0;  
  Double_t mu3 = (3.*mu1)-(2.*mu0);
  Double_t mu4 = (4.*mu1)-(3.*mu0);
  Double_t mu5 = (5.*mu1)-(4.*mu0);

  Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));  
  Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
  Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
  Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));

  Double_t lambda2 = lambda*lambda;
  Double_t lambda3 = lambda2*lambda;
  Double_t lambda4 = lambda3*lambda;
  Double_t lambda5 = lambda4*lambda;

  // k=0:
  arg = (x[0] - mu0)/sigma0;
  sum = TMath::Exp(-0.5*arg*arg)/sigma0;
  
  // k=1:
  arg = (x[0] - mu1)/sigma1;
  sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
  
  // k=2:
  arg = (x[0] - mu2)/sigma2;
  sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;

  // k=3:
  arg = (x[0] - mu3)/sigma3;
  sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;

  // k=4:
  arg = (x[0] - mu4)/sigma4;
  sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
  
  // k=5:
  arg = (x[0] - mu5)/sigma5;
  sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;

  return par[5]*sum;

};

inline Double_t gfKto6(Double_t *x, Double_t *par)
{

  Double_t lambda = par[0];  
  
  Double_t sum = 0.;
  Double_t arg = 0.;

  Double_t mu0 = par[1];
  Double_t mu1 = par[2];

  if (mu1 < mu0)
    return GIMMEABREAK;

  Double_t sigma0 = par[3];
  Double_t sigma1 = par[4];

  if (sigma1 < sigma0)
    return GIMMEABREAK;


  Double_t mu2 = (2.*mu1)-mu0;  
  Double_t mu3 = (3.*mu1)-(2.*mu0);
  Double_t mu4 = (4.*mu1)-(3.*mu0);
  Double_t mu5 = (5.*mu1)-(4.*mu0);
  Double_t mu6 = (6.*mu1)-(5.*mu0);

  Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));  
  Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
  Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
  Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
  Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));

  Double_t lambda2 = lambda*lambda;
  Double_t lambda3 = lambda2*lambda;
  Double_t lambda4 = lambda3*lambda;
  Double_t lambda5 = lambda4*lambda;
  Double_t lambda6 = lambda5*lambda;
  
  // k=0:
  arg = (x[0] - mu0)/sigma0;
  sum = TMath::Exp(-0.5*arg*arg)/sigma0;
  
  // k=1:
  arg = (x[0] - mu1)/sigma1;
  sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
  
  // k=2:
  arg = (x[0] - mu2)/sigma2;
  sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;

  // k=3:
  arg = (x[0] - mu3)/sigma3;
  sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;

  // k=4:
  arg = (x[0] - mu4)/sigma4;
  sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
  
  // k=5:
  arg = (x[0] - mu5)/sigma5;
  sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;

  // k=6:
  arg = (x[0] - mu6)/sigma6;
  sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
  
  return par[5]*sum;

};

inline Double_t gfKto7(Double_t *x, Double_t *par)
{

  Double_t lambda = par[0];  
  
  Double_t sum = 0.;
  Double_t arg = 0.;

  Double_t mu0 = par[1];
  Double_t mu1 = par[2];

  if (mu1 < mu0)
    return GIMMEABREAK;

  Double_t sigma0 = par[3];
  Double_t sigma1 = par[4];

  if (sigma1 < sigma0)
    return GIMMEABREAK;


  Double_t mu2 = (2.*mu1)-mu0;  
  Double_t mu3 = (3.*mu1)-(2.*mu0);
  Double_t mu4 = (4.*mu1)-(3.*mu0);
  Double_t mu5 = (5.*mu1)-(4.*mu0);
  Double_t mu6 = (6.*mu1)-(5.*mu0);
  Double_t mu7 = (7.*mu1)-(6.*mu0);
  
  Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));  
  Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
  Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
  Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
  Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
  Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0));

  Double_t lambda2 = lambda*lambda;
  Double_t lambda3 = lambda2*lambda;
  Double_t lambda4 = lambda3*lambda;
  Double_t lambda5 = lambda4*lambda;
  Double_t lambda6 = lambda5*lambda;
  Double_t lambda7 = lambda6*lambda;
  
  // k=0:
  arg = (x[0] - mu0)/sigma0;
  sum = TMath::Exp(-0.5*arg*arg)/sigma0;
  
  // k=1:
  arg = (x[0] - mu1)/sigma1;
  sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
  
  // k=2:
  arg = (x[0] - mu2)/sigma2;
  sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;

  // k=3:
  arg = (x[0] - mu3)/sigma3;
  sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;

  // k=4:
  arg = (x[0] - mu4)/sigma4;
  sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
  
  // k=5:
  arg = (x[0] - mu5)/sigma5;
  sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;

  // k=6:
  arg = (x[0] - mu6)/sigma6;
  sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
  
  // k=7:
  arg = (x[0] - mu7)/sigma7;
  sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
  
  return par[5]*sum;

};


inline Double_t gfKto8(Double_t *x, Double_t *par)
{

  Double_t lambda = par[0];  
  
  Double_t sum = 0.;
  Double_t arg = 0.;

  Double_t mu0 = par[1];
  Double_t mu1 = par[2];

  if (mu1 < mu0)
    return GIMMEABREAK;

  Double_t sigma0 = par[3];
  Double_t sigma1 = par[4];

  if (sigma1 < sigma0)
    return GIMMEABREAK;


  Double_t mu2 = (2.*mu1)-mu0;  
  Double_t mu3 = (3.*mu1)-(2.*mu0);
  Double_t mu4 = (4.*mu1)-(3.*mu0);
  Double_t mu5 = (5.*mu1)-(4.*mu0);
  Double_t mu6 = (6.*mu1)-(5.*mu0);
  Double_t mu7 = (7.*mu1)-(6.*mu0);
  Double_t mu8 = (8.*mu1)-(7.*mu0);    
  
  Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));  
  Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
  Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
  Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
  Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
  Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0));
  Double_t sigma8 = TMath::Sqrt((8.*sigma1*sigma1) - (7.*sigma0*sigma0));          

  Double_t lambda2 = lambda*lambda;
  Double_t lambda3 = lambda2*lambda;
  Double_t lambda4 = lambda3*lambda;
  Double_t lambda5 = lambda4*lambda;
  Double_t lambda6 = lambda5*lambda;
  Double_t lambda7 = lambda6*lambda;
  Double_t lambda8 = lambda7*lambda;
  
  // k=0:
  arg = (x[0] - mu0)/sigma0;
  sum = TMath::Exp(-0.5*arg*arg)/sigma0;
  
  // k=1:
  arg = (x[0] - mu1)/sigma1;
  sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
  
  // k=2:
  arg = (x[0] - mu2)/sigma2;
  sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;

  // k=3:
  arg = (x[0] - mu3)/sigma3;
  sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;

  // k=4:
  arg = (x[0] - mu4)/sigma4;
  sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
  
  // k=5:
  arg = (x[0] - mu5)/sigma5;
  sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;

  // k=6:
  arg = (x[0] - mu6)/sigma6;
  sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
  
  // k=7:
  arg = (x[0] - mu7)/sigma7;
  sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
  
  // k=8:
  arg = (x[0] - mu8)/sigma8;
  sum += 0.000024801587315*lambda8*TMath::Exp(-0.5*arg*arg)/sigma7;
  
  return par[5]*sum;

};

#endif  /* MARS_MCalibrationFits */