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/JPseudoExperiment.hh"
#include "JAstronomy/JAstronomyToolkit.hh"
 
#include "JGizmo/JGizmoToolkit.hh"
 
#include "Jeep/JContainer.hh"
#include "Jeep/JParser.hh"
 
 
namespace {
 
  /**
   * Auxiliary data structure for experiment.
   */
  struct experiment_type {
 
    JGIZMO::JRootObjectID Hd;         //!< data
    JGIZMO::JRootObjectID Hs;         //!< signal
    JGIZMO::JRootObjectID Hb;         //!< background
 
    /**
     * 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.Hd
                >> experiment.Hs
                >> experiment.Hb;
    }
 
 
    /**
     * 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.Hd << ' '
                 << experiment.Hs << ' '
                 << experiment.Hb;
    }
  };
 
 
  /**
   * 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;
  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;
  }
}