FittingFunctions.icc

Go to the documentation of this file.

/*  FITTING FUNCTIONS */

/* Different power-law-like fit functions to describe the measured flux */


enum modeltype {StandardICRC2015, StandardICRC2013, StandardMultiSmooth,
                InfillHard
               };


double SpectrumLgSmoothFitFunction2013(double* x, double* par)
{

  /*   ICRC 2013 model:

       0:  J0 (normalization);
       1:  log10(Eankle)    (energy in log10 at ankle position);
       2:  gamma1            (spectral index before ankle);
       3:  gamma2             (spectral index after ankle);
       4:  log10(Es)             (energy in log10 at the 1/2 suppression);
       5:  log10(Wc)           (suppression smoothness )
  */

  const double lgE = x[0];

  double lgJ;

  if (lgE < par[1])
    lgJ =  par[0] - par[2] * (lgE - par[1]);
  else {
    double lgJ_exp   =  -par[3] * lgE - log10(1.0 + exp((lgE - par[4]) / par[5]));
    double lgJ_ankle =  -par[3] * par[1] - log10(1.0 + exp((par[1] - par[4]) / par[5]));
    lgJ = par[0] + lgJ_exp - lgJ_ankle;
  }

  return pow(10, lgJ);
}


double SpectrumLgSmoothFitFunction2015(double* x, double* par)
{

  /*   ICRC 2015/2017 model:

      0:  J0 (normalization);
      1:  log10(Eankle)    (energy in log10 at ankle position);
      2:  gamma1            (spectral index before ankle);
      3:  gamma2             (spectral index after ankle);
      4:  log10(Es)             (energy in log10 at the 1/2 suppression);
      5:  Delta gamma     (increment of the spectral index beyond suppression)
  */


  const double lgE = x[0];

  double lgJ;

  if (lgE < par[1])
    lgJ =  par[0] - par[2] * (lgE - par[1]);
  else {
    double E = pow(10., lgE);
    double Ea = pow(10, par[1]);
    double Es = pow(10, par[4]);

    double lgJ_exp   =  log10(1.0 + pow(E / Es, par[5]));
    double lgJ_ankle =  log10(1.0 + pow(Ea / Es, par[5]));
    lgJ = par[0] - par[3] * (lgE - par[1]) + lgJ_ankle - lgJ_exp ;
  }

  return pow(10, lgJ);
}


double
FluxModelMultiSmoothBreak(double* x, double* par)
{
  /*  ICRC 2017? model for unfolding (not for spectrum fit):
      Smooth model with additional break as proposed by Oliver
      Should be used only for the unfolding to avoid non-diff features in correction factor

      0:  J0 (normalization);
      1:  log10(Eankle)    (energy in log10 at ankle position);
      2:  log10(Ex)        (energy in log10 at additional break);
      3:  log10(Es)        (energy in log10 at the 1/2 suppression);
      4:  gamma1           (spectral index before ankle);
      5:  gamma2           (spectral index after ankle but before x-break);
      6:  gamma3           (spectral index after x-break but before suppression);
      7:  gamma4           (spectral index after suppression);

  */

  const size_t n = 8;
  const size_t breaks = n / 2 - 1;
  const double lgE = x[0];

  double lgJ = par[0] - par[breaks + 1] * (lgE - 18.5);

  for (size_t i = 1; i <= breaks; ++i) {
    lgJ += TMath::Log10(1 + pow(10, par[breaks + i] * (lgE - par[i])));
    lgJ -= TMath::Log10(1 + pow(10, par[breaks + i + 1] * (lgE - par[i])));
  }

  return pow(10, lgJ);
}


double SpectrumInfillFitFunction(double* x, double* par)
{
  // Model: only upto ankle --> J0, Eank, gamma1, gamma2

  const double lgE = x[0];

  double lgJ;

  if (lgE < par[1])
    lgJ =  par[0] - par[2] * (lgE - par[1]);
  else
    lgJ =  par[0] - par[3] * (lgE - par[1]);

  return pow(10, lgJ);
}


// TVectorD versions of the above functions

TVectorD SpectrumLgSmoothFitFunction2013(TVectorD x, double* par)
{
  int ncol = x.GetNoElements();
  TVectorD J(ncol);

  for (int i = 0; i < ncol; ++i)
    J[i] = SpectrumLgSmoothFitFunction2013(&x[i], par);

  return J;
}


TVectorD SpectrumLgSmoothFitFunction2015(TVectorD x, double* par)
{
  int ncol = x.GetNoElements();
  TVectorD J(ncol);

  for (int i = 0; i < ncol; ++i)
    J[i] = SpectrumLgSmoothFitFunction2015(&x[i], par);

  return J;
}


TVectorD FluxModelMultiSmoothBreak(TVectorD x, double* par)
{
  int ncol = x.GetNoElements();
  TVectorD J(ncol);

  for (int i = 0; i < ncol; ++i)
    J[i] = FluxModelMultiSmoothBreak(&x[i], par);

  return J;
}


TVectorD SpectrumInfillFitFunction(TVectorD x, double* par)
{
  int ncol = x.GetNoElements();
  TVectorD J(ncol);

  for (int i = 0; i < ncol; ++i)
    J[i] = SpectrumInfillFitFunction(&x[i], par);

  return J;
}


// selects one of the implemented models according to the modeltype
double FluxModel(const modeltype mtype, double* x, double* par)
{
  switch (mtype) {
    case StandardICRC2013:
      return SpectrumLgSmoothFitFunction2013(x, par);
    case StandardICRC2015:
      return SpectrumLgSmoothFitFunction2015(x, par);
    case StandardMultiSmooth:
      return FluxModelMultiSmoothBreak(x, par);
    case InfillHard:
      return SpectrumInfillFitFunction(x, par);
  }

  return SpectrumLgSmoothFitFunction2013(x, par);
}


// change also this function to changes of the above function
// unfortunately, cannot be easily changed to template
TVectorD FluxModel(const modeltype mtype, TVectorD x, double* par)
{
  switch (mtype) {
    case StandardICRC2013:
      return SpectrumLgSmoothFitFunction2013(x, par);
    case StandardICRC2015:
      return SpectrumLgSmoothFitFunction2015(x, par);
    case StandardMultiSmooth:
      return FluxModelMultiSmoothBreak(x, par);
    case InfillHard:
      return SpectrumInfillFitFunction(x, par);
  }

  return SpectrumLgSmoothFitFunction2013(x, par);
}

double
(*FluxModels(const modeltype mtype))(double*, double*)
{
  switch (mtype) {
    case StandardICRC2013:
      return &SpectrumLgSmoothFitFunction2013;
    case StandardICRC2015:
      return &SpectrumLgSmoothFitFunction2015;
    case StandardMultiSmooth:
      return &FluxModelMultiSmoothBreak;
    case InfillHard:
      return &SpectrumInfillFitFunction;
  }

  return &SpectrumLgSmoothFitFunction2013;
}


TVectorD GetCorrectionFactor(const TVectorD x, const TMatrixD m, double* par, const modeltype mtype)
{

  int ncol = x.GetNoElements();
  TVectorD f0s = FluxModel(mtype, x, par);

  for (int i = 0; i < ncol; ++i)
    f0s[i] = f0s[i] * pow(10, x[i]);

  TVectorD f1s = m * f0s;
  f1s *= (x[1] - x[0]);

  TVectorD corr(ncol);
  for (int i = 0; i < ncol; ++i)
    corr[i] = f0s[i] / f1s[i];

  return corr;

}


TVectorD GetCorrectionFactor(TVectorD parsFit, const modeltype mtype)
{

  int ncol = vecLgE.GetNoElements();

  double* par = &parsFit[0];
  TVectorD f0s = FluxModel(mtype, vecLgE, par);

  for (int i = 0; i < ncol; ++i)
    f0s[i] = f0s[i] * pow(10, vecLgE[i]);

  TVectorD f1s = kmatrix * f0s;
  f1s *= (vecLgE[1] - vecLgE[0]);

  TVectorD corr(ncol);
  for (int i = 0; i < ncol; ++i)
    corr[i] = f0s[i] / f1s[i];

  return corr;

}