JPseudoExperiment.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/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"

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 JPseudoExperiment.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 125 of file JPseudoExperiment.cc.

{
  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();
}