ICRC2017/UnfoldSpectrum.cc

Go to the documentation of this file.

/* Authors: I. Valiño, G. Rodriguez, Alan Coleman, Alex Schulz
   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;

const bool HeraldAnalysis = false;

#include "UnfoldUtilities.icc"
#include "FittingFunctions.icc"
#include "Statistics.icc"


//  Format for parameters is (Regular, Infill)

/*  energy limits of the spectrum */
const double SpectrumMin[2] = {18.4, 17.5};
const double SpectrumMax[2] = {20.2, 19.7};
const double SpectrumBinSize[2] = {0.1, 0.1};

/* energy integration limits */
const double IntXmin[2] = {17.9, 16.9};
const double IntXmax[2] = {20.9, 19.9};
const int IntegrationSteps = 10; // Number of integration steps per bin (should be >~ 10)

const bool DoArray[2] = {true, true};  // Which spectra to do

const char* RegularDataFileName = "Events_SD1500_ICRC2017_Offline.dat";
const char* InfillDataFileName = "Events_SD750_ICRC2017_Offline.dat";
//onst char* InfillDataFileName = "Events_SD750_ICRC2017_Herald.dat";
//const char* RegularDataFileName = "Events_SD1500_ICRC2017_Herald.dat";

const bool ResolutionFromData = false;  //True = Resolution calculated from data; False = from sim
const bool Rebinning[2] = {true, true};

// Re-binning for SD1500 and SD750
const double SpectrumBins[15][15] = {{18.4, 18.5, 18.6, 18.7, 18.8, 18.9, 19.0, 19.1, 19.2, 19.3, 19.4, 19.5, 19.7, 19.9, 20.2},
                                    {17.5, 17.6, 17.7, 17.8, 17.9, 18.0, 18.1, 18.2, 18.3, 18.4, 18.5, 18.7, 18.9, 19.2, 19.7}};
const int BinNum[2] = {14, 14};

// see FittingFunctions.h
const modeltype modelfunctions[2] = {StandardMultiSmooth, InfillHard};
// const modeltype modelfunctions[2] = {StandardICRC2015, InfillHard};



// Number of fit parameters per array
const int NParArray[2] = {modelfunctions[0] == StandardMultiSmooth ? 8 : 6, 4};
const int Npar = 8; //Set this to the larger of the two values above

int iarray = 0;

double fParFitInit[Npar];
TVectorD parsFit(Npar);
TVectorD parsErrorsFit(Npar);

double exposure;

double fParFit[Npar], fParFitErrors[Npar];

TVectorD UnfoldCorrectionFactor;
TVectorD UnfoldCorrectionFactorMin;
TVectorD UnfoldCorrectionFactorMax;
TMatrixD UnfoldCorrectionFactorCoV;

TVectorD MinuitMinimization(double* start, bool computeCovStat);
void FitFCN(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t iflag);
double Likelihood(double* pars);
void PrintResults();

int main(void)
{

  gROOT->SetStyle("Plain");
  gStyle->SetCanvasDefH(700);
  gStyle->SetCanvasDefW(900);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetCanvasColor(0);

  gStyle->SetPadTopMargin(0.05);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPadLeftMargin(0.18);
  gStyle->SetPadRightMargin(0.05);
  gStyle->SetPadColor(0);
  gStyle->SetPadBorderMode(0);

  gStyle->SetFrameBorderMode(0);
  gStyle->SetTitleFillColor(0);
  gStyle->SetTitleBorderSize(1);
  gStyle->SetStatColor(0);
  gStyle->SetStatBorderSize(1);

  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);

  gStyle->SetPalette(1, 0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetTitleX(0.1f);
  gStyle->SetTitleW(0.8f);

  gStyle->SetTitleFont(42, "X");
  gStyle->SetTitleFont(42, "Y");
  gStyle->SetLabelFont(42, "X");
  gStyle->SetLabelFont(42, "Y");

  gStyle->SetTitleAlign(21);


  gStyle->SetMarkerStyle(8);
  gStyle->SetMarkerSize(1.4);

  gStyle->SetTitleXOffset(1.);
  gStyle->SetTitleYOffset(1.4);
  gStyle->SetTitleSize(0.06, "X");
  gStyle->SetTitleSize(0.06, "Y");

  gStyle->SetTextSize(20);
  gStyle->SetStripDecimals(0);

  gStyle->SetErrorX(0.);
  gStyle->cd();

  double arrayExposure[2];
  ifstream fOpenExposure("exposure");
  for (iarray = 0; iarray < 2; iarray++)
    fOpenExposure >> arrayExposure[iarray];
  fOpenExposure.close();

  for (iarray = 0; iarray < 2; iarray++) {

    if (!DoArray[iarray])
      continue;

    int hnbins = round((SpectrumMax[iarray] - SpectrumMin[iarray]) / SpectrumBinSize[iarray]);
    exposure = arrayExposure[iarray];

    TH1D* thEnergySpectrum = new TH1D(Form("thEnergySpectrum%d", iarray), "", hnbins, SpectrumMin[iarray], SpectrumMax[iarray]);

    TH1D* thEnergySpectrum_rebinned;

    if (Rebinning[iarray]) {
      std::cerr << "Re-binning active for " << Form("%s", iarray ? "Infill" : "Regular") << " array data...\n";
      thEnergySpectrum_rebinned = new TH1D(Form("thEnergySpectrum_rebinned%d", iarray), "", BinNum[iarray], SpectrumBins[iarray]);
    }

    rawAxis_lgE = thEnergySpectrum->GetXaxis();

    // Open file with Rec. events

    /*  variables */
    long int AugerId;
    long int EventId;
    int nst;
    double lgE, S1000, dS1000, theta, phi, dec, ra;

    ifstream evtFile;
    evtFile.open((iarray ? InfillDataFileName : RegularDataFileName), ios::in);

    char header[1024];
    evtFile.getline(header, 1024);

    std::cerr << "Reading " << Form("%s", iarray ? "Infill" : "Regular") << " array data...\n";

    while (!evtFile.eof()) {
      evtFile >> AugerId >> EventId >> lgE >> S1000 >> dS1000 >> nst >> theta >> phi >> dec >> ra;

      if (HeraldAnalysis) // energies given in EeV
        lgE = log10(lgE) + 18;

      if (!evtFile.eof()){
        thEnergySpectrum->Fill(lgE);
        if (Rebinning[iarray])
          thEnergySpectrum_rebinned->Fill(lgE);
      }
    }
    evtFile.close();

    //*****************************************************
    // 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
    //*****************************************************************
    if (Rebinning[iarray]) {
      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("RegularFunction", FluxModels(modelfunctions[iarray]), SpectrumMin[iarray], SpectrumMax[iarray], NParArray[iarray]);

    if (!iarray) {

      if (modelfunctions[iarray] == StandardICRC2015) {

        FitFunction->SetParameter(0, -19.);
        FitFunction->SetParameter(1, 18.7);
        FitFunction->SetParameter(2, 3.3);
        FitFunction->SetParameter(3, 2.8);
        FitFunction->SetParameter(4, 19.8);
        FitFunction->SetParameter(5, 0.15);

      } else if (modelfunctions[iarray] == StandardMultiSmooth) {

        FitFunction->SetParameter(0, -22.5);
        FitFunction->SetParameter(1, 18.7);
        FitFunction->SetParameter(2, 19.15);
        FitFunction->SetParameter(3, 19.65);
        FitFunction->SetParameter(4, 3.3);
        FitFunction->SetParameter(5, 2.6);
        FitFunction->SetParameter(6, 3.0);
        FitFunction->SetParameter(7, 5.0);
      }

    } else {

      FitFunction->SetParameter(0, -19);
      FitFunction->SetParameter(1, 18.7);
      FitFunction->SetParameter(2, 3.3);
      FitFunction->SetParameter(3, 2.8);
    }

    FitFunction->SetLineColor(1);
    FitFunction->SetLineStyle(9);
    FitFunction->SetLineWidth(1);

    tgEnergySpectrum->Fit(FitFunction, "q");

    for (int i = 0; i < NParArray[iarray]; 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[iarray] - IntXmin[iarray]) / thEnergySpectrum->GetBinWidth(1));

    int N = nbins * IntegrationSteps;

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

    TVectorD pCal = GetSDCalPars(iarray);

    double resolutionOffset[2] = {0, 0};
    kmatrix.ResizeTo(N, N);
    kmatrix = kResolutionMatrix(pCal, vecLgE, resolutionOffset, ResolutionFromData, iarray);

    double lgExp = log10(exposure);

    // Set Initial guess parameters from fit
    double start[NParArray[iarray]];
    start[0] = fParFitInit[0] + lgExp;  // Special treatment for norm
    for (int ipar = 1; ipar < NParArray[iarray]; ipar++)
      start[ipar] = fParFitInit[ipar];

    UnfoldCorrectionFactor.ResizeTo(N);
    UnfoldCorrectionFactor = MinuitMinimization(start, 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);

    TGraphAsymmErrors* tgUnfoldStat = new TGraphAsymmErrors(N, &vecLgE[0], &UnfoldCorrectionFactor[0], &dlgE[0], &dlgE[0], &unfoldstat[0], &unfoldstat[0]);
    tgUnfoldStat->Draw("a3");
    tgUnfoldStat->SetMinimum(0);
    tgUnfoldStat->SetMaximum(1.2);
    tgUnfoldStat->SetFillColor(kGray);

    TGraph* tgUnfold = new TGraph(N, &vecLgE[0], &UnfoldCorrectionFactor[0]);
    tgUnfold->Draw("l");
    tgUnfold->SetMarkerStyle(3);

    tgUnfold->SetFillColor(kBlack);
    tgUnfold->SetLineWidth(3);

    tgUnfoldStat->GetYaxis()->CenterTitle(kTRUE);
    tgUnfoldStat->GetXaxis()->CenterTitle(kTRUE);

    tgUnfoldStat->GetYaxis()->SetTitle("Correction factor:   J_{true}/J_{folded}");
    tgUnfoldStat->GetXaxis()->SetTitle("log_{10}(E/eV)");

    tgUnfoldStat->GetXaxis()->SetRangeUser(SpectrumMin[iarray], SpectrumMax[iarray]);
    tgUnfoldStat->GetYaxis()->SetRangeUser(0.7, 1.05);


    TGraphAsymmErrors* tgUnfoldSpectrum = new TGraphAsymmErrors();

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

    ofstream fOut(Form("%sSpectrum.dat", iarray ? "Infill" : "Regular"), 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);
    }

    if (Rebinning[iarray]){
      ofstream fOutRebinned(Form("%sSpectrum_rebinned.dat", iarray ? "Infill" : "Regular"), 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]);
        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();
    

    TCanvas* c2 = new TCanvas();
    c2->SetLogy();
    tgEnergySpectrum->SetMarkerStyle(22);
    tgEnergySpectrum->SetMarkerSize(2.2);
    tgEnergySpectrum->SetMarkerColor(1);
    tgEnergySpectrum->GetYaxis()->CenterTitle(kTRUE);
    tgEnergySpectrum->GetXaxis()->CenterTitle(kTRUE);
    tgEnergySpectrum->GetYaxis()->SetTitle("J(E)  [eV^{-1} km^{-2} sr^{-1} yr^{-1}]");
    tgEnergySpectrum->GetXaxis()->SetTitle("log_{10}(E/eV)");
    tgEnergySpectrum->Draw("ap");

    tgUnfoldSpectrum->SetMarkerStyle(20);
    tgUnfoldSpectrum->SetMarkerSize(1.8);
    tgUnfoldSpectrum->SetMarkerColor(2);
    tgUnfoldSpectrum->SetLineColor(2);
    tgUnfoldSpectrum->Draw("p");

    double hymin = FitFunction->Eval(SpectrumMax[iarray]) / 10.;
    double hymax = FitFunction->Eval(SpectrumMin[iarray]) * 10.;

    tgEnergySpectrum->GetXaxis()->SetRangeUser(SpectrumMin[iarray], SpectrumMax[iarray]);
    tgEnergySpectrum->GetYaxis()->SetRangeUser(hymin, hymax);

    FitFunction->SetLineColor(2);
    FitFunction->SetLineStyle(9);
    FitFunction->SetLineWidth(1);
    tgUnfoldSpectrum->Fit(FitFunction, "q");

    TLegend* leg = new TLegend(0.2, 0.3, 0.5, 0.45);
    leg->SetFillColor(0);
    leg->SetLineColor(0);
    leg->SetBorderSize(1);
    leg->SetTextFont(42);
    leg->SetTextSize(0.04);
    leg->AddEntry(tgEnergySpectrum, "measured spectrum", "p");
    leg->AddEntry(tgUnfoldSpectrum, "unfolded spectrum", "p");
    leg->Draw();



    TImage* img1 = TImage::Create();
    img1->FromPad(c1);
    img1->WriteImage(Form("%sCorrectionFactor.png", iarray ? "Infill" : "Regular"));

    TImage* img2 = TImage::Create();
    img2->FromPad(c2);
    img2->WriteImage(Form("%sSpectrum.png", iarray ? "Infill" : "Regular"));

    PrintResults();
  }

  return 0;

} //close main function


TVectorD MinuitMinimization(double* start, bool computeCovStat = true)
{

  const int N = NParArray[iarray];

  TMinuit gMinuit(N);
  double arglist[10];
  int ierflg = 0;
  arglist[0] = -1; // Minuit not verbose
  //gMinuit.mnexcm("SET  PRI", arglist, 1, ierflg);
  //gMinuit.mnexcm("SET  NOW", arglist, 1, ierflg);
  gMinuit.SetPrintLevel(-1);

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

  double step = 0.01;

  gMinuit.mnparm(0, "Norm", start[0], step, 0, 0, ierflg);
  gMinuit.mnparm(1, "lgEa", start[1], step, 0, 0, ierflg);

  if (!iarray) {

    if (modelfunctions[iarray] == StandardICRC2015) {

      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[iarray] == StandardMultiSmooth) {

      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 {

    gMinuit.mnparm(2, "g1", start[2], step, 0, 0, ierflg);
    gMinuit.mnparm(3, "g2", start[3], step, 0, 0, ierflg);
  }

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

  gMinuit.mnexcm("MINIMIZE", arglist, 2, ierflg);

  if (ierflg)
    cout << "\t Spectrum fit failed during Minuit." << endl;


  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, &parsFit[0], modelfunctions[iarray]);

  if (computeCovStat) {

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

  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[iarray], vecLgE, pars);

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

  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()
{
  cout << setprecision(3) << endl;

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

  if (!iarray) {

    if (modelfunctions[iarray] == StandardICRC2015) {

      cout << endl;
      cout << "\t Resulting parameters fitting the regular array flux model (old 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 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[iarray] == StandardMultiSmooth) {

      cout << endl;
      cout << "\t Resulting parameters fitting the ICRC 2017 folding 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 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 {

    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;

  }
}