Logo ROOT  
Reference Guide
StandardFrequentistDiscovery.C File Reference

Detailed Description

View in nbviewer Open in SWAN StandardFrequentistDiscovery

This is a standard demo that can be used with any ROOT file prepared in the standard way. You specify:

  • name for input ROOT file
  • name of workspace inside ROOT file that holds model and data
  • name of ModelConfig that specifies details for calculator tools
  • name of dataset

With default parameters the macro will attempt to run the standard hist2workspace example and read the ROOT file that it produces.

␛[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby␛[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
=== Using the following for ModelConfigNull ===
Observables: RooArgSet:: = (obs_x_channel1)
Parameters of Interest: RooArgSet:: = (SigXsecOverSM)
Nuisance Parameters: RooArgSet:: = (alpha_syst3,beta_syst2,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
Global Observables: RooArgSet:: = (nom_beta_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
PDF: RooProdPdf::model_channel1[ lumiConstraint * alpha_syst1Constraint * beta_syst2Constraint * alpha_syst3Constraint * gamma_stat_channel1_bin_0_constraint * gamma_stat_channel1_bin_1_constraint * channel1_model(obs_x_channel1) ] = 0.208662
1) 0x564ff6932ed0 RooRealVar:: SigXsecOverSM = 0 L(0 - 3) "SigXsecOverSM"
=== Using the following for ModelConfig ===
Observables: RooArgSet:: = (obs_x_channel1)
Parameters of Interest: RooArgSet:: = (SigXsecOverSM)
Nuisance Parameters: RooArgSet:: = (alpha_syst3,beta_syst2,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
Global Observables: RooArgSet:: = (nom_beta_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
PDF: RooProdPdf::model_channel1[ lumiConstraint * alpha_syst1Constraint * beta_syst2Constraint * alpha_syst3Constraint * gamma_stat_channel1_bin_0_constraint * gamma_stat_channel1_bin_1_constraint * channel1_model(obs_x_channel1) ] = 0.208662
1) 0x564ff69328b0 RooRealVar:: SigXsecOverSM = 1 L(0 - 3) "SigXsecOverSM"
[#0] PROGRESS:Generation -- Test Statistic on data: 2.00338
[#1] INFO:InputArguments -- Profiling conditional MLEs for Null.
[#1] INFO:InputArguments -- Using a ToyMCSampler. Now configuring for Null.
[#0] PROGRESS:Generation -- generated toys: 500 / 1000
[#1] INFO:InputArguments -- Profiling conditional MLEs for Alt.
[#1] INFO:InputArguments -- Using a ToyMCSampler. Now configuring for Alt.
[#0] PROGRESS:Generation -- generated toys: 500 / 1000
Results HypoTestCalculator_result:
- Null p-value = 0.515 +/- 0.0158043
- Significance = -0.0376083 +/- 0.0396435 sigma
- Number of Alt toys: 1000
- Number of Null toys: 1000
- Test statistic evaluated on data: 2.00338
- CL_b: 0.515 +/- 0.0158043
- CL_s+b: 0.767 +/- 0.0133683
- CL_s: 1.48932 +/- 0.0525612
total CPU time: 13.01
total real time: 13.0172
(double) 0.51500000
#include "TFile.h"
#include "TROOT.h"
#include "TH1F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TStopwatch.h"
#include "RooWorkspace.h"
#include "RooAbsData.h"
#include "RooRandom.h"
#include "RooRealSumPdf.h"
#include "TSystem.h"
#include <vector>
using namespace RooFit;
using namespace RooStats;
double StandardFrequentistDiscovery(const char *infile = "", const char *workspaceName = "channel1",
const char *modelConfigNameSB = "ModelConfig", const char *dataName = "obsData",
int toys = 1000, double poiValueForBackground = 0.0, double poiValueForSignal = 1.0)
// The workspace contains the model for s+b. The b model is "autogenerated"
// by copying s+b and setting the one parameter of interest to zero.
// To keep the script simple, multiple parameters of interest or different
// functional forms of the b model are not supported.
// for now, assume there is only one parameter of interest, and these are
// its values:
// -------------------------------------------------------
// First part is just to access a user-defined file
// or create the standard example file if it doesn't exist
const char *filename = "";
if (!strcmp(infile, "")) {
filename = "results/example_channel1_GammaExample_model.root";
bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
// if file does not exists generate with histfactory
if (!fileExist) {
#ifdef _WIN32
cout << "HistFactory file cannot be generated on Windows - exit" << endl;
return -1;
// Normally this would be run on the command line
cout << "will run standard hist2workspace example" << endl;
gROOT->ProcessLine(".! prepareHistFactory .");
gROOT->ProcessLine(".! hist2workspace config/example.xml");
cout << "\n\n---------------------" << endl;
cout << "Done creating example input" << endl;
cout << "---------------------\n\n" << endl;
} else
filename = infile;
// Try to open the file
TFile *file = TFile::Open(filename);
// if input file was specified byt not found, quit
if (!file) {
cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
return -1;
// -------------------------------------------------------
// Tutorial starts here
// -------------------------------------------------------
TStopwatch *mn_t = new TStopwatch;
// get the workspace out of the file
RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
if (!w) {
cout << "workspace not found" << endl;
return -1.0;
// get the modelConfig out of the file
ModelConfig *mc = (ModelConfig *)w->obj(modelConfigNameSB);
// get the data out of the file
RooAbsData *data = w->data(dataName);
// make sure ingredients are found
if (!data || !mc) {
cout << "data or ModelConfig was not found" << endl;
return -1.0;
// create null model
ModelConfig *mcNull = mc->Clone("ModelConfigNull");
// ----------------------------------------------------
// Configure a ProfileLikelihoodTestStat and a SimpleLikelihoodRatioTestStat
// to use simultaneously with ToyMCSampler
// ----------------------------------------------------
// configure the ToyMCImportanceSampler with two test statistics
ToyMCSampler toymcs(*plts, 50);
// Since this tool needs to throw toy MC the PDF needs to be
// extended or the tool needs to know how many entries in a dataset
// per pseudo experiment.
// In the 'number counting form' where the entries in the dataset
// are counts, and not values of discriminating variables, the
// datasets typically only have one entry and the PDF is not
// extended.
if (!mc->GetPdf()->canBeExtended()) {
if (data->numEntries() == 1) {
} else
cout << "Not sure what to do about this model" << endl;
// We can use PROOF to speed things along in parallel
// ProofConfig pc(*w, 2, "user@yourfavoriteproofcluster", false);
ProofConfig pc(*w, 2, "", false);
// toymcs.SetProofConfig(&pc); // enable proof
// instantiate the calculator
FrequentistCalculator freqCalc(*data, *mc, *mcNull, &toymcs);
freqCalc.SetToys(toys, toys); // null toys, alt toys
// Run the calculator and print result
HypoTestResult *freqCalcResult = freqCalc.GetHypoTest();
freqCalcResult->GetNullDistribution()->SetTitle("b only");
double pvalue = freqCalcResult->NullPValue();
// stop timing
cout << "total CPU time: " << mn_t->CpuTime() << endl;
cout << "total real time: " << mn_t->RealTime() << endl;
// plot
TCanvas *c1 = new TCanvas();
HypoTestPlot *plot = new HypoTestPlot(*freqCalcResult, 100, -0.49, 9.51);
// add chi2 to plot
int nPOI = 1;
TF1 *f = new TF1("f", TString::Format("1*ROOT::Math::chisquared_pdf(2*x,%d,0)", nPOI), 0, 20);
plot->AddTF1(f, TString::Format("#chi^{2}(2x,%d)", nPOI));
return pvalue;
#define f(i)
Definition: RSha256.hxx:104
@ kBlack
Definition: Rtypes.h:65
#define gROOT
Definition: TROOT.h:406
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
RooAbsArg * first() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:49
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
Bool_t canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition: RooAbsPdf.h:238
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:118
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:282
Does a frequentist hypothesis test.
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
HypoTestResult is a base class for results from hypothesis tests.
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
SamplingDistribution * GetNullDistribution(void) const
void Print(const Option_t *="") const
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
SamplingDistribution * GetAltDistribution(void) const
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
virtual void SetSnapshot(const RooArgSet &set)
Set parameter values for a particular hypothesis if using a common PDF by saving a snapshot in the wo...
virtual ModelConfig * Clone(const char *name="") const override
Definition: ModelConfig.h:54
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:237
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:234
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
virtual void SetVarName(const char *name)
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:46
void AddTF1(TF1 *f, const char *title=NULL, Option_t *drawOptions="SAME")
add a TF1
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:74
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
void Print(Option_t *opts=0) const
Print contents of the workspace.
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
The Canvas class.
Definition: TCanvas.h:23
1-Dim function class
Definition: TF1.h:213
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3997
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
Stopwatch class.
Definition: TStopwatch.h:28
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition: TString.cxx:2331
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition: TSystem.cxx:1294
return c1
Definition: legend1.C:41
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition: Asimov.h:19
static constexpr double pc
Definition: file.py:1
Sven Kreiss, Kyle Cranmer

Definition in file StandardFrequentistDiscovery.C.