JPseudoExperimentUpperLimit.cc

Go to the documentation of this file.

#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"
 
 
namespace {
 
  /**
   * Auxiliary data structure for experiment.
   */
  struct experiment_type {
 
    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.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.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 upper limit evaluations.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  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
  const double  MUMIN  =  1.0e-10;                  // minimal signal
  const double  MUMAX  =  1.0e+10;                  // maximal signal
  const size_t  N      =  numberOfTests * 5;        // statistics for probability calculation
 
  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);
 
  const JPseudoExperiment::h0_type mu_0(px, numberOfTests);
 
  for (const double x : mu_0) {
    h0.Fill(x);
  }
 
 
  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);
 
    const JAspera::fit_type result = aspera(true);
 
    DEBUG("pseudo experiment " << setw(3) << aspera.size() << ' ' << FIXED(12,5) << result.signal << ' ' << FIXED(12,5) << result.likelihood << endl); 
 
    double mu = 0.0;
 
    if (mu_0.getFractionBelow(result.signal) > 1.0 - Q) {
 
      double mumin = 1.0;
      double mumax = 1.0;
 
      {
        for ( ; ; mumax *= 2.0) {
 
          const double ts = aspera.getTestStatisticForUpperLimit(mumax);
 
          px.set(mumax);
 
          const double ps = px.getProbabilityForUpperLimit(mumax, ts, N);
 
          DEBUG("maximal value " << SCIENTIFIC(12,3) << mumax << ' ' << FIXED(12,5) << ts << ' ' << FIXED(12,5) << ps << endl); 
 
          if (ps < 1.0 - Q || mumax > MUMAX) {
            break;
          }
        }
      }
 
      if (mumax == 1.0) {
 
        mumin = 0.5*mumax;
 
        for ( ; ; mumin *= 0.5) {
 
          const double ts = aspera.getTestStatisticForUpperLimit(mumin);
 
          px.set(mumin);
 
          const double ps = px.getProbabilityForUpperLimit(mumin, ts, N);
 
          DEBUG("minimal value " << SCIENTIFIC(12,3) << mumin << ' ' << FIXED(12,5) << ts << ' ' << FIXED(12,5) << ps << endl); 
 
          if (ps > 1.0 - Q || mumin < MUMIN) {
            break;
          }
        }
 
        mumax = 2.0*mumin;
 
      } else {
 
        mumin = 0.5*mumax;
      }
 
      if      (mumin < MUMIN)
        mu = mumin;
      else if (mumax > MUMAX)
        mu = mumax;
      else {
 
        for ( ; ; ) {                               // binary search
 
          mu = 0.5 * (mumin + mumax);
 
          const double ts = aspera.getTestStatisticForUpperLimit(mu);
 
          px.set(mu);
 
          const double ps = px.getProbabilityForUpperLimit(mu, ts, N);
 
          DEBUG("current value " << SCIENTIFIC(12,3) << mu    << ' ' << FIXED(12,5) << ts << ' ' << FIXED(12,5) << ps << endl); 
 
          if (fabs(ps - (1.0 - Q)) <= precision) {
            break;
          }
        
          if (ps >= 1.0 - Q)
            mumin = mu;
          else
            mumax = mu;
        }
      }
    }
 
    mu_up.push_back(mu);
  }
  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();
}