EnergyFitter.h

Go to the documentation of this file.

#ifndef _EnergyFitterKG_EnergyFitter_h_
#define _EnergyFitterKG_EnergyFitter_h_

#include <fwk/VModule.h>


namespace fevt {
  class Eye;
}

namespace oBLAS {
  class lowerTriangularMatrix;
}

namespace FdProfileReconstructorKG {
  class CherenkovFluorescenceMatrix;
}

namespace FdProfileReconstructorKG {

  class EnergyFitter {

  public:
    bool Init();
    bool Run(fevt::Eye& eye, CherenkovFluorescenceMatrix* const chfl);
    bool Finish() const { return true; }

    bool ReFitProfile(const fevt::Eye& eye, CherenkovFluorescenceMatrix* const chfl);
    double GetXmax() const;
    double GetX0() const;
    double GetdEdXmax() const;
    double GetLambda() const;
    double GetEmEnergy() const;
    int GetVerbosity() const { return fVerbosity; }
    void SetVerbosity(const int iv) { fVerbosity = iv; }
    void SetYieldRefit(const bool what) { fYieldRefit = what; }

  private:
    enum EFitOption {
      eInitialize,
      eEnergyDeposit,
      eLight2Param,
      eLightAllParam,
      eLightEnergy
    };

    enum EMinimizationType {
      eMaximumLikelihood,
      eLeastSquares,
      eModifiedLeastSquares
    };

    enum EGHFitParameters {
      eEnergy = 0,  // or dEdXmax
      eXmax,
      eX0,
      eLambda,
      eNGHParameters
    };

    enum EParConstraints {
      eFree,
      eFixed,
      eConstrained
    };

    bool InitializeGHFit(const fevt::Eye&, const int iMode);
    bool FitProfile(const EFitOption opt);
    bool CalculateEnergy(fevt::Eye&);
    void FillProfilesAtAperture(fevt::Eye&);
    bool UpDateFCM();
    void CleanUp();
    void DumpCurrentParameters(const EFitOption) const;

    // MINUIT fit functions
    static void dEdXFitFunction(int&, double*, double&, double*, int);
    static void LightFitFunction(int&, double*, double&, double*, int);
    static void EnergyFitFunction(int&, double*, double&, double*, int);
    static double GaisserHillasChi2(const double, const double,
                                    const double, const double);
    static double GaisserHillasLogLike(const double, const double,
                                       const double, const double);

    // dEdX profile
    static unsigned int fnBins;  // number of bins
    static double* fdEdX;        // energy deposit
    static double* fvdEdX;       // its variance
    static double* fLight;       // light at diaphragm (measured)
    static double* fvLight;      // its variance from FdApertureLightFinder
    static double* fvBackGround; // baseline variance
    static double* fDepth;       // slant depth
    static double* fChFlDepth;

    // telescope properties (for Npe calculation)
    static double fDiaphragmArea;
    static double fGainVariance;
    static double fPhotonToPhotoElectron;

    // cherenkov fluorescence matrix
    CherenkovFluorescenceMatrix* fChFlPtr = nullptr;
    static const oBLAS::lowerTriangularMatrix* fChFl;

    // current GH parameters and their uncertainties
    double fX0 = 0;
    double fXmax = 0;
    double fdEdXmax = 0;
    double fLambda = 0;
    double fCovariance[eNGHParameters][eNGHParameters] = { { 0 } };
    double fChisqGH = 0;
    unsigned int fnDof = 0;
    double fEem = 0;

    // options
    static EParConstraints fFixX0;
    static EParConstraints fFixLambda;
    static EParConstraints fFixkUniv;
    static double fDefX0;
    static double fDefLambda;
    static double fVarX0;
    static double fVarLambda;

    static double fDefA1;
    static double fDefA2;
    static double fVarK;

    static EMinimizationType fMinimizationMethod;
    int fVerbosity = 0;
    bool fYieldRefit = false;

  };

}


#endif