LDFFinderKG/LikelihoodFunctions.h

Go to the documentation of this file.

#ifndef _LDFFinderKG_LikelihoodFunctions_h_
#define _LDFFinderKG_LikelihoodFunctions_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>

#include "LDF.h"


namespace LDFFinderKG {

  extern utl::CoordinateSystemPtr gBaryCS; // defined in LDFFinder.h


  // perpendicular distance from the axis
  inline
  double
  RPerp(const utl::Vector& axis, const utl::Vector& station)
  {
    const double scal = axis*station;
    return std::sqrt(station*station - scal*scal);
  }


  enum CoreType {
    eMC,
    eHybrid,
    eFit,
  };


  enum ParameterTreatment {
    eModeled = 0,
    eFitted,
    eFixed
  };


  struct StationFitData {
    unsigned int fId = 0;
    bool fIsSaturated = false;
    utl::Point fPos;
    double fSignal = 0;
    double fRecoveredSignal = 0;
    double fRecoveredSignalErr = 0;
    double fCTime = 0;
    double fT50 = 0;
    double fGPSTimeVariance = 0;
  };


  struct LDFFitConfig {
    bool fUseSilentStations = false;
    bool fUseSaturationRecovery = false;
    bool fFixCore = false;
    double fSilentMaxRadius = 0;
    double fSilentRadiusTransition = 0;
    double fSilentSignalThreshold = 0;
    double fRecoveryThreshold = 0;
    LDF fLDF;
    utl::wcd::ESignalVarianceModel fSignalVarianceModel = utl::wcd::eGAP2014_035;
    std::vector<ParameterTreatment> fLDFShapeTreatment;
    bool fUseAdditionalGPSTimeVariance = false;
  };


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

  public:
    LDFLikelihoodFunction(const LDFFitConfig& config,
                          const utl::Vector& showerAxis,
                          const std::vector<StationFitData>& data)
      : fConfig(config), fShowerAxis(showerAxis), fData(data) { }

    double
    operator()(const std::vector<double>& pars)
      const
    {
      // One comment on the following conversion factor called PoissonFactor:
      // The factor reflects the electron-gamma to muon ratio and gives rise to
      // criticism of the algorithm. But as long as we would like to include
      // "active" Zero stations, we have to assume a conversion factor. Otherwise we
      // cannot employ a maximum likelihood method! By definion a \chi^2 method gives
      // biased results (more strongly stated: in a mathematical sense they are
      // wrong) !
      // This factor has to be studied in detail by means of simulated AND real
      // data. See Haverah Park analysis for more details (Lapikens? or England?)

      const double showerSize = pars[0];

      const utl::Point core(pars[1], pars[2], 0, gBaryCS);
      const double cosTheta = fShowerAxis.GetCosTheta(gBaryCS);

      std::vector<double> shape(pars.begin()+3, pars.end());
      const std::vector<double> modeledShape =
        fConfig.fLDF.ShapeModel(cosTheta, showerSize);
      for (unsigned int i = 0, n = shape.size(); i < n; ++i)
        if (fConfig.fLDFShapeTreatment[i] == eModeled)
          shape[i] = modeledShape[i];

      const double sigmaFactor = utl::wcd::SignalUncertaintyFactor(fConfig.fSignalVarianceModel, cosTheta);
      const double poissonFactor = utl::wcd::PoissonFactor(sigmaFactor);

      double logLikelihood = 0;

      for (std::vector<StationFitData>::const_iterator sIt = fData.begin(), end = fData.end();
           sIt != end; ++sIt) {
        const utl::Vector coreToStation = sIt->fPos - core;
        const double rho = RPerp(fShowerAxis, coreToStation);

        const double signalPrediction = showerSize * fConfig.fLDF(rho, shape);  // TODO: add asymmetry
        const double particlePrediction = poissonFactor * signalPrediction;

        if (sIt->fSignal > 0) {  // candidate

          const double particlesInStation = poissonFactor * sIt->fSignal;
          const double sigma = poissonFactor *
            utl::wcd::SignalUncertainty(fConfig.fSignalVarianceModel, cosTheta, signalPrediction);

          if (sIt->fIsSaturated) {
            if (fConfig.fUseSaturationRecovery && sIt->fRecoveredSignal > 0 && sIt->fRecoveredSignalErr > 0) {
              // use recovered signal, add uncertainty of saturation recovery to signal uncertainty
              const double n = poissonFactor * sIt->fRecoveredSignal;
              const double dn = poissonFactor * sIt->fRecoveredSignalErr;
              const double s = std::sqrt(utl::Sqr(sigma) + utl::Sqr(dn));
              logLikelihood += utl::LogarithmOfNormalPDF(n, particlePrediction, s);
            } else {
              const double signal = (sIt->fRecoveredSignal > sIt->fSignal) ?
                (sIt->fRecoveredSignal - sIt->fRecoveredSignalErr) : sIt->fSignal;
              const double n = poissonFactor * signal;
              logLikelihood += utl::LogarithmOfNormalComplementCDF(n, particlePrediction, sigma);
            }
          } else {
            if (particlesInStation > 30) {

              // Stations with large signals: approximately gaussian probability
              logLikelihood += utl::LogarithmOfNormalPDF(particlesInStation, particlePrediction, sigma);

            } else {

              // Stations with low signal: Poisson probability
              logLikelihood += particlesInStation * std::log(particlePrediction) - particlePrediction;

              double logarithmicFactor = 0;
              const int stop = int(particlesInStation + 0.5);
              for (int j = 1; j <= stop; ++j)
                logarithmicFactor += std::log(double(j));

              logLikelihood -= logarithmicFactor;

            }
          }
        } else if (fConfig.fUseSilentStations) {  // silent stations

          logLikelihood += particlePrediction > 0.03 ?
            -particlePrediction +
            std::log(1 +                            // 0
                     particlePrediction * (1 +      // 1
                     0.5*particlePrediction * (1 +  // 2
                     (1./3)*particlePrediction      // 3
                     )))
            :
            -(1./24)*utl::Sqr(utl::Sqr(particlePrediction));

          // this threshold above is too large
          // it should be replaced with the 1 - trigger probability but depending on
          // whether new triggers are silenced or not...

        }

      } // loop over stations

      // max(LogLikelihood) = min(-LogLikelihood)
      return std::isnan(logLikelihood) ? std::numeric_limits<double>::max() : -logLikelihood;
    }

    double Up() const { return 0.5; }

  protected:
    const LDFFitConfig& fConfig;
    const utl::Vector& fShowerAxis;
    const std::vector<StationFitData>& fData;

  };


  class CurvatureChi2Function : public ROOT::Minuit2::FCNBase {

  public:
    CurvatureChi2Function(const utl::Point& core,
                          const std::vector<StationFitData>& data)
      : fCore(core), fData(data), fTimeVarModel(sdet::STimeVariance::GetInstance()) { }

    double
    operator()(const std::vector<double>& pars)
      const
    {
      const double& pu     = pars[0];
      const double& pv     = pars[1];
      const double& rCore  = pars[2];
      const double& cTCore = pars[3];

      const double pw = std::sqrt(1 + pu*pu + pv*pv);
      const double u = pu/pw;
      const double v = pv/pw;

      const double w = std::sqrt(1 - u*u - v*v);

      const utl::Vector axis(u,v,w, gBaryCS);

      double chi2 = 0;

      for (const auto& st : fData) {
        if (!st.fSignal)
          continue;

        const utl::Vector stationToCore = fCore - st.fPos;
        const utl::Vector stationToCenter = stationToCore + rCore*axis;
        const double rStation = stationToCenter.GetMag();

        const double cosThetaStation = stationToCenter.GetZ(gBaryCS) / rStation;

        const double sigma2 =
          fTimeVarModel.GetTimeSigma2(st.fSignal, st.fT50, cosThetaStation, RPerp(axis, stationToCore)) +
          st.fGPSTimeVariance;

        const double cdt = (st.fCTime - rStation) - (cTCore - rCore);
        chi2 += utl::Sqr(cdt) / sigma2;
      }

      return std::isnan(chi2) ? std::numeric_limits<double>::max() : chi2 / utl::Sqr(utl::kSpeedOfLight);
    }

    double Up() const { return 1; }

  protected:
    const utl::Point& fCore;
    const std::vector<StationFitData>& fData;
    const sdet::STimeVariance& fTimeVarModel;

  };


  class GlobalFitFunction : public ROOT::Minuit2::FCNBase {

  public:
    GlobalFitFunction(const LDFFitConfig& config, const std::vector<StationFitData>& data)
      : fCurvPart(fCore, data), fLDFPart(config, fShowerAxis, data) { }

    double
    operator()(const std::vector<double>& pars)
      const
    {
      const double u = pars[0];
      const double v = pars[1];
      const double w = std::sqrt(1 - u*u - v*v);

      fShowerAxis = utl::Vector(u, v, w, gBaryCS);

      fCore = utl::Point(pars[5], pars[6], 0, gBaryCS);

      const std::vector<double> curvPars(pars.begin(), pars.begin()+4);
      const std::vector<double> ldfPars(pars.begin()+4, pars.end());
      return 0.5*fCurvPart(curvPars) + fLDFPart(ldfPars);
    }

    double Up() const { return 0.5; }

  protected:
    mutable utl::Vector fShowerAxis;
    mutable utl::Point fCore;
    const CurvatureChi2Function fCurvPart;
    const LDFLikelihoodFunction fLDFPart;

  };


  class LDFChi2Function : public ROOT::Minuit2::FCNBase {

  public:
    LDFChi2Function(const LDFFitConfig& config,
                    const utl::Vector& showerAxis,
                    const std::vector<StationFitData>& data)
      : fConfig(config), fShowerAxis(showerAxis), fData(data) { }

    double
    operator()(const std::vector<double>& pars)
      const
    {
      const double showerSize = pars[0];

      const utl::Point core(pars[1], pars[2], 0, gBaryCS);

      const double cosTheta = fShowerAxis.GetCosTheta(gBaryCS);
      std::vector<double> shape(pars.begin()+3, pars.end());
      const std::vector<double> modeledShape =
        fConfig.fLDF.ShapeModel(cosTheta, showerSize);
      for (unsigned int i = 0, n = shape.size(); i < n; ++i)
        if (fConfig.fLDFShapeTreatment[i] == eModeled)
          shape[i] = modeledShape[i];

      double chi2 = 0;

      for (std::vector<StationFitData>::const_iterator sIt = fData.begin(), end = fData.end();
           sIt != end; ++sIt) {
        const double rho = RPerp(fShowerAxis, sIt->fPos - core);

        if (sIt->fSignal > 0) {  // candidate
          const double predictedSignal = showerSize * fConfig.fLDF(rho, shape);
          const bool useRecovery = fConfig.fUseSaturationRecovery &&
                                   sIt->fRecoveredSignal > 0 &&
                                   sIt->fRecoveredSignalErr > 0;

          if (sIt->fIsSaturated && !useRecovery)
            continue;

          const double signal = sIt->fIsSaturated ? sIt->fRecoveredSignal : sIt->fSignal;
          const double sigmaNotSat = utl::wcd::SignalUncertainty(fConfig.fSignalVarianceModel, cosTheta, predictedSignal);

          const double sigma = sIt->fIsSaturated ?
            std::sqrt(utl::Sqr(sIt->fRecoveredSignalErr) + utl::Sqr(sigmaNotSat)) : sigmaNotSat;

          // we do NOT add any contribution to the chi2, as it would be a bias.
          const double dev = (signal - predictedSignal) / sigma;
          chi2 += utl::Sqr(dev);
        } else if (fConfig.fUseSilentStations) {  // silent
          // Warning: the treatment of silent stations in the Chi2 is deprecated
          // and only kept for making comparisons with CDAS
          chi2 +=
            utl::Fermi(rho, fConfig.fSilentMaxRadius, fConfig.fSilentRadiusTransition) *
              showerSize * fConfig.fLDF(rho, shape) / fConfig.fSilentSignalThreshold;
        }
      } // loop over stations

      return std::isnan(chi2) ? std::numeric_limits<double>::max() : chi2;
    }

    double Up() const { return 1; }

  protected:
    const LDFFitConfig& fConfig;
    const utl::Vector& fShowerAxis;
    const std::vector<StationFitData>& fData;

  };


}


#endif