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 179 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(px.nuisance.signal.getDictionary()) <<"\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(px.nuisance.signal.getDictionary()) <<"\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();
}