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

Double_t SpectrumLipari(double* x, double* par)
{
 /*
  0: J0
  1: E12
  2: gamma1
  3: E23
  4: gamma2
  5: E34
  6: gamma3
  7: gamma4
  8: w12
  9: w23
  10: w34
 */
  double lgE = x[0];
  double w[3] = {par[8],par[9],par[10]};
  double J = par[0] - par[4] * (lgE - 18.5);
  double factor =  (par[4] - par[5]) * w[0];//
  factor *= log10( 1. + pow(10,(lgE - par[1]) / w[0] ) );
  J += factor;
  factor =  (par[5] - par[6]) * w[1];//
  factor *= log10( 1. + pow(10,(lgE - par[2]) / w[1] ) );
  J += factor;
  factor =  (par[6] - par[7]) * w[2];//
  factor *= log10( 1. + pow(10,(lgE - par[3]) / w[2] ) );

  J += factor;
  
  return pow(10,J);
}



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

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

  for (int i = 0; i < ncol; ++i)
    J[i] = SpectrumLipari(&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 if (model == "SpectrumLipari")
    return SpectrumLipari(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 if (model == "SpectrumLipari")
    return &SpectrumLipari;
  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;
}


TVectorD GetCorrectionFactor_Likelihood(const TVectorD x, const TMatrixD m, double* par,string model, TAxis *frawAxis_lgE) {

  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];
*/
  int N =  frawAxis_lgE->GetNbins();
  TVectorD mws(N);
  mws.Zero();
  TVectorD mws0(N);
  mws0.Zero();

  TVectorD corr_new(N);
  for (int i = 0; i < ncol; ++i) {
    double ilgE = vecLgE[i];
    int bin = frawAxis_lgE->FindBin(ilgE);

    if (bin < 1)
      continue;

    if (bin > N)
      break;

    //cout << vecLgE[i] << " " << bin << endl;
    mws[bin - 1] += f1s[i];   // expected number of events (delta in poisson distrib.)
    mws0[bin - 1] += f0s[i];   // expected 'true' number of events (delta in poisson distrib.)
  }

  for (int j = 0; j < N; ++j){
    corr_new[j] = mws0[j]/mws[j];
  }
 
  return corr_new;

}

TVectorD GetCorrectionFactorError_Likelihood(const TVectorD x, const TMatrixD m, int npar, TVectorD par, TVectorD parerr, TMatrixD parcov, string model, TAxis *fAxis_lgE) {
  //
  int Nbins = fAxis_lgE->GetNbins();
  TVectorD UnfoldCorrectionFactor_Likelihood(Nbins);
  UnfoldCorrectionFactor_Likelihood = GetCorrectionFactor_Likelihood(x, m, &par[0], model.c_str(), fAxis_lgE);

  //
  TVectorD UnfoldCorrectionFactorError_Likelihood(Nbins);
  //  
  for (int j = 0; j < Nbins; ++j){
    //
    double dunf_dp[npar];
    for (int k=0; k<npar; k++){
      //
      dunf_dp[k] = 0.;
      //
      TVectorD fpar;
      fpar.ResizeTo(npar);
      for (int l=0; l<npar; l++){
        fpar[l] = par[l];
      }
      double ff = 1.;
      fpar[k] = par[k] + parerr[k]*ff;
      //fpar[k] = par[k] *(1.+1e-2);
      //cout << k << " " << parerr[k]*ff << endl;
      TVectorD UnfoldCorrectionFactor_Likelihood_1(Nbins);
      UnfoldCorrectionFactor_Likelihood_1 = GetCorrectionFactor_Likelihood(x, m, &fpar[0], model.c_str(), fAxis_lgE);
      //
      if (parerr[k]>0)
        dunf_dp[k] = (UnfoldCorrectionFactor_Likelihood[j] - UnfoldCorrectionFactor_Likelihood_1[j]) / (fpar[k]-par[k]);
      //
    }
    //
    double unf_err = 0;
    for (int k1=0; k1<npar; k1++){
      for (int k2=0; k2<npar; k2++){
        unf_err = unf_err + dunf_dp[k1]*dunf_dp[k2] * parcov[k1][k2];
      }
    }
    if (unf_err>0.){
      unf_err = sqrt(unf_err);
    } else {
      unf_err = 0.;
    }
    //
    // VV
    // avoid to0 small errors in the unfolding coefficients
    // (it does not affect the results)
    double err_th = 0.15/1e2;
    if (unf_err<err_th) unf_err = err_th;
    UnfoldCorrectionFactorError_Likelihood[j] = unf_err;
    //
    //cout << j << " " << fAxis_lgE->GetBinCenter(j+1) << " " << UnfoldCorrectionFactor_Likelihood[j] << " " << UnfoldCorrectionFactorError_Likelihood[j] << "  -->  error in % = " << UnfoldCorrectionFactorError_Likelihood[j]/UnfoldCorrectionFactor_Likelihood[j]*1e2 << endl;
  }
  
  return UnfoldCorrectionFactorError_Likelihood;
  
}