JGen2.hh

Go to the documentation of this file.

#ifndef __JASTRONOMY__JGEN2__
#define __JASTRONOMY__JGEN2__
 
#include <vector>
#include <numeric>
 
#include "JAstronomy/JPseudoExperiment.hh"
#include "JAstronomy/JAspera.hh"
 
 
/**
 * \file
 *
 * Set of pseudo experments.
 * \author mdejong
 */
namespace JASTRONOMY {}
namespace JPP { using namespace JASTRONOMY; }
 
namespace JASTRONOMY {
 
  /**
   * Auxiliary data structure to fit signal strength using likelihood ratio for multiple pseudo experiments.
   */
  struct JGen2 :
    public std::vector<JPseudoExperiment>,
    public JPseudoExperiment_t
  {
    /**
     * Add remnant signal and background.
     */
    void add()
    {
      for (auto& i : static_cast<std::vector<JPseudoExperiment>&>(*this)) {
        i.add();
      }
    }
 
    /**
     * Configure lookup tables.
     *
     * \param  N                  number of bins
     */
    void configure(size_t N)
    {
      for (auto& i : static_cast<std::vector<JPseudoExperiment>&>(*this)) {
        i.configure(N);
      }
    }
 
 
    /**
     * Get total signal.
     *
     * \return           signal
     */
    double getSignal() const
    {
      return std::accumulate(this->begin(), this->end(), 0.0, [](const double y, const JPseudoExperiment& px) { return y + px.getSignal(); } );
    }
 
 
    /**
     * Get total background.
     *
     * \return           background
     */
    double getBackground() const
    {
      return std::accumulate(this->begin(), this->end(), 0.0, [](const double y, const JPseudoExperiment& px) { return y + px.getBackground(); } );
    }
 
 
    /**
     * 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 : static_cast<std::vector<JPseudoExperiment>&>(*this)) {
        i.set(fS, fB);
      }
    }
 
 
    /**
     * Get fit method.
     *
     * \return            fit
     */
    virtual JAspera& getAspera() override
    {
      return this->fit;
    }
 
 
    /**
     * Generate pseudo experiment and transfer S/N values to fit method.
     *
     * \param  out        output
     * \return            result
     */
    virtual stats_type run(JAspera& out) const
    {
      stats_type ps = { 0, 0 };
      
      for (const auto& px : static_cast<const std::vector<JPseudoExperiment>&>(*this)) {
        ps += px.run(out);
      }
 
      return ps;
    }
 
 
    /**
     * 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
    {
      stats_type ps = { 0, 0 };
 
      const double B = this->getBackground();
 
      for (size_t i = 0; i != nb; ++i) {
 
        const double X = gRandom->Uniform(0.0, B);
 
        double x = 0.0;
 
        for (const auto& px : static_cast<const std::vector<JPseudoExperiment>&>(*this)) {
 
          x += px.getBackground();
 
          if (x >= X) {
 
            ps += px.run(out, 1);
 
            break;
          }
        }
      }
 
      return ps;
    }
 
  protected:
    JAspera  fit;                                            //!< fit
  };
}
 
#endif