JAspera.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JASPERA__
#define __JASTRONOMY__JASPERA__
 
#include <vector>
#include <cmath>
 
 
/**
 * \file
 *
 * Per aspera ad astra.
 * \author mdejong
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  /**
   * Auxiliary data structure to fit signal strength using likelihood ratio.
   */
  struct JAspera :
    public std::vector<double>                      // N/S values
  {
 
    static constexpr double EPSILON = 1.0e-3;       //!< precision determination of signal strength
 
 
    /**
     * Result of fit.
     */
    struct fit_type {
      double likelihood  =  0.0;                    //!< likelihood
      double signal      =  0.0;                    //!< signal strength
    };
 
 
    /**
     * Put signal and background to list of pre-computed N/S values.
     *
     * \param  s                signal
     * \param  b                background
     */
    void put(const double s,
             const double b)
    {
      push_back(b/s);
 
      ws += s;
    }
 
 
    /**
     * Put signal and background to list of pre-computed N/S values.
     *
     * \param  n                data
     * \param  s                signal
     * \param  b                background
     */
    void put(const size_t n,
             const double s,
             const double b)
    {
      for (size_t i = 0; i != n; ++i) {
        push_back(b/s);
      }
 
      ws += s;
    }
 
 
    /**
     * Get likelihood for given signal strength.
     *
     * \param  p                signal strength
     * \return                  likelihood
     */
    double getLikelihood(const double p) const
    {
      double y = -p * ws;
 
      for (const auto& i : static_cast<const std::vector<double>&>(*this)) {
        y += log1p(p/i);
      }
 
      if (y > 0.0)
        return y;
      else
        return 0.0;
    }
 
 
    /**
     * Get derivative of likelihood for given signal strength.
     *
     * \param  p                signal strength
     * \return                  derivative of likelihood
     */
    double getDerivative(const double p) const
    {
      double y = -ws;
 
      for (const auto& i : static_cast<const std::vector<double>&>(*this)) {
        y += 1.0 / (p + i);
      }
 
      return y;
    }
 
 
    /**
     * Fit signal strength.
     *
     * \return                  result
     */
    fit_type operator()() const
    {
      using namespace std;
 
      double x1 = 0.0;                          // lower limit corresponds to no signal; less to under-fluctuation
      double x2 = (double) this->size() / ws;   // upper limit corresponds to no background (i.e. all N/S = 0)
      // bisection
      /*
      if (getDerivative(x1) > 0.0) {
 
        for ( ; ; ) {
 
          const double x = 0.5 * (x1 + x2);
 
          if (x2 - x1 <= EPSILON) {
            return { getLikelihood(x), x };
          }
 
          if (getDerivative(x) < 0.0)
            x2 = x;
          else
            x1 = x;
        }
      }
      */
      // Ridder's method
      double f1 = getDerivative(x1);
 
      if (f1 > 0.0) {
  
        double f2 = getDerivative(x2);
 
        while (x2 - x1 > EPSILON) {
 
          const double xm = 0.5 * (x1 + x2);
          const double fm = getDerivative(xm);
 
          const double s  = sqrt(fm*fm - f1*f2);
 
          if (s == 0.0) {
            break;
          }
 
          const double xn = xm + (xm - x1) * fm/s;
          const double fn = getDerivative(xn);
 
          if (signbit(fn) != signbit(fm)) {
 
            x1 = xm;
            f1 = fm;
            x2 = xn;
            f2 = fn;
 
          } else {
 
            if (signbit(fn)) {
 
              x2 = xn;
              f2 = fn;
 
            } else {
 
              x1 = xn;
              f1 = fn;
            }
          }
        }
 
        const double x = 0.5 * (x1 + x2);
 
        return { getLikelihood(x), x };
      }
      //
      return { 0.0, 0.0 };
    }
 
 
    /**
     * Get total signal strength.
     *
     * \return           signal strength
     */
    double getSignal() const
    {
      return ws;
    }
 
 
    /**
     * Set signal strength.
     *
     * \param  wS        signal strength
     */
    void setSignal(const double wS)
    {
      ws = wS;
    }
 
 
    /**
     * Add signal strength.
     *
     * \param  wS        signal strength
     */
    void addSignal(const double wS)
    {
      ws += wS;
    }
 
  protected:
 
    double ws  =  0.0;                                     //!< total signal strength
  };
}
 
#endif