Jpp 20.0.0-27-g39925593c-D
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JPseudoExperiment.cc File Reference

Test application for pseudo experiment generation and likelihood ratio evaluations. More...

#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"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Detailed Description

Test application for pseudo experiment generation and likelihood ratio evaluations.

Author
mdejong

Definition in file JPseudoExperiment.cc.

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 127 of file JPseudoExperiment.cc.

128{
129 using namespace std;
130 using namespace JPP;
131
132 typedef JAbstractHistogram<double> histogram_type;
133 typedef float data_type;
134
135 JPseudoExperiment px;
136
137 JContainer< vector<experiment_type> > setup;
138 string outputFile;
139 size_t numberOfTests;
140 double Fs;
141 double Fb;
142 double boost;
143 size_t M;
144 double SNR;
145 bool add;
146 histogram_type X;
147 data_type Q;
148 vector<double> P;
149 JRandom seed;
150 int debug;
151
152 try {
153
154 JParser<> zap;
155
156 zap['E'] = make_field(setup,
157 "douplets of signal and background histograms, "
158 << "each of which defined by <file name>:<histogram name>");
159 zap['o'] = make_field(outputFile, "output file with histograms") = "benchmark.root";
160 zap['n'] = make_field(numberOfTests, "number of tests");
161 zap['s'] = make_field(Fs, "signal strength");
162 zap['b'] = make_field(Fb, "background strength") = 1.0;
163 zap['B'] = make_field(boost, "boost livetime of experiment") = 1.0;
164 zap['M'] = make_field(M, "lookup table for CDFs") = 0;
165 zap['R'] = make_field(SNR, "signal-to-noise ratio") = 0.0;
166 zap['A'] = make_field(add, "add remnant signal and noise");
167 zap['N'] = make_field(px.nuisance,
168 "signal and background nuisances, "
169 << "each of which defined by <type> (values), \n"
170 << "\twhere <type> can be:" << get_keys(nuisance_helper)) = JPARSER::initialised();
171 zap['x'] = make_field(X, "x-axis likelihood histogram") = histogram_type(110, -10.0, +100.0);
172 zap['Q'] = make_field(Q, "minimal likelihood ratio") = 0.0;
173 zap['P'] = make_field(P, "maximal probability") = JPARSER::initialised();
174 zap['S'] = make_field(seed) = 0;
175 zap['d'] = make_field(debug) = 1;
176
177 zap(argc, argv);
178 }
179 catch(const exception& error) {
180 FATAL(error.what() << endl);
181 }
182
183
184 seed.set(gRandom);
185
186 JExperiment::setSNR(SNR);
187
188 for (const auto& i : setup) {
189
190 TObject* ps = getObject(i.Hs);
191 TObject* pb = getObject(i.Hb);
192
193 dynamic_cast<TH1*>(ps)->Scale(boost);
194 dynamic_cast<TH1*>(pb)->Scale(boost);
195
196 STATUS(printer("Signal:", ps) << endl);
197 STATUS(printer("Background:", pb) << endl);
198
199 px.add(ps, pb);
200 }
201
202 px.set(Fs, Fb);
203
204 if (M != 0) {
205 px.configure(M);
206 }
207
208 if (add) {
209 px.add();
210 }
211
212 STATUS("Signal: " << FIXED(10,3) << px.getSignal() << endl);
213 STATUS("Background: " << FIXED(10,3) << px.getBackground() << endl);
214 STATUS("Number of tests: " << setw(10) << numberOfTests << endl);
215
216
217 TFile out(outputFile.c_str(), "recreate");
218
219 out << *gRandom;
220
221 TH1D hl("hl", NULL, X.getNumberOfBins(), X.getLowerLimit(), X.getUpperLimit());
222 TH1D hn("hn", NULL, 100, -15.0, +15.0);
223
224
225 if (P.size() > 1) { // special case
226
227 JGarage<data_type, 100> storage;
228
229 double Ps = 0.0; // total probability
230 size_t nb = px.getBackground() + 1; // start value for number of background events
231 size_t mb = 0;
232
233 {
234 const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
235
236 for ( ; ; ++nb) {
237
238 const double p = poisson(nb, px.getBackground());
239 const size_t n = p * numberOfTests;
240
241 if (p > P[1]) {
242 continue;
243 }
244
245 if (n == 0) {
246 break;
247 }
248
249 if (mb == 0) {
250 mb = nb;
251
252 }
253 for (size_t i = 0; i != n; ++i) {
254 storage.put(px(nb).likelihood);
255 }
256
257 Ps += p;
258 }
259
260 const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
261
262 STATUS("Total time: " << setw(8) << (t1 - t0) / chrono::milliseconds(1) << " [ms]" << endl);
263 STATUS("Average time: " << setw(6) << (t1 - t0) / chrono::nanoseconds(numberOfTests) << " [ns]" << endl);
264 }
265
266 STATUS("Poisson: " << setw(3) << mb << ' ' << SCIENTIFIC(12,3) << Ps << ' ' << setw(8) << storage.size() << endl);
267
268 storage([&hl](const data_type x) { hl.Fill(x); });
269
270 double x = 0.0;
271 double p = 0.0;
272
273 {
274 const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
275
276 x = storage.getValue(P[0]/Ps);
277 p = storage.getProbability(x) * Ps;
278
279 const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
280
281 STATUS("Time to result: " << setw(6) << (t1 - t0) / chrono::nanoseconds(1) << " [ns]" << endl);
282 }
283
284 cout << "Minimal likelihood ratio: " << FIXED(9,5) << x << ' ' << SCIENTIFIC(12,3) << p << endl;
285
286 } else if (P.size() > 0 && Q > 0.0) {
287
288 const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
289
290 const double signal = discovery<float, 100>(px, Q, P[0], numberOfTests);
291
292 const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
293
294 STATUS("Total time: " << setw(8) << (t1 - t0) / chrono::milliseconds(1) << " [ms]" << endl);
295
296 cout << "Signal strength: " << FIXED(9,5) << signal << endl;
297
298 } else if (P.size() > 0 || Q > 0.0) {
299
300 JGarage<data_type, 100> storage;
301
302 {
303 const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
304
305 for (size_t i = 0; i != numberOfTests; ++i) {
306 storage.put(px().likelihood);
307 }
308
309 const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
310
311 STATUS("Total time: " << setw(8) << (t1 - t0) / chrono::milliseconds(1) << " [ms]" << endl);
312 STATUS("Average time: " << setw(6) << (t1 - t0) / chrono::nanoseconds(numberOfTests) << " [ns]" << endl);
313 }
314
315 storage([&hl](const data_type x) { hl.Fill(x); });
316
317 double x = 0.0;
318 double p = 0.0;
319
320 {
321 const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
322
323 if (P.size() == 1) { // determine minimal likelihood ratio given maximal probability
324 x = storage.getValue(P[0]);
325 p = storage.getProbability(x);
326 }
327
328 if (Q > 0.0) { // determine maximal probability given minimal likelihood ratio
329 p = storage.getProbability(Q);
330 x = storage.getValue(p);
331 }
332
333 const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
334
335 STATUS("Time to result: " << setw(6) << (t1 - t0) / chrono::microseconds(1) << " [us]" << endl);
336 }
337
338 if (P.size() == 1) { cout << "Minimal likelihood ratio: " << FIXED(9,5) << x << ' ' << SCIENTIFIC(12,3) << p << endl; }
339 if (Q > 0.0) { cout << "Maximal probability: " << SCIENTIFIC(12,3) << p << ' ' << FIXED(9,5) << x << endl; }
340
341 } else {
342
343 std::vector<JPseudoExperiment::result_type> storage(numberOfTests);
344
345 const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
346
347 px(storage);
348
349 const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
350
351 STATUS("Average time: " << setw(6) << (t1 - t0) / chrono::nanoseconds(storage.size()) << " [ns]" << endl);
352
353 for (const auto& result : storage) {
354
355 DEBUG("result: "
356 << setw(3) << result.ns << ' '
357 << setw(3) << result.nb << ' '
358 << FIXED(12,5) << result.likelihood << ' '
359 << FIXED(12,5) << result.signal << endl);
360
361 hl.Fill(result.likelihood);
362 hn.Fill(result.signal - result.ns);
363 }
364 }
365
366 out.Write();
367 out.Close();
368}
outputFile
string outputFile
Definition JDAQTimesliceSelector.cc:37
DEBUG
#define DEBUG(A)
Message macros.
Definition JMessage.hh:62
STATUS
#define STATUS(A)
Definition JMessage.hh:63
FATAL
#define FATAL(A)
Definition JMessage.hh:67
debug
int debug
debug level
Definition JSirene.cc:72
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JOSCPROB::JOscProbInterpolator
Template definition of a multi-dimensional oscillation probability interpolation table.
Definition JOscProbInterpolator.hh:71
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JGIZMO::getObject
TObject * getObject(const JRootObjectID &id)
Get first TObject with given identifier.
Definition JGizmoToolkit.hh:165
JLANG::get_keys
const array_type< JKey_t > & get_keys(const std::map< JKey_t, JValue_t, JComparator_t, JAllocator_t > &data)
Method to create array of keys of map.
Definition JVectorize.hh:139
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:45
JTOOLS::n
const int n
Definition JPolint.hh:791
FIXED
Auxiliary data structure for floating point format specification.
Definition JManip.hh:448
JASTRONOMY::JGarage::getProbability
double getProbability(const T value) const
Get maximal probability corresponding given minimal value.
Definition JGarage.hh:199
JASTRONOMY::JGarage::getValue
T getValue(const double P) const
Get minimal value corresponding given maximal probability.
Definition JGarage.hh:157
JASTRONOMY::JGarage::put
void put(const T value)
Add value to container.
Definition JGarage.hh:108
JASTRONOMY::JPseudoExperiment
Pseudo experiment using CDF for combined generation and likelihood evaluation.
Definition JPseudoExperiment.hh:341
JASTRONOMY::JPseudoExperiment::add
void add(const TObject *ps, const TObject *pb)
Add objects with PDFs of signal and background.
Definition JPseudoExperiment.hh:389
JASTRONOMY::JPseudoExperiment::nuisance
struct JASTRONOMY::JPseudoExperiment::parameters_type nuisance
JASTRONOMY::JPseudoExperiment::getSignal
double getSignal() const
Get total signal.
Definition JPseudoExperiment.hh:613
JASTRONOMY::JPseudoExperiment::set
virtual void set(const double fS, const double fB=1.0) override
Set scaling factors of signal and background strengths.
Definition JPseudoExperiment.hh:636
JASTRONOMY::JPseudoExperiment::getBackground
double getBackground() const
Get total background.
Definition JPseudoExperiment.hh:624
JASTRONOMY::JPseudoExperiment::configure
void configure(size_t N)
Configure lookup tables.
Definition JPseudoExperiment.hh:738
JCALIBRATE::data_type
Data structure for measured coincidence rates of all pairs of PMTs in optical module.
Definition JFitK40.hh:103
JMATH::JRandom
Template definition of random value generator.
Definition JMath/JRandom.hh:26
JPARSER::initialised
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition JParser.hh:68
SCIENTIFIC
Auxiliary data structure for floating point format specification.
Definition JManip.hh:488
Generated by doxygen 1.10.0