JGen2.cc File Reference

Test application for pseudo experiment generation and likelihood ratio evaluations. More...

#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"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Test application for pseudo experiment generation and likelihood ratio evaluations.

Author

mdejong

Definition in file JGen2.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 135 of file JGen2.cc.

{
  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;
  bool           add;
  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(nuisance_helper));
    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['A'] = make_field(add,           "add remnant signal and noise");
    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);
 
    JPseudoExperiment pi(pS, pB, ps, pb);
    pi.nuisance = i.nuisance;
    px.push_back(pi);
 
    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);
    STATUS(        "Nuisance:                  "<< pi.nuisance << endl);
 
  }
 
  px.set(Fs, Fb);
 
  if (M != 0) {
    px.configure(M);
  }
 
  if (add) {
    px.add();
  }
  
  STATUS("Signal:          " << FIXED(10,3) << px.getSignal()     << endl);
  STATUS("Background:      " << FIXED(10,3) << px.getBackground() << endl);
  STATUS("Number of tests: " << setw(10)    << numberOfTests << endl);
 
  if (px.getBackground() <= 0.0) {
    FATAL("No background." << endl);
  }
 
  if (debug >= debug_t) {
    for (size_t i = 0; i != px.size(); ++i) {
      cout << setw(2) << i << ' ' << px[i].nuisance.signal << ' ' << px[i].nuisance.background << 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 if (numberOfTests != 0) {
 
    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();
}