JGen2.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/JGen2.hh"
#include "JAstronomy/JGarage.hh"
#include "JAstronomy/JAstronomyToolkit.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 for generation
    JGIZMO::JRootObjectID HB;         //!< background for generation
    JGIZMO::JRootObjectID Hs;         //!< signal for evaluation
    JGIZMO::JRootObjectID Hb;         //!< background for evaluation
 
    JASTRONOMY::JPseudoExperiment::parameters_type nuisance;
 
    /**
     * 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
                >> experiment.Hs
                >> experiment.Hb
                >> experiment.nuisance;
    }
 
 
    /**
     * 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 << ' '
                 << experiment.Hs << ' '
                 << experiment.Hb << ' '
                 << experiment.nuisance;
    }
  };
 
 
  /**
   * 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;
 
  JContainer< vector<experiment_type> >  setup;
  string         outputFile;
  size_t         numberOfTests;
  double         Fs;
  double         Fb;
  size_t         M;
  double         SNR;
  histogram_type X;
  data_type      Q;
  double         P;
  JRandom        seed;
  int            debug;
 
  try { 
 
    JParser<> zap;
 
    JPseudoExperiment px;
 
    zap['E'] = make_field(setup,
                          "douplets of signal and background histograms and doublets of signal and background nuisances, \n"
                          << "\teach of which defined by <file name>:<histogram name> and <type> (values), respectively, \n"
                          << "\twhere <type> can be:" << get_keys(px.nuisance.signal.getDictionary()));
    zap['o'] = make_field(outputFile,    "output file with histograms")                     = "benchmark.root";
    zap['s'] = make_field(Fs,            "signal strength");
    zap['b'] = make_field(Fb,            "background strength")                             = 1.0;
    zap['n'] = make_field(numberOfTests, "number of tests");
    zap['M'] = make_field(M,             "lookup table for CDFs")                           = 0;
    zap['R'] = make_field(SNR,           "signal-to-noise ratio")                           = 0.0;
    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);
 
  JExperiment::setSNR(SNR);
 
 
  JGen2 px;
 
  for (const auto& i : setup) {
 
    const TObject* pS  =  getObject(i.HS);
    const TObject* pB  =  getObject(i.HB);
    const TObject* ps  =  getObject(i.Hs);
    const TObject* pb  =  getObject(i.Hb);
 
    STATUS(printer("Signal for generation:",     pS) << endl);
    STATUS(printer("Background for generation:", pB) << endl);
    STATUS(printer("Signal for evaluation:",     ps) << endl);
    STATUS(printer("Background for evaluation:", pb) << endl);
 
    JPseudoExperiment pi(pS, pB, ps, pb);
 
    pi.nuisance = i.nuisance;
 
    px.push_back(pi);
  }
 
  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 > 0.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, 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 > 0.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("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 > 0.0) {                // determine minimal likelihood ratio given maximal probability
        x = storage.getValue(P);
        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 > 0.0) { 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();
}