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/JAstronomyToolkit.hh"
 
#include "JLang/JVectorize.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JRandom.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;
  histogram_type X;
  data_type      Q;
  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['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")                             = 0.0;
    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);
 
 
  for (const auto& i : setup) {
 
    const TObject* ps  =  getObject(i.Hs);
    const TObject* pb  =  getObject(i.Hb);
 
    STATUS(printer("Signal:",     ps) << endl);
    STATUS(printer("Background:", pb) << endl);
 
    px.add(ps, pb);
  }
 
  px.set(Fs, Fb);
 
  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 > 0.0 || Q > 0.0) {
 
    JQuantile<data_type> storage(400.0);
 
    {
      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("Average time:   " << setw(6) << (t1 - t0) / chrono::nanoseconds(numberOfTests) << " [ns]" << endl);
    }
 
    for (const auto& result : storage) {
      for (const auto x : result.second) {
        hl.Fill(x);
      }
    }
 
    const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
    if (P > 0.0) {                // determine minimal likelihood ratio given maximal probability
 
      const double x = storage.getValue(P);
      const double p = storage.getProbability(x);
 
      cout << "Minimal likelihood ratio: " << FIXED(9,5) << x << ' ' << SCIENTIFIC(12,3) << p << endl;
    }
 
    if (Q > 0.0) {                // determine maximal probability given minimal likelihood ratio
 
      const double p = storage.getProbability(Q);
      const double x = storage.getValue(p);
 
      cout << "Maximal probability:      " << SCIENTIFIC(12,3) << p << ' ' << FIXED(9,5) << x << endl;
    }
 
    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);
 
  } 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();
}