PRD2020/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>
#include <TFeldmanCousins.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;
TMatrixD parsCov(20, 20);
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);

//VV
int NPar;

//VV
double fChi2_lik;
int fNdf;
TVectorD Chi2_lik(double* pars);
int IerFlg;

//VV
const double EeV = 1e18;

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"];
  //VV
  NPar = Npar; // in common on the top of the code
  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 (MatrixType == "Standard")
    exposure = mapParameters["Exposure"];
  if (MatrixType == "Declination")
    exposure = DirectionalExposure(mapParameters["AugerLatitude"]*TMath::DegToRad(), mapParameters["MinDeclination"] * TMath::DegToRad(), mapParameters["MaxDeclination"] * TMath::DegToRad());
  if (MatrixType == "Zenith")
    exposure = DirectionalExposureZenith(mapParameters["MinZenith"],mapParameters["MaxZenith"]);
  
  rawAxis_lgE = thEnergySpectrum->GetXaxis();


  // Open file with Rec. events
  vector <double> vLogESD;
  if (VersionInput == "ICRC2019")
    ReadEventListICRC2019(vLogESD,thEnergySpectrum, thEnergySpectrum_rebinned, DataFileName);
  else
    ReadEventList(vLogESD,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 (modelfunctions == "SpectrumLipari") {
    if ("no" != Verbose)
      cerr << "Performing initial fit to SpectrumLipari\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"]);
    
    FitFunction->SetParameter(8, mapParameters["par8"]);
    FitFunction->SetParLimits (8,mapParameters["par8"],mapParameters["par8"]);
    FitFunction->SetParameter(9, mapParameters["par9"]);
    FitFunction->SetParLimits (9,mapParameters["par9"],mapParameters["par9"]);
    FitFunction->SetParameter(10, mapParameters["par10"]);
    FitFunction->SetParLimits (10,mapParameters["par10"],mapParameters["par10"]);
  }


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

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

  // set the upper limit at the 90% CL
  // we can't do before because the fit of the raw spectrum is a chi2 and
  // the upper limits changes the initial conditions
  Nbins = thEnergySpectrum->GetNbinsX();
  for (int i = 1; i <= Nbins; ++i) {
    if (thEnergySpectrum->GetBinContent(i)==0){
      double xerr_lo, xerr_up, yerr_lo, yerr_up;
      xerr_lo = tgEnergySpectrum->GetErrorXlow(i -1);
      xerr_up = tgEnergySpectrum->GetErrorXhigh(i -1);
      yerr_lo = tgEnergySpectrum->GetErrorYlow(i -1);
      yerr_up = tgEnergySpectrum->GetErrorYhigh(i -1);
      double plow0, pup0, plow1, pup1;
      poisson_uncertainty( thEnergySpectrum->GetBinContent(i), plow0, pup0 );
      FC_UpperLimits(      thEnergySpectrum->GetBinContent(i), plow1, pup1 );
      double fac = pup1/pup0;
      tgEnergySpectrum->SetPointError(i - 1, xerr_lo, xerr_up, yerr_lo, yerr_up*fac);
      double dummy = raw_staterrup[i - 1];
      raw_staterrup[i - 1] = dummy*fac;
    }
  }
  Nbins_rebinned = thEnergySpectrum_rebinned->GetNbinsX();
  for (int i = 1; i <= Nbins_rebinned; ++i) {
   if (thEnergySpectrum_rebinned->GetBinContent(i)==0){
      double plow0, pup0, plow1, pup1;
      poisson_uncertainty( thEnergySpectrum_rebinned->GetBinContent(i), plow0, pup0 );
      FC_UpperLimits(      thEnergySpectrum_rebinned->GetBinContent(i), plow1, pup1 );
      double fac = pup1/pup0;
      double dummy = raw_staterrup_rebinned[i - 1];
      raw_staterrup_rebinned[i - 1] = dummy*fac;
    }
  }

  
  // ***************************************************
  // 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 (MatrixType == "Standard")
    kmatrix = kResolutionMatrix(pCal, vecLgE, resolutionOffset);
  if (MatrixType == "Declination")
    kmatrix = kResolutionMatrixDD(pCal, vecLgE, mapParameters["AugerLatitude"] * TMath::DegToRad(), mapParameters["MinDeclination"] * TMath::DegToRad(), mapParameters["MaxDeclination"] * TMath::DegToRad());
  if (MatrixType == "Zenith")
    kmatrix = kResolutionMatrixZenith(pCal, vecLgE, resolutionOffset);
  

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


  // *******************************************************************
  // VV
  // write the resulting spectrum in the output files (ascii and .root)
  //
  // *******************************************************************

  // output root file
  string OutPutRootFileName = OutPutFileName+".root";
  TFile *fRootOut = new TFile(OutPutRootFileName.c_str(),"RECREATE");

  // VV
  // function with the fitted parameters
  for (int i=0; i<Npar; i++){
    if (i==0) {
      FitFunction->SetParameter(i,parsFit[i]-log10(exposure));
      FitFunction->SetParError(i,parsErrorsFit[i]-log10(exposure));
    } else {
      FitFunction->SetParameter(i,parsFit[i]);      
      FitFunction->SetParError(i,parsErrorsFit[i]);
    }
  }

  //VV
  TGraphAsymmErrors* gre_JRaw = new TGraphAsymmErrors();
  TGraphAsymmErrors* gre_JRawE3 = new TGraphAsymmErrors();
  TGraphAsymmErrors* gre_J = new TGraphAsymmErrors();
  TGraphAsymmErrors* gre_JE3 = new TGraphAsymmErrors();
  TGraph* gr_unf = new TGraph();
  TGraphAsymmErrors* gre_unf = new TGraphAsymmErrors();
  
  // VV: get the unfolding correction
  // new calculation of unfolding coefficients
  Nbins = rawAxis_lgE->GetNbins();
  TVectorD UnfoldCorrectionFactor_Likelihood(Nbins);
  UnfoldCorrectionFactor_Likelihood = GetCorrectionFactor_Likelihood(vecLgE, kmatrix, &parsFit[0], modelfunctions.c_str(), rawAxis_lgE);
  TVectorD UnfoldCorrectionFactorError_Likelihood(Nbins);
  UnfoldCorrectionFactorError_Likelihood = GetCorrectionFactorError_Likelihood(vecLgE, kmatrix, NPar, parsFit, parsErrorsFit, parsCov, modelfunctions.c_str(), rawAxis_lgE);

  // output ascii file with the spectrum data points
  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) {

    //VV
    // unfolding coefficients. 
    double factorunfold = UnfoldCorrectionFactor_Likelihood[i];
    double statunc_factorunfold = UnfoldCorrectionFactorError_Likelihood[i];

    double nmod = raw_flux[i] * factorunfold;
    double lowerror = raw_staterrlow[i] * factorunfold;
    double uperror = raw_staterrup[i] * factorunfold;

    int Nevents = raw_nevents[i];
    char buffer2 [100];
    int n2 = sprintf (buffer2,"%4.2f  %3.2f  %6d  %6.1f  %5.4le  %5.4le  %5.4le  %4.3f  %4.3f", raw_lgEs[i], thEnergySpectrum->GetBinWidth(i) / 2., Nevents , exposure, nmod ,lowerror, uperror, factorunfold, statunc_factorunfold);
    fOut << buffer2 << endl;

    // VV for the output root file
    double En = pow(10.,raw_lgEs[i]);
    double JRaw = raw_flux[i];
    double JRaw_loerr = raw_staterrlow[i];
    double JRaw_uperr = raw_staterrup[i];
    double J = JRaw * factorunfold;
    double J_loerr = JRaw_loerr * factorunfold;
    double J_uperr = JRaw_uperr * factorunfold;
    double fac = pow(En,3);
    gre_JRaw->SetPoint(i, En, JRaw);
    gre_JRaw->SetPointError(i, 0., 0., JRaw_loerr, JRaw_uperr);
    gre_JRawE3->SetPoint(i, En, JRaw*fac);
    gre_JRawE3->SetPointError(i, 0., 0., JRaw_loerr*fac, JRaw_uperr*fac);
    gre_J->SetPoint(i, En, J);
    gre_J->SetPointError(i, 0., 0., J_loerr, J_uperr);
    gre_JE3->SetPoint(i, En, J*fac);
    gre_JE3->SetPointError(i, 0., 0., J_loerr*fac, J_uperr*fac);
    gr_unf->SetPoint(i, En, factorunfold);
    gre_unf->SetPoint(i, En, factorunfold);
    gre_unf->SetPointError(i, 0., 0., statunc_factorunfold, statunc_factorunfold);
    
  }

  // spectrum data points for the other choice of the bins
  // note: not further unfolding fit is done  

  //VV
  TGraphAsymmErrors* gre_JRaw_reb = new TGraphAsymmErrors();
  TGraphAsymmErrors* gre_JRawE3_reb = new TGraphAsymmErrors();
  TGraphAsymmErrors* gre_J_reb = new TGraphAsymmErrors();
  TGraphAsymmErrors* gre_JE3_reb = new TGraphAsymmErrors();
  TGraph* gr_unf_reb = new TGraph();
  TGraphAsymmErrors* gre_unf_reb = new TGraphAsymmErrors();

  // VV: get the unfolding correction
  TAxis *rawAxis_lgE_reb = thEnergySpectrum_rebinned->GetXaxis();
  int nbins_reb = rawAxis_lgE_reb->GetNbins();
  TVectorD UnfoldCorrectionFactor_Likelihood_reb(nbins_reb);
  UnfoldCorrectionFactor_Likelihood_reb = GetCorrectionFactor_Likelihood(vecLgE, kmatrix, &parsFit[0], modelfunctions.c_str(), rawAxis_lgE_reb);
  TVectorD UnfoldCorrectionFactorError_Likelihood_reb(nbins_reb);
  UnfoldCorrectionFactorError_Likelihood_reb = GetCorrectionFactorError_Likelihood(vecLgE, kmatrix, NPar, parsFit, parsErrorsFit, parsCov, modelfunctions.c_str(), rawAxis_lgE_reb);

 // output ascii file with the spectrum data points
  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    sigma[unfolding_correction]_stat" << endl;

  for (int i = 0; i < raw_lgEs_rebinned.GetNoElements(); ++i) {
    
    //VV
    // new calculation of unfolding coefficients. 
    double factorunfold_rebinned = UnfoldCorrectionFactor_Likelihood_reb[i];
    double statunc_factorunfold_rebinned = UnfoldCorrectionFactorError_Likelihood_reb[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;

    int Nevents = raw_nevents_rebinned[i];
    char buffer2 [100];
    int n2 = sprintf (buffer2,"%4.2f  %3.2f  %6d  %6.1f  %5.4le  %5.4le  %5.4le  %4.3f  %4.3f", raw_lgEs_rebinned[i], thEnergySpectrum_rebinned->GetBinWidth(i) / 2., Nevents , exposure, nmod_rebinned ,lowerror_rebinned, uperror_rebinned, factorunfold_rebinned, statunc_factorunfold_rebinned);
    fOutRebinned << buffer2 << endl;

    // VV for the output root file
    double En = pow(10.,raw_lgEs_rebinned[i]);
    double JRaw = raw_flux_rebinned[i];
    double JRaw_loerr = raw_staterrlow_rebinned[i];
    double JRaw_uperr = raw_staterrup_rebinned[i];
    double J = JRaw * factorunfold_rebinned;
    double J_loerr = JRaw_loerr * factorunfold_rebinned;
    double J_uperr = JRaw_uperr * factorunfold_rebinned;
    double fac = pow(En,3);
    gre_JRaw_reb->SetPoint(i, En, JRaw);
    gre_JRaw_reb->SetPointError(i, 0., 0., JRaw_loerr, JRaw_uperr);
    gre_JRawE3_reb->SetPoint(i, En, JRaw*fac);
    gre_JRawE3_reb->SetPointError(i, 0., 0., JRaw_loerr*fac, JRaw_uperr*fac);
    gre_J_reb->SetPoint(i, En, J);
    gre_J_reb->SetPointError(i, 0., 0., J_loerr, J_uperr);
    gre_JE3_reb->SetPoint(i, En, J*fac);
    gre_JE3_reb->SetPointError(i, 0., 0., J_loerr*fac, J_uperr*fac);
    gr_unf_reb->SetPoint(i, En, factorunfold_rebinned);
    gre_unf_reb->SetPoint(i, En, factorunfold_rebinned);
    gre_unf_reb->SetPointError(i, 0., 0., statunc_factorunfold_rebinned, statunc_factorunfold_rebinned);
    
  }
  
  fOutRebinned.close();
  fOut.close();
  PrintResults(parsFit, parsErrorsFit);

  

  // *******************************************************************
  //
  // Print compact migration matrix
  //
  // *******************************************************************
  if(PrintCompact=="yes")
    {
      cout<<"PRODUCING COMPACT MATRIX "<<endl;
      TF1* FitFunctionFitted;
      FitFunctionFitted = new TF1("FittingFunctionFitted", FluxModels(modelfunctions.c_str()), SpectrumMin, SpectrumMax, Npar);
      for (int u=0;u<Npar;u++)
        {
          FitFunctionFitted->SetParameter(u,parsFit[u]);
        }
      kResolutionMatrixZenith_J(pCal, vecLgE, resolutionOffset,FitFunctionFitted);
    }



  

  // **************************************************************************************************************************************

  /*
    In the following we calculate:

    1) the expected energy spectrum in very small energy bins using the fit function;
    2) the unfolding coefficients in the energy bins used for the spectrum but in small energy steps 
       (useful to show with a continuous function the behaviour with energy of the unfolding coefficients.
       
    Note that the the calculation of the expected energy spectrum is done applying a correction that accounts for the finite 
    energy bin size used for the data. This correction is negligible for energy bins smaller than 0.05. 
    The calculation is done assuming that all the energy bins used for the data have the same size. 

    Note also that the calculation is not valid of the rebinned spectrum. It is not extended to it since sometimes 
    the rebinned spectrum is defined using energy bins of different widths.     

   */

  TGraph* gr_ffitE3_NoBinSizeCorr = new TGraph();
  TGraphAsymmErrors* gre_ffitE3_NoBinSizeCorr = new TGraphAsymmErrors();
  TGraph* gr_ffitE3 = new TGraph();
  TGraphAsymmErrors* gre_ffitE3 = new TGraphAsymmErrors();
  TGraph* gr_ffitRawE3 = new TGraph();
  TGraphAsymmErrors* gre_ffitRawE3 = new TGraphAsymmErrors();
  TGraph* gr_f_unf = new TGraph();
  TGraphAsymmErrors* gre_f_unf = new TGraphAsymmErrors();
  int NbinsX = thEnergySpectrum->GetNbinsX();
  double xmin = thEnergySpectrum->GetBinLowEdge(1);
  double xmax = thEnergySpectrum->GetBinLowEdge(NbinsX) + thEnergySpectrum->GetBinWidth(NbinsX);
  int nnx = 1000;
  double step = (xmax-xmin)/nnx;
  step = 0.01;
  nnx = (xmax-xmin)/step;

  for (int i=0; i<=nnx; i++){
    //
    double logE = xmin + step*i;
    double En = pow(10.,logE);
    //
    // correction for the finite bin size --> fcorr_bin
    // to correctly visualize Jexp for finite energy bin size
    // important for 0.1 bins - negligible for bins < 0.05
    double flogE1 = logE - thEnergySpectrum->GetBinWidth(1)/2;
    int nl = 10;
    double step_flogE = thEnergySpectrum->GetBinWidth(1)/nl;
    double Jint = 0.;
    for (int l=0; l<nl; l++){
      double flogE = flogE1 + step_flogE*l + step_flogE/2.;
      double fDeltaE = (pow(10.,flogE+step_flogE/2.)/1e18 - pow(10.,flogE-step_flogE/2.)/1e18)*1e18;
      Jint = Jint + FitFunction->Eval(flogE) * fDeltaE;
    }
    double DeltaE = (  pow(10.,logE+thEnergySpectrum->GetBinWidth(1)/2.)/1e18
                     - pow(10.,logE-thEnergySpectrum->GetBinWidth(1)/2.)/1e18 ) * 1e18;
    double Jint_approx = FitFunction->Eval(logE) * DeltaE;
    double fcorr_bin = Jint / Jint_approx;
    //    
    //    
    // expected energy spectrum without correction for the finite size of the energy bins
    double Jexp_nocorr = FitFunction->Eval(logE);
    double Jexp = Jexp_nocorr * fcorr_bin;
    //
    // error on the expected energy spectrum
    // calculate the derivate of the function with respect to the parameters 
    double dfdpar_ffit[Npar];
    double x[1]={logE};
    for (int j=0; j<Npar; j++){
      dfdpar_ffit[j] = FitFunction->GradientPar(j,x,0.0001);
    }
    double ffit_err = 0.;
    for (int j1=0; j1<Npar; j1++){
      for (int j2=0; j2<Npar; j2++){
        ffit_err = ffit_err + dfdpar_ffit[j1]*dfdpar_ffit[j2] * parsCov[j1][j2];  
      }
    }
    ffit_err = sqrt(ffit_err);
    double Jexp_err = ffit_err;
    //
    // graphs for the output
    //
    double fac = pow(En,3);
    //
    // spectrum not corrected for the finite size of the energy bins
    // it coincides with the fit function
    gr_ffitE3_NoBinSizeCorr->SetPoint(i, En, Jexp_nocorr*fac);
    gre_ffitE3_NoBinSizeCorr->SetPoint(i, En, Jexp_nocorr*fac);
    gre_ffitE3_NoBinSizeCorr->SetPointError(i, 0., 0., Jexp_err*fac, Jexp_err*fac);
    //
    // spectrum corrected for the finite size of the energy bins
    gr_ffitE3->SetPoint(i, En, Jexp*fac);
    gre_ffitE3->SetPoint(i, En, Jexp*fac);
    gre_ffitE3->SetPointError(i, 0., 0., Jexp_err*fac, Jexp_err*fac);
    //

    // ***
    // unfolding coefficients calculated for the spectrum energy bin size but in the small steps
    // of the energy as defined for the loop. Useful to show the behaviour of the unfolding coefficients
    // when it is graphically superimpose to them
    //
    flogE1        = logE - thEnergySpectrum->GetBinWidth(1)/2;
    double flogE2 = logE + thEnergySpectrum->GetBinWidth(1)/2;
    int MM = 1;
    TAxis *thLgE_thEnergySpectrumBin = new TAxis(MM, flogE1, flogE2);
    //
    TVectorD UnfoldCorrectionFactor_Likelihood_thEnergySpectrum(MM);
    UnfoldCorrectionFactor_Likelihood_thEnergySpectrum =
      GetCorrectionFactor_Likelihood(vecLgE, kmatrix, &parsFit[0], modelfunctions.c_str(), thLgE_thEnergySpectrumBin);
    //
    TVectorD UnfoldCorrectionFactorError_Likelihood_thEnergySpectrum(MM);
    UnfoldCorrectionFactorError_Likelihood_thEnergySpectrum =
    GetCorrectionFactorError_Likelihood(vecLgE, kmatrix, NPar, parsFit, parsErrorsFit, parsCov, modelfunctions.c_str(), thLgE_thEnergySpectrumBin);
    //     
    double unf_thEnergySpectrum = UnfoldCorrectionFactor_Likelihood_thEnergySpectrum[0];
    double unf_err_thEnergySpectrum = UnfoldCorrectionFactorError_Likelihood_thEnergySpectrum[0];
    gr_f_unf->SetPoint(i, En, unf_thEnergySpectrum);
    gre_f_unf->SetPoint(i, En, unf_thEnergySpectrum);
    gre_f_unf->SetPointError(i, 0., 0., unf_err_thEnergySpectrum, unf_err_thEnergySpectrum);
    //
    // ***

    // calculate the expected raw energy spectrum
    fac = pow(En,3);
    double Jexp_raw = Jexp/unf_thEnergySpectrum;
    double Jexp_raw_err = Jexp_err/unf_thEnergySpectrum;
    gr_ffitRawE3->SetPoint(i, En, Jexp_raw*fac);          
    gre_ffitRawE3->SetPoint(i, En, Jexp_raw*fac);         
    gre_ffitRawE3->SetPointError(i, 0., 0., Jexp_raw_err*fac, Jexp_raw_err*fac);
    
  }
 
  // **************************************************************************************************************************************

  
  // tree with the events
  TTree *tree_ev = new TTree("tree_ev","events");
  double ESD = 0.;
  tree_ev->Branch("ESD",&ESD,"ESD/D");
  for (int l=0; l<vLogESD.size(); l++){
    double lgESD = vLogESD.at(l);
    ESD = pow(10.,lgESD);
    tree_ev->Fill();
    //cout << lgESD << endl;
  }
  
  // VV
  // parameters to be saved for the output
  // other info for the output root file
  // to be 
  TVectorD tExp(1);
  tExp(0) = exposure;
  TVectorD tPar(Npar), tParErr(Npar);
  for (int i=0; i<Npar; i++){
    double xpar, xpar_err;
    if (i==0) {
      xpar = pow(10., parsFit[0]) / exposure;
      xpar_err = xpar * log(10) * parsErrorsFit[0];
    }
    if (modelfunctions == "SpectrumLipari") {
      if (i==1 || i==2 || i==3) {
        xpar = pow(10., parsFit[i]);
        xpar_err = pow(10., parsFit[i])*log(10)*parsErrorsFit[i];
      }
      if (i>3){
        xpar = parsFit[i];
        xpar_err = parsErrorsFit[i];
      }
    }
    tPar(i) = xpar;
    tParErr(i) = xpar_err;
  }
        
  TVectorD tChi2fit(1);
  TVectorD tNdf(1);
  tChi2fit(0) = fChi2_lik;
  tNdf(0) = fNdf;
  TVectorD tIerFlg(1);
  tIerFlg(0) = IerFlg;


  // VV
  // output root file
  //
  tree_ev->Write();
  //
  thEnergySpectrum->Write();
  thEnergySpectrum_rebinned->Write();
  //
  gre_JRaw->Write("gre_JRaw");
  gre_JRawE3->Write("gre_JRawE3");
  gre_J->Write("gre_J");
  gre_JE3->Write("gre_JE3");
  gr_unf->Write("gr_unf");
  gre_unf->Write("gre_unf");
  //
  gre_JRaw_reb->Write("gre_JRaw_reb");
  gre_JRawE3_reb->Write("gre_JRawE3_reb");
  gre_J_reb->Write("gre_J_reb");
  gre_JE3_reb->Write("gre_JE3_reb");
  gr_unf_reb->Write("gr_unf_reb");
  gre_unf_reb->Write("gre_unf_reb");
  //
  gr_f_unf->Write("gr_f_unf");
  gre_f_unf->Write("gre_f_unf");
  gr_ffitE3->Write("gr_ffitE3");
  gre_ffitE3->Write("gre_ffitE3");
  gr_ffitRawE3->Write("gr_ffitRawE3");
  gre_ffitRawE3->Write("gre_ffitRawE3");
  gr_ffitE3_NoBinSizeCorr->Write("gr_ffitE3_NoBinSizeCorr");
  gre_ffitE3_NoBinSizeCorr->Write("gre_ffitE3_NoBinSizeCorr");
  //
  tPar.Write("par"); 
  tParErr.Write("par_err");
  tChi2fit.Write("Chi2");
  tNdf.Write("Ndf");
  tIerFlg.Write("IerFlg");
  //
  fRootOut->Close();
  //

  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);
  }
  else if (modelfunctions == "SpectrumLipari") {

    /* 
       VV
 
       Choice of the starting values:

       1) a first chi2 fit of the raw spectrum is done taking the starting values 
         stored in the configuration file 

       2) the values of the parameters output of the above fit are taken as 
         starting values for the unfolding fit

       The procedure to set the starting values of the unfolding fit is not good because, 
       at the ankle, the raw spectrum is different from the undolded one and a chi2 fit 
       is not good at the highest energies where one has to use a poissonian p.d.f..

       This generates several problems and ambiguities and sometimes, when the starting 
       values are cleary wrong, the fit does not converge. The problem emerges when many 
       parameters are fitted (like to include the new feature at 10 EeV) and when the statistics 
       is reduced (e.g. spectrum in different declination bands).

       For the above reasons, sometimes the starting values in the configuration 
       file are different and sometimes are hardcoded. The starting values below 
       are the ones that have been used for the PRD/PRL papers.

       For the future it is highly recommended to implement a poissonian fit of the raw 
       spectrum and to not use the output of the fit of the raw spectrum to set the 
       starting values for the fit of the unfolded spectrum. 
       
     */
    
    gMinuit.mnparm(0, "Norm", start[0], step, 0, 0, ierflg);
    gMinuit.mnparm(1, "lgE12", start[1], step, 0, 0, ierflg);
    gMinuit.mnparm(2, "lg23", start[2], step, 0, 0, ierflg);
    gMinuit.mnparm(3, "lg34", 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);
    gMinuit.mnparm(8, "w12", start[8], 0, start[8], start[8], ierflg);
    gMinuit.mnparm(9, "w23", start[9], 0, start[9], start[9], ierflg);
    gMinuit.mnparm(10, "w34", start[10], 0, start[10], start[10], ierflg);


    Double_t maxDecl = mapParameters["MaxDeclination"];
    Double_t minDecl = mapParameters["MinDeclination"];
    
    if (minDecl==-90. && maxDecl==-42.5){
      gMinuit.mnparm(0, "Norm", -13.539070, step, 0, 0, ierflg);
      gMinuit.mnparm(1, "lgE12", 18.711274, step, 0, 0, ierflg);
      gMinuit.mnparm(2, "lg23", 19.170022, step, 0, 0, ierflg);
      gMinuit.mnparm(3, "lg34", 19.660508, step, 0, 0, ierflg);
      gMinuit.mnparm(4, "g1", 3.313246, step, 0, 0, ierflg);
      gMinuit.mnparm(5, "g2", 2.570759, step, 0, 0, ierflg);
      gMinuit.mnparm(6, "g3", 3.070131, step, 0, 0, ierflg);
      gMinuit.mnparm(7, "g4", 4.886723, step, 0, 0, ierflg);
    }
    
    if (minDecl==-17.3 && maxDecl==25.){
      gMinuit.mnparm(1, "lgE12", log10(5.*EeV),  step, log10(3.*EeV),  log10(7.*EeV), ierflg);
      gMinuit.mnparm(2, "lgE23", log10(12.*EeV), step, log10(8.*EeV), log10(16.*EeV), ierflg);
      gMinuit.mnparm(3, "lgE34", log10(50.*EeV), step, log10(20.*EeV), log10(70.*EeV), ierflg);
      gMinuit.mnparm(4, "g1", 3.3, step, 0, 0, ierflg);
      gMinuit.mnparm(5, "g2", 2.5, step, 0, 0, ierflg);
      gMinuit.mnparm(6, "g3", 3.2, step, 0, 0, ierflg);
      gMinuit.mnparm(7, "g4", 5.0, step, 3., 7., ierflg);
    }

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

  gMinuit.mnexcm("MINIMIZE", arglist, 2, ierflg);
  //gMinuit.mnexcm("MIGRAD", arglist, 2, ierflg);
  IerFlg = 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];
  }

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

  if (computeCovStat) {

    // covariance matrix
    // attenzione
    // la matrice viene riempita shiftando nelle prime posizioni i valori dei parametri effettivamente fittati (non quelli fissi)
    // di conseguenza l'indice della matrice di covarianza non corrisponde a quello dei parametri
    // bisogna ritrovare la corrispondenza sapendo quali parametri sono stati fittati
    //
    // nota che l'ordine di riempimento viene preservato, ovvero:
    // il primo parametro libero che si trova nel vettore corrisponde al primo che si trova nella matrice, etc...

    // matrice come esce da minuit
    TMatrixD cov_minuit(N, N);
    gMinuit.mnemat(&cov_minuit[0][0], N);

    // matrice (da definire) la cui posizione degli elementi corrisponde a quelli del vettore dei parametri 
    TMatrixD cov(N, N);
    for (int j1=0; j1<N; j1++){
      for (int j2=0; j2<N; j2++){
        cov[j1][j2] = 0.;
      }
    }

    int NFreePars = gMinuit.GetNumFreePars();
    int iPar_cov[NFreePars];
    int iParFree = -1;

    for (int i=0; i<N; i++){
      //cout << i << " " <<  fitParameters[i] << " +/- " << fitParErrors[i] << "  -->  " << TMath::Abs(fitParErrors[i]/fitParameters[i]) << endl;
      if (TMath::Abs(fitParErrors[i]/fitParameters[i])>1e-40) { // fitted parameters. To be improved with a clever method without fitParameters at the denominator
        iParFree++;
        // corrispondenza tra l'indice di cov_minuit (iPar_cov) con quello della matrice cov (i) 
        iPar_cov[iParFree] = i;
        //cout << " --> " << i << "  --> " << iParFree << endl;
       }
    }
    //cout << endl;

    // loop only on the free parameters of parcov_dummy
    for (int j1=0; j1<NFreePars; j1++){
      for (int j2=0; j2<NFreePars; j2++){
        int i1 = iPar_cov[j1];
        int i2 = iPar_cov[j2];
        cov[i1][i2] = cov_minuit[j1][j2];
      }
    }


    UnfoldCorrectionFactorCoV.ResizeTo(ncol, ncol);
    UnfoldCorrectionFactorCoV = propagate_covariance(GetCorrectionFactor, parsFit, cov, modelfunctions.c_str());
    for (int i1=0; i1<N; i1++){
      for (int i2=0; i2<N; i2++){
        parsCov[i1][i2] = cov[i1][i2];
        if (i1==i2){
          // check
          //cout << i1 << " " << i2 << " " <<  sqrt(cov[i1][i2]) << " " << fitParErrors[i1] <<  endl;
          if (TMath::Abs(fitParErrors[i1])>0. && TMath::Abs(sqrt(cov[i1][i2])/fitParErrors[i1]-1.)>1e-20 )
            cout << i1 << " " << i2 << " ERROR --> " << TMath::Abs(sqrt(cov[i1][i2])/fitParErrors[i1]) << endl;
        }
      }
    }

    // print the correlation coefficient
    printf("\n");
    TMatrixD rho(N, N);
    for (int j1=0; j1<N; j1++){
      for (int j2=0; j2<N; j2++){
        rho[j1][j2] = 0.;
        if (fitParErrors[j1]>0. && fitParErrors[j2]>0.)
          rho[j1][j2] = parsCov[j1][j2]/fitParErrors[j1]/fitParErrors[j2];
      }
    }

    printf("%5.0d  ",0);
    for (int i2=0; i2<N; i2++)
      printf("%5.0d  ",i2+1);
    printf("\n");
    for (int i1=0; i1<N; i1++){
      printf("%5.0d  ",i1+1);
      for (int i2=0; i2<N; i2++){
        printf("%5.2f  ",rho[i1][i2]);
      }
      printf("\n");
    }
    
    
  }


  //
  // it calculates also the deviance --> fChi2_lik
  double ff = Likelihood(fitParameters); 
  fNdf = raw_nevents.GetNoElements() - gMinuit.GetNumFreePars();
  double pvalue = TMath::Prob(fChi2_lik,(int) fNdf);
  double chi2prob = pvalue*1e2;
  double nsigma = sqrt(2.) * TMath::ErfInverse(1.-pvalue);
  cout << endl;
  printf(" -ln(L) = %f \n",ff);
  printf(" Chi2 = %f \n",fChi2_lik);
  printf(" Ndf = %d  \n",(int) fNdf);
  printf(" Chi2/Ndf = %f    Chi2 prob = %f %%   -   p-value %le \n", fChi2_lik/fNdf, chi2prob, pvalue  );
  printf(" n x sigma = %f \n", nsigma);
  cout << endl;

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


  // deviance
  TVectorD Dev(N);
  for (int i = 0; i < N; ++i){
    double fN = raw_nevents[i];
    double fNexp =  mws[i];
    if (raw_nevents[i]>0){
      Dev[i] = 2. * ( (fNexp - fN * log(fNexp)) - (fN - fN * log(fN)) );
    } else {
      Dev[i] = 2. * ( (fNexp - fN * log(fNexp)) - 0.                  );
    }      
  }
  fChi2_lik = Dev.Sum();
  
  // 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;
  }
  else if (modelfunctions == "SpectrumLipari") {
    cout << endl;
    // VV
    cout << "\n" << endl;
    printf(" J0       = (%le +/- %le) \n",J0,errorJ0);
    printf(" E12      = (%3.1f +/- %3.1f) x 1e18 eV  \n", pow(10., parsFit[1])/EeV,     pow(10., parsFit[1])*log(10)*parsErrorsFit[1]/EeV); 
    printf(" E23      = (%3.1f +/- %3.1f) x 1e18 eV  \n", pow(10., parsFit[2])/EeV,     pow(10., parsFit[2])*log(10)*parsErrorsFit[2]/EeV);
    printf(" E34      = (%3.0f +/- %3.0f) x 1e18 eV  \n", pow(10., parsFit[3])/EeV,     pow(10., parsFit[3])*log(10)*parsErrorsFit[3]/EeV);
    printf(" gamma1   = %3.2f +/- %3.2f  \n",parsFit[4],parsErrorsFit[4]); 
    printf(" gamma2   = %3.2f +/- %3.2f  \n",parsFit[5],parsErrorsFit[5]); 
    printf(" gamma3   = %3.2f +/- %3.2f  \n",parsFit[6],parsErrorsFit[6]); 
    printf(" gamma4   = %3.1f +/- %3.1f  \n",parsFit[7],parsErrorsFit[7]); 
    printf(" w12        = %3.2f +/- %3.2f  \n",parsFit[8],parsErrorsFit[8]); 
    printf(" w23        = %3.2f +/- %3.2f  \n",parsFit[9],parsErrorsFit[9]); 
    printf(" w34        = %3.2f +/- %3.2f  \n",parsFit[10],parsErrorsFit[10]); 
    cout << "\n" << endl;
  }

}