RdGeoCeLDFFitter/LikelihoodFunction.h

Go to the documentation of this file.

#ifndef _RdGeoCeLDFFitter_LikelihoodFunction_h_
#define _RdGeoCeLDFFitter_LikelihoodFunction_h_

#include <iostream>
#include <cmath>
#include <vector>
#include <limits>

#include <Minuit2/FCNBase.h>

#include <utl/CoordinateSystemPtr.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/Math.h>
#include <utl/PhysicalFunctions.h>


namespace RdGeoCeLDFFitter {

  double getAtmosphere(double h);
  double getDensity(double h);
  double getVerticalHeight(double at);

  struct FitConfig {
    bool useStationsWithoutSignal;
    bool fitCore;

    FitConfig() :
        useStationsWithoutSignal(true),
        fitCore(true)
    {
    }
  };

  struct StationFitData {
    double x;
    double y;
    double f;
    double ferror;
    double f_vB;
    double ferror_vB;
    double f_vvB;
    double ferror_vvB;
    double distance;
    bool hasSignal;

    StationFitData() :
        x(0),
        y(0),
        f(0),
        ferror(0),
        f_vB(0),
        ferror_vB(0),
        f_vvB(0),
        ferror_vvB(0),
        distance(0),
        hasSignal(false)
    {
    }
  };

  struct EventFitData {
    utl::Vector fShowerAxis;
    utl::CoordinateSystemPtr fLocalCS;
    double sinalpha;
    double zenith;

    EventFitData() :
        fShowerAxis(),
        fLocalCS(),
        sinalpha(1),
        zenith(0)
    {
    }

  };

  class LDFLikelihoodFunction : public ROOT::Minuit2::FCNBase {

   public:
    LDFLikelihoodFunction(const FitConfig ldfconfig, const EventFitData eventData,
                          const std::vector<StationFitData>& stationData);

    double operator()(const std::vector<double>& pars) const override;

    double GetChi2FullEnergyFluence(const std::vector<double>& pars) const;

    double Up() const override
    {
      return fTheErrorDef;  // needs to be one because the neg. log. likelihood is defined as -2*log(Likelihood)
    }

    void SetErrorDef(double def) override
    {
      fTheErrorDef = def;
    }

   protected:

   private:
    double GetChi2LDFModel1(const double x, const double y, const double Erad, const double dxmax,
                            const double coreX, const double coreY, const double f_vB,
                            const double ferror_vB, const double f_vvB, const double ferror_vvB,
                            const double zenith, const double sinalpha) const;



    const FitConfig fFitConfig;

    const EventFitData fEventData;
    const std::vector<StationFitData>& fStationData;

    utl::Vector fMagneticFieldVector;

    double fTheErrorDef;
  };

}

#endif