JRealExperiment.cc

Go to the documentation of this file.

#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 "JLang/JVectorize.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JTools/JAbstractHistogram.hh" //for get_keys
#include "JROOT/JRandom.hh"
 
#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
 
 
namespace {
 
  /**
   * Auxiliary data structure for experiment.
   */
  struct experiment_type {
 
    JGIZMO::JRootObjectID HS;         //!< signal for generation
    JGIZMO::JRootObjectID HB;         //!< background for generation
    JGIZMO::JRootObjectID Hs;         //!< signal
    JGIZMO::JRootObjectID Hb;         //!< background
    JGIZMO::JRootObjectID Hd;         //!< data
 
    JASTRONOMY::JPseudoExperiment::parameters_type nuisance;
 
    /**
     * Read experiment from input stream.
     *
     * \param  in             input stream
     * \param  experiment     experiment
     * \return                input stream
     */
    friend inline std::istream& operator>>(std::istream& in, experiment_type& experiment)
    {
      return in >> experiment.HS
                >> experiment.HB
                >> experiment.Hs
                >> experiment.Hb
                >> experiment.Hd        
                >> experiment.nuisance;
    }
 
 
    /**
     * Write experiment to output stream.
     *
     * \param  out            output stream
     * \param  experiment     experiment
     * \return                output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const experiment_type& experiment)
    {
      return out << experiment.HS << ' '
                 << experiment.HB << ' '
                 << experiment.Hs << ' '
                 << experiment.Hb << ' '
                 << experiment.Hd << ' '
                 << experiment.nuisance;
 
    }
  };
 
 
  /**
   * Auxiliary data structure for printing.
   */
  struct printer {
    /**
     * Constructor.
     *
     * \param  title          title
     * \param  ps             pointer to object
     */
    printer(const char* const title,
            const TObject*    ps) :
      title(title),
      ps(ps)
    {}
 
    /**
     * Write printer to output stream.
     *
     * \param  out            output stream
     * \param  printer        printer
     * \return                output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const printer& printer)
    {
      using namespace std;
      using namespace JPP;
 
      out << setw(16) << left << printer.title << right;
 
      if (printer.ps != NULL) {
 
        out << ' ' << setw(16) << left << printer.ps->GetName() << right;
 
        const TH1* h1 = dynamic_cast<const TH1*>(printer.ps);
 
        if (h1 != NULL) {
          out << ' ' << FIXED(10,3) << h1->GetSumOfWeights();
        }
      }
 
      return out;
    }
 
  private:
    const char* const title;
    const TObject*    ps;
  };
}
 
 
/**
 * \file
 * Test application for pseudo experiment generation and likelihood ratio evaluations.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  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();
    rx.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;
  }
}