Logo ROOT  
Reference Guide
StandardTestStatDistributionDemo.C File Reference

Detailed Description

StandardTestStatDistributionDemo.C.

View in nbviewer Open in SWAN This simple script plots the sampling distribution of the profile likelihood ratio test statistic based on the input Model File. To do this one needs to specify the value of the parameter of interest that will be used for evaluating the test statistic and the value of the parameters used for generating the toy data. In this case, it uses the upper-limit estimated from the ProfileLikleihoodCalculator, which assumes the asymptotic chi-square distribution for -2 log profile likelihood ratio. Thus, the script is handy for checking to see if the asymptotic approximations are valid. To aid, that comparison, the script overlays a chi-square distribution as well. The most common parameter of interest is a parameter proportional to the signal rate, and often that has a lower-limit of 0, which breaks the standard chi-square distribution. Thus the script allows the parameter to be negative so that the overlay chi-square is the correct asymptotic distribution.

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 ModelConfig ===
Observables: RooArgSet:: = (obs_x_channel1,weightVar,channelCat)
Parameters of Interest: RooArgSet:: = (SigXsecOverSM)
Nuisance Parameters: RooArgSet:: = (alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
Global Observables: RooArgSet:: = (nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
PDF: RooSimultaneous::simPdf[ indexCat=channelCat channel1=model_channel1 ] = 0.174888
[#1] INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Minization -- createNLL: caching constraint set under name CONSTR_OF_PDF_simPdf_FOR_OBS_channelCat:obs_x_channel1 with 4 entries
[#1] INFO:Minization -- Including the following constraint terms in minimization: (alpha_syst2Constraint,alpha_syst3Constraint,gamma_stat_channel1_bin_0_constraint,gamma_stat_channel1_bin_1_constraint)
[#1] INFO:Minization -- The following global observables have been defined: (nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE
[#0] PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL - using Minuit / Migrad with strategy 1
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_simPdf_obsData_with_constr) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minization -- RooMinimizer::optimizeConst: activating const optimization
RooAbsTestStatistic::initSimMode: creating slave calculator #0 for state channel1 (2 dataset entries)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:Fitting -- RooAbsTestStatistic::initSimMode: created 1 slave calculators.
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
[#1] INFO:Minization -- The following expressions have been identified as constant and will be precalculated and cached: (signal_channel1_nominal,background1_channel1_nominal,background2_channel1_nominal)
[#1] INFO:Minization -- The following expressions will be evaluated in cache-and-track mode: (mc_stat_channel1)
[#1] INFO:Minization --
RooFitResult: minimized FCN value: -1033.78, estimated distance to minimum: 2.39183e-08
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
SigXsecOverSM 1.1153e+00 +/- 5.86e-01
alpha_syst2 -8.9016e-03 +/- 9.82e-01
alpha_syst3 1.7918e-02 +/- 9.48e-01
gamma_stat_channel1_bin_0 9.9955e-01 +/- 4.93e-02
gamma_stat_channel1_bin_1 1.0036e+00 +/- 8.01e-02
Warning: lower value for SigXsecOverSM is at limit 0
--------------------------------------
Will generate sampling distribution at SigXsecOverSM = 2.32744
1) 0x556d1ab5baf0 RooRealVar:: SigXsecOverSM = 2.32744 +/- 0.586082 L(-3 - 3) "SigXsecOverSM"
2) 0x556d1ab58e80 RooRealVar:: alpha_syst2 = 0.710796 +/- 0.982182 L(-5 - 5) "alpha_syst2"
3) 0x556d1ab5efc0 RooRealVar:: alpha_syst3 = 0.260222 +/- 0.947589 L(-5 - 5) "alpha_syst3"
4) 0x556d187fb030 RooRealVar:: gamma_stat_channel1_bin_0 = 1.03668 +/- 0.0493173 L(0 - 1.25) "gamma_stat_channel1_bin_0"
5) 0x556d18622270 RooRealVar:: gamma_stat_channel1_bin_1 = 1.05157 +/- 0.080065 L(0 - 1.5) "gamma_stat_channel1_bin_1"
[#0] PROGRESS:Generation -- generated toys: 500 / 1000
#include "TFile.h"
#include "TROOT.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TSystem.h"
#include "TF1.h"
#include "TSystem.h"
#include "RooWorkspace.h"
#include "RooAbsData.h"
using namespace RooFit;
using namespace RooStats;
bool useProof = false; // flag to control whether to use Proof
int nworkers = 0; // number of workers (default use all available cores)
// -------------------------------------------------------
// The actual macro
void StandardTestStatDistributionDemo(const char *infile = "", const char *workspaceName = "combined",
const char *modelConfigName = "ModelConfig", const char *dataName = "obsData")
{
// the number of toy MC used to generate the distribution
int nToyMC = 1000;
// The parameter below is needed for asymptotic distribution to be chi-square,
// but set to false if your model is not numerically stable if mu<0
bool allowNegativeMu = true;
// -------------------------------------------------------
// 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_combined_GaussExample_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;
#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;
}
// -------------------------------------------------------
// Now get the data and workspace
// get the workspace out of the file
RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
if (!w) {
cout << "workspace not found" << endl;
return;
}
// get the modelConfig out of the file
ModelConfig *mc = (ModelConfig *)w->obj(modelConfigName);
// get the modelConfig 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;
}
mc->Print();
// -------------------------------------------------------
// Now find the upper limit based on the asymptotic results
ProfileLikelihoodCalculator plc(*data, *mc);
LikelihoodInterval *interval = plc.GetInterval();
double plcUpperLimit = interval->UpperLimit(*firstPOI);
delete interval;
cout << "\n\n--------------------------------------" << endl;
cout << "Will generate sampling distribution at " << firstPOI->GetName() << " = " << plcUpperLimit << endl;
int nPOI = mc->GetParametersOfInterest()->getSize();
if (nPOI > 1) {
cout << "not sure what to do with other parameters of interest, but here are their values" << endl;
}
// -------------------------------------------------------
// create the test stat sampler
// to avoid effects from boundary and simplify asymptotic comparison, set min=-max
if (allowNegativeMu)
firstPOI->setMin(-1 * firstPOI->getMax());
// temporary RooArgSet
RooArgSet poi;
// create and configure the ToyMCSampler
ToyMCSampler sampler(ts, nToyMC);
sampler.SetPdf(*mc->GetPdf());
sampler.SetObservables(*mc->GetObservables());
if (!mc->GetPdf()->canBeExtended() && (data->numEntries() == 1)) {
cout << "tell it to use 1 event" << endl;
sampler.SetNEventsPerToy(1);
}
firstPOI->setVal(plcUpperLimit); // set POI value for generation
sampler.SetParametersForTestStat(*mc->GetParametersOfInterest()); // set POI value for evaluation
if (useProof) {
ProofConfig pc(*w, nworkers, "", false);
sampler.SetProofConfig(&pc); // enable proof
}
firstPOI->setVal(plcUpperLimit);
RooArgSet allParameters;
allParameters.add(*mc->GetParametersOfInterest());
allParameters.add(*mc->GetNuisanceParameters());
allParameters.Print("v");
SamplingDistribution *sampDist = sampler.GetSamplingDistribution(allParameters);
plot.AddSamplingDistribution(sampDist);
plot.GetTH1F(sampDist)->GetYaxis()->SetTitle(
Form("f(-log #lambda(#mu=%.2f) | #mu=%.2f)", plcUpperLimit, plcUpperLimit));
plot.SetAxisTitle(Form("-log #lambda(#mu=%.2f)", plcUpperLimit));
TCanvas *c1 = new TCanvas("c1");
c1->SetLogy();
plot.Draw();
double min = plot.GetTH1F(sampDist)->GetXaxis()->GetXmin();
double max = plot.GetTH1F(sampDist)->GetXaxis()->GetXmax();
TF1 *f = new TF1("f", Form("2*ROOT::Math::chisquared_pdf(2*x,%d,0)", nPOI), min, max);
f->Draw("same");
c1->SaveAs("standard_test_stat_distribution.pdf");
}
Author
Kyle Cranmer

Definition in file StandardTestStatDistributionDemo.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
RooStats::ModelConfig::GetObservables
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
Definition: ModelConfig.h:261
RooWorkspace.h
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:216
TH1F.h
PointSetInterval.h
f
#define f(i)
Definition: RSha256.hxx:122
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:176
RooRealVar::setMin
void setMin(const char *name, Double_t value)
Set minimum of name range to given value.
Definition: RooRealVar.cxx:420
RooAbsData
Definition: RooAbsData.h:46
RooStats::SamplingDistribution
Definition: SamplingDistribution.h:39
RooAbsRealLValue::getMax
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
Definition: RooAbsRealLValue.h:89
RooStats::ProfileLikelihoodCalculator
Definition: ProfileLikelihoodCalculator.h:28
Form
char * Form(const char *fmt,...)
ConfidenceBelt.h
SamplingDistPlot.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
ToyMCSampler.h
TGeant4Unit::pc
static constexpr double pc
Definition: TGeant4SystemOfUnits.h:136
RooStats::SamplingDistPlot
Definition: SamplingDistPlot.h:38
RooStats::ToyMCSampler::SetParametersForTestStat
virtual void SetParametersForTestStat(const RooArgSet &nullpoi)
Definition: ToyMCSampler.h:163
RooStats::ToyMCSampler::GetSamplingDistribution
virtual SamplingDistribution * GetSamplingDistribution(RooArgSet &paramPoint)
Definition: ToyMCSampler.cxx:303
RooArgSet::add
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:88
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
RooStats::SamplingDistPlot::GetTH1F
TH1F * GetTH1F(const SamplingDistribution *samplDist=NULL)
Returns the TH1F associated with the give SamplingDistribution.
Definition: SamplingDistPlot.cxx:577
TFile.h
TROOT.h
RooStats::ModelConfig::Print
virtual void Print(Option_t *option="") const override
overload the print method
Definition: ModelConfig.cxx:119
RooStats::ToyMCSampler::SetNEventsPerToy
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Definition: ToyMCSampler.h:157
LikelihoodInterval.h
TAxis::GetXmin
Double_t GetXmin() const
Definition: TAxis.h:133
TSystem.h
TH1::GetYaxis
TAxis * GetYaxis()
Definition: TH1.h:318
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
RooStats::ProofConfig
Definition: ProofConfig.h:52
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
RooStats::ProfileLikelihoodCalculator::GetInterval
virtual LikelihoodInterval * GetInterval() const
Return a likelihood interval.
Definition: ProfileLikelihoodCalculator.cxx:223
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
TNamed::SetTitle
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
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
FeldmanCousins.h
RooStats::ToyMCSampler::SetObservables
virtual void SetObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:181
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
RooStats::ToyMCSampler::SetProofConfig
void SetProofConfig(ProofConfig *pc=NULL)
Definition: ToyMCSampler.h:241
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
Definition: Asimov.h:19
file
Definition: file.py:1
RooStats::ToyMCSampler
Definition: ToyMCSampler.h:79
RooStats::ModelConfig::GetGlobalObservables
const RooArgSet * GetGlobalObservables() const
get RooArgSet for global observables (return NULL if not existing)
Definition: ModelConfig.h:267
ProfileLikelihoodCalculator.h
RooStats::SamplingDistPlot::SetAxisTitle
void SetAxisTitle(char *varName)
Definition: SamplingDistPlot.h:85
RooAbsCollection::Print
virtual void Print(Option_t *options=0) const
This method must be overridden when a class wants to print itself.
Definition: RooAbsCollection.h:199
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:53
RooStats::ProfileLikelihoodTestStat
Definition: ProfileLikelihoodTestStat.h:38
RooStats::ModelConfig::GetNuisanceParameters
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:252
TAxis::GetXmax
Double_t GetXmax() const
Definition: TAxis.h:134
TF1
1-Dim function class
Definition: TF1.h:212
RooRealVar
Definition: RooRealVar.h:36
TH1::GetXaxis
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:317
RooStats::ModelConfig
Definition: ModelConfig.h:36
RooStats::ToyMCSampler::SetPdf
virtual void SetPdf(RooAbsPdf &pdf)
Definition: ToyMCSampler.h:168
RooStats::SamplingDistPlot::AddSamplingDistribution
Double_t AddSamplingDistribution(const SamplingDistribution *samplingDist, Option_t *drawOptions="NORMALIZE HIST")
adds the sampling distribution and returns the scale factor
Definition: SamplingDistPlot.cxx:104
RooStats::LikelihoodInterval
Definition: LikelihoodInterval.h:34
RooStats::ToyMCSampler::SetGlobalObservables
virtual void SetGlobalObservables(const RooArgSet &o)
Definition: ToyMCSampler.h:183
RooAbsPdf::canBeExtended
Bool_t canBeExtended() const
Definition: RooAbsPdf.h:236
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooArgSet
Definition: RooArgSet.h:28
gROOT
#define gROOT
Definition: TROOT.h:406
RooStats::LikelihoodInterval::UpperLimit
Double_t UpperLimit(const RooRealVar &param)
return the upper bound of the interval on a given parameter
Definition: LikelihoodInterval.h:69
ProfileLikelihoodTestStat.h
c1
return c1
Definition: legend1.C:41