JAstronomyToolkit.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JASTRONOMYTOOLKIT__
#define __JASTRONOMY__JASTRONOMYTOOLKIT__
 
#include "JAstronomy/JPseudoExperiment.hh"
#include "JAstronomy/JGarage.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
 
  /**
   * Auxiliary data structure for discovery calculations.
   */
  template<class T, const size_t factor>
  struct JDiscovery :
    public JGarage<T, factor>
  {
    using JGarage<T, factor>::getProbability;
 
    /**
     * Constructor.
     *
     * \param  px             pseudo experiment
     */
    JDiscovery(JPseudoExperiment_t& px) :
      px(px)
    {}
 
 
    /**
     * Get signal strength for which given likelihood is exceeded with given probability.
     *
     * \param  Q              likelihood
     * \param  P              probability [0,1>
     * \param  numberOfTests  number of tests
     * \param  precision      precision of signal strength
     * \return                signal strength
     */
    inline double operator()(const double Q,
                             const double P,
                             const size_t numberOfTests,
                             const double precision = 1.0e-3)
    {
      double xmin = 1.0;
      double xmax = 2.0;
 
      while (getProbability(xmin, Q, numberOfTests) > P) {      // set lower limit
        xmax  = xmin;
        xmin *= 0.5;
      }
 
      while (getProbability(xmax, Q, numberOfTests) < P) {      // set upper limit
        xmin  = xmax;
        xmax *= 2.0;
      }
 
      while (xmax - xmin > precision) {                         // bisection
 
        const double x = 0.5 * (xmin + xmax);
 
        if (getProbability(x, Q, numberOfTests) < P)
          xmin = x;
        else
          xmax = x;
      }
 
      return 0.5 * (xmin + xmax);
    }
 
  protected:
    /**
     * Run pseudo experiments to get probability corresponding to given minimal likelihood.
     *
     * \param  fS             signal strength
     * \param  Q              likelihood
     * \param  numberOfTests  number of tests
     * \return                probability
     */
    inline double getProbability(const double fS,
                                 const double Q,
                                 const size_t numberOfTests)
    {
      px.set(fS);
 
      this->clear();
 
      for (size_t i = 0; i != numberOfTests; ++i) {
        this->put(px().likelihood);
      }
 
      return this->getProbability(Q);
    }
    
    JPseudoExperiment_t& px;
  };
 
 
  /**
   * Run pseudo experiments to determine signal strength for which given likelihood is exceeded with given probability.
   *
   * \param  px             pseudo experiment
   * \param  Q              likelihood
   * \param  P              probability [0,1>
   * \param  numberOfTests  number of tests
   * \param  precision      precision of signal strength
   * \return                signal strength
   */
  template<class T, const size_t factor>
  inline double discovery(JPseudoExperiment_t& px,
                          const double Q,
                          const double P,
                          const size_t numberOfTests,
                          const double precision = 1.0e-3)
  {
    JDiscovery<T, factor> discovery(px);
 
    return discovery(Q, P, numberOfTests, precision);
  }
}
 
#endif