JPseudoExperimentUpperLimit.cc File Reference

Test application for pseudo experiment generation and upper limit 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/JGarage.hh"
#include "JAstronomy/JAstronomyToolkit.hh"
#include "JLang/JVectorize.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JRandom.hh"
#include "JMath/JMathSupportkit.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 upper limit evaluations.

Author

mdejong

Definition in file JPseudoExperimentUpperLimit.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 127 of file JPseudoExperimentUpperLimit.cc.

{
  using namespace std;
  using namespace JPP;
 
  typedef JAbstractHistogram<double>   histogram_type;
 
  JContainer< vector<experiment_type> >  setup;
  string         outputFile;
  size_t         numberOfTests;
  double         boost;
  size_t         M;
  double         SNR;
  histogram_type X;
  JRandom        seed;
  double         precision;
  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['B'] = make_field(boost,         "boost livetime of experiment")                    = 1.0;
    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 signal histogram")                         = histogram_type(100, -1.0e3, +1.0e3);
    zap['p'] = make_field(precision,     "precision")                                       = 1.0e-4;
    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);
 
  const double  Q      =  0.9;                      // probability limit
 
  JPseudoExperiment px;
 
  for (const auto& i : setup) {
 
    TObject* ps  =  getObject(i.Hs);
    TObject* pb  =  getObject(i.Hb);
 
    dynamic_cast<TH1*>(ps)->Scale(boost);
    dynamic_cast<TH1*>(pb)->Scale(boost);
 
    STATUS(printer("Signal:",     ps) << endl);
    STATUS(printer("Background:", pb) << endl);
 
    px.add(ps, pb);
  }
 
  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 h0("h0", NULL, X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit());
  TH1D h1("h1", NULL, X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit());
 
 
  // generate background only pseudo experiments and store fitted signals
 
  px.set(0.0);
 
  for (size_t i = 0; i != numberOfTests; ++i) {
 
    JAspera aspera;
 
    px.run(aspera);
 
    h0.Fill(aspera(true).signal);
  }
 
 
  vector<double> mu_up;
 
  for (size_t i = 0; i != numberOfTests; ++i) {
 
    if ((i*10) < numberOfTests || (i*10)%numberOfTests == 0) {
      STATUS(setw(6) << i << '\r'); DEBUG(endl);
    }
 
    JAspera aspera;
 
    px.set(0.0);
 
    px.run(aspera);
 
    mu_up.push_back(px.getSignalStrengthForUpperLimit(aspera, Q, numberOfTests));
  }
  STATUS(endl);
 
  for (const double x : mu_up) {
    h1.Fill(x);
  }
 
  sort(mu_up.begin(), mu_up.end());
 
  cout << "<mu> " << FIXED(7,3) << mu_up[mu_up.size() / 2] << endl;
  
  out.Write();
  out.Close();
}