SdReconstruction/LDFFinderKG/LDFFinder.h

Go to the documentation of this file.

#ifndef _LDFFinderKG_LDFFinder_h_
#define _LDFFinderKG_LDFFinder_h_

#include <fwk/VModule.h>
#include <utl/Point.h>
#include <utl/Vector.h>
#include <utl/TimeStamp.h>
#include <utl/CovarianceMatrix.h>
#include <sdet/ArrayType.h>

#include <cmath>
#include <sstream>

#include "LDF.h"
#include "EnergyConversion.h"
#include "LikelihoodFunctions.h"


namespace evt {
  class Event;
}

namespace sevt {
  class SEvent;
}

namespace LDFFinderKG {

  utl::CoordinateSystemPtr gBaryCS; // is used by LikelihoodFunctions, too


  class LDFFinder : public fwk::VModule {
  public:
    fwk::VModule::ResultFlag Init() override;
    fwk::VModule::ResultFlag Run(evt::Event& event) override;
    fwk::VModule::ResultFlag Finish() override { return eSuccess; }

    enum InfoLevel {
      eNone = 0,
      eFinal = 1,
      eIntermediate = 2,
      eObscure = 3,
      eMinuit = 4
    };

    enum LDFParameter {
      eShowerSize = 0,
      eCoreX = 1,
      eCoreY = 2,
      eShapeParameters = 3
    };

    struct LDFFitResult {

      LDFFitResult(const LDFFitConfig& config) :
        fPar(eShapeParameters + config.fLDF.GetNShapeParameters()),
        fCov(fPar.size())
      { }

      std::vector<double> fPar;
      utl::CovarianceMatrix fCov;

      double fChi2 = 0;
      int fNdof = 0;
      double fLogLikelihood = 0;
      double fRecStage = 0;

      // for convenience

      utl::Point GetCore() const
      { return utl::Point(fPar[eCoreX], fPar[eCoreY], 0, gBaryCS); }

      std::vector<double> GetLDFShapeParameters() const
      { return std::vector<double>(fPar.begin() + eShapeParameters, fPar.end()); }

    };

    enum CurvatureParameter {
      eU = 0, // u = sin(theta)*cos(phi)
      eV = 1, // v = sin(theta)*sin(phi)
      eRCore = 2,
      eCTCore = 3
    };

    struct CurvFitResult {

      std::vector<double> fPar = std::vector<double>(4);
      utl::CovarianceMatrix fCov = utl::CovarianceMatrix(4);
      double fChi2 = 0;
      int fNdof = 0;

      // for convenience

      double
      GetTheta()
        const
      {
        const double u = fPar[eU];
        const double v = fPar[eV];
        const double w = std::sqrt(1 - u*u - v*v);
        return std::acos(w);
      }

      utl::Vector
      GetAxis()
        const
      {
        const double u = fPar[eU];
        const double v = fPar[eV];
        const double w = std::sqrt(1 - u*u - v*v);
        return utl::Vector(u, v, w, gBaryCS);
      }

    };

  private:
    // LDF part
    bool FixBeta(const LDFFitResult& result,
                 const utl::Vector& showerAxis,
                 const std::vector<StationFitData>& data) const;

    bool FixGamma(const LDFFitResult& result,
                  const utl::Vector& showerAxis,
                  const std::vector<StationFitData>& data) const;

    void FitLDFSimplified(LDFFitResult& result,
                          LDFFitConfig& config,
                          const utl::Vector& showerAxis,
                          const std::vector<StationFitData>& data) const;

    bool FitLDF(LDFFitResult& result,
                LDFFitConfig& config,
                const utl::Vector& showerAxis,
                const std::vector<StationFitData>& data) const;

    // Curvature part
    double ParameterizedRc(const double cosTheta, const double showerSize) const;

    double FitCurvature(CurvFitResult& result,
                        const utl::Point& core,
                        const std::vector<StationFitData>& data) const;

    // other
    std::vector<StationFitData> MakeStationFitData(const sevt::SEvent& sEvent) const;

    bool FixCore(utl::Point& core,
                 utl::TimeStamp& coreTime,
                 utl::Vector& showerAxis,
                 const evt::Event& event) const;

    void SetRecData(evt::Event& event,
                    const LDFFitResult& lresult,
                    const CurvFitResult& cresult) const;

    void OutputResults(const evt::Event& event) const;

    bool FitGlobal(CurvFitResult& cres, LDFFitResult& lres,
                   const LDFFitConfig& config,
                   const std::vector<StationFitData>& data) const;
    std::vector<int> CleanSilentStation(const LDFFitResult& ldfres,
                                        const LDFFitConfig& config,
                                        const utl::Vector& showerAxis,
                                        std::vector<StationFitData>& stationFitData) const;
    // vars
    int fInfoLevel = 0;
    sdet::Array::Type fArrayType = sdet::Array::e1500;
    double fMaxTheta = 0;
    double fMaxBaryToCoreDistance = 0;
    double fMaxChi2 = 0;
    bool fUseMaxLikelihood = false;
    bool fUseSilentStationsFromOtherGrids = false;
    CoreType fCoreType;
    bool fFixedAxis = false;
    bool fFixedCurvature = false;
    bool fGlobalFit = false;
    utl::Point fBarycenter;
    utl::TimeStamp fBaryTime;
    LDFFitConfig fLDFFitConfig;
    EnergyConversion fEnergyConversion;
    int fRoptN = 0;
    std::vector<int> fBadSilents;

    struct Stage0 {
      int fMinNumberOfStations = 0;
      bool fUseCorePosition = false;
    } fStage0;

    struct Stage1 {
      int fMinNumberOfStations = 0;
    } fStage1;

    struct Stage2 {
      int fMinNumberRelaxBeta = 0;
      std::vector<double> fFixBetaRange;
      std::vector<double> fFixBetaLeverArmRange;
      int fMinNumberRelaxGamma = 0;
      std::vector<double> fFixGammaRange;
      std::vector<double> fFixGammaLeverArmRange;
    } fStage2;

    struct Stage3 {
      int fMinNumberOfStations = 0;
      int fMinNumberForFullCurvatureFit = 0;
      double fMaxAxisAngleDifference = 0;
      bool fOutliersRejection = false;
    } fStage3;

    struct Silents {
      int fMinNumberOfStations = 0;
      double fMaxLDFValue = 0;
      double fMaxDistance = 0;
    } fSilentsVars;

    REGISTER_MODULE("LDFFinderKG", LDFFinder);

  };

}


#endif