ROOT   6.21/01 Reference Guide
HybridInstructional.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_roostats
3 /// \notebook -js
4 /// Example demonstrating usage of HybridCalcultor
5 ///
6 /// A hypothesis testing example based on number counting
7 /// with background uncertainty.
8 ///
9 /// NOTE: This example must be run with the ACLIC (the + option ) due to the
10 /// new class that is defined.
11 ///
12 /// This example:
13 /// - demonstrates the usage of the HybridCalcultor (Part 4-6)
14 /// - demonstrates the numerical integration of RooFit (Part 2)
15 /// - validates the RooStats against an example with a known analytic answer
16 /// - demonstrates usage of different test statistics
17 /// - explains subtle choices in the prior used for hybrid methods
18 /// - demonstrates usage of different priors for the nuisance parameters
19 /// - demonstrates usage of PROOF
20 ///
21 /// The basic setup here is that a main measurement has observed x events with an
22 /// expectation of s+b. One can choose an ad hoc prior for the uncertainty on b,
23 /// or try to base it on an auxiliary measurement. In this case, the auxiliary
24 /// measurement (aka control measurement, sideband) is another counting experiment
25 /// with measurement y and expectation tau*b. With an 'original prior' on b,
26 /// called \f$\eta(b)\f$ then one can obtain a posterior from the auxiliary measurement
27 /// \f$\pi(b) = \eta(b) * Pois(y|tau*b).\f$ This is a principled choice for a prior
28 /// on b in the main measurement of x, which can then be treated in a hybrid
29 /// Bayesian/Frequentist way. Additionally, one can try to treat the two
30 /// measurements simultaneously, which is detailed in Part 6 of the tutorial.
31 ///
32 /// This tutorial is related to the FourBin.C tutorial in the modeling, but
33 /// focuses on hypothesis testing instead of interval estimation.
34 ///
35 /// More background on this 'prototype problem' can be found in the
36 /// following papers:
37 ///
38 /// - Evaluation of three methods for calculating statistical significance
39 /// when incorporating a systematic uncertainty into a test of the
40 /// background-only hypothesis for a Poisson process
41 /// Authors: Robert D. Cousins, James T. Linnemann, Jordan Tucker
42 /// http://arxiv.org/abs/physics/0702156
43 /// NIM A 595 (2008) 480--501
44 ///
45 /// - Statistical Challenges for Searches for New Physics at the LHC
46 /// Author: Kyle Cranmer
47 /// http://arxiv.org/abs/physics/0511028
48 ///
49 /// - Measures of Significance in HEP and Astrophysics
50 /// Authors J. T. Linnemann
51 /// http://arxiv.org/abs/physics/0312059
52 ///
53 /// \macro_image
54 /// \macro_output
55 /// \macro_code
56 ///
57 /// \authors Kyle Cranmer, Wouter Verkerke, and Sven Kreiss
58
59 #include "RooGlobalFunc.h"
60 #include "RooRealVar.h"
61 #include "RooProdPdf.h"
62 #include "RooWorkspace.h"
63 #include "RooDataSet.h"
64 #include "TCanvas.h"
65 #include "TStopwatch.h"
66 #include "TH1.h"
67 #include "RooPlot.h"
68 #include "RooMsgService.h"
69
71
73 #include "RooStats/ToyMCSampler.h"
74 #include "RooStats/HypoTestPlot.h"
75
78
80
81 using namespace RooFit;
82 using namespace RooStats;
83
84 // ----------------------------------
85 // A New Test Statistic Class for this example.
86 // It simply returns the sum of the values in a particular
87 // column of a dataset.
88 // You can ignore this class and focus on the macro below
89 class BinCountTestStat : public TestStatistic {
90 public:
91  BinCountTestStat(void) : fColumnName("tmp") {}
92  BinCountTestStat(string columnName) : fColumnName(columnName) {}
93
94  virtual Double_t Evaluate(RooAbsData &data, RooArgSet & /*nullPOI*/)
95  {
96  // This is the main method in the interface
97  Double_t value = 0.0;
98  for (int i = 0; i < data.numEntries(); i++) {
99  value += data.get(i)->getRealValue(fColumnName.c_str());
100  }
101  return value;
102  }
103  virtual const TString GetVarName() const { return fColumnName; }
104
105 private:
106  string fColumnName;
107
108 protected:
109  ClassDef(BinCountTestStat, 1)
110 };
111
112 ClassImp(BinCountTestStat)
113
114  // ----------------------------------
115  // The Actual Tutorial Macro
116
117  void HybridInstructional()
118 {
119
120  // This tutorial has 6 parts
122  // Setup
123  // 1. Make the model for the 'prototype problem'
124  // Special cases
125  // 2. Use RooFit's direct integration to get p-value & significance
126  // 3. Use RooStats analytic solution for this problem
127  // RooStats HybridCalculator -- can be generalized
128  // 4. RooStats ToyMC version of 2. & 3.
129  // 5. RooStats ToyMC with an equivalent test statistic
130  // 6. RooStats ToyMC with simultaneous control & main measurement
131
132  // It takes ~4 min without PROOF and ~2 min with PROOF on 4 cores.
133  // Of course, everything looks nicer with more toys, which takes longer.
134
135  TStopwatch t;
136  t.Start();
137  TCanvas *c = new TCanvas;
138  c->Divide(2, 2);
139
140  // ----------------------------------------------------
141  // P A R T 1 : D I R E C T I N T E G R A T I O N
142  // ====================================================
143  // Make model for prototype on/off problem
144  // $Pois(x | s+b) * Pois(y | tau b )$
145  // for Z_Gamma, use uniform prior on b.
146  RooWorkspace *w = new RooWorkspace("w");
147  w->factory("Poisson::px(x[150,0,500],sum::splusb(s[0,0,100],b[100,0,300]))");
148  w->factory("Poisson::py(y[100,0,500],prod::taub(tau[1.],b))");
149  w->factory("PROD::model(px,py)");
150  w->factory("Uniform::prior_b(b)");
151
152  // We will control the output level in a few places to avoid
153  // verbose progress messages. We start by keeping track
154  // of the current threshold on messages.
156
157  // Use PROOF-lite on multi-core machines
158  ProofConfig *pc = NULL;
159  // uncomment below if you want to use PROOF
160  // ~~~
161  // pc = new ProofConfig(*w, 4, "workers=4", kFALSE); // machine with 4 cores
162  // pc = new ProofConfig(*w, 2, "workers=2", kFALSE); // machine with 2 cores
163  // ~~~
164
165  // ----------------------------------------------------
166  // P A R T 2 : D I R E C T I N T E G R A T I O N
167  // ====================================================
168  // This is not the 'RooStats' way, but in this case the distribution
169  // of the test statistic is simply x and can be calculated directly
170  // from the PDF using RooFit's built-in integration.
171  // Note, this does not generalize to situations in which the test statistic
172  // depends on many events (rows in a dataset).
173
174  // construct the Bayesian-averaged model (eg. a projection pdf)
175  // $p'(x|s) = \int db p(x|s+b) * [ p(y|b) * prior(b) ]$
176  w->factory("PROJ::averagedModel(PROD::foo(px|b,py,prior_b),b)");
177
178  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR); // lower message level
179  // plot it, red is averaged model, green is b known exactly, blue is s+b av model
180  RooPlot *frame = w->var("x")->frame(Range(50, 230));
181  w->pdf("averagedModel")->plotOn(frame, LineColor(kRed));
182  w->pdf("px")->plotOn(frame, LineColor(kGreen));
183  w->var("s")->setVal(50.);
184  w->pdf("averagedModel")->plotOn(frame, LineColor(kBlue));
185  c->cd(1);
186  frame->Draw();
187  w->var("s")->setVal(0.);
188
189  // compare analytic calculation of Z_Bi
190  // with the numerical RooFit implementation of Z_Gamma
191  // for an example with x = 150, y = 100
192
193  // numeric RooFit Z_Gamma
194  w->var("y")->setVal(100);
195  w->var("x")->setVal(150);
196  RooAbsReal *cdf = w->pdf("averagedModel")->createCdf(*w->var("x"));
197  cdf->getVal(); // get ugly print messages out of the way
198  cout << "-----------------------------------------" << endl;
199  cout << "Part 2" << endl;
200  cout << "Hybrid p-value from direct integration = " << 1 - cdf->getVal() << endl;
201  cout << "Z_Gamma Significance = " << PValueToSignificance(1 - cdf->getVal()) << endl;
202  RooMsgService::instance().setGlobalKillBelow(msglevel); // set it back
203
204  // ---------------------------------------------
205  // P A R T 3 : A N A L Y T I C R E S U L T
206  // =============================================
207  // In this special case, the integrals are known analytically
208  // and they are implemented in RooStats::NumberCountingUtils
209
210  // analytic Z_Bi
211  double p_Bi = NumberCountingUtils::BinomialWithTauObsP(150, 100, 1);
212  double Z_Bi = NumberCountingUtils::BinomialWithTauObsZ(150, 100, 1);
213  cout << "-----------------------------------------" << endl;
214  cout << "Part 3" << endl;
215  std::cout << "Z_Bi p-value (analytic): " << p_Bi << std::endl;
216  std::cout << "Z_Bi significance (analytic): " << Z_Bi << std::endl;
217  t.Stop();
218  t.Print();
219  t.Reset();
220  t.Start();
221
222  // -------------------------------------------------------------
223  // P A R T 4 : U S I N G H Y B R I D C A L C U L A T O R
224  // =============================================================
225  // Now we demonstrate the RooStats HybridCalculator.
226  //
227  // Like all RooStats calculators it needs the data and a ModelConfig
228  // for the relevant hypotheses. Since we are doing hypothesis testing
229  // we need a ModelConfig for the null (background only) and the alternate
230  // (signal+background) hypotheses. We also need to specify the PDF,
231  // the parameters of interest, and the observables. Furthermore, since
232  // the parameter of interest is floating, we need to specify which values
233  // of the parameter corresponds to the null and alternate (eg. s=0 and s=50)
234  //
235  // define some sets of variables obs={x} and poi={s}
236  // note here, x is the only observable in the main measurement
237  // and y is treated as a separate measurement, which is used
238  // to produce the prior that will be used in this calculation
239  // to randomize the nuisance parameters.
240  w->defineSet("obs", "x");
241  w->defineSet("poi", "s");
242
243  // create a toy dataset with the x=150
244  RooDataSet *data = new RooDataSet("d", "d", *w->set("obs"));
246
247  // Part 3a : Setup ModelConfigs
248  // -------------------------------------------------------
249  // create the null (background-only) ModelConfig with s=0
250  ModelConfig b_model("B_model", w);
251  b_model.SetPdf(*w->pdf("px"));
252  b_model.SetObservables(*w->set("obs"));
253  b_model.SetParametersOfInterest(*w->set("poi"));
254  w->var("s")->setVal(0.0); // important!
255  b_model.SetSnapshot(*w->set("poi"));
256
257  // create the alternate (signal+background) ModelConfig with s=50
258  ModelConfig sb_model("S+B_model", w);
259  sb_model.SetPdf(*w->pdf("px"));
260  sb_model.SetObservables(*w->set("obs"));
261  sb_model.SetParametersOfInterest(*w->set("poi"));
262  w->var("s")->setVal(50.0); // important!
263  sb_model.SetSnapshot(*w->set("poi"));
264
265  // Part 3b : Choose Test Statistic
266  // ----------------------------------
267  // To make an equivalent calculation we need to use x as the test
268  // statistic. This is not a built-in test statistic in RooStats
269  // so we define it above. The new class inherits from the
270  // RooStats::TestStatistic interface, and simply returns the value
271  // of x in the dataset.
272
273  BinCountTestStat binCount("x");
274
275  // Part 3c : Define Prior used to randomize nuisance parameters
276  // -------------------------------------------------------------
277  //
278  // The prior used for the hybrid calculator is the posterior
279  // from the auxiliary measurement y. The model for the aux.
280  // measurement is Pois(y|tau*b), thus the likelihood function
281  // is proportional to (has the form of) a Gamma distribution.
282  // if the 'original prior' $\eta(b)$ is uniform, then from
283  // Bayes's theorem we have the posterior:
284  // $\pi(b) = Pois(y|tau*b) * \eta(b)$
285  // If $\eta(b)$ is flat, then we arrive at a Gamma distribution.
286  // Since RooFit will normalize the PDF we can actually supply
287  // $py=Pois(y,tau*b)$ that will be equivalent to multiplying by a uniform.
288  //
289  // Alternatively, we could explicitly use a gamma distribution:
290  // w->factory("Gamma::gamma(b,sum::temp(y,1),1,0)");
291  //
292  // or we can use some other ad hoc prior that do not naturally
293  // follow from the known form of the auxiliary measurement.
294  // The common choice is the equivalent Gaussian:
295  w->factory("Gaussian::gauss_prior(b,y, expr::sqrty('sqrt(y)',y))");
296  // this corresponds to the "Z_N" calculation.
297  //
298  // or one could use the analogous log-normal prior
299  w->factory("Lognormal::lognorm_prior(b,y, expr::kappa('1+1./sqrt(y)',y))");
300  //
301  // Ideally, the HybridCalculator would be able to inspect the full
302  // model $Pois(x | s+b) * Pois(y | tau b )$ and be given the original
303  // prior $\eta(b)$ to form $\pi(b) = Pois(y|tau*b) * \eta(b)$.
304  // This is not yet implemented because in the general case
305  // it is not easy to identify the terms in the PDF that correspond
306  // to the auxiliary measurement. So for now, it must be set
307  // explicitly with:
308  // - ForcePriorNuisanceNull()
309  // - ForcePriorNuisanceAlt()
310  // the name "ForcePriorNuisance" was chosen because we anticipate
311  // this to be auto-detected, but will leave the option open
312  // to force to a different prior for the nuisance parameters.
313
314  // Part 3d : Construct and configure the HybridCalculator
315  // -------------------------------------------------------
316
317  HybridCalculator hc1(*data, sb_model, b_model);
318  ToyMCSampler *toymcs1 = (ToyMCSampler *)hc1.GetTestStatSampler();
319  toymcs1->SetNEventsPerToy(1); // because the model is in number counting form
320  toymcs1->SetTestStatistic(&binCount); // set the test statistic
321  hc1.SetToys(20000, 1000);
322  hc1.ForcePriorNuisanceAlt(*w->pdf("py"));
323  hc1.ForcePriorNuisanceNull(*w->pdf("py"));
324  // if you wanted to use the ad hoc Gaussian prior instead
325  // ~~~
326  // hc1.ForcePriorNuisanceAlt(*w->pdf("gauss_prior"));
327  // hc1.ForcePriorNuisanceNull(*w->pdf("gauss_prior"));
328  // ~~~
329  // if you wanted to use the ad hoc log-normal prior instead
330  // ~~~
331  // hc1.ForcePriorNuisanceAlt(*w->pdf("lognorm_prior"));
332  // hc1.ForcePriorNuisanceNull(*w->pdf("lognorm_prior"));
333  // ~~~
334
335  // enable proof
336  // NOTE: This test statistic is defined in this macro, and is not
337  // working with PROOF currently. Luckily test stat is fast to evaluate.
338  // if(pc) toymcs1->SetProofConfig(pc);
339
340  // these lines save current msg level and then kill any messages below ERROR
342  // Get the result
343  HypoTestResult *r1 = hc1.GetHypoTest();
344  RooMsgService::instance().setGlobalKillBelow(msglevel); // set it back
345  cout << "-----------------------------------------" << endl;
346  cout << "Part 4" << endl;
347  r1->Print();
348  t.Stop();
349  t.Print();
350  t.Reset();
351  t.Start();
352
353  c->cd(2);
354  HypoTestPlot *p1 = new HypoTestPlot(*r1, 30); // 30 bins, TS is discrete
355  p1->Draw();
356
357  // -------------------------------------------------------------------------
358  // # P A R T 5 : U S I N G H Y B R I D C A L C U L A T O R
359  // # W I T H A N A L T E R N A T I V E T E S T S T A T I S T I C
360  //
361  // A likelihood ratio test statistics should be 1-to-1 with the count x
362  // when the value of b is fixed in the likelihood. This is implemented
363  // by the SimpleLikelihoodRatioTestStat
364
365  SimpleLikelihoodRatioTestStat slrts(*b_model.GetPdf(), *sb_model.GetPdf());
366  slrts.SetNullParameters(*b_model.GetSnapshot());
367  slrts.SetAltParameters(*sb_model.GetSnapshot());
368
369  // HYBRID CALCULATOR
370  HybridCalculator hc2(*data, sb_model, b_model);
371  ToyMCSampler *toymcs2 = (ToyMCSampler *)hc2.GetTestStatSampler();
372  toymcs2->SetNEventsPerToy(1);
373  toymcs2->SetTestStatistic(&slrts);
374  hc2.SetToys(20000, 1000);
375  hc2.ForcePriorNuisanceAlt(*w->pdf("py"));
376  hc2.ForcePriorNuisanceNull(*w->pdf("py"));
377  // if you wanted to use the ad hoc Gaussian prior instead
378  // ~~~
379  // hc2.ForcePriorNuisanceAlt(*w->pdf("gauss_prior"));
380  // hc2.ForcePriorNuisanceNull(*w->pdf("gauss_prior"));
381  // ~~~
382  // if you wanted to use the ad hoc log-normal prior instead
383  // ~~~
384  // hc2.ForcePriorNuisanceAlt(*w->pdf("lognorm_prior"));
385  // hc2.ForcePriorNuisanceNull(*w->pdf("lognorm_prior"));
386  // ~~~
387
388  // enable proof
389  if (pc)
390  toymcs2->SetProofConfig(pc);
391
392  // these lines save current msg level and then kill any messages below ERROR
394  // Get the result
395  HypoTestResult *r2 = hc2.GetHypoTest();
396  cout << "-----------------------------------------" << endl;
397  cout << "Part 5" << endl;
398  r2->Print();
399  t.Stop();
400  t.Print();
401  t.Reset();
402  t.Start();
404
405  c->cd(3);
406  HypoTestPlot *p2 = new HypoTestPlot(*r2, 30); // 30 bins
407  p2->Draw();
408
409  // -----------------------------------------------------------------------------
410  // # P A R T 6 : U S I N G H Y B R I D C A L C U L A T O R W I T H A N A L T E R N A T I V E T E S T
411  // # S T A T I S T I C A N D S I M U L T A N E O U S M O D E L
412  //
413  // If one wants to use a test statistic in which the nuisance parameters
414  // are profiled (in one way or another), then the PDF must constrain b.
415  // Otherwise any observation x can always be explained with s=0 and b=x/tau.
416  //
417  // In this case, one is really thinking about the problem in a
418  // different way. They are considering x,y simultaneously.
419  // and the PDF should be $Pois(x | s+b) * Pois(y | tau b )$
420  // and the set 'obs' should be {x,y}.
421
422  w->defineSet("obsXY", "x,y");
423
424  // create a toy dataset with the x=150, y=100
425  w->var("x")->setVal(150.);
426  w->var("y")->setVal(100.);
427  RooDataSet *dataXY = new RooDataSet("dXY", "dXY", *w->set("obsXY"));
429
430  // now we need new model configs, with PDF="model"
431  ModelConfig b_modelXY("B_modelXY", w);
432  b_modelXY.SetPdf(*w->pdf("model")); // IMPORTANT
433  b_modelXY.SetObservables(*w->set("obsXY"));
434  b_modelXY.SetParametersOfInterest(*w->set("poi"));
435  w->var("s")->setVal(0.0); // IMPORTANT
436  b_modelXY.SetSnapshot(*w->set("poi"));
437
438  // create the alternate (signal+background) ModelConfig with s=50
439  ModelConfig sb_modelXY("S+B_modelXY", w);
440  sb_modelXY.SetPdf(*w->pdf("model")); // IMPORTANT
441  sb_modelXY.SetObservables(*w->set("obsXY"));
442  sb_modelXY.SetParametersOfInterest(*w->set("poi"));
443  w->var("s")->setVal(50.0); // IMPORTANT
444  sb_modelXY.SetSnapshot(*w->set("poi"));
445
446  // without this print, their can be a crash when using PROOF. Strange.
447  // w->Print();
448
449  // Test statistics like the profile likelihood ratio
450  // (or the ratio of profiled likelihoods (Tevatron) or the MLE for s)
451  // will now work, since the nuisance parameter b is constrained by y.
452  // ratio of alt and null likelihoods with background yield profiled.
453  //
454  // NOTE: These are slower because they have to run fits for each toy
455
456  // Tevatron-style Ratio of profiled likelihoods
457  // $Q_Tev = -log L(s=0,\hat\hat{b})/L(s=50,\hat\hat{b})$
458  RatioOfProfiledLikelihoodsTestStat ropl(*b_modelXY.GetPdf(), *sb_modelXY.GetPdf(), sb_modelXY.GetSnapshot());
459  ropl.SetSubtractMLE(false);
460
461  // profile likelihood where alternate is best fit value of signal yield
462  // $\lambda(0) = -log L(s=0,\hat\hat{b})/L(\hat{s},\hat{b})$
463  ProfileLikelihoodTestStat profll(*b_modelXY.GetPdf());
464
465  // just use the maximum likelihood estimate of signal yield
466  // $MLE = \hat{s}$
467  MaxLikelihoodEstimateTestStat mlets(*sb_modelXY.GetPdf(), *w->var("s"));
468
469  // However, it is less clear how to justify the prior used in randomizing
470  // the nuisance parameters (since that is a property of the ensemble,
471  // and y is a property of each toy pseudo experiment. In that case,
472  // one probably wants to consider a different y0 which will be held
473  // constant and the prior $\pi(b) = Pois(y0 | tau b) * \eta(b)$.
474  w->factory("y0[100]");
475  w->factory("Gamma::gamma_y0(b,sum::temp0(y0,1),1,0)");
476  w->factory("Gaussian::gauss_prior_y0(b,y0, expr::sqrty0('sqrt(y0)',y0))");
477
478  // HYBRID CALCULATOR
479  HybridCalculator hc3(*dataXY, sb_modelXY, b_modelXY);
480  ToyMCSampler *toymcs3 = (ToyMCSampler *)hc3.GetTestStatSampler();
481  toymcs3->SetNEventsPerToy(1);
482  toymcs3->SetTestStatistic(&slrts);
483  hc3.SetToys(30000, 1000);
484  hc3.ForcePriorNuisanceAlt(*w->pdf("gamma_y0"));
485  hc3.ForcePriorNuisanceNull(*w->pdf("gamma_y0"));
486  // if you wanted to use the ad hoc Gaussian prior instead
487  // ~~~{.cpp}
488  // hc3.ForcePriorNuisanceAlt(*w->pdf("gauss_prior_y0"));
489  // hc3.ForcePriorNuisanceNull(*w->pdf("gauss_prior_y0"));
490  // ~~~
491
492  // choose fit-based test statistic
493  toymcs3->SetTestStatistic(&profll);
494  // toymcs3->SetTestStatistic(&ropl);
495  // toymcs3->SetTestStatistic(&mlets);
496
497  // enable proof
498  if (pc)
499  toymcs3->SetProofConfig(pc);
500
501  // these lines save current msg level and then kill any messages below ERROR
503  // Get the result
504  HypoTestResult *r3 = hc3.GetHypoTest();
505  cout << "-----------------------------------------" << endl;
506  cout << "Part 6" << endl;
507  r3->Print();
508  t.Stop();
509  t.Print();
510  t.Reset();
511  t.Start();
513
514  c->cd(4);
516  HypoTestPlot *p3 = new HypoTestPlot(*r3, 50); // 50 bins
517  p3->Draw();
518
519  c->SaveAs("zbi.pdf");
520
521  // -----------------------------------------
522  // OUTPUT W/O PROOF (2.66 GHz Intel Core i7)
523  // =========================================
524
525  // -----------------------------------------
526  // Part 2
527  // Hybrid p-value from direct integration = 0.00094165
528  // Z_Gamma Significance = 3.10804
529  // -----------------------------------------
530  //
531  // Part 3
532  // Z_Bi p-value (analytic): 0.00094165
533  // Z_Bi significance (analytic): 3.10804
534  // Real time 0:00:00, CP time 0.610
535
536  // -----------------------------------------
537  // Part 4
538  // Results HybridCalculator_result:
539  // - Null p-value = 0.00115 +/- 0.000228984
540  // - Significance = 3.04848 sigma
541  // - Number of S+B toys: 1000
542  // - Number of B toys: 20000
543  // - Test statistic evaluated on data: 150
544  // - CL_b: 0.99885 +/- 0.000239654
545  // - CL_s+b: 0.476 +/- 0.0157932
546  // - CL_s: 0.476548 +/- 0.0158118
547  // Real time 0:00:07, CP time 7.620
548
549  // -----------------------------------------
550  // Part 5
551  // Results HybridCalculator_result:
552  // - Null p-value = 0.0009 +/- 0.000206057
553  // - Significance = 3.12139 sigma
554  // - Number of S+B toys: 1000
555  // - Number of B toys: 20000
556  // - Test statistic evaluated on data: 10.8198
557  // - CL_b: 0.9991 +/- 0.000212037
558  // - CL_s+b: 0.465 +/- 0.0157726
559  // - CL_s: 0.465419 +/- 0.0157871
560  // Real time 0:00:34, CP time 34.360
561
562  // -----------------------------------------
563  // Part 6
564  // Results HybridCalculator_result:
565  // - Null p-value = 0.000666667 +/- 0.000149021
566  // - Significance = 3.20871 sigma
567  // - Number of S+B toys: 1000
568  // - Number of B toys: 30000
569  // - Test statistic evaluated on data: 5.03388
570  // - CL_b: 0.999333 +/- 0.000149021
571  // - CL_s+b: 0.511 +/- 0.0158076
572  // - CL_s: 0.511341 +/- 0.0158183
573  // Real time 0:05:06, CP time 306.330
574
575  // ---------------------------------------------------------
576  // OUTPUT w/ PROOF (2.66 GHz Intel Core i7, 4 virtual cores)
577  // =========================================================
578
579  // -----------------------------------------
580  // Part 5
581  // Results HybridCalculator_result:
582  // - Null p-value = 0.00075 +/- 0.000173124
583  // - Significance = 3.17468 sigma
584  // - Number of S+B toys: 1000
585  // - Number of B toys: 20000
586  // - Test statistic evaluated on data: 10.8198
587  // - CL_b: 0.99925 +/- 0.000193577
588  // - CL_s+b: 0.454 +/- 0.0157443
589  // - CL_s: 0.454341 +/- 0.0157564
590  // Real time 0:00:16, CP time 0.990
591
592  // -----------------------------------------
593  // Part 6
594  // Results HybridCalculator_result:
595  // - Null p-value = 0.0007 +/- 0.000152699
596  // - Significance = 3.19465 sigma
597  // - Number of S+B toys: 1000
598  // - Number of B toys: 30000
599  // - Test statistic evaluated on data: 5.03388
600  // - CL_b: 0.9993 +/- 0.000152699
601  // - CL_s+b: 0.518 +/- 0.0158011
602  // - CL_s: 0.518363 +/- 0.0158124
603  // Real time 0:01:25, CP time 0.580
604
605  // ----------------------------------
606  // Comparison
607  // ----------------------------------
608  // LEPStatToolsForLHC
609  // https://plone4.fnal.gov:4430/P0/phystat/packages/0703002
610  // Uses Gaussian prior
611  // CL_b = 6.218476e-04, Significance = 3.228665 sigma
612  //
613  // ----------------------------------
614  // Comparison
615  // ----------------------------------
616  // Asymptotic
617  // From the value of the profile likelihood ratio (5.0338)
618  // The significance can be estimated using Wilks's theorem
619  // significance = sqrt(2*profileLR) = 3.1729 sigma
620 }
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:46
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:219
RooCmdArg LineColor(Color_t color)
virtual const RooArgSet * get() const
Definition: RooAbsData.h:87
void Print(const Option_t *="") const
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
Definition: Rtypes.h:64
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
Double_t PValueToSignificance(Double_t pvalue)
returns one-sided significance corresponding to a p-value
Definition: RooStatsUtils.h:51
HypoTestResult is a base class for results from hypothesis tests.
RooFit::MsgLevel globalKillBelow() const
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Definition: ToyMCSampler.h:149
MsgLevel
Verbosity level for RooMsgService::StreamConfig in RooMsgService.
Definition: RooGlobalFunc.h:65
Basic string class.
Definition: TString.h:131
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
Definition: Rtypes.h:64
static RooMsgService & instance()
Return reference to singleton instance.
virtual RooPlot * plotOn(RooPlot *frame, const RooCmdArg &arg1=RooCmdArg::none(), const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none(), const RooCmdArg &arg9=RooCmdArg::none(), const RooCmdArg &arg10=RooCmdArg::none()) const
Definition: RooAbsPdf.h:118
Double_t BinomialWithTauObsP(Double_t nObs, Double_t bExp, Double_t tau)
RooAbsReal * createCdf(const RooArgSet &iset, const RooArgSet &nset=RooArgSet())
Create a cumulative distribution function of this p.d.f in terms of the observables listed in iset...
Definition: RooAbsPdf.cxx:3370
virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i)
Definition: ToyMCSampler.h:184
RooCmdArg Range(const char *rangeName, Bool_t adjustNorm=kTRUE)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
#define ClassDef(name, id)
Definition: Rtypes.h:326
void SetNullParameters(const RooArgSet &nullParameters)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Double_t BinomialWithTauObsZ(Double_t nObs, Double_t bExp, Double_t tau)
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:252
static constexpr double pc
Bool_t defineSet(const char *name, const RooArgSet &aset, Bool_t importMissing=kFALSE)
Define a named RooArgSet with given constituents.
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:89
void SetProofConfig(ProofConfig *pc=NULL)
Definition: ToyMCSampler.h:233
void setGlobalKillBelow(RooFit::MsgLevel level)
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:44
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:71
RooDataSet is a container class to hold unbinned data.
Definition: RooDataSet.h:31
TestStatistic that returns the ratio of profiled likelihoods.
Same purpose as HybridCalculatorOriginal, but different implementation.
RooPlot * frame(const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create a new RooPlot on the heap with a drawing frame initialized for this object, but no plot contents.
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:44
The Canvas class.
Definition: TCanvas.h:31
Namespace for the RooStats classes.
Definition: Asimov.h:20
#define ClassImp(name)
Definition: Rtypes.h:365
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:59
TestStatistic class that returns -log(L[null] / L[alt]) where L is the likelihood.
RooRealVar * var(const char *name) const
Retrieve real-valued variable (RooRealVar) with given name. A null pointer is returned if not found...
MaxLikelihoodEstimateTestStat: TestStatistic that returns maximum likelihood estimate of a specified ...
RooFactoryWSTool & factory()
Return instance to factory tool.
Double_t getRealValue(const char *name, Double_t defVal=0, Bool_t verbose=kFALSE) const
Get value of a RooAbsReal stored in set with given name.
Definition: RooArgSet.cxx:472
const RooArgSet * set(const char *name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
#define c(i)
Definition: RSha256.hxx:101
void Reset()
Definition: TStopwatch.h:52
Definition: Rtypes.h:64
TestStatistic is an interface class to provide a facility for construction test statistics distributi...
Definition: TestStatistic.h:31
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:306
virtual void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: RooPlot.cxx:712
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the &#39;data&#39; argset, to the data set...
Stopwatch class.
Definition: TStopwatch.h:28