Logo ROOT  
Reference Guide
StandardFrequentistDiscovery.C File Reference

Detailed Description

StandardFrequentistDiscovery.

View in nbviewer Open in SWAN 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.

RooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby
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
Snapshot:
1) 0x55be771db1f0 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
Snapshot:
1) 0x55be771daba0 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: 15.71
total real time: 15.7004
(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;
#endif
// 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;
mn_t->Start();
// 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) {
w->Print();
cout << "data or ModelConfig was not found" << endl;
return -1.0;
}
firstPOI->setVal(poiValueForSignal);
// create null model
ModelConfig *mcNull = mc->Clone("ModelConfigNull");
firstPOI->setVal(poiValueForBackground);
// ----------------------------------------------------
// Configure a ProfileLikelihoodTestStat and a SimpleLikelihoodRatioTestStat
// to use simultaneously with ToyMCSampler
plts->SetOneSidedDiscovery(true);
plts->SetVarName("q_{0}/2");
// ----------------------------------------------------
// 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) {
toymcs.SetNEventsPerToy(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");
freqCalcResult->GetAltDistribution()->SetTitle("s+b");
freqCalcResult->Print();
double pvalue = freqCalcResult->NullPValue();
// stop timing
mn_t->Stop();
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);
plot->SetLogYaxis(true);
// 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);
f->SetLineColor(kBlack);
f->SetLineStyle(7);
plot->AddTF1(f, TString::Format("#chi^{2}(2x,%d)", nPOI));
plot->Draw();
c1->SaveAs("standard_discovery_output.pdf");
return pvalue;
}
Authors
Sven Kreiss, Kyle Cranmer

Definition in file StandardFrequentistDiscovery.C.

RooWorkspace::data
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1368
RooWorkspace.h
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:216
TH1F.h
f
#define f(i)
Definition: RSha256.hxx:122
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:176
RooAbsData
Definition: RooAbsData.h:46
RooNumIntConfig.h
RooStats::HypoTestResult::GetNullDistribution
SamplingDistribution * GetNullDistribution(void) const
Definition: HypoTestResult.h:82
RooStats::HypoTestResult::Print
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: HypoTestResult.cxx:339
RooStats::HypoTestResult::NullPValue
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
Definition: HypoTestResult.h:57
TStopwatch.h
TFile::Open
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:3946
TGeant4Unit::pc
static constexpr double pc
Definition: TGeant4SystemOfUnits.h:136
RooStats::ModelConfig::Clone
virtual ModelConfig * Clone(const char *name="") const override
clone
Definition: ModelConfig.h:66
TString::Format
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:2311
TCanvas.h
TSystem::AccessPathName
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
HypoTestResult.h
ToyMCImportanceSampler.h
TFile.h
RooStats::HypoTestPlot
Definition: HypoTestPlot.h:28
TROOT.h
TStopwatch::RealTime
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
RooArgSet::snapshot
RooArgSet * snapshot(bool deepCopy=true) const
Use RooAbsCollection::snapshot(), but return as RooArgSet.
Definition: RooArgSet.h:134
kBlack
@ kBlack
Definition: Rtypes.h:65
RooStats::FrequentistCalculator
Definition: FrequentistCalculator.h:31
LikelihoodInterval.h
TSystem.h
FrequentistCalculator.h
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
ModelConfig.h
RooFit
Definition: RooCFunction1Binding.h:29
SimpleLikelihoodRatioTestStat.h
RooStats::ProofConfig
Definition: ProofConfig.h:52
RooRandom.h
RooWorkspace::obj
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
Definition: RooWorkspace.cxx:2106
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
TStopwatch::CpuTime
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
RooStats::HypoTestResult::GetAltDistribution
SamplingDistribution * GetAltDistribution(void) const
Definition: HypoTestResult.h:83
TFile
Definition: TFile.h:54
RooStats::ModelConfig::GetPdf
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
Definition: ModelConfig.h:246
SamplingDistribution.h
gSystem
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
RooStats::ModelConfig::GetParametersOfInterest
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
Definition: ModelConfig.h:249
RooWorkspace
Definition: RooWorkspace.h:43
RooAbsData.h
TF1.h
RooStats::SamplingDistPlot::Draw
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: SamplingDistPlot.cxx:297
TCanvas
Definition: TCanvas.h:23
RooStats::ProfileLikelihoodTestStat::SetOneSidedDiscovery
void SetOneSidedDiscovery(Bool_t flag=true)
Definition: ProfileLikelihoodTestStat.h:100
RooStats
Definition: Asimov.h:19
file
Definition: file.py:1
RooStats::ToyMCSampler
Definition: ToyMCSampler.h:79
TStopwatch
Definition: TStopwatch.h:28
RooStats::ProfileLikelihoodTestStat::SetVarName
virtual void SetVarName(const char *name)
Definition: ProfileLikelihoodTestStat.h:149
ProfileLikelihoodCalculator.h
RooStats::HypoTestResult
Definition: HypoTestResult.h:28
RooStats::ProfileLikelihoodTestStat
Definition: ProfileLikelihoodTestStat.h:38
LikelihoodIntervalPlot.h
TF1
1-Dim function class
Definition: TF1.h:212
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
RooRealSumPdf.h
RooRealVar
Definition: RooRealVar.h:36
RooStats::ModelConfig
Definition: ModelConfig.h:36
HypoTestPlot.h
RooStats::ModelConfig::SetSnapshot
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...
Definition: ModelConfig.cxx:209
RooAbsPdf::canBeExtended
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
RooStats::SamplingDistPlot::AddTF1
void AddTF1(TF1 *f, const char *title=NULL, Option_t *drawOptions="SAME")
add a TF1
Definition: SamplingDistPlot.cxx:235
RooArgSet
Definition: RooArgSet.h:28
gROOT
#define gROOT
Definition: TROOT.h:406
ProfileLikelihoodTestStat.h
c1
return c1
Definition: legend1.C:41
RooStats::SamplingDistPlot::SetLogYaxis
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
Definition: SamplingDistPlot.h:105