FdEnergyDepositFinderKG/ProfileFitter.h

Go to the documentation of this file.

#ifndef _FdEnergyDepositFinderKG_ProfileFitter_h_
#define _FdEnergyDepositFinderKG_ProfileFitter_h_

#include "ProfileFitFunctions.h"
#include "GHShapeParameters.h"
#include "CFMatrixCalculator.h"

#include <fdet/Telescope.h>

#include <utl/Function.h>

#include <cmath>
#include <vector>


class TMinuit;
namespace utl {
  class TabulatedFunctionErrors;
}

namespace evt {
  class GaisserHillas4Parameter;
}


namespace FdEnergyDepositFinderKG {

  class CFMatrixCalculator;


  class ProfileFitter {

  public:
    enum EFitMode {
      eNMax,
      eIntegral
    };

  public:
    ProfileFitter();
    void Init();

    void SetVerbosity(const int v) { fVerbosity = v;};

    void SetFunctionType(const EFitFunctionType type) { fFunctionType = type; }

    void SetFitMode(const EFitMode mode) { fFitMode = mode; }

    void FillGHParameters(evt::GaisserHillas4Parameter& ghPars) const;

    bool Fit();

    void SetProfileData(const utl::TabulatedFunctionErrors&);
    void SetLightFluxData(const CFMatrixCalculator& cfMatrixData) { SetLightFluxData(cfMatrixData, cfMatrixData); }
    void SetLightFluxData(const CFMatrixCalculator& cfMatrixData, const CFMatrixCalculator& cfMatrixDataDense);
    void SetStartParameters(const evt::GaisserHillas4Parameter&);
    void SetUnivConstrained(const bool constrained, const utl::Function& func, const double ksigma);
    void SetShapePars(const GHShapeParameters&);
    void ResetStartParameters() { fStartParameters.clear(); }

    void GetEnergy(double& energy, double& energyError) const;

    double GetShapeConstraint1Chi2() const { return fShapeConstraintChi2[0]; }
    double GetShapeConstraint2Chi2() const { return fShapeConstraintChi2[1]; }
    double GetShapeConstraint3Chi2() const { return fShapeConstraintChi2[2]; }

    unsigned int GetNDof() const { return fDepth.size() - (fFunctionType == eGH2dEdXChi2 ? 2 : 4); }

    void SetAafCorrection(const int aaf = 1, const bool useNoiseBins = true) { fAafCorrection = aaf; fUseNoiseBins = useNoiseBins; }
    double GetLHChi2() const { return fLHChi2; }

  private:
    bool GHFit();
    void InitMinuit(TMinuit&) const;
    bool Minimize(TMinuit&);

    double CorrCoeff(const EGHFunctionPar i, const EGHFunctionPar j) const
    { return fCovariance[i][j] / std::sqrt(fCovariance[i][i]*fCovariance[j][j]); }

    static void GHFitFunction(int&,double* const,double&, double* const, const int);

    static double GaisserHillasLogLike(const double* const par);
    static double GaisserHillasLogLikeConvoluted(const double* const par);
    static double GaisserHillasChi2(const double* const par);

    static void PrepareTimeConvolution(const fdet::Telescope& detTel);
    static void DoTimeConvolution(const std::vector<double>& originalTrace,
                                  std::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);

    static std::vector<std::vector<double>> fErrFuncFactorsBuf;

    int fVerbosity;
    static EFitFunctionType fFunctionType;
    static std::vector<double> fDepth;
    static std::vector<double> fSize;
    static std::vector<double> fSizeVariance;
    static bool fIsConstrained;

    static std::string fGaisserHillasType;

    static bool fIsUnivConstrained;
    static utl::Function fUnivFunction;
    static double fVarK;

    static bool fEMGLConstrained;
    static double fTauOverSigma;

    static int fAafCorrection;
    static bool fUseNoiseBins;

    static const CFMatrixCalculator* fCFMatrixData;
    static const CFMatrixCalculator* fCFMatrixDataDense;

    // GaisserHillas shape parameters, constraints etc
    static GHShapeParameters fGHShapePars;

    // fit result
    std::vector<double> fFitParameters;
    std::vector<double> fStartParameters;
    double fCovariance[eNGHpars][eNGHpars];
    double fChi2;
    double fLHChi2;
    int fStatus;
    static unsigned int fNdof;
    static EFitMode fFitMode;

    static double fShapeConstraintChi2[3];

  };

}


#endif