ICRC2019/UnfoldSpectrum.cc

Go to the documentation of this file.

/* Authors: I. Valiño, G. Rodriguez, Alan Coleman, Alex Schulz, F. Fenu
   based on an original code of H. Dembinski */
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <vector>
#include <cstring>
#include <sstream>
#include <TFile.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TF1.h>
#include "Math/WrappedTF1.h"
#include "Math/Integrator.h"
#include <TH1D.h>
#include <TNtupleD.h>
#include <TMath.h>
#include <TMinuit.h>
#include <TGraphAsymmErrors.h>
#include <TGraphErrors.h>
#include <TGraph.h>
#include <TAxis.h>
#include <TMatrixD.h>
#include <TVectorD.h>
#include <TApplication.h>
#include <TMinuit.h>
#include <TImage.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TLegend.h>

using namespace std;
TVectorD raw_lgEs, raw_nevents, raw_flux, raw_staterrlow, raw_staterrup;
TVectorD raw_lgEs_rebinned, raw_nevents_rebinned, raw_flux_rebinned, raw_staterrlow_rebinned, raw_staterrup_rebinned;
TAxis* thLgE;
TVectorD vecLgE;
TMatrixD kmatrix;
TAxis* rawAxis_lgE;
bool HeraldAnalysis;
double InverseSdCalibration(const double, TVectorD );
double ScaledErrorFunction(const double,  const std::vector<double>&);
void TestFormat(string, int);

#include "BiasFunctions.hh"
#include "ResolutionFunctions.hh"
#include "TriggerFunctions.hh"
#include "UnfoldUtilities.hh"
#include "FittingFunctions.hh"
#include "Statistics.hh"

string modelfunctions;
double IntXmin;
double IntXmax;
int IntegrationSteps;
double exposure;
TVectorD UnfoldCorrectionFactor;
TVectorD UnfoldCorrectionFactorMin;
TVectorD UnfoldCorrectionFactorMax;
TMatrixD UnfoldCorrectionFactorCoV;
TVectorD parsFit;
TVectorD parsErrorsFit;
TVectorD MinuitMinimization(double*, Int_t, bool);
void FitFCN(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t iflag);
double Likelihood(double* pars);
void PrintResults(TVectorD, TVectorD);

int main(int argc, char* argv[]) {
  if (2 != argc) {
    PrintUsage(argc);
    return 0;
  }
  LoadParam(argv[1]);
  string DataFileName = SDFile;
  modelfunctions = Model;
  Double_t SpectrumMin = mapParameters["SpectrumMin"];
  Double_t SpectrumMax = mapParameters["SpectrumMax"];
  Double_t SpectrumBinSize = mapParameters["SpectrumBinSize"];
  const Int_t BinNum = mapParameters["NumBins"];
  double SpectrumBins[BinNum + 1];
  const int Npar = mapParameters["NumParam"];
  double fParFitInit[Npar];
  parsFit.ResizeTo(Npar);
  parsErrorsFit.ResizeTo(Npar);
  IntXmin = mapParameters["MinIntegration"];
  IntXmax = mapParameters["MaxIntegration"];
  IntegrationSteps = mapParameters["IntegrationSteps"];

  if (CodeUsed == "Offline")
    HeraldAnalysis = false;
  else
    HeraldAnalysis = true;

  for (Int_t i = 0; i < BinNum + 1; i++) {
    ostringstream convert;
    convert << i + 1;
    string nameParameter = "r";
    nameParameter += convert.str();
    SpectrumBins[i] = mapParameters[nameParameter.c_str()];
  }

  cout << "Input File: " << DataFileName << endl;
  cout << "Output File: " << OutPutFileName << endl;
  cout << "Output File Rebinned: " << OutPutFileNameRebinned << endl;

  int hnbins = round((SpectrumMax - SpectrumMin) / SpectrumBinSize);
  TH1D* thEnergySpectrum = new TH1D("thEnergySpectrum", "", hnbins, SpectrumMin, SpectrumMax);
  TH1D* thEnergySpectrum_rebinned = new TH1D("thEnergySpectrum_rebinned", "", BinNum, SpectrumBins); //Check! index of SpectrumBins
  if (DeclinationDepentMatrix == "no")
    exposure = mapParameters["Exposure"];
  if (DeclinationDepentMatrix == "yes")
    exposure = DirectionalExposure(mapParameters["AugerLatitude"]*TMath::DegToRad(), mapParameters["MinDeclination"] * TMath::DegToRad(), mapParameters["MaxDeclination"] * TMath::DegToRad());
  rawAxis_lgE = thEnergySpectrum->GetXaxis();

  // Open file with Rec. events
  if (VersionInput == "ICRC2019")
    ReadEventListICRC2019(thEnergySpectrum, thEnergySpectrum_rebinned, DataFileName);
  else
    ReadEventList(thEnergySpectrum, thEnergySpectrum_rebinned, DataFileName);

  if ("no" != Verbose)
    cerr << "\nRead in " << thEnergySpectrum->GetEntries() << " events\n" << endl;

  // *****************************************************
  // RAW SPECTRUM
  // *****************************************************
  int Nbins = thEnergySpectrum->GetNbinsX();
  raw_lgEs.ResizeTo(Nbins), raw_nevents.ResizeTo(Nbins), raw_flux.ResizeTo(Nbins), raw_staterrlow.ResizeTo(Nbins), raw_staterrup.ResizeTo(Nbins);

  TGraphAsymmErrors* tgEnergySpectrum = new TGraphAsymmErrors();
  for (int i = 1; i <= Nbins; ++i) {
    double binW = thEnergySpectrum->GetBinWidth(i);
    double we = pow(10, thEnergySpectrum->GetBinCenter(i) + binW / 2) - pow(10, thEnergySpectrum->GetBinCenter(i) - binW / 2);
    double factor = exposure * we;
    double lgE = thEnergySpectrum->GetBinCenter(i);
    double content = thEnergySpectrum->GetBinContent(i) / factor;
    tgEnergySpectrum->SetPoint(i - 1, lgE, content);
    double plow, pup;
    poisson_uncertainty(thEnergySpectrum->GetBinContent(i), plow, pup);
    tgEnergySpectrum->SetPointError(i - 1, 0., 0., plow / factor, pup / factor);
    raw_lgEs[i - 1]    = lgE;
    raw_nevents[i - 1] = thEnergySpectrum->GetBinContent(i);
    raw_flux[i - 1]    = content;
    raw_staterrlow[i - 1] = plow / factor;
    raw_staterrup[i - 1] = pup / factor;
  }

  // *****************************************************************
  // Rebinned raw spectrum
  // *****************************************************************
  int Nbins_rebinned = thEnergySpectrum_rebinned->GetNbinsX();
  raw_lgEs_rebinned.ResizeTo(Nbins_rebinned), raw_nevents_rebinned.ResizeTo(Nbins_rebinned), raw_flux_rebinned.ResizeTo(Nbins_rebinned), raw_staterrlow_rebinned.ResizeTo(Nbins_rebinned), raw_staterrup_rebinned.ResizeTo(Nbins_rebinned);

  for (int i = 1; i <= Nbins_rebinned; ++i) {

    double binW_rebinned = thEnergySpectrum_rebinned->GetBinWidth(i);
    double we_rebinned = pow(10, thEnergySpectrum_rebinned->GetBinCenter(i) + binW_rebinned / 2) - pow(10, thEnergySpectrum_rebinned->GetBinCenter(i) - binW_rebinned / 2);
    double factor_rebinned = exposure * we_rebinned;
    double lgE_rebinned = thEnergySpectrum_rebinned->GetBinCenter(i);
    double content_rebinned = thEnergySpectrum_rebinned->GetBinContent(i) / factor_rebinned;

    double plow_rebinned, pup_rebinned;
    poisson_uncertainty(thEnergySpectrum_rebinned->GetBinContent(i), plow_rebinned, pup_rebinned);
    raw_lgEs_rebinned[i - 1]    = lgE_rebinned;
    raw_nevents_rebinned[i - 1] = thEnergySpectrum_rebinned->GetBinContent(i);
    raw_flux_rebinned[i - 1]    = content_rebinned;
    raw_staterrlow_rebinned[i - 1] = plow_rebinned / factor_rebinned;
    raw_staterrup_rebinned[i - 1] = pup_rebinned / factor_rebinned;
  }

  /*  Model parameters:   */
  /*  0: normalization; 1:  log10(energy) at ankle position; 2: spectral index before ankle; 3 spectra index after ankle;
      4: log10(energy) at the 1/2 suppression; 5: increment of the spectral index beyond suppression */

  TF1* FitFunction;
  FitFunction = new TF1("FittingFunction", FluxModels(modelfunctions.c_str()), SpectrumMin, SpectrumMax, Npar);


  if (modelfunctions == "StandardMultiSmooth") {
    if ("no" != Verbose)
      cerr << "Performing initial fit to StandardMultiSmooth\n";
    FitFunction->SetParameter(0, mapParameters["par0"]);
    FitFunction->SetParameter(1, mapParameters["par1"]);
    FitFunction->SetParameter(2, mapParameters["par2"]);
    FitFunction->SetParameter(3, mapParameters["par3"]);
    FitFunction->SetParameter(4, mapParameters["par4"]);
    FitFunction->SetParameter(5, mapParameters["par5"]);
    FitFunction->SetParameter(6, mapParameters["par6"]);
    FitFunction->SetParameter(7, mapParameters["par7"]);
  }
  if (modelfunctions == "StandardICRC2015") {
    if ("no" != Verbose)
      cerr << "Performing initial fit to StandardICRC2015\n";
    FitFunction->SetParameter(0, mapParameters["par0"]);
    FitFunction->SetParameter(1, mapParameters["par1"]);
    FitFunction->SetParameter(2, mapParameters["par2"]);
    FitFunction->SetParameter(3, mapParameters["par3"]);
    FitFunction->SetParameter(4, mapParameters["par4"]);
    FitFunction->SetParameter(5, mapParameters["par5"]);
  }
  if (modelfunctions == "Spectrum2015_2017SmoothAnkle") {
    if ("no" != Verbose)
      cerr << "Performing initial fit to Spectrum2015_2017SmoothAnkle\n";
    FitFunction->SetParameter(0, mapParameters["par0"]);
    FitFunction->SetParameter(1, mapParameters["par1"]);
    FitFunction->SetParameter(2, mapParameters["par2"]);
    FitFunction->SetParameter(3, mapParameters["par3"]);
    FitFunction->SetParameter(4, mapParameters["par4"]);
    FitFunction->SetParameter(5, mapParameters["par5"]);
  }
  if (modelfunctions == "InfillHard") {
    if ("no" != Verbose)
      cerr << "Performing initial fit to InfillHard\n";
    FitFunction->SetParameter(0, mapParameters["par0"]);
    FitFunction->SetParameter(1, mapParameters["par1"]);
    FitFunction->SetParameter(2, mapParameters["par2"]);
    FitFunction->SetParameter(3, mapParameters["par3"]);
  }
  if (modelfunctions == "SpectrumInfillICRC2019") {
    if ("no" != Verbose)
      cerr << "Performing initial fit to SpectrumInfillICRC2019\n";
    FitFunction->SetParameter(0, mapParameters["par0"]);
    FitFunction->SetParameter(1, mapParameters["par1"]);
    FitFunction->SetParameter(2, mapParameters["par2"]);
    FitFunction->SetParameter(3, mapParameters["par3"]);
    FitFunction->SetParameter(4, mapParameters["par4"]);
    FitFunction->SetParameter(5, mapParameters["par5"]);
    FitFunction->SetParameter(6, mapParameters["par6"]);

    if (mapParameters["MaxZenith"] > 40) {
      FitFunction->FixParameter(1, mapParameters["par1"]);  //lgEbreak
      FitFunction->FixParameter(3, mapParameters["par3"]);  //g0
      FitFunction->FixParameter(6, mapParameters["par6"]);  //dg
    }
  }

  if ("no" == Verbose)
    tgEnergySpectrum->Fit(FitFunction, "q");
  else
    tgEnergySpectrum->Fit(FitFunction);

  for (int i = 0; i < Npar; i++) {
    fParFitInit[i] = FitFunction->GetParameter(i);
  }

  // ***************************************************
  // UNFOLDING starts ....
  // ***************************************************
  /*  Numerical integration made using a Rienmann sum */
  //do unfolding, x = lgE, y = count

  int nbins = round((IntXmax - IntXmin) / thEnergySpectrum->GetBinWidth(1));
  int N = nbins * IntegrationSteps;

  thLgE = new TAxis(N, IntXmin, IntXmax);
  vecLgE.ResizeTo(N);
  for (int i = 1; i <= N; ++i)
    vecLgE[i - 1] = thLgE->GetBinCenter(i);

  TVectorD pCal = GetSDCalPars();

  double resolutionOffset[2] = {0, 0};
  kmatrix.ResizeTo(N, N);
  if (DeclinationDepentMatrix == "no")
    kmatrix = kResolutionMatrix(pCal, vecLgE, resolutionOffset);
  if (DeclinationDepentMatrix == "yes")
    kmatrix = kResolutionMatrixDD(pCal, vecLgE, mapParameters["AugerLatitude"] * TMath::DegToRad(), mapParameters["MinDeclination"] * TMath::DegToRad(), mapParameters["MaxDeclination"] * TMath::DegToRad());

  double lgExp = log10(exposure);

  // Set Initial guess parameters from fit
  double start[Npar];
  start[0] = fParFitInit[0] + lgExp;  // Special treatment for norm

  for (int ipar = 1; ipar < Npar; ipar++)
    start[ipar] = fParFitInit[ipar];

  UnfoldCorrectionFactor.ResizeTo(N);
  UnfoldCorrectionFactor = MinuitMinimization(start, Npar, true);


  // statistics errors due to unfolding.
  TVectorD unfoldstat(N);
  TVectorD dlgE(N);

  for (int i = 0; i < N; ++i) {

    unfoldstat[i] = sqrt(UnfoldCorrectionFactorCoV[i][i]);
    dlgE[i] = thLgE->GetBinCenter(1) / 2.;
  }
  TGraph* tgUnfoldStatValue = new TGraph(N, &vecLgE[0], &unfoldstat[0]);

  TCanvas* c1 = new TCanvas();
  c1->SetGridx(1);
  c1->SetGridy(1);

  TGraph* tgUnfold = new TGraph(N, &vecLgE[0], &UnfoldCorrectionFactor[0]);
  TGraphAsymmErrors* tgUnfoldSpectrum = new TGraphAsymmErrors();

  /*  Writing the resulting spectrum in an output file */

  ofstream fOut(OutPutFileName.c_str(), ios::out);
  fOut << "# lg(E/eV)  lgE_hw  raw_counts  exposure[km2 sr yr]  J [km-2 sr-1 yr-1 eV-1]  lower_sigma[J]_stat  upper_sigma[J]_stat   unfolding_correction    sigma[unfolding_correction]_stat" << endl;

  for (int i = 0; i < raw_lgEs.GetNoElements(); ++i) {

    double factorunfold = tgUnfold->Eval(raw_lgEs[i]);
    double statunc_factorunfold =   tgUnfoldStatValue->Eval(raw_lgEs[i]);
    double nmod = raw_flux[i] * factorunfold;
    double lowerror = raw_staterrlow[i] * factorunfold;
    double uperror = raw_staterrup[i] * factorunfold;

    fOut << raw_lgEs[i] << "  " << thEnergySpectrum->GetBinWidth(i) / 2. << "  " << raw_nevents[i] << "  " << exposure << "   " << nmod << "  " << lowerror << "  " << uperror << " " << factorunfold << " " << statunc_factorunfold << endl;

    tgUnfoldSpectrum->SetPoint(i, raw_lgEs[i], nmod);
    tgUnfoldSpectrum->SetPointError(i, 0., 0., lowerror, uperror);
  }

  ofstream fOutRebinned(OutPutFileNameRebinned.c_str(), ios::out);
  fOutRebinned << "# lg(E/eV)  lgE_hw  raw_counts  exposure[km2 sr yr]  J [km-2 sr-1 yr-1 eV-1]  lower_sigma[J]_stat  upper_sigma[J]_stat   unfolding_correction" << endl;
  for (int i = 0; i < raw_lgEs_rebinned.GetNoElements(); ++i) {
    double factorunfold_rebinned = tgUnfold->Eval(raw_lgEs_rebinned[i]);// factor is calculated on constant binning! This affects only raw.
    double statunc_factorunfold_rebinned =   tgUnfoldStatValue->Eval(raw_lgEs_rebinned[i]);
    double nmod_rebinned = raw_flux_rebinned[i] * factorunfold_rebinned;
    double lowerror_rebinned = raw_staterrlow_rebinned[i] * factorunfold_rebinned;
    double uperror_rebinned = raw_staterrup_rebinned[i] * factorunfold_rebinned;
    fOutRebinned << raw_lgEs_rebinned[i] << "  " << thEnergySpectrum_rebinned->GetBinWidth(i + 1) / 2. << "  " << raw_nevents_rebinned[i] << "  " << exposure << "   " << nmod_rebinned << "  " << lowerror_rebinned << "  " << uperror_rebinned << " " << factorunfold_rebinned << " " << statunc_factorunfold_rebinned << endl;
  }
  fOutRebinned.close();
  fOut.close();
  PrintResults(parsFit, parsErrorsFit);
  return 0;
} //close main function


TVectorD MinuitMinimization(double* start, Int_t Npar_Minim, bool computeCovStat = true) {
  const int N = Npar_Minim;
  TMinuit gMinuit(N);
  double arglist[10];
  int ierflg = 0;
  arglist[0] = -1; // Minuit not verbose
  gMinuit.SetPrintLevel(-1);

  arglist[0] = 0.5;  // Set error param for negative-log-likelihood
  gMinuit.mnexcm("SET  ERR", arglist, 1, ierflg);

  double step = 0.001;

  if (modelfunctions == "StandardICRC2015") {
    gMinuit.mnparm(0, "Norm", start[0], step, 0, 0, ierflg);
    gMinuit.mnparm(1, "lgEa", start[1], step, 0, 0, ierflg);
    gMinuit.mnparm(2, "g1", start[2], step, 0, 0, ierflg);
    gMinuit.mnparm(3, "g2", start[3], step, 0, 0, ierflg);
    gMinuit.mnparm(4, "lgEs", start[4], step, 0, 0, ierflg);
    gMinuit.mnparm(5, "Dgamma", start[5], step, 0, 0, ierflg);
  } else if (modelfunctions == "StandardMultiSmooth") {
    gMinuit.mnparm(0, "Norm", start[0], step, 0, 0, ierflg);
    gMinuit.mnparm(1, "lgEa", start[1], step, 0, 0, ierflg);
    gMinuit.mnparm(2, "lgEx", start[2], step, 0, 0, ierflg);
    gMinuit.mnparm(3, "lgEs", start[3], step, 0, 0, ierflg);
    gMinuit.mnparm(4, "g1", start[4], step, 0, 0, ierflg);
    gMinuit.mnparm(5, "g2", start[5], step, 0, 0, ierflg);
    gMinuit.mnparm(6, "g3", start[6], step, 0, 0, ierflg);
    gMinuit.mnparm(7, "g4", start[7], step, 0, 0, ierflg);
  } else if (modelfunctions == "InfillHard") {
    gMinuit.mnparm(0, "Norm", start[0], step, 0, 0, ierflg);
    gMinuit.mnparm(1, "lgEa", start[1], step, 0, 0, ierflg);
    gMinuit.mnparm(2, "g1", start[2], step, 0, 0, ierflg);
    gMinuit.mnparm(3, "g2", start[3], step, 0, 0, ierflg);
  } else if (modelfunctions == "Spectrum2015_2017SmoothAnkle") {
    gMinuit.mnparm(0, "Norm", start[0], step, 0, 0, ierflg);
    gMinuit.mnparm(1, "lgEa", start[1], step, 0, 0, ierflg);
    gMinuit.mnparm(2, "g1", start[2], step, 0, 0, ierflg);
    gMinuit.mnparm(3, "g2", start[3], step, 0, 0, ierflg);
    gMinuit.mnparm(4, "lgEs", start[4], step, 0, 0, ierflg);
    gMinuit.mnparm(5, "Dgamma", start[5], step, 0, 0, ierflg);
  } else if (modelfunctions == "SpectrumInfillICRC2019") {
    const bool is55 = (mapParameters["MaxZenith"] > 40);
    gMinuit.mnparm(0, "Norm", start[0], step, 0, 0, ierflg);
    gMinuit.mnparm(1, "lgEbreak", start[1], is55 ? 0 : step, 0, 0, ierflg);
    gMinuit.mnparm(2, "lgEankle", start[2], step, 0, 0, ierflg);
    gMinuit.mnparm(3, "g0", start[3], is55 ? 0 : step, 0, 0, ierflg);
    gMinuit.mnparm(4, "g1", start[4], step, 0, 0, ierflg);
    gMinuit.mnparm(5, "g2", start[5], step, 0, 0, ierflg);
    gMinuit.mnparm(6, "Dgamma", start[6], is55 ? 0 : step, 0, 0, ierflg);
  }

  gMinuit.SetFCN(FitFCN);
  // Now ready for minimization step
  arglist[0] = 1000;
  arglist[1] = 1.;

  gMinuit.mnexcm("MINIMIZE", arglist, 2, ierflg);
  if ("no" != Verbose)
    cout << "Minuit error flag: " << ierflg << endl;
  if (ierflg) {
    cout << "\t Spectrum fit failed during Minuit." << endl;
    cout << "\t Change initial conditions or change unfolding function" << endl;
    throw std::exception();
  }

  double fitParameters[N];
  double fitParErrors[N];
  for (int i = 0; i < N; i++)
    gMinuit.GetParameter(i, fitParameters[i], fitParErrors[i]);

  parsFit.ResizeTo(N);

  for (int i = 0; i < N; i++) {
    parsFit[i] = fitParameters[i];
    parsErrorsFit[i] = fitParErrors[i];
  }

  int ncol = vecLgE.GetNoElements();
  TVectorD CorrFactor(ncol);
  CorrFactor = GetCorrectionFactor(vecLgE, kmatrix, &fitParameters[0], modelfunctions.c_str());


  if (computeCovStat) {
    TMatrixD cov(N, N);
    gMinuit.mnemat(&cov[0][0], N);
    UnfoldCorrectionFactorCoV.ResizeTo(ncol, ncol);
    UnfoldCorrectionFactorCoV = propagate_covariance(GetCorrectionFactor, parsFit, cov, modelfunctions.c_str());
  }

  return CorrFactor;
}

void FitFCN(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t iflag) {
  f = Likelihood(par);
}


double Likelihood(double* pars) {

  int ncol = vecLgE.GetNoElements();
  TVectorD mwraws0 = FluxModel(modelfunctions.c_str(), vecLgE, pars); // dN/dE

  for (int i = 0; i < ncol; ++i)
    mwraws0[i] = mwraws0[i] * pow(10., vecLgE[i]) * log(10);//  dN/dlogE

  TVectorD mwraws = kmatrix * mwraws0;
  mwraws *= thLgE->GetBinWidth(1);

  int N = raw_nevents.GetNoElements();
  TVectorD mws(N);
  mws.Zero();

  for (int i = 0; i < ncol; ++i) {

    double ilgE = vecLgE[i];

    int bin = rawAxis_lgE->FindBin(ilgE);

    if (bin < 1)
      continue;

    if (bin > N)
      break;

    mws[bin - 1] += mwraws[i] * thLgE->GetBinWidth(1);   // expected number of events (delta in poisson distrib.)
  }

  TVectorD like(N);
  for (int i = 0; i < N; ++i)
    like[i] = mws[i] - raw_nevents[i] * log(mws[i]);

  // neg-log-likelihood fit assuming Poisson fluctuations
  return like.Sum();

}

void PrintResults(TVectorD parsFit, TVectorD parsErrorsFit) {
  cout << setprecision(5) << endl;

  double J0 = pow(10., parsFit[0]) / exposure;
  double errorJ0 = J0 * log(10) * parsErrorsFit[0];

  if (modelfunctions == "StandardICRC2015") {
    cout << endl;
    cout << "\t Resulting parameters fitting the regular array flux model (broken PL + smooth suppression):  \n" << endl;
    cout << "\t Flux normalization:   " << J0 << "  +- " << errorJ0 << endl;
    cout << "\t Eankle:               " << pow(10., parsFit[1]) << "  +- " << pow(10., parsFit[1])*log(10)*parsErrorsFit[1] << "  eV " << endl;
    cout << "\t Esuppression:         " << pow(10., parsFit[4]) << "  +- " << pow(10., parsFit[4])*log(10)*parsErrorsFit[4] << "  eV " << endl;
    cout << "\t gamma1                " << parsFit[2] << "  +- " << parsErrorsFit[2] << endl;
    cout << "\t gamma2                " << parsFit[3] << "  +- " << parsErrorsFit[3] << endl;
    cout << "\t Delta_gamma           " << parsFit[5] << "  +- " << parsErrorsFit[5] << endl;
    cout << "\n" << endl;
  } else if (modelfunctions == "StandardMultiSmooth") {
    cout << endl;
    cout << "\t Resulting parameters fitting the ICRC 2017 folding model (3 breaks PL):  \n" << endl;
    cout << "\t Flux normalization:   " << J0 << "  +- " << errorJ0 << endl;
    cout << "\t Eankle:               " << pow(10., parsFit[1]) << "  +- " << pow(10., parsFit[1])*log(10)*parsErrorsFit[1] << "  eV " << endl;
    cout << "\t Ex:               " << pow(10., parsFit[2]) << "  +- " << pow(10., parsFit[2])*log(10)*parsErrorsFit[2] << "  eV " << endl;
    cout << "\t Esuppression:         " << pow(10., parsFit[3]) << "  +- " << pow(10., parsFit[3])*log(10)*parsErrorsFit[3] << "  eV " << endl;
    cout << "\t gamma1                " << parsFit[4] << "  +- " << parsErrorsFit[4] << endl;
    cout << "\t gamma2                " << parsFit[5] << "  +- " << parsErrorsFit[5] << endl;
    cout << "\t gamma3                " << parsFit[6] << "  +- " << parsErrorsFit[6] << endl;
    cout << "\t gamma4                " << parsFit[7] << "  +- " << parsErrorsFit[7] << endl;
    cout << "\n" << endl;
  } else if (modelfunctions == "InfillHard") {
    cout << endl;
    cout << "\t Resulting parameters fitting the Infill array flux model:  \n" << endl;
    cout << "\t Flux normalization:   " << J0 << "  +- " << errorJ0 << endl;
    cout << "\t Eankle:               " << pow(10., parsFit[1]) << "  +- " << pow(10., parsFit[1])*log(10)*parsErrorsFit[1] << "  eV " << endl;
    cout << "\t gamma1                " << parsFit[2] << "  +- " << parsErrorsFit[2] << endl;
    cout << "\t gamma2                " << parsFit[3] << "  +- " << parsErrorsFit[3] << endl;
    cout << "\n" << endl;
  } else if (modelfunctions == "Spectrum2015_2017SmoothAnkle") {
    cout << endl;
    cout << "\t Resulting parameters fitting the regular array flux model (Smooth ankle + smooth suppression):  \n" << endl;
    cout << "\t Flux normalization:   " << J0 << "  +- " << errorJ0 << endl;
    cout << "\t Eankle:               " << pow(10., parsFit[1]) << "  +- " << pow(10., parsFit[1])*log(10)*parsErrorsFit[1] << "  eV " << endl;
    cout << "\t Esuppression:         " << pow(10., parsFit[4]) << "  +- " << pow(10., parsFit[4])*log(10)*parsErrorsFit[4] << "  eV " << endl;
    cout << "\t gamma1                " << parsFit[2] << "  +- " << parsErrorsFit[2] << endl;
    cout << "\t gamma2                " << parsFit[3] << "  +- " << parsErrorsFit[3] << endl;
    cout << "\t Delta_gamma           " << parsFit[5] << "  +- " << parsErrorsFit[5] << endl;
    cout << "\n" << endl;
  } else if (modelfunctions == "SpectrumInfillICRC2019") {
    cout << endl;
    cout << "\t Resulting parameters fitting the ICRC 2019 infill function:  \n" << endl;
    cout << "\t Flux normalization:   " << J0 << "  +- " << errorJ0 << endl;
    cout << "\t Ebreak:               " << pow(10., parsFit[1]) << "  +- " << pow(10., parsFit[1])*log(10)*parsErrorsFit[1] << "  eV " << endl;
    cout << "\t Eankle:               " << pow(10., parsFit[2]) << "  +- " << pow(10., parsFit[2])*log(10)*parsErrorsFit[2] << "  eV " << endl;
    cout << "\t gamma0                " << parsFit[3] << "  +- " << parsErrorsFit[3] << endl;
    cout << "\t gamma1                " << parsFit[4] << "  +- " << parsErrorsFit[4] << endl;
    cout << "\t gamma2                " << parsFit[5] << "  +- " << parsErrorsFit[5] << endl;
    cout << "\t dGamma                " << parsFit[6] << "  +- " << parsErrorsFit[6] << endl;
    cout << "\n" << endl;
  }

}