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/JPseudoExperiment.hh"
#include "JAstronomy/JAstronomyToolkit.hh"
#include "JGizmo/JGizmoToolkit.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 124 of file JRealExperiment.cc.

{
  using namespace std;
  using namespace JPP;
 
  JContainer< vector<experiment_type> > setup;
  size_t         numberOfTests;
  double         SNR;
  int            debug;
 
  try { 
 
    JParser<> zap;
 
    zap['E'] = make_field(setup,
                          "triplets of data, signal and background histograms, "\
                          "each of which defined by <file name>:<histogram name>");
    zap['n'] = make_field(numberOfTests, "number of tests for upper limit calculation")     = 0;
    zap['R'] = make_field(SNR,           "signal-to-noise ratio")                           = 0.0;
    zap['d'] = make_field(debug)       = 1;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  JExperiment::setSNR(SNR);
 
  JRealExperiment   rx;
  JPseudoExperiment px;
 
  for (const auto& i : setup) {
 
    const TObject* pd  =  getObject(i.Hd);
    const TObject* ps  =  getObject(i.Hs);
    const TObject* pb  =  getObject(i.Hb);
 
    STATUS(printer("Data:",       pd) << endl);
    STATUS(printer("Signal:",     ps) << endl);
    STATUS(printer("Background:", pb) << endl);
 
    rx.add(pd, ps, pb);
    px.add(ps, pb);
  }
 
  
  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 (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;
  }
}