EnergyFitter.cc

Go to the documentation of this file.

#include "EnergyFitter.h"
#include "diagonalMatrix.h"
#include "lowerTriangularMatrix.h"
#include "upperTriangularMatrix.h"
#include "matrixInversionExceptions.h"
#include "CherenkovFluorescenceMatrix.h"

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/GaisserHillas4Parameter.h>

#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/Pixel.h>
#include <fevt/FEvent.h>
#include <fevt/EyeRecData.h>

#include <fwk/CentralConfig.h>

#include <utl/Reader.h>
#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/TabulatedFunctionErrors.h>

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Camera.h>
#include <fdet/Telescope.h>
#include <fdet/Pixel.h>

#include <string>
#include <sstream>
#include <iostream>
#include <iomanip>

#include <TMinuit.h>
#include <TMath.h>

using namespace utl;
using namespace evt;
using namespace fevt;
using namespace fwk;
using namespace det;
using namespace fdet;
using namespace oBLAS;
using namespace std;


namespace FdProfileReconstructorKG {

  // static stuff for MINUIT

  EnergyFitter::EMinimizationType EnergyFitter::fMinimizationMethod = eModifiedLeastSquares;
  EnergyFitter::EParConstraints EnergyFitter::fFixX0 = eFree;
  EnergyFitter::EParConstraints EnergyFitter::fFixLambda = eFree;
  EnergyFitter::EParConstraints EnergyFitter::fFixkUniv = eFree;
  unsigned int EnergyFitter::fnBins = 0;
  double* EnergyFitter::fdEdX = nullptr;
  double* EnergyFitter::fvdEdX = nullptr;
  double* EnergyFitter::fLight = nullptr;
  double* EnergyFitter::fvLight = nullptr;
  double* EnergyFitter::fvBackGround = nullptr;
  double* EnergyFitter::fDepth = nullptr;
  double* EnergyFitter::fChFlDepth = nullptr;
  double EnergyFitter::fDefX0 = 0;
  double EnergyFitter::fVarX0 = 0;
  double EnergyFitter::fDefLambda = 0;
  double EnergyFitter::fVarLambda = 0;

  double EnergyFitter::fDefA1 = 0;
  double EnergyFitter::fDefA2 = 0;
  double EnergyFitter::fVarK = 0;

  const oBLAS::lowerTriangularMatrix* EnergyFitter::fChFl = nullptr;
  double EnergyFitter::fDiaphragmArea = 0;
  double EnergyFitter::fGainVariance = 0;
  double EnergyFitter::fPhotonToPhotoElectron = 0;


  bool
  EnergyFitter::Init()
  {
    Branch topB = fwk::CentralConfig::GetInstance()->GetTopBranch("FdProfileReconstructor");
    Branch eCalcB = topB.GetChild("profileFitting");

    Branch x0Branch = eCalcB.GetChild("X0");
    auto fitStrategy = x0Branch.GetChild("strategy").Get<string>();
    fFixX0 = (fitStrategy == "constrained") ? eConstrained : ((fitStrategy == "free") ? eFree : eFixed);
    x0Branch.GetChild("mean").GetData(fDefX0);
    fDefX0 /= g/cm2;
    const auto x0RMS = x0Branch.GetChild("sigma").Get<double>();
    fVarX0 = Sqr(x0RMS / (g/cm2));

    Branch lambdaBranch = eCalcB.GetChild("lambda");
    fitStrategy = lambdaBranch.GetChild("strategy").Get<string>();
    fFixLambda = (fitStrategy == "constrained") ? eConstrained : ((fitStrategy == "free") ? eFree : eFixed);
    lambdaBranch.GetChild("mean").GetData(fDefLambda);
    fDefLambda /= g/cm2;
    const double lambdaRMS = lambdaBranch.GetChild("sigma").Get<double>();
    fVarLambda = Sqr(lambdaRMS / (g/cm2));

    Branch kUnivBranch = eCalcB.GetChild("kUniv");
    fitStrategy = kUnivBranch.GetChild("strategy").Get<string>();
    fFixkUniv = (fitStrategy == "constrained") ? eConstrained : ((fitStrategy == "free") ? eFree : eFixed);
    kUnivBranch.GetChild("A1").GetData(fDefA1);
    fDefA1 /= g/cm2;
    kUnivBranch.GetChild("A2").GetData(fDefA2);
    fDefA2 /= g/cm2;
    const auto kRMS = kUnivBranch.GetChild("sigma").Get<double>();
    fVarK = Sqr(kRMS / (g/cm2));

    eCalcB.GetChild("verbosity").GetData(fVerbosity);

    const auto miniMeth = eCalcB.GetChild("minimizationMethod").Get<string>();
    if (miniMeth == "maximumLikelihood")
      fMinimizationMethod = eMaximumLikelihood;
    else if (miniMeth == "leastSquares")
      fMinimizationMethod = eLeastSquares;
    else if (miniMeth == "modifiedLeastSquares")
      fMinimizationMethod = eModifiedLeastSquares;

    if (fVerbosity > -1) {
      std::ostringstream info;
      info << "\n\t X0 is ";
      switch (fFixX0) {
      case eFree:
        info << "free\n";
        break;
      case eFixed:
        info << "fixed at " << fDefX0 << " g/cm2\n";
        break;
      default:
        info << "constrained at " << fDefX0 << " +/- " << sqrt(fVarX0) << " g/cm2\n";
        break;
      }

      info << "\t lambda is ";
      switch (fFixLambda) {
      case eFree:
        info << "free\n";
        break;
      case eFixed:
        info << "fixed at " << fDefLambda << " g/cm2\n";
        break;
      default:
        info << "constrained at " << fDefLambda << " +/- " << sqrt(fVarLambda) << " g/cm2\n";
        break;
      }

      info << "\t kUniv is ";
      switch (fFixkUniv) {
      case eFree:
        info << "free\n";
        break;
      case eFixed:
        info << "fixed at " << fDefA1 << " g/cm2, " << fDefA2 << " g/cm2\n";
        break;
      default:
        info << "constrained at (" << fDefA1 << " " << fDefA2 << ") +/- " << sqrt(fVarK) << " g/cm2\n";
        break;
      }

      info << "\t minimization method: " << miniMeth;

      INFO(info);
    }

    return true;
  }


  // determines energy for each eye of "event"
  bool
  EnergyFitter::Run(fevt::Eye& eye, CherenkovFluorescenceMatrix* const chfl)
  {
    fChFlPtr = chfl;

    if (InitializeGHFit(eye, 0) &&  // starting parameters
        FitProfile(eEnergyDeposit) &&  // fit two-parameter GH to dEdX profile
        FitProfile(eLightAllParam)) {  // fit four-parameter GH to light profile

      FillProfilesAtAperture(eye);

      if (CalculateEnergy(eye)) {
        CleanUp();
        return true;
      }
    }
    CleanUp();
    return false;
  }


  bool
  EnergyFitter::InitializeGHFit(const fevt::Eye& eye, const int iMode)
  {
    // Cleanup() ?
    fdEdX = nullptr;
    fvdEdX = nullptr;
    fDepth = nullptr;
    fLight = nullptr;
    fvLight = nullptr;
    fvBackGround = nullptr;
    fChFl = nullptr;
    fChFlDepth = nullptr;
    fChisqGH = 0;
    fnDof = 0;

    const auto* depthPtr = fChFlPtr->GetDepth();
    if (!depthPtr || depthPtr->size() < 1)
      return false;

    const auto& depth = *depthPtr;

    fChFlDepth = new double[depth.size()];
    for (unsigned int i = 0, n = depth.size(); i < n; ++i)
      fChFlDepth[i] = depth[i];

    // fill light profile
    if (!eye.GetRecData().HasLightFlux())
      return false;

    const auto& lightFlux = eye.GetRecData().GetLightFlux(fevt::FdConstants::eTotal);
    const auto& bgFlux = eye.GetRecData().GetLightFlux(fevt::FdConstants::eBackground);

    fnBins = lightFlux.GetNPoints();
    fLight = new double[fnBins];
    fvLight = new double[fnBins];
    fvBackGround = new double[fnBins];
    for (unsigned int i = 0; i < fnBins; ++i) {
      fLight[i] = lightFlux.GetY(i);
      fvLight[i] = Sqr(lightFlux.GetYErr(i));
      fvBackGround[i] = Sqr(bgFlux.GetYErr(i));
    }

    // fill dEdX profile
    const auto& recShower = eye.GetRecData().GetFRecShower();
    if (!recShower.HasEnergyDeposit())
      return false;

    const auto& dedxProf = recShower.GetEnergyDeposit();

    if (fnBins != dedxProf.GetNPoints()) {
      ostringstream err;
      err << "Inconsistent binning dedx/light: " << fnBins
          << ' ' << dedxProf.GetNPoints();
      ERROR(err);
      return false;
    }

    fdEdX = new double[fnBins];
    fvdEdX = new double[fnBins];
    fDepth = new double[fnBins];

    // work with decent units
    const double kXconvert = g/cm2;
    const double kEconvert = PeV/kXconvert;
    const double kEconvert2 = Sqr(kEconvert);

    for (unsigned int i = 0; i < fnBins; ++i) {
      fdEdX[i] = dedxProf.GetY(i) / kEconvert;
      fvdEdX[i] = dedxProf.GetYErr(i) * dedxProf.GetYErr(i) / kEconvert2;
      fDepth[i] = dedxProf.GetX(i) / kXconvert;
    }

    const auto& detFD = det::Detector::GetInstance().GetFDetector();
    const double referenceLambda = detFD.GetReferenceLambda();

    for (auto tIt = eye.TelescopesBegin(ComponentSelector::eHasData),
         end = eye.TelescopesEnd(ComponentSelector::eHasData);
         tIt != end; ++tIt) {
      double avgPEfactor = 0;
      int countPix = 0;
      for (auto pIt = tIt->PixelsBegin(ComponentSelector::eHasData),
           end = tIt->PixelsEnd(ComponentSelector::eHasData);
           pIt != end; ++pIt) {
        const auto& detPixel = detFD.GetPixel(*pIt);
        avgPEfactor += detPixel.GetDiaPhoton2PEFactor(referenceLambda);
        ++countPix;
      }

      if (countPix)
        avgPEfactor /= countPix;

      const auto& detTel = detFD.GetTelescope(*tIt);
      fDiaphragmArea = detTel.GetDiaphragmArea();
      fGainVariance = detTel.GetCamera().GetGainVariance();
      fPhotonToPhotoElectron = avgPEfactor;
      break;
    }

    if (iMode > 0)
      return true;

    // combine to kAve bins for estimate of dEdXmax and Xmax
    const int kAve = 10;
    int nAv = int(double(fnBins) / kAve);
    if (nAv < 1)
      nAv = 1;

    bool firstTime = true;
    double xMax = -1;
    double dEdXmax = -1;
    for (unsigned int i = 0; i < fnBins - nAv; ++i) {
      double eSum = 0;
      double wSum = 0;
      double xSum = 0;
      for (unsigned int j = i; j < i + nAv; ++j) {
        const double w = 1 / fvdEdX[j];
        eSum += fdEdX[j] * w;
        wSum += w;
        xSum += fDepth[j];
      }

      const double eM = eSum / wSum;
      const double xM = xSum / nAv;
      if (firstTime || eM > dEdXmax) {
        firstTime = false;
        dEdXmax = eM;
        xMax = xM;
      }
    }

    if (xMax < 0) {
      if (fVerbosity > -1)
        ERROR("InitializeGHFit - no initial values found! "
              "Shouldn't happen! This is only a sanity check!");
      return false;
    } else {
      fXmax = xMax;
      fdEdXmax = dEdXmax;
      fX0 = fDefX0;
      fLambda = fDefLambda;

      if (fVerbosity > -1)
        DumpCurrentParameters(eInitialize);

      return true;
    }
  }


  // recalculates contributions at aperture and dEdX

  void
  EnergyFitter::FillProfilesAtAperture(fevt::Eye& eye)
  {
    const double kXconvert = g/cm2;
    const double kEconvert = PeV/kXconvert;

    // a) light at aperture
    const unsigned int mSiz = fChFl->size1();
    ublas::vector<double> dEdX(mSiz);

    for (unsigned int i = 0; i < mSiz; ++i) {
      const double x = fChFlDepth[i] / kXconvert;
      dEdX(i) = GaisserHillas(x, fX0, fXmax, fdEdXmax, fLambda) * kEconvert;
    }

    const auto vFluo = ublas::prod(*fChFlPtr->GetFluorescenceMatrix(), dEdX);
    const auto vChDir = ublas::prod(*fChFlPtr->GetDirectCherenkovMatrix(), dEdX);
    const auto vChMScat = ublas::prod(*fChFlPtr->GetMieScatteredCherenkovMatrix(), dEdX);
    const auto vChRScat = ublas::prod(*fChFlPtr->GetRayleighScatteredCherenkovMatrix(), dEdX);

    auto& eyeRecData = eye.GetRecData();
    const auto& lightFlux = eyeRecData.GetLightFlux(fevt::FdConstants::eTotal);
    auto& dirCher = eyeRecData.GetLightFlux(fevt::FdConstants::eCherDirect);
    dirCher.Clear();
    auto& scatMCher = eyeRecData.GetLightFlux(fevt::FdConstants::eCherMieScattered);
    scatMCher.Clear();
    auto& scatRCher = eyeRecData.GetLightFlux(fevt::FdConstants::eCherRayleighScattered);
    scatRCher.Clear();
    auto& cherMS = eyeRecData.GetLightFlux(fevt::FdConstants::eCherMultScattered);
    cherMS.Clear();
    auto& fluorTot = eyeRecData.GetLightFlux(fevt::FdConstants::eFluorTotal);
    fluorTot.Clear();
    auto& fluorDir = eyeRecData.GetLightFlux(fevt::FdConstants::eFluorDirect);
    fluorDir.Clear();
    auto& fluorMS = eyeRecData.GetLightFlux(fevt::FdConstants::eFluorMultScattered);
    fluorMS.Clear();

    const auto& fluorMSFactors = fChFlPtr->GetFluorescenceMSFactor();
    const auto& cherMSFactors = fChFlPtr->GetCherenkovMSFactor();

    for (unsigned int i = mSiz - fnBins, j = 0; i < mSiz; ++i, ++j) {
      const double time = lightFlux.GetX(j);
      const double fl = vFluo(i);
      fluorTot.PushBack(time, 0, fl, 0);
      const double fluoMSfac = fluorMSFactors[i];
      const double scattFluo = fl * (1 - 1/fluoMSfac);
      fluorDir.PushBack(time, 0, fl - scattFluo, 0);
      fluorMS.PushBack(time, 0, scattFluo, 0);
      const double cd = vChDir(i);
      const double cherMSfac = cherMSFactors[i];
      const double scattDirCher = cd * (1 - 1/cherMSfac);
      dirCher.PushBack(time, 0, cd - scattDirCher, 0);
      const double cms = vChMScat(i);
      const double scattMieCher = cms * (1 - 1/cherMSfac);
      scatMCher.PushBack(time, 0, cms - scattMieCher, 0);
      const double crs = vChRScat(i);
      const double scattRayCher = crs * (1 - 1/cherMSfac);
      scatRCher.PushBack(time, 0, crs - scattRayCher, 0);
      cherMS.PushBack(time, 0, scattDirCher + scattRayCher + scattMieCher, 0);
    }

    // b) dEdX

    // add GH light to detected light vector outside FOV
    ublas::vector<double> n(mSiz);
    diagonalMatrix Vn(mSiz);
    for (unsigned int i = 0; i < mSiz - fnBins; ++i) {
      n(i) = vFluo(i) + vChDir(i) + vChMScat(i) + vChRScat(i);
      Vn(i, i) = 0;
    }

    for (unsigned int i = mSiz - fnBins, j = 0; i < mSiz; ++i, ++j) {
      n(i) = lightFlux.GetY(j);
      const double en = lightFlux.GetYErr(j);
      Vn(i, i) = Sqr(en);
    }

    auto& CF = *fChFl;
    const lowerTriangularMatrix* CFinvPtr = nullptr;
    try {
      CFinvPtr = CF.GetInverse();
    } catch (SingularMatrixException& s) {
      ERROR("SingularMatrixException for CF");
      return;
    }
    const auto& CFinv = *CFinvPtr;

    // calculate covariance matrix

    if (fVerbosity > 1)
      cerr << "\n    recalculating dEdX\n        covariance matrix ... " << endl;

    // calculate diagonal elements of  CF*V*(CF^T)^(-1)

    vector<double> covdEdX(mSiz);
    for (unsigned int i = 0; i < mSiz; ++i) {
      double sum = 0;
      for (unsigned int k = i; k < mSiz; ++k)
        sum += Sqr(CFinv(k, i)) * Vn(i, i);
      covdEdX[i] = sum;
    }
    const auto dEdX2 = ublas::prod(CFinv, n);

    auto& dedxProf = eyeRecData.GetFRecShower().GetEnergyDeposit();
    auto& longProf = eyeRecData.GetFRecShower().GetLongitudinalProfile();
    const double Ecut = eyeRecData.GetFRecShower().GetEnergyCutoff();

    for (unsigned int i = 0, j = mSiz - fnBins; i < fnBins; ++i, ++j) {
      double& eDep = dedxProf.GetY(i);
      double& eDepErr = dedxProf.GetYErr(i);
      double& Ne = longProf.GetY(i);
      double& NeErr = longProf.GetYErr(i);
      eDep = dEdX2(j);
      eDepErr = sqrt(covdEdX[j]);
      const double dEdXperElectron = utl::EnergyDeposit(fDepth[i], fXmax, Ecut);
      Ne = eDep / dEdXperElectron;
      NeErr = eDepErr / dEdXperElectron;
    }
  }


  // simple Gaisser-Hillas MINUIT fit-function
  void
  EnergyFitter::dEdXFitFunction(int& /*npar*/, double* /*gin*/, double& f,
                                double* par, int /*iflag*/)
  {
    double chisq = 0;

    const double dEdXmax= par[0];
    const double Xmax = par[1];
    const double X0 = par[2];
    const double lambda = par[3];

    if (X0 > Xmax) {
      f = fnBins * 1e2;
      return;
    }

    for (unsigned int i = 0; i < fnBins; ++i) {
      const double X = fDepth[i];
      const double dEdX = fdEdX[i];
      const double vdEdX = fvdEdX[i];
      const double dedx = GaisserHillas(X, X0, Xmax, dEdXmax, lambda);
      const double resid = dEdX - dedx;
      chisq += Sqr(resid) / vdEdX;
    }

    // add constraints
    if (fFixX0 == eConstrained)
      chisq += Sqr(X0 - fDefX0) / fVarX0;
    if (fFixLambda == eConstrained)
      chisq += Sqr(lambda - fDefLambda) / fVarLambda;

    f = chisq;
  }


  void
  EnergyFitter::LightFitFunction(int& /*npar*/, double* /*gin*/, double& f,
                                 double* par, int /*iflag*/)
  {
    const double dEdXmax = par[0];
    const double Xmax = par[1];
    const double X0 = par[2];
    const double lambda = par[3];

    if (X0 > Xmax)
      f = fnBins * 1e2;
    else
      f = (fMinimizationMethod == eMaximumLikelihood) ?
        GaisserHillasLogLike(dEdXmax, Xmax, X0, lambda) :
        GaisserHillasChi2(dEdXmax, Xmax, X0, lambda);
  }


  void
  EnergyFitter::EnergyFitFunction(int& /*npar*/, double* /*gin*/, double& f,
                                  double* par, int /*iflag*/)
  {
    const double Eem = par[0];
    const double Xmax = par[1];
    const double X0 = par[2];
    const double lambda = par[3];

    if (X0 > Xmax) {
      f = fnBins * 1e2;
      return;
    }

    const GaisserHillas4Parameter::ShapeParameter xPar(gh::eX0, X0);
    const GaisserHillas4Parameter::ShapeParameter lPar(gh::eLambda, lambda);
    const GaisserHillas4Parameter gh4(1., Xmax, xPar, lPar);

    const double dEdXmax = Eem / gh4.GetIntegral();

    f = (fMinimizationMethod == eMaximumLikelihood) ?
      GaisserHillasLogLike(dEdXmax, Xmax, X0, lambda) :
      GaisserHillasChi2(dEdXmax, Xmax, X0, lambda);
  }


  // calculates -2*log(L) of data wrt. forward folded Gaisser-Hillas function
  double
  EnergyFitter::GaisserHillasLogLike(const double dEdXmax,
                                     const double Xmax,
                                     const double X0,
                                     const double lambda)
  {
    double logLikeSum = 0;

    // conversion factors
    const double kXconvert = g/cm2;
    const double kEconvert = PeV/kXconvert;
    const double kFluxToPe = fDiaphragmArea * fPhotonToPhotoElectron;
    const double kFluxToPe2 = kFluxToPe * kFluxToPe;

    // threshold for gaussian approximation
    const double maxNpe = 30;
    // maximum n for folding of Gauss and Poisson
    const int maxSum = 100;

    // the matrix...
    const auto& m = *fChFl;
    const unsigned int mSiz = m.size1();
    const unsigned int offSet = mSiz - fnBins;

    // vectors with GH-values and light at aperture values
    vector<double> dedx(mSiz);
    vector<double> lightSum(fnBins);

    for (unsigned int i = 0; i < mSiz; ++i) {
      const double X = fChFlDepth[i] / kXconvert;
      if (X > X0)
        dedx[i] = GaisserHillas(X, X0, Xmax, dEdXmax, lambda);
      else
        dedx[i] = 0;

      if (i >= offSet) {
        double sum = 0;
        for (unsigned int k = 0; k <= i; ++k)
          sum += dedx[k] * m(i, k);
        lightSum[i - offSet] = sum*kEconvert;
      }
    }

    double nPeObsTot = 0;
    double nPeExpTot = 0;

    // sum up log-likelihoods
    for (unsigned int i = 0; i < fnBins; ++i) {

      // expected and observed photo electrons
      const double nPeBackground = fvBackGround[i] * kFluxToPe2 / (1 + fGainVariance);
      const double nPeSignalExp = lightSum[i]*kFluxToPe;
      const double nPeExp = nPeSignalExp + nPeBackground;
      const double nPeObs = fLight[i] * kFluxToPe + nPeBackground;

      nPeObsTot += nPeObs;
      nPeExpTot += nPeExp;

      double logLike;
      if (nPeObs < maxNpe) {

        const double gaussVar = nPeExp * fGainVariance;

        const double lnNpe = log(nPeExp);

        // do the Poisson/Gauss folding
        double sum = 0;
        for (int k = 0; k < maxSum; ++k) {
          const double gaussArg = -0.5*Sqr(k - nPeObs) / gaussVar;
          const double lnFaculty = LogGamma(k + 1);
          const double expArg = gaussArg - nPeExp + k*lnNpe - lnFaculty;
          const double value = exp(expArg);
          sum += value;
        }
        logLike = log(sum) - 0.5*log(gaussVar);

      } else {

        // Gaussian approximation
        const double gaussVar = nPeExp * (1 + fGainVariance);
        logLike = -0.5 * Sqr(nPeExp - nPeObs) / gaussVar - 0.5 * log(gaussVar);

      }

      logLikeSum += logLike;

    }

    // extended logLike: total number of nPe's
    const double gaussVar = nPeExpTot * (1 + fGainVariance);
    logLikeSum +=-0.5 * Sqr(nPeExpTot - nPeObsTot) / gaussVar - 0.5 * log(gaussVar);

    // add constraints
    if (fFixX0 == eConstrained)
      logLikeSum += -0.5 * log(fVarX0) - 0.5 * Sqr(X0 - fDefX0) / fVarX0;
    if (fFixLambda == eConstrained)
      logLikeSum += -0.5 * log(fVarLambda) - 0.5 * Sqr(lambda - fDefLambda) / fVarLambda;

    // New constraints

    if (fFixkUniv == eConstrained) {
      const double myEta = (Xmax - X0) / lambda;
      const double myEcal = lambda * dEdXmax * pow(kE / myEta, myEta) * TMath::Gamma(myEta + 1);
      const double fkUniv = myEcal / dEdXmax;  // convert Ecal in PeV
      const double fkUnivTrue = fDefA1 + fDefA2 * log10(myEcal * 1e15);

      logLikeSum += -0.5 * log(fVarK) - 0.5 * Sqr(fkUniv - fkUnivTrue) / fVarK;
    }

    return -2 * logLikeSum;
  }


  // calculates chi^2 of data to Gaisser-Hillas function by forward folding of dEdX profile
  double
  EnergyFitter::GaisserHillasChi2(const double dEdXmax,
                                  const double Xmax,
                                  const double X0,
                                  const double lambda)
  {
    double chisq = 0;

    const double kXconvert = g/cm2;
    const double kEconvert = PeV/kXconvert;

    const auto& m = *fChFl;
    const unsigned int mSiz = m.size1();
    const unsigned int offSet = mSiz - fnBins;

    vector<double> dedx(mSiz);

    for (unsigned int i = 0; i < mSiz; ++i) {
      const double X = fChFlDepth[i] / kXconvert;
      if (X > X0)
        dedx[i] = GaisserHillas(X, X0, Xmax, dEdXmax, lambda);
      else
        dedx[i] = 0;
    }

    for (unsigned int i = 0; i < fnBins; ++i) {

      const unsigned int j = offSet + i;

      double lightSum = 0;
      for (unsigned int k = 0; k <= j; ++k)
        lightSum += dedx[k] * m(j, k);

      lightSum *= kEconvert;

      const double resid = lightSum - fLight[i];

      double variance;
      if (fMinimizationMethod == eModifiedLeastSquares)
        variance = fvLight[i];
      else {
        const double vLight = lightSum * (1 + fGainVariance) / (fDiaphragmArea * fPhotonToPhotoElectron);
        const double bgVariance = fvBackGround[i];
        variance = vLight+bgVariance;
      }

      chisq += Sqr(resid) / variance;

    }

    // add constraints
    if (fFixX0 == eConstrained)
      chisq += Sqr(X0 - fDefX0) / fVarX0;
    if (fFixLambda == eConstrained)
      chisq += Sqr(lambda - fDefLambda) / fVarLambda;

    // New constraints
    if (fFixkUniv == eConstrained) {
      double myEta = (Xmax -X0) / lambda;
      double myEcal = lambda * dEdXmax * pow(kE / myEta, myEta) * TMath::Gamma(myEta + 1);
      double fkUniv = myEcal / dEdXmax;  // convert Ecal in PeV
      double fkUnivTrue = fDefA1 + fDefA2 * log10(myEcal * 1e15);

      chisq += Sqr(fkUniv - fkUnivTrue) / fVarK;
    }

    return chisq;
  }


  bool
  EnergyFitter::FitProfile(const EFitOption option)
  {
    // initialize MINUIT

    const int npar = eNGHParameters;  // number of parameters
    double arglist[2] = { 0 };
    int ierflag = 0;

    TMinuit minuit(npar);

    // set appropriate fit function
    switch (option) {
    case eEnergyDeposit:
      minuit.SetFCN(EnergyFitter::dEdXFitFunction);
      break;
    case eLightEnergy:
      minuit.SetFCN(EnergyFitter::EnergyFitFunction);
      break;
    default:
      minuit.SetFCN(EnergyFitter::LightFitFunction);
      break;
    }

    // update fluorescence cherenkov matrix after first fit
    if (option == eLightAllParam && !UpDateFCM())
      return false;

    fChFl = fChFlPtr->GetCherenkovFluorescenceMatrix();

    // MINUIT verbosity
    if (fVerbosity == 2)
      minuit.SetPrintLevel(0);
    else if (fVerbosity >= 3)
      minuit.SetPrintLevel(fVerbosity);
    else {
      minuit.SetPrintLevel(-1);
      minuit.mnexcm("SET NOW", arglist, 0, ierflag);
    }
    arglist[0] = 1;
    minuit.mnexcm("SET ERR", arglist, 1, ierflag);

    // set starting parameters
    if (option == eLightEnergy)
      minuit.mnparm(eEnergy,"Eem", fEem, 0.1*fabs(fEem), 0, 0, ierflag);
    else
      minuit.mnparm(eEnergy, "dEdXmax", fdEdXmax, 0.1*fabs(fdEdXmax), 0, 0, ierflag);
    minuit.mnparm(eXmax, "Xmax", fXmax, 20, 0, 0, ierflag);
    minuit.mnparm(eX0, "X0", fX0, 10, -1000, fDepth[0], ierflag);
    minuit.mnparm(eLambda, "lambda", fLambda, 5, 10, 150, ierflag);

    if (fFixLambda == eFixed || option == eEnergyDeposit)
      minuit.FixParameter(eLambda);
    if (fFixX0 == eFixed || option == eEnergyDeposit)
      minuit.FixParameter(eX0);

    // Perform minimization
    arglist[0] = 1000;
    arglist[1] = 1;
    minuit.mnexcm("MIGRAD", arglist, 2, ierflag);

    if (ierflag) {
      if (fVerbosity > -1)
        cerr << "            MIGRAD failed in dEdX fit " << ierflag << endl;
      return false;
    }

    double amin = 0;
    double edm = 0;
    double errdef = 0;
    int nvpar = 0;
    int nparx = 0;
    int icstat = 0;
    minuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);

    // store chi2, parameters at minimum and covariance

    fnDof = fnBins - 2;
    if (fFixX0 == eFree)
      fnDof -= 1;
    if (fFixLambda == eFree)
      fnDof -= 1;

    double p[eNGHParameters] = { 0 };
    if (option != eLightEnergy)
      minuit.GetParameter(eEnergy, fdEdXmax, p[eEnergy]);
    else
      minuit.GetParameter(eEnergy, fEem, p[eEnergy]);

    minuit.GetParameter(eXmax, fXmax, p[eXmax]);
    minuit.GetParameter(eX0, fX0, p[eX0]);
    minuit.GetParameter(eLambda, fLambda, p[eLambda]);
    minuit.mnemat(&fCovariance[0][0], eNGHParameters);

    if (fMinimizationMethod == eMaximumLikelihood) {
      fChisqGH = GaisserHillasChi2(fdEdXmax, fXmax, fX0, fLambda);
      //PrintLightFluxes();
    } else
      fChisqGH = amin;

    if (fVerbosity > -1 )
      DumpCurrentParameters(option);

    return true;
  }


  // calculates total energy and its error, returns true if successful
  bool
  EnergyFitter::CalculateEnergy(fevt::Eye& eye)
  {
    // create and set GH parameters
    auto& recShower = eye.GetRecData().GetFRecShower();
    if (!recShower.HasGHParameters()) {
      evt::GaisserHillas4Parameter ghTmp;
      recShower.MakeGHParameters(ghTmp);
    }

    constexpr double gcm2 = g/cm2;
    auto& gh = recShower.GetGHParameters();
    gh.SetXMax(fXmax*gcm2, sqrt(fCovariance[eXmax][eXmax])*gcm2);
    gh.SetNMax(fdEdXmax*PeV/gcm2, sqrt(fCovariance[eEnergy][eEnergy])*PeV/gcm2, true);
    gh.SetChiSquare(fChisqGH, fnDof);

    auto* gh4 = dynamic_cast<GaisserHillas4Parameter*>(&gh);
    if (gh4) {
      gh4->SetShapeParameter(gh::eX0, fX0*gcm2, sqrt(fCovariance[eX0][eX0])*gcm2);
      gh4->SetShapeParameter(gh::eLambda, fLambda*gcm2, sqrt(fCovariance[eLambda][eLambda])*gcm2);
    }

    // fit energy to obtain realistic errors

    fEem = GetEmEnergy() / PeV;
    if (!fEem) {
      if (fVerbosity > -1)
        ERROR("Error calculating gamma function");
      recShower.SetEmEnergy(0, 0);
      recShower.SetTotalEnergy(0, 0);
      return false;
    }

    if (FitProfile(eLightEnergy)) {
      // electromagnetic energy
      fEem *= PeV;
      const double EemStatVar = fCovariance[eEnergy][eEnergy] * PeV*PeV;
      recShower.SetEmEnergy(fEem, sqrt(EemStatVar), ShowerFRecData::eStatistical);
      return true;
    }

    return false;
  }


  // calculates Gaisser-Hillas integral for current parameters (a.k.a. calorimetric energy) Offline units are returned
  double
  EnergyFitter::GetEmEnergy()
    const
  {
    const GaisserHillas4Parameter::ShapeParameter xPar(gh::eX0, fX0);
    const GaisserHillas4Parameter::ShapeParameter lPar(gh::eLambda, fLambda);
    const GaisserHillas4Parameter gh4(fdEdXmax, fXmax, xPar, lPar);
    return gh4.GetIntegral()*PeV;
  }


  // updates Cherenkov-Fluorescence matrix given new fXmax value
  bool
  EnergyFitter::UpDateFCM()
  {
    const double kMaxXmaxDiff = 20;
    if (std::fabs(fChFlPtr->GetXmax()/(g/cm2) - fXmax) > kMaxXmaxDiff &&
        !fYieldRefit)
      fChFlPtr->UpdateXmax(fXmax*(g/cm2));

    return true;
  }


  void
  EnergyFitter::DumpCurrentParameters(const EFitOption option)
    const
  {
    switch (option) {
    case eInitialize:
      cout << "                  X0 [g/cm^2]  lambda [g/cm^2]  "
              "Xm [g/cm^2]  dEdXm [PeV/(g/cm^2)]  chi2/ndf    \n"
              "       ---------------------------------------------"
              "------------------------------------------\n"
              "    ==> init   :";
      break;
    case eEnergyDeposit:
      cout << "    ==> prefit :";
      break;
    case eLightAllParam:
      cout << "    ==> final  :";
      break;
    default:
      return;
    }
    cout << setw(12) << fX0
         << setw(15) << fLambda
         << setw(15) << fXmax
         << setw(17) << fdEdXmax
         << setw(12) << int(fChisqGH) << "/"
         << fnDof << endl;
  }


  // steering routine for profile refit
  bool
  EnergyFitter::ReFitProfile(const fevt::Eye& eye, CherenkovFluorescenceMatrix* const chfl)
  {
    fChFlPtr = chfl;

    if (InitializeGHFit(eye, 1) && FitProfile(EnergyFitter::eLightAllParam)) {
      CleanUp();
      return true;
    }
    CleanUp();
    return false;
  }


  void
  EnergyFitter::CleanUp()
  {
    delete[] fdEdX;
    fdEdX = nullptr;

    delete[] fvdEdX;
    fvdEdX = nullptr;

    delete[] fLight;
    fLight = nullptr;

    delete[] fvLight;
    fvLight = nullptr;

    delete[] fvBackGround;
    fvBackGround = nullptr;

    delete[] fvdEdX;
    fvdEdX = nullptr;

    delete[] fDepth;
    fDepth = nullptr;

    delete[] fChFlDepth;
    fChFlDepth = nullptr;
  }


  // returns xmax in Offline units
  double
  EnergyFitter::GetXmax()
    const
  {
    return fXmax*g/cm2;
  }


  // returns dEdXmax in Offline units
  double
  EnergyFitter::GetdEdXmax()
    const
  {
    return fdEdXmax*PeV/g*cm2;
  }


  // returns X0 in Offline units
  double
  EnergyFitter::GetX0()
    const
  {
    return fX0*g/cm2;
  }


  // returns lambda in Offline units
  double
  EnergyFitter::GetLambda()
    const
  {
    return fLambda*g/cm2;
  }

}