FdEnergyDepositFinderKG/ProfileFitter.cc

Go to the documentation of this file.

#include "ProfileFitter.h"
#include "CFMatrixCalculator.h"
#include <utl/TabulatedFunctionErrors.h>
#include <evt/GaisserHillas4Parameter.h>
#include <utl/ErrorLogger.h>
#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/Reader.h>
#include <utl/AugerException.h>
#include <fwk/CentralConfig.h>
#include <evt/GaisserHillas4Parameter.h>
#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <fdet/Camera.h>
#include <fdet/Eye.h>

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

#include <iostream>
#include <sstream>
#include <algorithm>
#include <limits>

using namespace std;
using namespace utl;
using namespace FdEnergyDepositFinderKG;

using evt::GaisserHillas4Parameter;



// Defining a function for managing a non-gaussian exponential tail for the L parameter in the USP reconstruction

inline  double 
GaussExp(const double x, const double mu, const double v, const double tau_sigma)
{
  TF1 *fGaussExp_L = new TF1("fGaussExp_L","[0] * [3]/2. * exp( [3]/2. * ( 2.*[1] + [3]*[2]*[2] - 2.*(x+[4]) ) ) * TMath::Erfc(  ( [1] + [3]*[2]*[2] - (x+[4]) ) / (sqrt(2.)*[2]) )",10.,2000.); 
  const double grammage =  g/(cm*cm);
  const double Lc_e = mu / grammage;
  const double Lc_sigma_e = sqrt(v) / grammage;
  const double Lc_tau = tau_sigma * Lc_sigma_e;
  
  fGaussExp_L->SetParameters(1., Lc_e, Lc_sigma_e , 1/Lc_tau, 0.);
  const double Lc_shift = fGaussExp_L->GetMaximumX() - mu / grammage;
  fGaussExp_L->SetParameter(4,Lc_shift);
  double ff = fGaussExp_L->Eval( x / grammage);
  
  fGaussExp_L->Delete(); 
  return ff;
}




unsigned int ProfileFitter::fNdof = 0;
vector<double> ProfileFitter::fDepth;
vector<double> ProfileFitter::fSize;
vector<double> ProfileFitter::fSizeVariance;
EFitFunctionType ProfileFitter::fFunctionType = eGH4dEdXChi2;
bool ProfileFitter::fIsConstrained = false;

std::string ProfileFitter::fGaisserHillasType = "";

bool ProfileFitter::fIsUnivConstrained = false;
utl::Function ProfileFitter::fUnivFunction;
double ProfileFitter::fVarK = 0;

bool ProfileFitter::fEMGLConstrained = false;
double ProfileFitter::fTauOverSigma = 0;


const CFMatrixCalculator* ProfileFitter::fCFMatrixData = nullptr;
const CFMatrixCalculator* ProfileFitter::fCFMatrixDataDense = nullptr;
GHShapeParameters ProfileFitter::fGHShapePars;
ProfileFitter::EFitMode ProfileFitter::fFitMode = eNMax;

double ProfileFitter::fShapeConstraintChi2[3] = { 0, 0, 0 };

int ProfileFitter::fAafCorrection = 0;
bool ProfileFitter::fUseNoiseBins = false;
vector<vector<double>> ProfileFitter::fErrFuncFactorsBuf;


ProfileFitter::ProfileFitter() :
  fVerbosity(0)
{
  for (int i = 0; i < 3; ++i)
    fShapeConstraintChi2[i] = 0;
}


void
ProfileFitter::Init()
{
  Branch topBranch =
    fwk::CentralConfig::GetInstance()->GetTopBranch("FdEnergyDepositFinder");

  fGHShapePars = GHShapeParameters(topBranch);

  fGaisserHillasType = topBranch.GetChild("gaisserHillasType").Get<string>();

  Branch kUnivBranch = topBranch.GetChild("kUnivConstrained");

  if (kUnivBranch.GetChild("constrained").Get<bool>()) {
    fIsUnivConstrained = true;
    kUnivBranch.GetChild("function").GetData(fUnivFunction);
    const double kRMS = kUnivBranch.GetChild("ksigma").Get<double>();
    fVarK = pow(kRMS, 2); // / (g/cm2), 2);
  }
 

// Introducing a non-gaussian exponential tail for the L parameter in the USP reconstruction
 
 Branch EMGLConstraintBranch = topBranch.GetChild("ExponentiallyModifiedGaussianLConstraint");

  if (EMGLConstraintBranch.GetChild("activateEMGLConstraint").Get<bool>()) {
    fEMGLConstrained = true;
    fTauOverSigma = EMGLConstraintBranch.GetChild("tauOverSigma").Get<double>();
    if (fGaisserHillasType != "eUSP")
      WARNING("The exponentially modified gaussian constraint is defined only for USP Gaisser-Hillas function (L parameter). "
              "For the other Gaisser-Hillas functions the constraint is always gaussian"); 
  }



  topBranch.GetChild("antiAliasingFilterCorrection").GetData(fAafCorrection);
  topBranch.GetChild("useNoiseBins").GetData(fUseNoiseBins);

  if (topBranch.GetChild("kUnivConstrained").GetChild("constrained").Get<bool>() &&
      fGaisserHillasType == "eUSP") {
    WARNING("k-constraint used in USP Gaisser-Hillas function - k is a "
            "function of L and R and should be already constrained"
            " through them.");
  }

  ostringstream info;
  info << "\n"
          "\t anti-aliasing filter correction is " << (fAafCorrection ? "on" : "off") << "\n"
          "\t use noise bins is " << (fUseNoiseBins ? "on" : "off") << '\n';
  for (int index = 0; index < 2; ++index) {
    const EOneTwo i = EOneTwo(index);
    info << "\t " << fGHShapePars.Name(i) << " is constrained at (E=1e19eV) "
         << fGHShapePars.Mean(i, 1e19*eV)/fGHShapePars.Unit(i) << "+/-" << fGHShapePars.Sigma(i, 1e19*eV)/fGHShapePars.Unit(i) << " "
            "[" << fGHShapePars.Min(i)/fGHShapePars.Unit(i) << ", " << fGHShapePars.Max(i)/fGHShapePars.Unit(i) << "] "
         << fGHShapePars.UnitName(i) << '\n';
  }
  INFO(info);
  for (int i = 0; i < 3; ++i)
    fShapeConstraintChi2[i] = 0;
}


void
ProfileFitter::SetLightFluxData(const CFMatrixCalculator& cfMatrixData, const CFMatrixCalculator& cfMatrixDataDense)
{
  fCFMatrixData = &cfMatrixData;
  fCFMatrixDataDense = &cfMatrixDataDense;

  //prepare anti-aliasing filter convolution parameters for telescopes in FOV
  if (fAafCorrection) {
    fErrFuncFactorsBuf.clear();
    for (auto telIter = (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataBegin() : fCFMatrixDataDense->AllTelDataBegin());
         telIter != (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataEnd() : fCFMatrixDataDense->AllTelDataEnd());
         ++telIter) {
      if (telIter->GetType() == TelescopeData::eInsideFOV) {
        const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
        const unsigned int eyeId = fCFMatrixData->GetEyeId();
        const unsigned int telId = telIter->GetId();
        const fdet::Telescope& detTel = detFD.GetEye(eyeId).GetTelescope(telId);
        PrepareTimeConvolution(detTel);
      }
    }
    //copy zeta pixel light fractions to noise data bins
    if (fUseNoiseBins && fAafCorrection == 2) {
      auto telIterNoNoise = fCFMatrixDataDense->AllTelDataBegin();
      for (auto telIter = fCFMatrixDataDense->NoiseTelDataBegin(); telIter != fCFMatrixDataDense->NoiseTelDataEnd(); ++telIter) {
        if (telIter->GetType() == TelescopeData::eInsideFOV) {
          while (telIterNoNoise->GetType() == TelescopeData::eOutsideFOV) {
            ++telIterNoNoise;
          }
          vector<TelescopeDataBin>::const_iterator binIterNoNoise = telIterNoNoise->TelDataBinsBegin();
          for (auto binIter = telIter->TelDataBinsBegin(); binIter !=  telIter->TelDataBinsEnd(); ++binIter)  {
            if (binIter->fMinDepth != 0 || binIter->fMaxDepth != 0) { //not a noise bin
              auto& zp = binIter->GetZetaPixels();
              for (unsigned int z = 0, n = zp.size(); z < n; ++z) {
                zp[z].fLightFraction = binIterNoNoise->GetZetaPixels()[z].fLightFraction;
              }
              ++binIterNoNoise;
            }
          }
          ++telIterNoNoise;
        }
      }
    }
  }
}


void
ProfileFitter::SetProfileData(const TabulatedFunctionErrors& dEdXProfile)
{
  fDepth.clear();
  fSize.clear();
  fSizeVariance.clear();

  for (unsigned int i = 0, n = dEdXProfile.GetNPoints(); i < n; ++i) {
    fDepth.push_back(dEdXProfile.GetX(i));
    fSize.push_back(dEdXProfile.GetY(i));
    fSizeVariance.push_back(pow(dEdXProfile.GetYErr(i), 2));
  }
}


void
ProfileFitter::SetShapePars(const GHShapeParameters &pars)
{
  fGHShapePars = pars;
}


void
ProfileFitter::SetUnivConstrained(const bool constrained,
                                  const utl::Function& func,
                                  const double ksigma)
{
  fIsUnivConstrained = constrained;
  fUnivFunction = func;
  fVarK = pow(ksigma, 2);
}


bool
ProfileFitter::Fit()
{
  if (fSizeVariance.empty()) {
    ERROR("no data to fit!");
    return false;
  }

  fIsConstrained = true;

  switch (fFunctionType) {
  case eUndefined:
    return false;
  case eGH2dEdXChi2:
  case eGH4dEdXChi2:
  case eGH4LightLogLike:
    return GHFit();
  }

  return false;
}


bool
ProfileFitter::GHFit()
{
  for (int i = 0; i < 3; ++i)
    fShapeConstraintChi2[i] = 0;

  const unsigned int nPar = eNGHpars;
  TMinuit theMinuit(nPar);
  theMinuit.SetFCN(ProfileFitter::GHFitFunction);
  InitMinuit(theMinuit);

  if (fStartParameters.size() < eNGHpars) {
    ERROR("too few start parameters for fit");
    return false;
  }

  int ierflag = 0;
  theMinuit.mnparm(eGHnMax, "nMax", fStartParameters[eGHnMax],
                   fabs(fStartParameters[eGHnMax])*5e-2, 0, 0, ierflag);
  theMinuit.mnparm(eGHXmax, "xMax", fStartParameters[eGHXmax],
                   10*g/cm2, 0, 0, ierflag);
  theMinuit.mnparm(eGHShape1, fGHShapePars.Name(eOne),
                   fStartParameters[eGHShape1], fGHShapePars.Step(eOne),
                   fGHShapePars.Min(eOne), fGHShapePars.Max(eOne), ierflag);
  theMinuit.mnparm(eGHShape2, fGHShapePars.Name(eTwo),
                   fStartParameters[eGHShape2], fGHShapePars.Step(eTwo),
                   fGHShapePars.Min(eTwo), fGHShapePars.Max(eTwo), ierflag);
  if (fFunctionType == eGH2dEdXChi2) {
    theMinuit.FixParameter(eGHShape1);
    theMinuit.FixParameter(eGHShape2);
  }

  bool retVal = Minimize(theMinuit);

  if (fFunctionType == eGH2dEdXChi2)
    fNdof -= 2;

  return retVal;
}


void
ProfileFitter::GHFitFunction(int& /*npar*/, double* const /*gin*/,
                             double& value, double* const par,
                             const int /*iFlag*/)
{
  double fitParameters[eNGHpars];
  for (int i = 0; i < eNGHpars; ++i)
    fitParameters[i] = par[i];

  // nMax from integral
  if (fFitMode == eIntegral) {
    const GaisserHillas4Parameter::ShapeParameter xPar(fGHShapePars.Type(eOne), par[eGHShape1]);
    const GaisserHillas4Parameter::ShapeParameter lPar(fGHShapePars.Type(eTwo), par[eGHShape2]);
    const GaisserHillas4Parameter gh4(1, par[eGHXmax], xPar, lPar, fGHShapePars.Type());
    fitParameters[eGHnMax] = par[eGHIntegral] / gh4.GetIntegral();
  }

  if (fAafCorrection ||
      fUseNoiseBins ||
      fCFMatrixDataDense->GetCFMatrix().GetSize() != fCFMatrixData->GetCFMatrix().GetSize())
    value = fFunctionType == eGH4LightLogLike ?
      GaisserHillasLogLikeConvoluted(fitParameters) : GaisserHillasChi2(fitParameters);
  else
    value = fFunctionType == eGH4LightLogLike ?
      GaisserHillasLogLike(fitParameters) : GaisserHillasChi2(fitParameters);
}


// calculates -2*log(L) of data wrt.
// forward folded Gaisser-Hillas function
double
ProfileFitter::GaisserHillasLogLike(const double* const par)
{
  double logLikeSum = 0;
  fNdof = 0;

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

  const GaisserHillas4Parameter::ShapeParameter xPar(fGHShapePars.Type(eOne), par[eGHShape1]);
  const GaisserHillas4Parameter::ShapeParameter lPar(fGHShapePars.Type(eTwo), par[eGHShape2]);
  const GaisserHillas4Parameter ghFunction(par[eGHnMax], par[eGHXmax], xPar, lPar, fGHShapePars.Type());

  if (std::isnan(par[eGHXmax]) || std::isnan(par[eGHnMax]) ||
      std::isnan(par[eGHShape1]) ||  std::isnan(par[eGHShape2]))
    return std::numeric_limits<double>::max();

  LowerTriangularMatrix cfMatrix = fCFMatrixData->GetCFMatrix();
  const unsigned int matrixSize = cfMatrix.GetSize();

  // vectors with GH-values and light at aperture values
  ColumnVector energyDeposit;
  energyDeposit.SetSize(matrixSize);

  int i = 0;
  for (auto telIter = fCFMatrixData->AllTelDataBegin(); telIter != fCFMatrixData->AllTelDataEnd(); ++telIter) {
    for (auto binIter = telIter->TelDataBinsBegin(); binIter !=  telIter->TelDataBinsEnd(); ++binIter) {
      energyDeposit(i) = ghFunction.Eval(binIter->GetMeanDepth());
      ++i;
    }
  }

  ColumnVector predictedLightFlux = cfMatrix*energyDeposit;

  double nPeObsTot = 0;
  double nPeExpTot = 0;

  double gainVariance = 0;

  // sum up log-likelihoods
  i = 0;
  for (auto telIter = fCFMatrixData->AllTelDataBegin(); telIter != fCFMatrixData->AllTelDataEnd(); ++telIter) {

    if (telIter->GetType() == TelescopeData::eInsideFOV) {

      // conversion factors for this telescope
      const double kFluxToPe = telIter->GetDiaphragmArea() * telIter->GetPhotonToPhotoElectron();
      const double kFluxToPe2 = kFluxToPe * kFluxToPe;
      gainVariance = telIter->GetGainVariance();
      
      for (auto binIter = telIter->TelDataBinsBegin(); binIter != telIter->TelDataBinsEnd(); ++binIter) {
        
        const double lightFlux = binIter->fSignal;
        const double bgVariance = binIter->fBackgroundVariance;

        // expected and observed photo electrons
        const double nPeBackground = bgVariance*kFluxToPe2 / (1 + gainVariance);
        const double nPeSignalExp = predictedLightFlux(i) * kFluxToPe;
        const double nPeExp = nPeSignalExp + nPeBackground;
        const double nPeObs = lightFlux*kFluxToPe + nPeBackground;

        nPeObsTot += nPeObs;
        nPeExpTot += nPeExp;

        double logLike;
        if (nPeObs < maxNpe) {

          const double gaussVar = nPeExp*gainVariance;

          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*pow(k - nPeObs, 2) / 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 + gainVariance);
          logLike = -0.5*pow(nPeExp - nPeObs, 2) / gaussVar - 0.5*log(gaussVar);
        }

        logLikeSum += logLike;
        ++fNdof;
        ++i;
      }
    } else
      i += telIter->GetTelescopeDataBins().size();

  }

  // extended logLike: total number of nPe's
  const double gaussVar = nPeExpTot*(1 + gainVariance);
  logLikeSum += -0.5*pow(nPeExpTot - nPeObsTot, 2) / gaussVar - 0.5*log(gaussVar);
  --fNdof;

  // add constraints
  if (fIsConstrained || fIsUnivConstrained) {

    const double Xmax = par[eGHXmax];
    const double dEdXmax = par[eGHnMax];
    const double X0 = ghFunction.GetShapeParameter(evt::gh::eX0);
    const double lambda = ghFunction.GetShapeParameter(evt::gh::eLambda);
    const double myEta = (Xmax -X0)/lambda;
    const double myEcal = lambda * dEdXmax * pow(kE/myEta,myEta) * TMath::Gamma(myEta+1);

    if (fIsConstrained) {
      const double p[2] = {par[eGHShape1], par[eGHShape2]};
      for (int index = 0; index < 2; ++index) {
        
        
        const EOneTwo i = EOneTwo(index);
        const double v = fGHShapePars.Variance(i, myEcal);
        const double mu = fGHShapePars.Mean(i, myEcal);

        if (fEMGLConstrained==true && fGHShapePars.Name(i)=="L_{USP}"){ 
          
          double ff = GaussExp(p[i], mu, v , fTauOverSigma);
          
          logLikeSum += (log(ff));
          fShapeConstraintChi2[index] = -2.*(log(ff));
          
        } else {
          
          logLikeSum += (-0.5*log(v) - 0.5*pow(p[i]-mu, 2)/v);
          fShapeConstraintChi2[index] = -2.*(-0.5*log(v) - 0.5*pow(p[i]-mu, 2)/v);
          
        }
      }
    }

    if (fIsUnivConstrained) {
      const double fkUniv = myEcal/dEdXmax;
      const double fkUnivTrue = fUnivFunction(myEcal);
      logLikeSum +=(-0.5*log(fVarK) - 0.5*pow(fkUniv - fkUnivTrue, 2)/fVarK);
      fShapeConstraintChi2[2] = -2.*(-0.5*log(fVarK)-0.5*pow(fkUniv-fkUnivTrue ,2)/fVarK);
    }
  }

  return -2 * logLikeSum;
}


// version with forward folded anti-aliasing filter convolution
double
ProfileFitter::GaisserHillasLogLikeConvoluted(const double* const par)
{
  double logLikeSum = 0;
  fNdof = 0;

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

  const GaisserHillas4Parameter::ShapeParameter xPar(fGHShapePars.Type(eOne), par[eGHShape1]);
  const GaisserHillas4Parameter::ShapeParameter lPar(fGHShapePars.Type(eTwo), par[eGHShape2]);
  const GaisserHillas4Parameter ghFunction(par[eGHnMax], par[eGHXmax], xPar, lPar, fGHShapePars.Type());

  if (std::isnan(par[eGHXmax]) || std::isnan(par[eGHnMax]) ||
      std::isnan(par[eGHShape1]) || std::isnan(par[eGHShape2]))
    return std::numeric_limits<double>::max();

  LowerTriangularMatrix cfMatrix = fCFMatrixDataDense->GetCFMatrix();
  const unsigned int matrixSize = cfMatrix.GetSize();

  // vectors with GH-values and light at aperture values
  ColumnVector energyDeposit;
  energyDeposit.SetSize(matrixSize);

  int i = 0;
  for (auto telIter = fCFMatrixDataDense->AllTelDataBegin(); telIter != fCFMatrixDataDense->AllTelDataEnd(); ++telIter) {
    for (auto binIter = telIter->TelDataBinsBegin(); binIter !=  telIter->TelDataBinsEnd(); ++binIter) {
      energyDeposit(i) = ghFunction.Eval(binIter->GetMeanDepth());
      ++i;
    }
  }

  ColumnVector predictedLightFlux = cfMatrix*energyDeposit;

  vector<double> predLF;
  vector<int> telBinSizes;
  vector<int> telBinSizesPredictedLightFlux;
  // include noise bins and sum dense bins
  i = 0;
  for (auto telIter = (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataBegin() : fCFMatrixDataDense->AllTelDataBegin());
       telIter != (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataEnd() : fCFMatrixDataDense->AllTelDataEnd());
       ++telIter) {

    int lastParentBin = -1;
    int telBinSize = 0;
    int k = 0;
    for (auto binIter = telIter->TelDataBinsBegin(); binIter != telIter->TelDataBinsEnd(); ++binIter) {
      if (telIter->GetType() == TelescopeData::eInsideFOV) {
        if (binIter->fParentBin == lastParentBin) {
          if (binIter->fMinDepth != 0 || binIter->fMaxDepth != 0) {
            predLF.back() += predictedLightFlux(i);
            ++i;
            ++k;
          } //nothing for noise bins
        } else { //new real bin in aperture light
          if (binIter->fMinDepth != 0 || binIter->fMaxDepth != 0) {
            predLF.push_back(predictedLightFlux(i));
            ++i;
            ++k;
          } else { //noise bin
            predLF.push_back(0);
          }
          lastParentBin = binIter->fParentBin;
          ++telBinSize;
        }
      } else { //outside FOV bins are "not dense"
        predLF.push_back(predictedLightFlux(i));
        ++telBinSize;
        ++i;
        ++k;
      }
    }
    telBinSizes.push_back(telBinSize);
    telBinSizesPredictedLightFlux.push_back(k);
  }

  // prepare expected photo electrons

  double gainVariance = 0;

  i = 0;
  // vector<double> nPeExpArray(predictedLightFlux.GetSize(),0);
  vector<double> nPeExpArray(predLF.size(), 0);
  for (auto telIter = (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataBegin() : fCFMatrixDataDense->AllTelDataBegin());
       telIter != (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataEnd() : fCFMatrixDataDense->AllTelDataEnd());
       ++telIter) {

    if (telIter->GetType() == TelescopeData::eInsideFOV) {
      const double kFluxToPe = telIter->GetDiaphragmArea() * telIter->GetPhotonToPhotoElectron();
      //const double kFluxToPe2 = kFluxToPe * kFluxToPe;
      gainVariance = telIter->GetGainVariance();

      int lastParentBin = -1;
      for (auto binIter = telIter->TelDataBinsBegin(); binIter !=  telIter->TelDataBinsEnd(); ++binIter) {

        if (binIter->fParentBin != lastParentBin) { //sum only for real aperture light bins
          // const double bgVariance = binIter->fBackgroundVariance;

          // expected and observed photo electrons
          //const double nPeBackground = bgVariance*kFluxToPe2/(1.+gainVariance);
          //const double nPeSignalExp = predictedLightFlux(i)*kFluxToPe;
          const double nPeSignalExp = predLF[i]*kFluxToPe;
          const double nPeExp = nPeSignalExp; // +nPeBackground;
          nPeExpArray[i] = nPeExp;
          ++i;
          lastParentBin = binIter->fParentBin;
        }

      }
    } else
      i += telIter->GetTelescopeDataBins().size(); // for outside of FOV the same as in telBinSizes
  }

  // do anti-aliasing filter convolution
  i = 0;
  int k = 0;
  int t = 0;
  if (fAafCorrection) {
    unsigned int telIndex = 0;
    for (auto telIter = (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataBegin() : fCFMatrixDataDense->AllTelDataBegin());
         telIter != (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataEnd() : fCFMatrixDataDense->AllTelDataEnd());
         ++telIter) {

      unsigned int telBinSize = telBinSizes[t];

      if (telIter->GetType() == TelescopeData::eInsideFOV) {
        //const double kFluxToPe = telIter->GetDiaphragmArea() * telIter->GetPhotonToPhotoElectron();
        //const double kFluxToPe2 = kFluxToPe * kFluxToPe;
        gainVariance = telIter->GetGainVariance();

        const unsigned int nSafety = fErrFuncFactorsBuf[telIndex].size() - 1;

        const unsigned int convTraceSize = telBinSize;

        const unsigned int origTraceSize = convTraceSize + 2*nSafety;
        vector<double> originalTrace(origTraceSize, 0.);
        for (unsigned int j = 0; j < convTraceSize; j++) {
          originalTrace[j+nSafety] = nPeExpArray[i+j];
        }

        vector<double> convolutedTrace(convTraceSize, 0.);

        DoTimeConvolution(originalTrace, convolutedTrace, nSafety, telIndex, *telIter, predictedLightFlux, predLF, k, i);

        for (unsigned int j = 0; j < convTraceSize; ++j) {
          nPeExpArray[i+j] = convolutedTrace[j];
        }
        ++telIndex;
      }
      i += telBinSize;
      k += telBinSizesPredictedLightFlux[t];
      ++t;
    }
  }

  double nPeObsTot = 0;
  double nPeExpTot = 0;

  // sum up log-likelihoods
  i = 0;
  for (auto telIter = (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataBegin() : fCFMatrixDataDense->AllTelDataBegin());
       telIter != (fUseNoiseBins ? fCFMatrixDataDense->NoiseTelDataEnd() : fCFMatrixDataDense->AllTelDataEnd());
       ++telIter) {

    if (telIter->GetType() == TelescopeData::eInsideFOV) {

      // conversion factors for this telescope
      const double kFluxToPe = telIter->GetDiaphragmArea() * telIter->GetPhotonToPhotoElectron();
      const double kFluxToPe2 = kFluxToPe * kFluxToPe;
      gainVariance = telIter->GetGainVariance();

      int lastParentBin = -1;

      for (auto binIter = telIter->TelDataBinsBegin(); binIter != telIter->TelDataBinsEnd(); ++binIter) {

        if (binIter->fParentBin != lastParentBin) { //sum only for real aperture light bins

          const double lightFlux = binIter->fSignal;
          const double bgVariance = binIter->fBackgroundVariance;

          // expected and observed photo electrons
          const double nPeBackground = bgVariance * kFluxToPe2 / (1 + gainVariance);
          const double nPeExp = nPeExpArray[i] + nPeBackground;
          const double nPeObs = lightFlux * kFluxToPe + nPeBackground;

          nPeObsTot += nPeObs;
          nPeExpTot += nPeExp;

          double logLike;
          if (nPeObs < maxNpe) {

            const double gaussVar = nPeExp*gainVariance;

            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 * pow(k - nPeObs, 2) / 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 + gainVariance);
            logLike = -0.5*pow(nPeExp - nPeObs, 2) / gaussVar - 0.5*log(gaussVar);
          }

          logLikeSum += logLike;
          ++fNdof;
          ++i;
          lastParentBin = binIter->fParentBin;
        }

      }

    } else
      i += telIter->GetTelescopeDataBins().size(); // for outside of FOV the same as in telBinSizes

  }

  // extended logLike: total number of nPe's
  const double gaussVar = nPeExpTot * (1 + gainVariance);
  logLikeSum += -0.5*pow(nPeExpTot - nPeObsTot, 2) / gaussVar - 0.5*log(gaussVar);
  --fNdof;

  // add constraints
  if (fIsConstrained || fIsUnivConstrained) {

    const double Xmax = par[eGHXmax];
    const double dEdXmax = par[eGHnMax];
    const double X0 = ghFunction.GetShapeParameter(evt::gh::eX0);
    const double lambda = ghFunction.GetShapeParameter(evt::gh::eLambda);
    const double myEta = (Xmax -X0)/lambda;
    const double myEcal = lambda * dEdXmax * pow(kE/myEta,myEta) * TMath::Gamma(myEta+1);

    if (fIsConstrained) {
      const double p[2] = {par[eGHShape1], par[eGHShape2]};
      for (int index = 0; index < 2; ++index) {
        const EOneTwo i = EOneTwo(index);
        const double v = fGHShapePars.Variance(i, myEcal);
        const double mu = fGHShapePars.Mean(i, myEcal);
        

       if (fEMGLConstrained==true && fGHShapePars.Name(i)=="L_{USP}"){ 

          double ff = GaussExp(p[i], mu, v , fTauOverSigma);
          
          logLikeSum += (log(ff));
          fShapeConstraintChi2[index] = -2.*(log(ff));
                  
        } else {
          
          logLikeSum += (-0.5*log(v) - 0.5*pow(p[i]-mu, 2)/v);
          fShapeConstraintChi2[index] = -2.*(-0.5*log(v) - 0.5*pow(p[i]-mu, 2)/v);
          
        }
      }
    }

    if (fIsUnivConstrained) {
      const double fkUniv = myEcal/dEdXmax;
      const double fkUnivTrue = fUnivFunction(myEcal);
      logLikeSum += -0.5*log(fVarK) - 0.5*pow(fkUniv - fkUnivTrue, 2) / fVarK;
      fShapeConstraintChi2[2] = log(fVarK) + pow(fkUniv - fkUnivTrue, 2) / fVarK;
    }
  }

  return -2*logLikeSum;
}


// calculates chi2 of dE/dX profile wrt.
// Gaisser-Hillas function
double
ProfileFitter::GaisserHillasChi2(const double* par)
{
  double chi2 = 0;
  const GaisserHillas4Parameter::ShapeParameter xPar(fGHShapePars.Type(eOne), par[eGHShape1]);
  const GaisserHillas4Parameter::ShapeParameter lPar(fGHShapePars.Type(eTwo), par[eGHShape2]);
  const GaisserHillas4Parameter ghFunction(par[eGHnMax], par[eGHXmax], xPar, lPar, fGHShapePars.Type());
  fNdof = 0;
  if (std::isnan(par[eGHXmax]) || std::isnan(par[eGHnMax]) ||
      std::isnan(par[eGHShape1]) ||  std::isnan(par[eGHShape2]))
    return std::numeric_limits<double>::max();

  for (unsigned int i = 0; i < fDepth.size(); ++i) {

    const double y = fSize[i];

    const double x = fDepth[i];
    const double fx = ghFunction.Eval(x);
    const double vy = fSizeVariance[i];

    chi2 += pow(y - fx, 2) / vy;
    ++fNdof;
  }

  // constraints a la FdProfileReconstructor
  if (fIsConstrained || fIsUnivConstrained) {

    const double Xmax = par[eGHXmax];
    const double dEdXmax = par[eGHnMax];
    const double X0 = ghFunction.GetShapeParameter(evt::gh::eX0);
    const double lambda = ghFunction.GetShapeParameter(evt::gh::eLambda);
    const double myEta = (Xmax - X0) / lambda;
    const double myEcal = lambda * dEdXmax * pow(kE / myEta, myEta) * TMath::Gamma(myEta+1);

    if (fIsConstrained) {

      const double p[2] = { par[eGHShape1], par[eGHShape2] };
      for (int index = 0; index < 2; ++index) {
        const EOneTwo i = EOneTwo(index);
        const double v = fGHShapePars.Variance(i, myEcal);
        const double mu = fGHShapePars.Mean(i, myEcal);
        
        if (fEMGLConstrained==true && fGHShapePars.Name(i)=="L_{USP}"){ 
          
          double ff = GaussExp(p[i], mu, v , fTauOverSigma);
          
          chi2 += -2.*log(ff);
          fShapeConstraintChi2[index] = -2.*log(ff);
                  
        } else {
          
          const double addchi2 = pow(p[i] - mu, 2) / v;
          chi2 += addchi2;
          fShapeConstraintChi2[index] = addchi2;

        }
      }
    }

    if (fIsUnivConstrained) {
      const double fkUniv = myEcal/dEdXmax;
      const double fkUnivTrue = fUnivFunction(myEcal);
      const double addchi2 = pow(fkUniv - fkUnivTrue, 2) / fVarK;
      chi2 += addchi2;
      fShapeConstraintChi2[2] = addchi2;
    }
  }

  return chi2;
}


bool
ProfileFitter::Minimize(TMinuit& theMinuit)
{
  int ierflag = 0;
  double arglist[2] = { 1000, 1 };
  theMinuit.mnexcm("MINIMIZE", arglist, 2, ierflag);

  if (ierflag) {
    if (fVerbosity > -1)
      cerr << " MINIMIZE failed " << ierflag << endl;
    return false;
  }

  theMinuit.mnexcm("HESSE", arglist, 2, ierflag);

  double amin, edm, errdef;
  int nvpar, nparx, icstat;
  theMinuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
  fStatus = icstat;

  if (fStatus != 3) {
    // from TMinuit.cxx:
    //           matrix:  0= not calculated at all
    //                    1= approximation only, not accurate
    //                    2= full matrix, but forced positive-definite
    //                    3= full accurate covariance matrix
    if (fVerbosity > -1)
      cerr << " ProfileFitter::Minimize(): status from mnstat "
           << fStatus << endl;
    return false;
  }

  if (nparx != eNGHpars) {
    cerr << " ProfileFitter::Minimize(): parameter mismatch??? " << endl;
    return false;
  }

  fFitParameters.resize(nparx);
  double dummy = 0;
  for (int i = 0; i < nparx; ++i)
    theMinuit.GetParameter(i, fFitParameters[i], dummy);
  theMinuit.mnemat(&fCovariance[0][0], eNGHpars);

  if (fFunctionType != eGH4LightLogLike) {
    fChi2 = amin;

    //VN: to be able to use fFitModes and fill fShapeConstrainChi2
    int dummy1 = 0;
    double dummy2 = 0;
    double dummy3 = 0;
    int dummy4 = 0;
    GHFitFunction(dummy1, &dummy2, dummy3, &fFitParameters.front(), dummy4);
    //GaisserHillasChi2(&fFitParameters.front()); // in order to fill fShapeConstrainChi2 !

    fLHChi2 = amin;
  } else {
    //VN: to be able to use fFitModes and fill fShapeConstrainChi2
    int dummy1 = 0;
    double dummy2 = 0;
    double dummy3 = 0;
    int dummy4 = 0;
    // simple chi2 for logLike
    fFunctionType = eGH4dEdXChi2;
    GHFitFunction(dummy1, &dummy2, fChi2, &fFitParameters.front(), dummy4);
    //fChi2 = GaisserHillasChi2(&fFitParameters.front());

    //VN: light fit for fShapeConstrainChi2
    fFunctionType = eGH4LightLogLike;
    GHFitFunction(dummy1, &dummy2, dummy3, &fFitParameters.front(), dummy4);
    fLHChi2 = amin;
  }

  return true;
}


void
ProfileFitter::InitMinuit(TMinuit& theMinuit)
  const
{
  int ierflag;
  if (fVerbosity == 2)
    theMinuit.SetPrintLevel(0);
  else if (fVerbosity >= 3)
    theMinuit.SetPrintLevel(fVerbosity);
  else {
    theMinuit.SetPrintLevel(-1);
    theMinuit.mnexcm("SET NOW", 0 ,0,ierflag);
  }

  double arglist = 1;
  theMinuit.mnexcm("SET ERR", &arglist, 1, ierflag);
  arglist = 2;
  theMinuit.mnexcm("SET STRAT", &arglist, 1, ierflag);
}


void
ProfileFitter::SetStartParameters(const evt::GaisserHillas4Parameter& gh4)
{
  fStartParameters.resize(eNGHpars);
  if (fFitMode == eIntegral)
    fStartParameters[eGHnMax] = gh4.GetIntegral();
  else
    fStartParameters[eGHnMax] = gh4.GetNMax();
  fStartParameters[eGHXmax] = gh4.GetXMax();
  fStartParameters[eGHShape1] = gh4.GetShapeParameter(fGHShapePars.Type(eOne));
  fStartParameters[eGHShape2] = gh4.GetShapeParameter(fGHShapePars.Type(eTwo));
}


void
ProfileFitter::FillGHParameters(GaisserHillas4Parameter& ghPars)
  const
{
  if (ghPars.GetFunctionType() != fGHShapePars.Type()) {
    using namespace evt::gh;
    ostringstream errMsg;
    throw DoesNotComputeException("GH type mismatch: " +
                                  GetFunctionTypeName(ghPars.GetFunctionType()) +
                                  "!=" +
                                  GetFunctionTypeName(fGHShapePars.Type()));
  } else if (fFitParameters.size() != eNGHpars) {
    ERROR("FillGHParameters() without data!");
    ghPars.SetXMax(0, 0);
    ghPars.SetNMax(0, 0, true);
    ghPars.SetShapeParameter(fGHShapePars.Type(eOne), 0, 0);
    ghPars.SetShapeParameter(fGHShapePars.Type(eTwo), 0, 0);
    ghPars.SetChiSquare(0, 0);
  } else {
    ghPars.SetXMax(fFitParameters[eGHXmax], sqrt(fCovariance[eGHXmax][eGHXmax]));
    ghPars.SetNMax(fFitParameters[eGHnMax], sqrt(fCovariance[eGHnMax][eGHnMax]), true);
    ghPars.SetShapeParameter(fGHShapePars.Type(eOne), fFitParameters[eGHShape1], sqrt(fCovariance[eGHShape1][eGHShape1]));
    ghPars.SetShapeParameter(fGHShapePars.Type(eTwo), fFitParameters[eGHShape2], sqrt(fCovariance[eGHShape2][eGHShape2]));
    ghPars.SetChiSquare(fChi2, fNdof);

    ghPars.SetNMaxXMaxCorrelation(CorrCoeff(eGHnMax, eGHXmax));
    ghPars.SetCorrelationNMaxShapeParameter(fGHShapePars.Type(eTwo), CorrCoeff(eGHnMax, eGHShape2));
    ghPars.SetCorrelationNMaxShapeParameter(fGHShapePars.Type(eOne), CorrCoeff(eGHnMax, eGHShape1));
    ghPars.SetCorrelationXMaxShapeParameter(fGHShapePars.Type(eTwo), CorrCoeff(eGHXmax, eGHShape2));
    ghPars.SetCorrelationXMaxShapeParameter(fGHShapePars.Type(eOne), CorrCoeff(eGHXmax, eGHShape1));
    ghPars.SetCorrelationShapeParameters(CorrCoeff(eGHShape2, eGHShape1));
  }
}


void
ProfileFitter::GetEnergy(double& energy, double& energyError)
  const
{
  if (fFitParameters.size() != eNGHpars) {
    ERROR(" cannot calculate energy without fit!");
    energy = 0;
    energyError = 0;
    return;
  }

  if (fFitMode == eIntegral) {
    energy = fFitParameters[eGHIntegral];
    energyError = sqrt(fCovariance[eGHIntegral][eGHIntegral]);
  } else {
    GaisserHillas4Parameter ghPars(fGHShapePars.Type());
    FillGHParameters(ghPars);
    energy = ghPars.GetIntegral();
    energyError = energy * ghPars.GetIntegralError();
  }
}


void
ProfileFitter::PrepareTimeConvolution(const fdet::Telescope& detTel)
{
  const fdet::Camera& detCamera = detTel.GetCamera();
  const double timeBinSize = detCamera.GetFADCBinSize();

  // see IEEE Trans.Nucl.Sci.50:1208-1213,2003
  const double cutoffFrequency = detCamera.GetCutoffFrequency();
  const double tau = sqrt(log(2.)) / (kTwoPi * cutoffFrequency);

  vector<double> fErrFuncFactors;
  fErrFuncFactors.clear();

  const double epsilon = 1e-6;
  const int maxIter = 50;
  int iter = 0;

  double sumE = 0;
  double sumESquare = 0;

  while (sumE < 1 - epsilon && iter < maxIter) {

    const double t1 = (iter - 0.5) * timeBinSize;
    const double t2 = (iter + 0.5) * timeBinSize;

    const double E1 = 0.5 * erf(t1 / (sqrt(2.)*tau));
    const double E2 = 0.5 * erf(t2 / (sqrt(2.)*tau));
    const double E = E2 - E1;

    fErrFuncFactors.push_back(E);

    if (iter == 0) {
      sumESquare += E*E;
      sumE += E;
    } else {
      sumESquare += 2*E*E;
      sumE += 2*E;
    }

    ++iter;
  }

  if (iter+1 >= maxIter || fErrFuncFactors.empty()) {
    ostringstream msg;
    msg << " Error filling fErrFuncFactors !!! "
        << iter << ' ' << fErrFuncFactors.size();
    ERROR(msg);
    throw AugerException(msg.str());
  }

  fErrFuncFactorsBuf.push_back(fErrFuncFactors);
}


void
ProfileFitter::DoTimeConvolution(const vector<double>& originalTrace,
                                 vector<double>& convolutedTrace,
                                 unsigned int timeOffset,
                                 unsigned int ErrFuncFactorsIndex,
                                 const TelescopeData& telData,
                                 const ColumnVector& predictedLightFlux,
                                 const std::vector<double>& predLF,
                                 const int startBinDense,
                                 const int startBin)
{
  const int convolutedTraceSize = convolutedTrace.size();
  const int originialTraceSize  = originalTrace.size();
  vector<double> fErrFuncFactors;

  const std::vector<TelescopeDataBin>& telDataBins = telData.GetTelescopeDataBins();

  if (ErrFuncFactorsIndex < fErrFuncFactorsBuf.size()) {
    fErrFuncFactors = fErrFuncFactorsBuf[ErrFuncFactorsIndex];
  } else {
    ostringstream msg;
    msg << " ErrFuncFactorIndex out of range ";
    ERROR(msg);
    throw AugerException(msg.str());
  }
  const int nConv = int(fErrFuncFactors.size());

  // check trace sizes...
  if (originialTraceSize-convolutedTraceSize < 2*(nConv - 1)) {
    // that should never happend !!!!
    ostringstream msg;
    msg << " Trace size mismatch!!! "
        << convolutedTraceSize << ' ' << originialTraceSize << ' ' << fErrFuncFactors.size();
    ERROR(msg);
    throw AugerException(msg.str());
  }

  //prepare zeta pixel signals
  vector<vector<pair<int, double>>> telDataBinZeta;
  int lastParentBin = -1;
  int i = 0;
  if (fAafCorrection == 2) {
    for (unsigned int t = 0; t < telDataBins.size(); ++t) {
      if (telDataBins[t].fParentBin == lastParentBin) {
        if (telDataBins[t].fMinDepth != 0 || telDataBins[t].fMaxDepth != 0) {
          bool noZeta = bool(telDataBinZeta.back().size());
          for (unsigned int z = 0; z < telDataBins[t].GetZetaPixels().size(); ++z) {
            if (noZeta) {
              telDataBinZeta.back()[z].second += telDataBins[t].GetZetaPixels()[z].fLightFraction * predictedLightFlux(i+startBinDense);
            } else {
              telDataBinZeta.back().push_back(make_pair(telDataBins[t].GetZetaPixels()[z].fPixelId, telDataBins[t].GetZetaPixels()[z].fLightFraction * predictedLightFlux(i+startBinDense)));
            }
          }
          ++i;
        } //nothing for noise bins
      } else { //new real bin in aperture light
        vector<pair<int, double>> zetaLF;
        if (telDataBins[t].fMinDepth != 0 || telDataBins[t].fMaxDepth != 0) {  //not a noise bin
          for (unsigned int z = 0; z < telDataBins[t].GetZetaPixels().size(); ++z) {
            zetaLF.push_back(make_pair(telDataBins[t].GetZetaPixels()[z].fPixelId, telDataBins[t].GetZetaPixels()[z].fLightFraction * predictedLightFlux(i+startBinDense)));
          }
          ++i;
        }
        telDataBinZeta.push_back(zetaLF);
        lastParentBin = telDataBins[t].fParentBin;
      }
    }
    for (unsigned int t = 0, nt = telDataBinZeta.size(); t < nt; ++t) {
      for (unsigned int z = 0, nz = telDataBinZeta[t].size(); z < nz; ++z) {
        if (predLF[startBin+t] > 0)
          telDataBinZeta[t][z].second /= predLF[startBin + t];
        else
          telDataBinZeta[t][z].second = 0;
      }
    }
  }

  const int tOff = timeOffset;
  for (int i = -nConv + 1; i < nConv; ++i) {
    const double w = fErrFuncFactors[abs(i)];
    for (int tConv = 0; tConv < convolutedTraceSize; ++tConv) {
      const int tOrig = tConv + tOff + i;
      // zeta weight factors
      double wzeta = 1;
      if (fAafCorrection == 2 && tOrig >= tOff && tOrig < (convolutedTraceSize + tOff)) {
        wzeta = 0;
        const int tO = tOrig - tOff;
        const int tC = tConv;
        for (unsigned int zO = 0, nO = telDataBinZeta[tO].size(); zO < nO; ++zO) {
          for (unsigned int zC = 0, nC = telDataBinZeta[tC].size(); zC < nC; ++zC) {
            if (telDataBinZeta[tO][zO].first == telDataBinZeta[tC][zC].first) {
              wzeta += telDataBinZeta[tO][zO].second; // relative to the total ldf light fraction -> same zeta pixels => wzeta == 1
              break;
            }
          }
        }
      }
      convolutedTrace[tConv] += wzeta * w * originalTrace[tOrig];
    }
  }
}