JPseudoExperimentCondNLLR.cc File Reference

Application for generating conditional Likelihood histograms of part of a dataset compared to the complete dataset. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <chrono>
#include <vector>
#include <map>
#include <exception>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"
#include "TRandom3.h"
#include "JAstronomy/JPseudoExperiment.hh"
#include "JAstronomy/JAspera.hh"
#include "JAstronomy/JAstronomyToolkit.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "JTools/JAbstractHistogram.hh"
#include "JROOT/JRandom.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

Application for generating conditional Likelihood histograms of part of a dataset compared to the complete dataset.

Used to generate trial factors with JTrialFactorsFromCondNLLR.cc

Author

mdejong fvazquezdesola

Definition in file JPseudoExperimentCondNLLR.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 188 of file JPseudoExperimentCondNLLR.cc.

{
  using namespace std;
  using namespace JPP;
 
  typedef JAbstractHistogram<double> histogram_type;
 
  JContainer< vector<experiment_type> >  setup_vec1;
  JContainer< vector<experiment_type> >  setup_vec2;
  string              outputFile;
  size_t              numberOfTests;
  double              Fs;
  double              Fb;
  size_t              M;
  double              SNR;
  histogram_type      X;
  JRandom             seed;
  int                 debug;
 
  try { 
 
    JParser<> zap;
    JPseudoExperiment px;
 
    zap['E'] = make_field(setup_vec1,
                          "inputs for 1st data period: signal and background histograms for generation and evaluation,\n"
                          << "\texposure boost, and signal and background nuisances, provided as\n"
                          << "\t'<file>:<hgen_s name> <file>:<hgen_b name> <file>:<heval_s name> <file>:<heval_b name> <boost> <type_s> (values) <type_b> (values)' \n"
                          << "\twhere <type> can be:" << get_keys(nuisance_helper) <<"\n."
                          << "\t Can be called repeatedly to add multiple datasets to this period");
    zap['F'] = make_field(setup_vec2,
                          "inputs for 2nd data period: signal and background histograms for generation and evaluation,\n"
                          << "\texposure boost, and signal and background nuisances, provided as\n"
                          << "\t'<file>:<hgen_s name> <file>:<hgen_b name> <file>:<heval_s name> <file>:<heval_b name> <boost> <type_s> (values) <type_b> (values)' \n"
                          << "\twhere <type> can be:" << get_keys(nuisance_helper) <<"\n."
                          << "\t Can be called repeatedly to add multiple datasets to this period");
    zap['o'] = make_field(outputFile,       "output file with histograms")                  = "CondNLLR.root";
    zap['s'] = make_field(Fs,               "signal strength");
    zap['b'] = make_field(Fb,               "background strength")                          = 1.0;
    zap['M'] = make_field(M,                "lookup table for CDFs")                        = 0;
    zap['R'] = make_field(SNR,              "signal-to-noise ratio")                        = 0.0;
    zap['n'] = make_field(numberOfTests,    "number of tests / PEs");
    zap['x'] = make_field(X,                "x-axis for likelihood histogram")              = histogram_type(200, 0, +20.0);
    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);
 
  cout << "**Setting up 1st period**" << endl;
  experiment_group pxg_1(setup_vec1, Fs, Fb, M);
  cout << "**Setting up 2nd period**" << endl;
  experiment_group pxg_2(setup_vec2, Fs, Fb, M);
 
  TFile out(outputFile.c_str(), "recreate");
 
  TH2D h2d_nllr("h2d_nllr", "Negative Log-Likelihood Ratio; NLLR 1st period; NLLR total",
                X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(),
                X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit());
  TH2D h2d_nllr_constr("h2d_nllr_constr", "Negative Log-Likelihood Ratio (filtered); NLLR 1st period; NLLR total",
                       X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit(),
                       X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit());
 
  size_t nexp_passfilter = 0;
 
  for (size_t nexp =0; nexp != numberOfTests; ++nexp) {
    if (nexp%100000 == 0) {
      cout << "Done " << nexp << " pseudo-experiments out of " << numberOfTests << endl;
    }
 
    JAspera aspera;
    JAspera::fit_type result_1;
    JAspera::fit_type result_A;
 
    for (const auto& px : pxg_1.px_vec) { px.run(aspera); }
    result_1 = aspera();
 
    for (const auto& px : pxg_2.px_vec) { px.run(aspera); }
    result_A = aspera();
 
    auto NLLR_1 = result_1.likelihood;
    auto NLLR_A = result_A.likelihood;
 
    if (NLLR_1 < 0) {
      NLLR_1 = 0;
    }
 
    if (NLLR_A < 0) {
      NLLR_A = 0;
    }
 
    h2d_nllr.Fill(NLLR_1, NLLR_A, 1./numberOfTests);
 
    if( result_1.signal * pxg_1.signal_weight > 1) {
      nexp_passfilter++;
      h2d_nllr_constr.Fill(NLLR_1, NLLR_A, 1./numberOfTests);
    }
  }
 
  cout << nexp_passfilter << " PEs pass the first period filter, out of " << numberOfTests << endl;
 
  h2d_nllr.Write();
  h2d_nllr_constr.Write();
 
  out.Write();
  out.Close();
}