ICRC2019/FittingFunctions.hh

Go to the documentation of this file.

/*  FITTING FUNCTIONS */

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


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_t Spectrum2015_2017SmoothAnkle(Double_t *v, Double_t *par) {
  Double_t J;
  //par[0] J0
  //par[1] E12
  //par[2] gamma1
  //par[3] gamma2
  //par[4] E23
  //par[5] Delta gamma
  double J0 = pow(10, par[0]);
  double E12 = pow(10, par[1]);
  double E23 = pow(10, par[4]);

  double E = pow(10, v[0]);

  J = pow(E / E12, -par[2]) * ( 1 + pow(E / E12, par[2]) ) / (1 + pow(E / E12, par[3] ) ) * (1. + pow(E12 / E23, par[5]) ) / (1 + pow(E / E23, par[5]));

  J *= J0;
  return J;
}


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);
}

double SpectrumInfillICRC2019(double* x, double* par) {
  /*
  0: J0 (normalization at break)
  1: log10(E0)
  2: log10(Eankle)
  3: gamma0
  4: gamma1
  5: gamma2
  6: dGamma0
  */

  const double lgNorm = par[0];
  const double lgEbreak = par[1];
  const double lgEa = par[2];
  const double gamma0 = par[3];
  const double gamma1 = par[4];
  const double gamma2 = par[5];
  const double rate = par[6];

  const double lgE = x[0];

  //Maybe one day we can actually fit curve
  //For now, set to infinity
  const double lgEs = 1000000;//log10(39.e18);
  const double dGamma = 2.5;

  double lgJ;

  if (lgE < lgEa) {
    const double factor = pow(10, lgE) / pow(10, lgEbreak);
    lgJ =  lgNorm - gamma0 * (lgE - lgEbreak) + (gamma0 - gamma1) / rate * log10(1 + pow(factor, rate));
  } else {
    const double factor = pow(10, lgEa) / pow(10, lgEbreak);
    double E = pow(10., lgE);
    double Ea = pow(10, lgEa);
    double Es = pow(10, lgEs);
    double lgJ_exp   =  log10(1.0 + pow(E / Es, dGamma));
    double lgJ_ankle =  log10(1.0 + pow(Ea / Es, dGamma));
    lgJ =  lgNorm - gamma0 * (lgEa - lgEbreak) + (gamma0 - gamma1) / rate * log10(1 + pow(factor, rate)) - gamma2 * (lgE - lgEa) + lgJ_ankle - lgJ_exp;
  }

  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  Spectrum2015_2017SmoothAnkle(TVectorD x, double* par) {
  int ncol = x.GetNoElements();
  TVectorD J(ncol);

  for (int i = 0; i < ncol; ++i)
    J[i] = Spectrum2015_2017SmoothAnkle(&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;
}


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

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

  return J;
}



// selects one of the implemented models according to the modeltype
TVectorD FluxModel(string model, TVectorD x, double* par) {
  if (model == "StandardICRC2013")
    return SpectrumLgSmoothFitFunction2013(x, par);
  else if (model == "StandardICRC2015")
    return SpectrumLgSmoothFitFunction2015(x, par);
  else if (model == "StandardMultiSmooth")
    return FluxModelMultiSmoothBreak(x, par);
  else if (model == "InfillHard")
    return SpectrumInfillFitFunction(x, par);
  else if (model == "Spectrum2015_2017SmoothAnkle")
    return Spectrum2015_2017SmoothAnkle(x, par);
  else if (model == "SpectrumInfillICRC2019")
    return SpectrumInfillICRC2019(x, par);
  else {
    cerr << endl;
    cerr << model << ": no valid unfolding function chosen " << endl;
    throw std::exception();
  }

  return SpectrumLgSmoothFitFunction2013(x, par);
}

double
(*FluxModels(string model))(double*, double*) {

  if (model == "StandardICRC2013")
    return &SpectrumLgSmoothFitFunction2013;
  else if (model == "StandardICRC2015")
    return &SpectrumLgSmoothFitFunction2015;
  else if (model == "StandardMultiSmooth")
    return &FluxModelMultiSmoothBreak;
  else if (model == "InfillHard")
    return &SpectrumInfillFitFunction;
  else if (model == "Spectrum2015_2017SmoothAnkle")
    return &Spectrum2015_2017SmoothAnkle;
  else if (model == "SpectrumInfillICRC2019")
    return &SpectrumInfillICRC2019;
  else {
    cerr << endl;
    cerr << model << ": no valid unfolding function chosen " << endl;
    throw std::exception();
  }

  return &SpectrumLgSmoothFitFunction2013;

}


TVectorD GetCorrectionFactor(const TVectorD x, const TMatrixD m, double* par, string model) {
  int ncol = x.GetNoElements();
  TVectorD f0s = FluxModel(model, 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, string model) {
  int ncol = vecLgE.GetNoElements();

  double* par = &parsFit[0];
  TVectorD f0s = FluxModel(model, 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;
}