ScintillatorLDFFinderKG/LikelihoodFunctions.h

Go to the documentation of this file.

#ifndef _ScintillatorLDFFinderKG_LikelihoodFunctions_h_
#define _ScintillatorLDFFinderKG_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 ScintillatorLDFFinderKG {

  extern utl::CoordinateSystemPtr gBaryCS; // defined in ScintillatorLDFFinder.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);
  }

  inline
  double
  DistanceCut(const double showerSize, const double cosTheta)
  {
    const double s = std::log10(showerSize);
    const double sec = 1. / cosTheta;
  
    const double f0 = 303.242;
    const double f1 = -355.686;
    const double f2 = 171.659;
    const double f3 = 586.443;
    const double f4 = 1309.454;
    const double f5 = -544.919;
    const double f6 = 156.465;
    const double f7 = -558.186;
    const double f8 = 236.191;
    
    const double distcut =  f0 + f1*sec + f2*sec*sec 
                         + (f3 + f4*sec + f5*sec*sec)*s 
                         + (f6 + f7*sec + f8*sec*sec)*s*s;

    return distcut;
  }

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


  struct StationFitData {
    unsigned int fId = 0;
    bool fIsSaturated = false;
    utl::Point fPos;
    double fSignal = 0;
  };


  struct LDFFitConfig {
    // add use silent
    bool fFixCore;
    LDF fLDF;
    utl::ssd::ESignalVarianceModel fSignalVarianceModel;
    std::vector<ParameterTreatment> fLDFShapeTreatment;
    bool fCorePropagationEnabled;
    double deltaR;
  };


  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
    {

      /*
      The concept of Poisson Factor is also used for scintillator signals in order to obtain Poisson p.d.f.s for small signals.
      The conversion from signal to effective number of particles is given by:
      
      Neff = p*S,
      
      where the factor p is derived in such a way that the observed signal uncertainties are formulated in terms of an effective particle number 
      under a Poisson assumption (see SD reconstrucion manual GAP2005_035).
      
      */
      
      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::ssd::SignalUncertaintyFactor(fConfig.fSignalVarianceModel, cosTheta);
      const double poissonFactor = utl::ssd::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);
        const double particlePrediction = poissonFactor * signalPrediction;

        if (sIt->fSignal > 0) {

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

          if (sIt->fIsSaturated) // ToDo: treat saturated SSDs
            continue;

          if (particlePrediction > 30) {
            // large signals: approximately gaussian probability
            logLikelihood += utl::LogarithmOfNormalPDF(particlesInScintillator, particlePrediction, sigma);
          } else {
            // low signals: Poisson probability
            logLikelihood += particlesInScintillator * std::log(particlePrediction) - particlePrediction;
            double logarithmicFactor = 0;
            const int stop = int(particlesInScintillator + 0.5);
            for (int j = 1; j <= stop; ++j)
              logarithmicFactor += std::log(double(j));

            logLikelihood -= logarithmicFactor;

          }
        } // if Signal > 0
      } // loop over stations

      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 LDFLikelihoodFunctionTN : public ROOT::Minuit2::FCNBase {

  public:
    LDFLikelihoodFunctionTN(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
    {

      /*
        Usage of a truncated normal distribution for describing SSD signals. 
        As SSD is triggered by the WCD, SSD signals can be low as to integral 
        of baseline (~0 MIPs). A truncated normal distribution was proposed 
        as a p.d.f of the SSD signals and a corresponding likelihood function 
        is derived under the assumption that such p.d.f should work for both 
        small and large signals. No Poisson factor is needed. 
      */
      
      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 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 signal = sIt->fSignal;
        const double signalPrediction = showerSize * fConfig.fLDF(rho, shape);

        /*
          Zero and negative signals should be properly handled in the likelihood.
          We keep this like that until signal uncertainty and likelihood are 
          well understood for small signals. Until then, a distance cut is applied 
          to reject low or zero signals according to a parameterization in 
          lgS and sec theta. 
        */

        if (signal > 0) {

          double sigma = utl::ssd::SignalUncertainty(fConfig.fSignalVarianceModel, cosTheta, signalPrediction);
          /*
            If desired, the core uncertainties deltaR can be propagated via Gaussian error propagation as sigmaR
            to the signal uncertainty.
          */
          const bool corePropagationEnabled = fConfig.fCorePropagationEnabled;
          if (corePropagationEnabled){
            const double deltaR = fConfig.deltaR;
            const double sigmaR = deltaR * showerSize * fConfig.fLDF.FirstDerivative(rho,shape);
            sigma = std::sqrt(utl::Sqr(sigma) + utl::Sqr(sigmaR));
          }
          
          /* 
            Current LDF parameterization was determined using no saturated SSDs.
            This would need further studies in the future. 
          */
          if (sIt->fIsSaturated) 
            continue;

          logLikelihood += - 0.5*std::log(2*M_PI*utl::Sqr(sigma))
                           - 0.5*utl::Sqr((signal - signalPrediction) / sigma)
                           - std::log(0.5*(1 + std::erfc(signalPrediction / (std::sqrt(2)*sigma))));
                      
        } // if Signal > 0
      } // loop over stations

      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 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);
          if (sIt->fIsSaturated)
            continue;

          const double signal = sIt->fSignal;
          const double sigma = utl::ssd::SignalUncertainty(fConfig.fSignalVarianceModel, cosTheta, predictedSignal);

          // 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);
        } // signal > 0
      } // 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