JPseudoExperiment.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <chrono>
#include <vector>
#include <map>
#include <exception>
#include <algorithm>
 
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TRandom3.h"
 
#include "JAstronomy/JPseudoExperiment.hh"
#include "JAstronomy/JGarage.hh"
#include "JAstronomy/JAstronomyToolkit.hh"
 
#include "JLang/JVectorize.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JRandom.hh"
#include "JMath/JMathSupportkit.hh"
 
#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
 
 
namespace {
 
  /**
   * Auxiliary data structure for experiment.
   */
  struct experiment_type {
 
    JGIZMO::JRootObjectID Hs;         //!< signal
    JGIZMO::JRootObjectID Hb;         //!< background
 
    /**
     * Read experiment from input stream.
     *
     * \param  in             input stream
     * \param  experiment     experiment
     * \return                input stream
     */
    friend inline std::istream& operator>>(std::istream& in, experiment_type& experiment)
    {
      return in >> experiment.Hs
                >> experiment.Hb;
    }
 
 
    /**
     * Write experiment to output stream.
     *
     * \param  out            output stream
     * \param  experiment     experiment
     * \return                output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const experiment_type& experiment)
    {
      return out << experiment.Hs << ' '
                 << experiment.Hb;
    }
  };
 
 
  /**
   * Auxiliary data structure for printing.
   */
  struct printer {
    /**
     * Constructor.
     *
     * \param  title          title
     * \param  ps             pointer to object
     */
    printer(const char* const title,
            const TObject*    ps) :
      title(title),
      ps(ps)
    {}
 
    /**
     * Write printer to output stream.
     *
     * \param  out            output stream
     * \param  printer        printer
     * \return                output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const printer& printer)
    {
      using namespace std;
      using namespace JPP;
 
      out << setw(16) << left << printer.title << right;
 
      if (printer.ps != NULL) {
 
        out << ' ' << setw(16) << left << printer.ps->GetName() << right;
 
        const TH1* h1 = dynamic_cast<const TH1*>(printer.ps);
 
        if (h1 != NULL) {
          out << ' ' << FIXED(10,3) << h1->GetSumOfWeights();
        }
      }
 
      return out;
    }
 
  private:
    const char* const title;
    const TObject*    ps;
  };
}
 
 
/**
 * \file
 * Test application for pseudo experiment generation and likelihood ratio evaluations.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  typedef JAbstractHistogram<double>   histogram_type;
  typedef float                        data_type;
 
  JPseudoExperiment px;
 
  JContainer< vector<experiment_type> >  setup;
  string         outputFile;
  size_t         numberOfTests;
  double         Fs;
  double         Fb;
  double         boost;
  size_t         M;
  double         SNR;
  histogram_type X;
  data_type      Q;
  vector<double> P;
  JRandom        seed;
  int            debug;
 
  try { 
 
    JParser<> zap;
 
    zap['E'] = make_field(setup,
                          "douplets of signal and background histograms, "
                          << "each of which defined by <file name>:<histogram name>");
    zap['o'] = make_field(outputFile,    "output file with histograms")                     = "benchmark.root";
    zap['n'] = make_field(numberOfTests, "number of tests");
    zap['s'] = make_field(Fs,            "signal strength");
    zap['b'] = make_field(Fb,            "background strength")                             = 1.0;
    zap['B'] = make_field(boost,         "boost livetime of experiment")                    = 1.0;
    zap['M'] = make_field(M,             "lookup table for CDFs")                           = 0;
    zap['R'] = make_field(SNR,           "signal-to-noise ratio")                           = 0.0;
    zap['N'] = make_field(px.nuisance,
                          "signal and background nuisances, "
                          << "each of which defined by <type> (values), \n"
                          << "\twhere <type> can be:" 
                          << get_keys(px.nuisance.signal.getDictionary()))                  = JPARSER::initialised();
    zap['x'] = make_field(X,             "x-axis likelihood histogram")                     = histogram_type(110, -10.0, +100.0);
    zap['Q'] = make_field(Q,             "minimal likelihood ratio")                        = 0.0;
    zap['P'] = make_field(P,             "maximal probability")                             = JPARSER::initialised();
    zap['S'] = make_field(seed)        = 0;
    zap['d'] = make_field(debug)       = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  seed.set(gRandom);
 
  JExperiment::setSNR(SNR);
 
  for (const auto& i : setup) {
 
    TObject* ps  =  getObject(i.Hs);
    TObject* pb  =  getObject(i.Hb);
 
    dynamic_cast<TH1*>(ps)->Scale(boost);
    dynamic_cast<TH1*>(pb)->Scale(boost);
 
    STATUS(printer("Signal:",     ps) << endl);
    STATUS(printer("Background:", pb) << endl);
 
    px.add(ps, pb);
  }
 
  px.set(Fs, Fb);
 
  if (M != 0) {
    px.configure(M);
  }
 
  STATUS("Signal:          " << FIXED(10,3) << px.getSignal()     << endl);
  STATUS("Background:      " << FIXED(10,3) << px.getBackground() << endl);
  STATUS("Number of tests: " << setw(10)    << numberOfTests      << endl);
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  out << *gRandom;
 
  TH1D hl("hl", NULL, X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit());
  TH1D hn("hn", NULL, 100, -15.0, +15.0);
 
 
  if (P.size() > 1) {                                     // special case
 
    JGarage<data_type, 100> storage;
 
    double Ps = 0.0;                                      // total probability
    size_t nb = px.getBackground() + 1;                   // start value for number of background events
    size_t mb = 0;
 
    {
      const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
      for ( ; ; ++nb) {
 
        const double p = poisson(nb, px.getBackground());
        const size_t n = p * numberOfTests;
 
        if (p > P[1]) {
          continue;
        }
 
        if (n == 0) {
          break;
        }
 
        if (mb == 0) {
          mb = nb;
 
        }
        for (size_t i = 0; i != n; ++i) {
          storage.put(px(nb).likelihood);
        }
 
        Ps += p;
      }
 
      const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
      STATUS("Total time:   "   << setw(8) << (t1 - t0) / chrono::milliseconds(1)            << " [ms]" << endl);
      STATUS("Average time:   " << setw(6) << (t1 - t0) / chrono::nanoseconds(numberOfTests) << " [ns]" << endl);
    }
 
    STATUS("Poisson: " << setw(3) << mb << ' ' << SCIENTIFIC(12,3) << Ps << ' ' << setw(8) << storage.size() << endl);
 
    storage([&hl](const data_type x) { hl.Fill(x); });
 
    double x = 0.0;
    double p = 0.0;
 
    {
      const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
      x = storage.getValue(P[0]/Ps);
      p = storage.getProbability(x) * Ps;
 
      const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
      STATUS("Time to result: " << setw(6) << (t1 - t0) / chrono::nanoseconds(1) << " [ns]" << endl);
    }
 
    cout << "Minimal likelihood ratio: " << FIXED(9,5) << x << ' ' << SCIENTIFIC(12,3) << p << endl;
 
  } else if (P.size() > 0 && Q > 0.0) {
 
    const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
    const double signal = discovery<float, 100>(px, Q, P[0], numberOfTests);
 
    const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
    STATUS("Total time:   "   << setw(8) << (t1 - t0) / chrono::milliseconds(1)            << " [ms]" << endl);
 
    cout << "Signal strength: " << FIXED(9,5) << signal << endl;
 
  } else if (P.size() > 0 || Q > 0.0) {
 
    JGarage<data_type, 100> storage;
 
    {
      const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
      for (size_t i = 0; i != numberOfTests; ++i) {
        storage.put(px().likelihood);
      }
 
      const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
      STATUS("Total time:   "   << setw(8) << (t1 - t0) / chrono::milliseconds(1)            << " [ms]" << endl);
      STATUS("Average time:   " << setw(6) << (t1 - t0) / chrono::nanoseconds(numberOfTests) << " [ns]" << endl);
    }
 
    storage([&hl](const data_type x) { hl.Fill(x); });
 
    double x = 0.0;
    double p = 0.0;
 
    {
      const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
      if (P.size() == 1) {          // determine minimal likelihood ratio given maximal probability
        x = storage.getValue(P[0]);
        p = storage.getProbability(x);
      }
 
      if (Q > 0.0) {                // determine maximal probability given minimal likelihood ratio
        p = storage.getProbability(Q);
        x = storage.getValue(p);
      }
 
      const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
      STATUS("Time to result: " << setw(6) << (t1 - t0) / chrono::microseconds(1) << " [us]" << endl);
    }
 
    if (P.size() == 1) { cout << "Minimal likelihood ratio: " << FIXED(9,5) << x << ' ' << SCIENTIFIC(12,3) << p << endl; }
    if (Q > 0.0)       { cout << "Maximal probability:      " << SCIENTIFIC(12,3) << p << ' ' << FIXED(9,5) << x << endl; }
 
  } else {
 
    std::vector<JPseudoExperiment::result_type> storage(numberOfTests);
 
    const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
    px(storage);
      
    const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
    STATUS("Average time:   " << setw(6) << (t1 - t0) / chrono::nanoseconds(storage.size()) << " [ns]" << endl);
 
    for (const auto& result : storage) {
 
      DEBUG("result: "
            << setw(3)     << result.ns          << ' '
            << setw(3)     << result.nb          << ' '
            << FIXED(12,5) << result.likelihood  << ' '
            << FIXED(12,5) << result.signal      << endl);
 
      hl.Fill(result.likelihood);
      hn.Fill(result.signal - result.ns);
    }
  }
 
  out.Write();
  out.Close();
}