JRealExperiment.cc File Reference

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

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <exception>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TRandom3.h"
#include "JAstronomy/JRealExperiment.hh"
#include "JAstronomy/JGen2.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 JRealExperiment.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 139 of file JRealExperiment.cc.

{
  using namespace std;
  using namespace JPP;
 
  JContainer< vector<experiment_type> > setup;
  size_t         numberOfTests;
  size_t         M;
  double         SNR;
  bool           add;
  int            debug;
  double         testsignal;
 
  try { 
 
    JParser<> zap;
 
    JPseudoExperiment px;
 
    zap['E'] = make_field(setup,
                          "quintuplets of signal and background histograms (generation then evaluation) and data, \n"
                          << "\teach of which defined by <file name>:<histogram name> and <type> (values), respectively, \n"
                          << "\twhere <type> can be:" << get_keys(nuisance_helper));
    zap['n'] = make_field(numberOfTests, "number of tests for upper limit calculation")     = 0;
    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['s'] = make_field(testsignal,           "signal to test")                           = -1.;
    zap['d'] = make_field(debug)       = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  JExperiment::setSNR(SNR);
 
  JRealExperiment   rx;
  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);
    const TObject* pd  =  getObject(i.Hd);
 
    JPseudoExperiment pi(pS, pB, ps, pb);
 
    pi.nuisance = i.nuisance;
 
    px.push_back(pi);
    rx.add(pd, ps, pb);
 
    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(printer("Data:                     ", pd) << endl);
    STATUS(        "Nuisance:                  " << pi.nuisance << endl);
  }
 
  if (M != 0) {
    px.configure(M);
  }
 
  if (add) {
    px.add();
  }
 
  const JRealExperiment::fit_type result = rx();
 
  cout << "signal:     " << FIXED(9,3) << rx.getSignal()        << "/" << setw(6) << rx.size() << endl;
  cout << "derivative: " << FIXED(9,3) << rx.getDerivative(0.0) << endl;
  if (testsignal >= 0) {
    cout << "check TS: " << FIXED(12,5) << rx.getLikelihood(testsignal) << ' ' << FIXED(12,5) << testsignal << endl;
  }
 
  if (debug >= debug_t) {
 
    size_t n = 0;
 
    for (double x : rx) {
 
      cout << SCIENTIFIC(12,3) << 1.0/x << ((n + 1)%10 == 0 ? "\n" : " ");
 
      n += 1;
    }
    cout << endl;
  }
 
  cout << "result: "
       << FIXED(12,5) << result.likelihood  << ' '
       << FIXED(12,5) << result.signal      << endl;
 
  if (numberOfTests != 0) {
 
    const double Q  = 0.9;  // probability limit
 
    const double mu = px.getSignalStrengthForUpperLimit(rx, Q, numberOfTests);
 
    cout << "upper limit: " << FIXED(12,5) << mu << endl;
  }
}