JPseudoExperiment.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JPSEUDOEXPERIMENT__
#define __JASTRONOMY__JPSEUDOEXPERIMENT__
 
#include <istream>
#include <ostream>
#include <vector>
#include <algorithm>
#include <limits>
 
#include "TObject.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TRandom3.h"
 
#include "JAstronomy/JExperiment.hh"
#include "JAstronomy/JAspera.hh"
#include "JAstronomy/JNuisance.hh"
 
 
/**
 * \file
 *
 * Pseudo experiment.
 * \author mdejong
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  
  /**
   * Auxiliary interface for pseudo experiment.
   */
  struct JPseudoExperiment_t :
    public JExperiment
  {
    /**
     * Statistics of pseudo experiment.
     */
    struct stats_type {
      /**
       * Addition operator.
       *
       * \param  px              statistics of pseudo experiment
       * \return                 this statistics
       */
      stats_type& operator+=(const stats_type& px)
      {
        this->ns += px.ns;
        this->nb += px.nb;
 
        return *this;
      }
 
      size_t ns  =  0;                              //!< number of generated signal events
      size_t nb  =  0;                              //!< number of generated background events
    };
 
 
    typedef JAspera::fit_type  fit_type;            //!< fit type
 
 
    /**
     * Combined result of pseudo experiment and fit.
     */
    struct result_type :
      public stats_type,
      public fit_type
    {};
 
 
    /**
     * Set scaling factors of signal and background strengths.
     *
     * \param  fS        signal strength
     * \param  fB        background strength
     */
    virtual void set(const double fS, const double fB = 1.0) = 0;
 
 
    /**
     * Get fit method.
     *
     * \return            fit
     */
    virtual JAspera& getAspera() = 0;
 
 
    /**
     * Generate pseudo experiment and transfer S/N values to fit method.
     *
     * \param  out        output
     * \return            result
     */
    virtual stats_type run(JAspera& out) const = 0;
 
 
    /**
     * Generate background only pseudo experiment and transfer S/N values to fit method.
     *
     * \param  out        output
     * \param  nb         number of background events
     * \return            result
     */
    virtual stats_type run(JAspera& out, const size_t nb) const = 0;
 
 
    /**
     * Generate pseudo experiment and fit signal strength.
     *
     * \return            result
     */
    inline result_type operator()()
    {
      JAspera& aspera = getAspera();
 
      // reset
 
      aspera.clear();
      aspera.setSignal(0.0);
 
      return { run(aspera), aspera() };
    }
 
 
    /**
     * Generate background only pseudo experiment and fit signal strength.
     *
     * \param  nb         number of background events
     * \return            result
     */
    inline result_type operator()(const size_t nb)
    {
      JAspera& aspera = getAspera();
 
      // reset
 
      aspera.clear();
      aspera.setSignal(0.0);
 
      return { run(aspera, nb), aspera() };
    }
 
 
    /**
     * Run pseudo experiments using given storage.
     *
     * \param  storage             storage
     */
    inline void operator()(std::vector<result_type>& storage)
    {
      for (auto& i : storage) {
        i = (*this)();
      }
    }
 
 
    /**
     * Run pseudo experiments using given storage.
     *
     * \param  pm                  pointer to data member of result
     * \param  storage             storage
     */
    template<class T, class JValue_t> 
    inline void operator()(JValue_t result_type::*pm, std::vector<T>& storage)
    {
      for (auto& i : storage) {
        i = (*this)().*pm;
      }
    }
 
 
    /**
     * Run pseudo experiments using given storage.
     *
     * \param  pm                  pointer to data member of result
     * \param  storage             storage
     */
    template<class T, class JValue_t> 
    inline void operator()(JValue_t JAspera::fit_type::*pm, std::vector<T>& storage)
    {
      (*this)(static_cast<JValue_t result_type::*>(pm), storage);
    }
 
 
    static constexpr double MINIMAL_SIGNAL_STRENGTH  =  1.0e-10;          // minimal signal strength for upper limit
    static constexpr double MAXIMAL_SIGNAL_STRENGTH  =  1.0e+10;          // maximal signal strength for upper limit
 
 
    /**
     * Get signal strength given result of experiment and probability of upper limit.
     *
     * \param  aspera         result of experiment
     * \param  Q              probability
     * \param  numberOfTests  number of tests
     * \param  precision      precision
     * \return                signal strength
     */
    double getSignalStrengthForUpperLimit(const JAspera& aspera,
                                          const double   Q,
                                          const size_t   numberOfTests,
                                          const double   precision = 1.0e-4)
    {
      using namespace std;
 
      const JAspera::fit_type result = aspera(true);
 
      size_t n = 0;
 
      this->set(0.0);
 
      for (size_t i = 0; i != numberOfTests; ++i) {
 
        JAspera aspera;
 
        this->run(aspera);
 
        if (aspera(true).signal < result.signal) {
          n += 1;
        }
      }
 
      if (n > (1.0 - Q) * numberOfTests) {
 
        double mumin = 1.0;
        double mumax = 1.0;
 
        {
          for ( ; ; mumax *= 2.0) {
 
            const double ts = aspera.getTestStatisticForUpperLimit(mumax);
 
            this->set(mumax);
 
            const double ps = this->getProbabilityForUpperLimit(mumax, ts, numberOfTests);
 
            if (ps < 1.0 - Q || mumax > MAXIMAL_SIGNAL_STRENGTH) {
              break;
            }
          }
        }
 
        if (mumax == 1.0) {
 
          mumin = 0.5*mumax;
 
          for ( ; ; mumin *= 0.5) {
 
            const double ts = aspera.getTestStatisticForUpperLimit(mumin);
 
            this->set(mumin);
 
            const double ps = this->getProbabilityForUpperLimit(mumin, ts, numberOfTests);
 
            if (ps > 1.0 - Q || mumin < MINIMAL_SIGNAL_STRENGTH) {
              break;
            }
          }
 
          mumax = 2.0*mumin;
 
        } else {
 
          mumin = 0.5*mumax;
        }
 
        if      (mumin < MINIMAL_SIGNAL_STRENGTH)
          return mumin;
        else if (mumax > MAXIMAL_SIGNAL_STRENGTH)
          return mumax;
        else {
 
          for ( ; ; ) {                               // binary search
 
            const double mu = 0.5 * (mumin + mumax);
 
            const double ts = aspera.getTestStatisticForUpperLimit(mu);
 
            this->set(mu);
 
            const double ps = this->getProbabilityForUpperLimit(mu, ts, numberOfTests);
 
            if (fabs(ps - (1.0 - Q)) <= precision) {
              return mu;
            }
 
            if (ps >= 1.0 - Q)
              mumin = mu;
            else
              mumax = mu;
 
            // we have slight overcoverage and can't reach requested CL precisely
 
            if (ps <= 1.0 - Q && mumax - mumin < precision) {
              return mu; 
            }
          }
        }
      }
 
      return 0.0;
    }
 
  protected:
    /**
     * Get probability for given pseudo experiment and signal strength to exceed minimal test statistic for upper limit.
     *
     * \param  ps             signal strength
     * \param  ts             test statistic
     * \param  nx             number of pseudo experiments
     */
    inline double getProbabilityForUpperLimit(const double ps,
                                              const double ts,
                                              const size_t nx) const
    {
      size_t ns = 0;
 
      for (size_t i = 0; i != nx; ++i) {
 
        JAspera aspera;
 
        this->run(aspera);
 
        if (ps <= std::numeric_limits<double>::min()) {
 
          if (aspera().signal <= ps) {
            ns += 1;
          }
 
        } else if (aspera.getTestStatisticForUpperLimit(ps) > ts) {
          ns += 1;
        }
      }
 
      return (double) ns / (double) nx;
    }
  };
 
 
  /**
   * Pseudo experiment using CDF for combined generation and likelihood evaluation.
   */
  struct JPseudoExperiment :
    public JPseudoExperiment_t
  {
    using JPseudoExperiment_t::operator();
 
    /**
     * Default constructor.
     */
    JPseudoExperiment()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  hs         histogram with PDF of signal
     * \param  hb         histogram with PDF of background
     */
    template<class H_t>
    JPseudoExperiment(const H_t&   hs,
                      const H_t&   hb)
    {
      add(hs, hb);
    }
 
 
    /**
     * Constructor.
     *
     * \param  hS         histogram with PDF for generation of signal
     * \param  hB         histogram with PDF for generation of background
     * \param  hs         histogram with PDF for evaluation of signal
     * \param  hb         histogram with PDF for evaluation of background
     */
    template<class H_t>
    JPseudoExperiment(const H_t&   hS,
                      const H_t&   hB,
                      const H_t&   hs,
                      const H_t&   hb)
    {
      add(hS, hB, hs, hb);
    }
 
 
    /**
     * Add objects with PDFs of signal and background.
     *
     * \param  ps         pointer to object with PDF of signal
     * \param  pb         pointer to object with PDF of background
     */
    void add(const TObject* ps,
             const TObject* pb)
    {
      if (add<TH3>(ps, pb)) { return; }
      if (add<TH2>(ps, pb)) { return; }
      if (add<TH1>(ps, pb)) { return; }
    }
 
 
    /**
     * Add objects with PDFs of signal and background.
     *
     * \param  pS         pointer to object with PDF for generation of signal
     * \param  pB         pointer to object with PDF for generation of background
     * \param  ps         pointer to object with PDF for evaluation of signal
     * \param  pb         pointer to object with PDF for evaluation of background
     */
    void add(const TObject* pS,
             const TObject* pB,
             const TObject* ps,
             const TObject* pb)
    {
      if (add<TH3>(pS, pB, ps, pb)) { return; }
      if (add<TH2>(pS, pB, ps, pb)) { return; }
      if (add<TH1>(pS, pB, ps, pb)) { return; }
    }
 
 
    /**
     * Add histograms with PDFs of signal and background.
     *
     * \param  hs         histogram with PDF of signal
     * \param  hb         histogram with PDF of background
     */
    void add(const TH1&   hs,
             const TH1&   hb)
    {
      for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
        add(hs.GetBinContent(ix),
            hb.GetBinContent(ix));
      }
    }
 
 
    /**
     * Add histograms with PDFs of signal and background.
     *
     * \param  hS         histogram with PDF for generation of signal
     * \param  hB         histogram with PDF for generation of background
     * \param  hs         histogram with PDF for evaluation of signal
     * \param  hb         histogram with PDF for evaluation of background
     */
    void add(const TH1&   hS,
             const TH1&   hB,
             const TH1&   hs,
             const TH1&   hb)
    {
      for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
        add(hS.GetBinContent(ix),
            hB.GetBinContent(ix),
            hs.GetBinContent(ix),
            hb.GetBinContent(ix));
      }
    }
 
 
    /**
     * Add histograms with PDFs of signal and background.
     *
     * \param  hs         histogram with PDF of signal
     * \param  hb         histogram with PDF of background
     */
    void add(const TH2&   hs,
             const TH2&   hb)
    {
      for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
        for (Int_t iy = 1; iy <= hs.GetYaxis()->GetNbins(); ++iy) {
          add(hs.GetBinContent(ix, iy),
              hb.GetBinContent(ix, iy));
        }
      }
    }
 
 
    /**
     * Add histograms with PDFs of signal and background.
     *
     * \param  hS         histogram with PDF for generation of signal
     * \param  hB         histogram with PDF for generation of background
     * \param  hs         histogram with PDF for evaluation of signal
     * \param  hb         histogram with PDF for evaluation of background
     */
    void add(const TH2&   hS,
             const TH2&   hB,
             const TH2&   hs,
             const TH2&   hb)
    {
      for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
        for (Int_t iy = 1; iy <= hs.GetYaxis()->GetNbins(); ++iy) {
          add(hS.GetBinContent(ix, iy),
              hB.GetBinContent(ix, iy),
              hs.GetBinContent(ix, iy),
              hb.GetBinContent(ix, iy));
        }
      }
    }
 
 
    /**
     * Add histograms with PDFs of signal and background.
     *
     * \param  hs         histogram with PDF of signal
     * \param  hb         histogram with PDF of background
     */
    void add(const TH3&   hs,
             const TH3&   hb)
    {
      for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
        for (Int_t iy = 1; iy <= hs.GetYaxis()->GetNbins(); ++iy) {
          for (Int_t iz = 1; iz <= hs.GetZaxis()->GetNbins(); ++iz) {
            add(hs.GetBinContent(ix, iy, iz),
                hb.GetBinContent(ix, iy, iz));
          }
        }
      }
    }
 
 
    /**
     * Add histograms with PDFs of signal and background.
     *
     * \param  hS         histogram with PDF for generation of signal
     * \param  hB         histogram with PDF for generation of background
     * \param  hs         histogram with PDF for evaluation of signal
     * \param  hb         histogram with PDF for evaluation of background
     */
    void add(const TH3&   hS,
             const TH3&   hB,
             const TH3&   hs,
             const TH3&   hb)
    {
      for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
        for (Int_t iy = 1; iy <= hs.GetYaxis()->GetNbins(); ++iy) {
          for (Int_t iz = 1; iz <= hs.GetZaxis()->GetNbins(); ++iz) {
            add(hS.GetBinContent(ix, iy, iz),
                hB.GetBinContent(ix, iy, iz),
                hs.GetBinContent(ix, iy, iz),
                hb.GetBinContent(ix, iy, iz));
          }
        }
      }
    }
 
 
    /**
     * Add signal and background.
     *
     * \param  s          signal
     * \param  b          background
     */
    void add(const double s,
             const double b)
    {
      if (check(s, b)) {
 
        cs.put(s);
        cb.put(b);
 
        aspera.put(s, b);
 
      } else {
 
        remnant.add(s, b, s, b);
      }
    }
 
 
    /**
     * Add signal and background.
     *
     * \param  S          signal for generation
     * \param  B          background for generation
     * \param  s          signal for evaluation
     * \param  b          background for evaluation
     */
    void add(const double S,
             const double B,
             const double s,
             const double b)
    {
      if (//check(S, B) &&
          check(s, b)) {
 
        cs.put(S);
        cb.put(B);
 
        aspera.put(s, b);
 
      } else {
 
        remnant.add(S, B, s, b);
      }
    }
 
 
    /**
     * Add remnant signal and background.
     */
    void add()
    {
      if (remnant.n != 0) {
 
        add(remnant.S, remnant.B, remnant.s, remnant.b);
 
        remnant.reset();
      }
    }
 
 
    /**
     * Get total signal.
     *
     * \return           signal
     */
    double getSignal() const
    {
      return cs.back() * this->fs;
    }
 
 
    /**
     * Get total background.
     *
     * \return           background
     */
    double getBackground() const
    {
      return cb.back() * this->fb;
    }
 
 
    /**
     * Set scaling factors of signal and background strengths.
     *
     * \param  fS        signal strength
     * \param  fB        background strength
     */
    virtual void set(const double fS, const double fB = 1.0) override
    {
      for (auto& i : aspera) {
        i *= fB / this->fb;
      }
      
      this->fs = fS;
      this->fb = fB;
    }
 
 
    /**
     * Get fit method.
     *
     * \return            fit
     */
    virtual JAspera& getAspera() override
    {
      return this->fit;
    }
 
 
    /**
     * Generate pseudo experiment and transfer values to fit method.
     *
     * \param  out        output
     * \return            result
     */
    virtual stats_type run(JAspera& out) const
    {
      out.addSignal(aspera.getSignal());
 
      const size_t ns = gRandom->Poisson(getSignal()     * nuisance.signal    ->get());           // number of signal events
      const size_t nb = gRandom->Poisson(getBackground() * nuisance.background->get());           // number of background events
 
      for (size_t i = ns; i != 0; --i) { out.push_back(aspera[cs.get_index(gRandom->Rndm())]); }  // store distributed signal events
      for (size_t i = nb; i != 0; --i) { out.push_back(aspera[cb.get_index(gRandom->Rndm())]); }  // store distributed background events
 
      return { ns, nb };
    }
 
 
    /**
     * Generate background only pseudo experiment and transfer S/N values to fit method.
     *
     * \param  out        output
     * \param  nb         number of background events
     * \return            result
     */
    virtual stats_type run(JAspera& out, const size_t nb) const
    {
      out.addSignal(aspera.getSignal());
 
      for (size_t i = nb; i != 0; --i) { out.push_back(aspera[cb.get_index(gRandom->Rndm())]); }  // store distributed background events
 
      return { 0, nb };
    }
 
 
    /**
     * Nuisance parameters.
     */
    struct parameters_type {
 
      /**
       * Read parameters from input stream.
       *
       * \param  in                 input stream
       * \param  parameters         parameters
       * \return                    input stream
       */
      friend inline std::istream& operator>>(std::istream& in, parameters_type& parameters)
      {
        return in >> parameters.signal
                  >> parameters.background;
      }
 
 
      /**
       * Write parameters to output stream.
       *
       * \param  out                output stream
       * \param  parameters         parameters
       * \return                    output stream
       */
      friend inline std::ostream& operator<<(std::ostream& out, const parameters_type& parameters)
      {
        return out << parameters.signal     << ' '
                   << parameters.background;
      }
 
      nuisance_type  signal     = std::make_shared<JNuisanceFixed>();
      nuisance_type  background = std::make_shared<JNuisanceFixed>();
 
    } nuisance;
 
 
    /**
     * Configure lookup tables.
     *
     * \param  N                  number of bins
     */
    void configure(size_t N)
    {
      cs.configure(N);
      cb.configure(N);
    }
 
  protected:
    /**
     * Add objects with PDF of signal and background.
     *
     * \param  ps         pointer to object with PDF of signal
     * \param  pb         pointer to object with PDF of background
     * \return            true if added; else false
     */
    template<class H_t>
    bool add(const TObject* ps,
             const TObject* pb)
    {
      if (dynamic_cast<const H_t*>(ps) != NULL &&
          dynamic_cast<const H_t*>(pb) != NULL) {
 
        const H_t& hs = dynamic_cast<const H_t&>(*ps);
        const H_t& hb = dynamic_cast<const H_t&>(*pb);
 
        if (check(hs, hb)) {
 
          add(hs, hb);
 
          return true;
        }
      }
 
      return false;
    }
 
 
    /**
     * Add objects with PDF of signal and background.
     *
     * \param  pS         pointer to object with PDF for generation of signal
     * \param  pB         pointer to object with PDF for generation of background
     * \param  ps         pointer to object with PDF for evaluation of signal
     * \param  pb         pointer to object with PDF for evaluation of background
     * \return            true if added; else false
     */
    template<class H_t>
    bool add(const TObject* pS,
             const TObject* pB,
             const TObject* ps,
             const TObject* pb)
    {
      if (dynamic_cast<const H_t*>(pS) != NULL &&
          dynamic_cast<const H_t*>(pB) != NULL &&
          dynamic_cast<const H_t*>(ps) != NULL &&
          dynamic_cast<const H_t*>(pb) != NULL) {
 
        const H_t& hS = dynamic_cast<const H_t&>(*pS);
        const H_t& hB = dynamic_cast<const H_t&>(*pB);
        const H_t& hs = dynamic_cast<const H_t&>(*ps);
        const H_t& hb = dynamic_cast<const H_t&>(*pb);
 
        if (check(hS, hB) && check(hB, hs) && check(hs, hb)) {
 
          add(hS, hB, hs, hb);
 
          return true;
        }
      }
 
      return false;
    }
 
 
    /**
     * Auxiliary data structure for CDF.
     */
    struct cdf_type :
      public std::vector<double>
    {
      /**
       * Configure lookup table to substitute binary search.
       *
       * \param  N                  number of bins
       */
      void configure(size_t N)
      {
        if (this->size() > 1) {
 
          index.resize(N);
 
          size_t p = 0;
 
          for (size_t i = 0; i != N; ++i) {
 
            const double x = this->back() * (double) i / (double) N;
 
            while ((*this)[p+1] < x) {
              ++p;
            }
 
            index[i] = p;
          }
        }
      }
 
      /**
       * Get index corresponding to given random value.
       *
       * \param  rv                 random value [0,1>
       * \param  option             option to force use of binary search
       * \return                    index
       */
      inline size_t get_index(const double rv, const bool option = false) const
      {
        const double value = this->back() * rv;
 
        if (option || index.empty()) {
 
          size_t first = 0;
          size_t count = this->size();
 
          for (size_t i, step; count != 0; ) {
 
            step = count / 2;
            i    = first + step;
 
            if ((*this)[i] < value) {
              first  = ++i;
              count -= step + 1;
            } else {
              count  = step;
            }
          }
 
          return first;
 
        } else {
 
          size_t i = index[rv * index.size()];
 
          while (i != this->size() && (*this)[i] < value) {
            ++i;
          }
 
          return i;
        }
      }
 
 
      /**
       * Put given value.
       *
       * \param  x                  value
       */
      void put(const double x)
      {
        if (this->empty())
          this->push_back(x);
        else
          this->push_back(this->back() + x);
      }
 
    private:
      std::vector<size_t> index;                             //!< lookup table
    };
 
 
    cdf_type cs;                                             //!< CDF of signal
    cdf_type cb;                                             //!< CDF of background
 
    JAspera  aspera;                                         //!< pre-computed N/S values
 
    double   fs  =  1.0;                                     //!< scaling factor signal strength
    double   fb  =  1.0;                                     //!< scaling factor background strength
 
    JAspera  fit;                                            //!< fit
 
    struct remnant_type {
 
      void add(const double S, const double B, const double s, const double b)
      {
        this->n += 1;
        this->S += S;
        this->B += B;
        this->s += s;
        this->b += b;
      }
 
      void reset()
      {
        this->n = 0;
        this->S = 0.0;
        this->B = 0.0;
        this->s = 0.0;
        this->b = 0.0;
      }
 
      size_t n = 0;
      double S = 0.0;
      double B = 0.0;
      double s = 0.0;
      double b = 0.0;
 
    } remnant;
  };
}
 
#endif