Logo ROOT  
Reference Guide
StandardHypoTestDemo.C File Reference

Detailed Description

View in nbviewer Open in SWAN Standard tutorial macro for hypothesis test (for computing the discovery significance) using all RooStats hypothesis tests calculators and test statistics.

Usage:

root>.L StandardHypoTestDemo.C
root> StandardHypoTestDemo("fileName","workspace name","S+B modelconfig name","B model name","data set
name",calculator type, test statistic type, number of toys)
type = 0 Freq calculator
type = 1 Hybrid calculator
type = 2 Asymptotic calculator
type = 3 Asymptotic calculator using nominal Asimov data sets (not using fitted parameter values but nominal ones)
testStatType = 0 LEP
= 1 Tevatron
= 2 Profile Likelihood
= 3 Profile Likelihood one sided (i.e. = 0 if mu_hat < 0)
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
RooWorkspace(combined) combined contents
variables
---------
(Lumi,SigXsecOverSM,alpha_syst1,alpha_syst2,alpha_syst3,binWidth_obs_x_channel1_0,binWidth_obs_x_channel1_1,binWidth_obs_x_channel1_2,channelCat,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1,nominalLumi,obs_x_channel1,weightVar)
p.d.f.s
-------
RooGaussian::alpha_syst1Constraint[ x=alpha_syst1 mean=nom_alpha_syst1 sigma=1 ] = 1
RooGaussian::alpha_syst2Constraint[ x=alpha_syst2 mean=nom_alpha_syst2 sigma=1 ] = 1
RooGaussian::alpha_syst3Constraint[ x=alpha_syst3 mean=nom_alpha_syst3 sigma=1 ] = 1
RooRealSumPdf::channel1_model[ binWidth_obs_x_channel1_0 * L_x_signal_channel1_overallSyst_x_Exp + binWidth_obs_x_channel1_1 * L_x_background1_channel1_overallSyst_x_StatUncert + binWidth_obs_x_channel1_2 * L_x_background2_channel1_overallSyst_x_StatUncert ] = 220
RooPoisson::gamma_stat_channel1_bin_0_constraint[ x=nom_gamma_stat_channel1_bin_0 mean=gamma_stat_channel1_bin_0_poisMean ] = 0.019943
RooPoisson::gamma_stat_channel1_bin_1_constraint[ x=nom_gamma_stat_channel1_bin_1 mean=gamma_stat_channel1_bin_1_poisMean ] = 0.039861
RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.1 ] = 1
RooProdPdf::model_channel1[ lumiConstraint * alpha_syst1Constraint * alpha_syst2Constraint * alpha_syst3Constraint * gamma_stat_channel1_bin_0_constraint * gamma_stat_channel1_bin_1_constraint * channel1_model(obs_x_channel1) ] = 0.174888
RooSimultaneous::simPdf[ indexCat=channelCat channel1=model_channel1 ] = 0.174888
functions
--------
RooProduct::L_x_background1_channel1_overallSyst_x_StatUncert[ Lumi * background1_channel1_overallSyst_x_StatUncert ] = 0
RooProduct::L_x_background2_channel1_overallSyst_x_StatUncert[ Lumi * background2_channel1_overallSyst_x_StatUncert ] = 100
RooProduct::L_x_signal_channel1_overallSyst_x_Exp[ Lumi * signal_channel1_overallSyst_x_Exp ] = 10
RooStats::HistFactory::FlexibleInterpVar::background1_channel1_epsilon[ paramList=(alpha_syst2) ] = 1
RooHistFunc::background1_channel1_nominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 0
RooProduct::background1_channel1_overallSyst_x_Exp[ background1_channel1_nominal * background1_channel1_epsilon ] = 0
RooProduct::background1_channel1_overallSyst_x_StatUncert[ mc_stat_channel1 * background1_channel1_overallSyst_x_Exp ] = 0
RooStats::HistFactory::FlexibleInterpVar::background2_channel1_epsilon[ paramList=(alpha_syst3) ] = 1
RooHistFunc::background2_channel1_nominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 100
RooProduct::background2_channel1_overallSyst_x_Exp[ background2_channel1_nominal * background2_channel1_epsilon ] = 100
RooProduct::background2_channel1_overallSyst_x_StatUncert[ mc_stat_channel1 * background2_channel1_overallSyst_x_Exp ] = 100
RooProduct::gamma_stat_channel1_bin_0_poisMean[ gamma_stat_channel1_bin_0 * gamma_stat_channel1_bin_0_tau ] = 400
RooProduct::gamma_stat_channel1_bin_1_poisMean[ gamma_stat_channel1_bin_1 * gamma_stat_channel1_bin_1_tau ] = 100
ParamHistFunc::mc_stat_channel1[ ] = 1
RooStats::HistFactory::FlexibleInterpVar::signal_channel1_epsilon[ paramList=(alpha_syst1) ] = 1
RooHistFunc::signal_channel1_nominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 10
RooProduct::signal_channel1_overallNorm_x_sigma_epsilon[ SigXsecOverSM * signal_channel1_epsilon ] = 1
RooProduct::signal_channel1_overallSyst_x_Exp[ signal_channel1_nominal * signal_channel1_overallNorm_x_sigma_epsilon ] = 10
datasets
--------
RooDataSet::asimovData(obs_x_channel1,weightVar,channelCat)
RooDataSet::obsData(channelCat,obs_x_channel1)
embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::signal_channel1nominalDHist(obs_x_channel1)
RooDataHist::background1_channel1nominalDHist(obs_x_channel1)
RooDataHist::background2_channel1nominalDHist(obs_x_channel1)
parameter snapshots
-------------------
NominalParamValues = (nom_alpha_syst2=0[C],nom_alpha_syst3=0[C],nom_gamma_stat_channel1_bin_0=400[C],nom_gamma_stat_channel1_bin_1=100[C],weightVar=0,obs_x_channel1=1.75,Lumi=1[C],nominalLumi=1[C],alpha_syst1=0[C],nom_alpha_syst1=0[C],alpha_syst2=0,alpha_syst3=0,gamma_stat_channel1_bin_0=1 +/- 0.05,gamma_stat_channel1_bin_1=1 +/- 0.1,SigXsecOverSM=1,binWidth_obs_x_channel1_0=2[C],binWidth_obs_x_channel1_1=2[C],binWidth_obs_x_channel1_2=2[C])
named sets
----------
ModelConfig_GlobalObservables:(nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
ModelConfig_NuisParams:(alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
ModelConfig_Observables:(obs_x_channel1,weightVar,channelCat)
ModelConfig_POI:(SigXsecOverSM)
globalObservables:(nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
observables:(obs_x_channel1,weightVar,channelCat)
generic objects
---------------
RooStats::ModelConfig::ModelConfig
=== Using the following for ModelConfigB_only ===
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
Snapshot:
1) 0x55c7d75a5c40 RooRealVar:: SigXsecOverSM = 0 L(0 - 3) "SigXsecOverSM"
=== 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
Snapshot:
1) 0x55c7d75a6430 RooRealVar:: SigXsecOverSM = 1 L(0 - 3) "SigXsecOverSM"
[#0] PROGRESS:Generation -- Test Statistic on data: 1.77404
[#1] INFO:InputArguments -- Profiling conditional MLEs for Null.
[#1] INFO:InputArguments -- Using a ToyMCSampler. Now configuring for Null.
[#0] PROGRESS:Generation -- generated toys: 500 / 5000
[#0] PROGRESS:Generation -- generated toys: 1000 / 5000
[#0] PROGRESS:Generation -- generated toys: 1500 / 5000
[#0] PROGRESS:Generation -- generated toys: 2000 / 5000
[#0] PROGRESS:Generation -- generated toys: 2500 / 5000
[#0] PROGRESS:Generation -- generated toys: 3000 / 5000
[#0] PROGRESS:Generation -- generated toys: 3500 / 5000
[#0] PROGRESS:Generation -- generated toys: 4000 / 5000
[#0] PROGRESS:Generation -- generated toys: 4500 / 5000
[#1] INFO:InputArguments -- Profiling conditional MLEs for Alt.
[#1] INFO:InputArguments -- Using a ToyMCSampler. Now configuring for Alt.
[#0] PROGRESS:Generation -- generated toys: 500 / 1250
[#0] PROGRESS:Generation -- generated toys: 1000 / 1250
Results HypoTestCalculator_result:
- Null p-value = 0.0304 +/- 0.002428
- Significance = 1.87495 +/- 0.0352942 sigma
- Number of Alt toys: 1250
- Number of Null toys: 5000
- Test statistic evaluated on data: 1.77404
- CL_b: 0.0304 +/- 0.002428
- CL_s+b: 0.4328 +/- 0.0140138
- CL_s: 14.2368 +/- 1.22696
Expected p -value and significance at -2 sigma = 0.839 significance -0.990356 sigma
Expected p -value and significance at -1 sigma = 0.2238 significance 0.759422 sigma
Expected p -value and significance at 0 sigma = 0.0434 significance 1.71252 sigma
Expected p -value and significance at 1 sigma = 0.0028 significance 2.77033 sigma
Expected p -value and significance at 2 sigma = 0.0002 significance 3.54008 sigma
#include "TFile.h"
#include "RooWorkspace.h"
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooRandom.h"
#include "TGraphErrors.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TSystem.h"
#include "TROOT.h"
#include <cassert>
using namespace RooFit;
using namespace RooStats;
struct HypoTestOptions {
bool noSystematics = false; // force all systematics to be off (i.e. set all nuisance parameters as constant)
double nToysRatio = 4; // ratio Ntoys Null/ntoys ALT
double poiValue = -1; // change poi snapshot value for S+B model (needed for expected p0 values)
int printLevel = 0;
bool generateBinned = false; // for binned generation
bool useProof = false; // use Proof
bool enableDetailedOutput = false; // for detailed output
};
HypoTestOptions optHT;
void StandardHypoTestDemo(const char *infile = "", const char *workspaceName = "combined",
const char *modelSBName = "ModelConfig", const char *modelBName = "",
const char *dataName = "obsData", int calcType = 0, /* 0 freq 1 hybrid, 2 asymptotic */
int testStatType = 3, /* 0 LEP, 1 TeV, 2 LHC, 3 LHC - one sided*/
int ntoys = 5000, bool useNC = false, const char *nuisPriorName = 0)
{
bool noSystematics = optHT.noSystematics;
double nToysRatio = optHT.nToysRatio; // ratio Ntoys Null/ntoys ALT
double poiValue = optHT.poiValue; // change poi snapshot value for S+B model (needed for expected p0 values)
int printLevel = optHT.printLevel;
bool generateBinned = optHT.generateBinned; // for binned generation
bool useProof = optHT.useProof; // use Proof
bool enableDetOutput = optHT.enableDetailedOutput;
// Other Parameter to pass in tutorial
// apart from standard for filename, ws, modelconfig and data
// type = 0 Freq calculator
// type = 1 Hybrid calculator
// type = 2 Asymptotic calculator
// testStatType = 0 LEP
// = 1 Tevatron
// = 2 Profile Likelihood
// = 3 Profile Likelihood one sided (i.e. = 0 if mu < mu_hat)
// ntoys: number of toys to use
// useNumberCounting: set to true when using number counting events
// nuisPriorName: name of prior for the nuisance. This is often expressed as constraint term in the global model
// It is needed only when using the HybridCalculator (type=1)
// If not given by default the prior pdf from ModelConfig is used.
// extra options are available as global parameters of the macro. They major ones are:
// generateBinned generate binned data sets for toys (default is false) - be careful not to activate with
// a too large (>=3) number of observables
// nToyRatio ratio of S+B/B toys (default is 2)
// printLevel
// disable - can cause some problems
// ToyMCSampler::SetAlwaysUseMultiGen(true);
SimpleLikelihoodRatioTestStat::SetAlwaysReuseNLL(true);
ProfileLikelihoodTestStat::SetAlwaysReuseNLL(true);
RatioOfProfiledLikelihoodsTestStat::SetAlwaysReuseNLL(true);
// RooRandom::randomGenerator()->SetSeed(0);
// to change minimizers
// ~~~{.bash}
// ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);
// ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
// ROOT::Math::MinimizerOptions::SetDefaultTolerance(1);
// ~~~
// -------------------------------------------------------
// 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 but not found, quit
if (!file) {
cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
return;
}
// -------------------------------------------------------
// Tutorial starts here
// -------------------------------------------------------
// get the workspace out of the file
RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
if (!w) {
cout << "workspace not found" << endl;
return;
}
w->Print();
// get the modelConfig out of the file
ModelConfig *sbModel = (ModelConfig *)w->obj(modelSBName);
// get the modelConfig out of the file
RooAbsData *data = w->data(dataName);
// make sure ingredients are found
if (!data || !sbModel) {
w->Print();
cout << "data or ModelConfig was not found" << endl;
return;
}
// make b model
ModelConfig *bModel = (ModelConfig *)w->obj(modelBName);
// case of no systematics
// remove nuisance parameters from model
if (noSystematics) {
const RooArgSet *nuisPar = sbModel->GetNuisanceParameters();
if (nuisPar && nuisPar->getSize() > 0) {
std::cout << "StandardHypoTestInvDemo"
<< " - Switch off all systematics by setting them constant to their initial values" << std::endl;
}
if (bModel) {
const RooArgSet *bnuisPar = bModel->GetNuisanceParameters();
if (bnuisPar)
}
}
if (!bModel) {
Info("StandardHypoTestInvDemo", "The background model %s does not exist", modelBName);
Info("StandardHypoTestInvDemo", "Copy it from ModelConfig %s and set POI to zero", modelSBName);
bModel = (ModelConfig *)sbModel->Clone();
bModel->SetName(TString(modelSBName) + TString("B_only"));
RooRealVar *var = dynamic_cast<RooRealVar *>(bModel->GetParametersOfInterest()->first());
if (!var)
return;
double oldval = var->getVal();
var->setVal(0);
// bModel->SetSnapshot( RooArgSet(*var, *w->var("lumi")) );
bModel->SetSnapshot(RooArgSet(*var));
var->setVal(oldval);
}
if (!sbModel->GetSnapshot() || poiValue > 0) {
Info("StandardHypoTestDemo", "Model %s has no snapshot - make one using model poi", modelSBName);
RooRealVar *var = dynamic_cast<RooRealVar *>(sbModel->GetParametersOfInterest()->first());
if (!var)
return;
double oldval = var->getVal();
if (poiValue > 0)
var->setVal(poiValue);
// sbModel->SetSnapshot( RooArgSet(*var, *w->var("lumi") ) );
sbModel->SetSnapshot(RooArgSet(*var));
if (poiValue > 0)
var->setVal(oldval);
// sbModel->SetSnapshot( *sbModel->GetParametersOfInterest() );
}
// part 1, hypothesis testing
// null parameters must include snapshot of poi plus the nuisance values
RooArgSet nullParams(*bModel->GetSnapshot());
if (bModel->GetNuisanceParameters())
nullParams.add(*bModel->GetNuisanceParameters());
slrts->SetNullParameters(nullParams);
RooArgSet altParams(*sbModel->GetSnapshot());
if (sbModel->GetNuisanceParameters())
altParams.add(*sbModel->GetNuisanceParameters());
slrts->SetAltParameters(altParams);
new RatioOfProfiledLikelihoodsTestStat(*bModel->GetPdf(), *sbModel->GetPdf(), sbModel->GetSnapshot());
ropl->SetSubtractMLE(false);
if (testStatType == 3)
profll->SetPrintLevel(printLevel);
if (enableDetOutput) {
}
/* profll.SetReuseNLL(mOptimize);*/
/* slrts.SetReuseNLL(mOptimize);*/
/* ropl.SetReuseNLL(mOptimize);*/
AsymptoticCalculator::SetPrintLevel(printLevel);
// note here Null is B and Alt is S+B
if (calcType == 0)
hypoCalc = new FrequentistCalculator(*data, *sbModel, *bModel);
else if (calcType == 1)
hypoCalc = new HybridCalculator(*data, *sbModel, *bModel);
else if (calcType == 2)
hypoCalc = new AsymptoticCalculator(*data, *sbModel, *bModel);
if (calcType == 0) {
((FrequentistCalculator *)hypoCalc)->SetToys(ntoys, ntoys / nToysRatio);
if (enableDetOutput)
((FrequentistCalculator *)hypoCalc)->StoreFitInfo(true);
}
if (calcType == 1) {
((HybridCalculator *)hypoCalc)->SetToys(ntoys, ntoys / nToysRatio);
// n. a. yetif (enableDetOutput) ((HybridCalculator*) hypoCalc)->StoreFitInfo(true);
}
if (calcType == 2) {
if (testStatType == 3)
((AsymptoticCalculator *)hypoCalc)->SetOneSidedDiscovery(true);
if (testStatType != 2 && testStatType != 3)
Warning("StandardHypoTestDemo",
"Only the PL test statistic can be used with AsymptoticCalculator - use by default a two-sided PL");
}
// check for nuisance prior pdf in case of nuisance parameters
if (calcType == 1 && (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters())) {
RooAbsPdf *nuisPdf = 0;
if (nuisPriorName)
nuisPdf = w->pdf(nuisPriorName);
// use prior defined first in bModel (then in SbModel)
if (!nuisPdf) {
Info("StandardHypoTestDemo",
"No nuisance pdf given for the HybridCalculator - try to deduce pdf from the model");
if (bModel->GetPdf() && bModel->GetObservables())
nuisPdf = RooStats::MakeNuisancePdf(*bModel, "nuisancePdf_bmodel");
else
nuisPdf = RooStats::MakeNuisancePdf(*sbModel, "nuisancePdf_sbmodel");
}
if (!nuisPdf) {
if (bModel->GetPriorPdf()) {
nuisPdf = bModel->GetPriorPdf();
Info("StandardHypoTestDemo",
"No nuisance pdf given - try to use %s that is defined as a prior pdf in the B model",
nuisPdf->GetName());
} else {
Error("StandardHypoTestDemo", "Cannot run Hybrid calculator because no prior on the nuisance parameter is "
"specified or can be derived");
return;
}
}
assert(nuisPdf);
Info("StandardHypoTestDemo", "Using as nuisance Pdf ... ");
nuisPdf->Print();
const RooArgSet *nuisParams =
RooArgSet *np = nuisPdf->getObservables(*nuisParams);
if (np->getSize() == 0) {
Warning("StandardHypoTestDemo",
"Prior nuisance does not depend on nuisance parameters. They will be smeared in their full range");
}
delete np;
((HybridCalculator *)hypoCalc)->ForcePriorNuisanceAlt(*nuisPdf);
((HybridCalculator *)hypoCalc)->ForcePriorNuisanceNull(*nuisPdf);
}
/* hypoCalc->ForcePriorNuisanceAlt(*sbModel->GetPriorPdf());*/
/* hypoCalc->ForcePriorNuisanceNull(*bModel->GetPriorPdf());*/
ToyMCSampler *sampler = (ToyMCSampler *)hypoCalc->GetTestStatSampler();
if (sampler && (calcType == 0 || calcType == 1)) {
// look if pdf is number counting or extended
if (sbModel->GetPdf()->canBeExtended()) {
if (useNC)
Warning("StandardHypoTestDemo", "Pdf is extended: but number counting flag is set: ignore it ");
} else {
// for not extended pdf
if (!useNC) {
int nEvents = data->numEntries();
Info("StandardHypoTestDemo",
"Pdf is not extended: number of events to generate taken from observed data set is %d", nEvents);
sampler->SetNEventsPerToy(nEvents);
} else {
Info("StandardHypoTestDemo", "using a number counting pdf");
sampler->SetNEventsPerToy(1);
}
}
if (data->isWeighted() && !generateBinned) {
Info("StandardHypoTestDemo", "Data set is weighted, nentries = %d and sum of weights = %8.1f but toy "
"generation is unbinned - it would be faster to set generateBinned to true\n",
data->numEntries(), data->sumEntries());
}
if (generateBinned)
sampler->SetGenerateBinned(generateBinned);
// use PROOF
if (useProof) {
ProofConfig pc(*w, 0, "", kFALSE);
sampler->SetProofConfig(&pc); // enable proof
}
// set the test statistic
if (testStatType == 0)
sampler->SetTestStatistic(slrts);
if (testStatType == 1)
sampler->SetTestStatistic(ropl);
if (testStatType == 2 || testStatType == 3)
sampler->SetTestStatistic(profll);
}
HypoTestResult *htr = hypoCalc->GetHypoTest();
htr->SetBackgroundAsAlt(false);
htr->Print(); // how to get meaningful CLs at this point?
delete sampler;
delete slrts;
delete ropl;
delete profll;
if (calcType != 2) {
HypoTestPlot *plot = new HypoTestPlot(*htr, 100);
plot->SetLogYaxis(true);
plot->Draw();
} else {
std::cout << "Asymptotic results " << std::endl;
}
// look at expected significances
// found median of S+B distribution
if (calcType != 2) {
HypoTestResult htExp("Expected Result");
htExp.Append(htr);
// find quantiles in alt (S+B) distribution
double p[5];
double q[5];
for (int i = 0; i < 5; ++i) {
double sig = -2 + i;
p[i] = ROOT::Math::normal_cdf(sig, 1);
}
std::vector<double> values = altDist->GetSamplingDistribution();
TMath::Quantiles(values.size(), 5, &values[0], q, p, false);
for (int i = 0; i < 5; ++i) {
htExp.SetTestStatisticData(q[i]);
double sig = -2 + i;
std::cout << " Expected p -value and significance at " << sig << " sigma = " << htExp.NullPValue()
<< " significance " << htExp.Significance() << " sigma " << std::endl;
}
} else {
// case of asymptotic calculator
for (int i = 0; i < 5; ++i) {
double sig = -2 + i;
// sigma is inverted here
double pval = AsymptoticCalculator::GetExpectedPValues(htr->NullPValue(), htr->AlternatePValue(), -sig, false);
std::cout << " Expected p -value and significance at " << sig << " sigma = " << pval << " significance "
<< ROOT::Math::normal_quantile_c(pval, 1) << " sigma " << std::endl;
}
}
// write result in a file in case of toys
bool writeResult = (calcType != 2);
if (enableDetOutput) {
writeResult = true;
Info("StandardHypoTestDemo", "Detailed output will be written in output result file");
}
if (htr != NULL && writeResult) {
// write to a file the results
const char *calcTypeName = (calcType == 0) ? "Freq" : (calcType == 1) ? "Hybr" : "Asym";
TString resultFileName = TString::Format("%s_HypoTest_ts%d_", calcTypeName, testStatType);
// strip the / from the filename
TString name = infile;
name.Replace(0, name.Last('/') + 1, "");
resultFileName += name;
TFile *fileOut = new TFile(resultFileName, "RECREATE");
htr->Write();
Info("StandardHypoTestDemo", "HypoTestResult has been written in the file %s", resultFileName.Data());
fileOut->Close();
}
}
Author
Lorenzo Moneta

Definition in file StandardHypoTestDemo.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:249
RooWorkspace.h
RooStats::SimpleLikelihoodRatioTestStat::EnableDetailedOutput
virtual void EnableDetailedOutput(bool e=true)
Definition: SimpleLikelihoodRatioTestStat.h:155
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:257
RooStats::RatioOfProfiledLikelihoodsTestStat::SetSubtractMLE
void SetSubtractMLE(bool subtract)
Definition: RatioOfProfiledLikelihoodsTestStat.h:133
TLine.h
RooStats::HybridCalculator
Same purpose as HybridCalculatorOriginal, but different implementation.
Definition: HybridCalculator.h:22
Warning
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition: TError.cxx:231
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:187
TNamed::SetName
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
RooAbsData
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:49
RooStats::SamplingDistribution
This class simply holds a sampling distribution of some test statistic.
Definition: SamplingDistribution.h:28
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
TString::Data
const char * Data() const
Definition: TString.h:369
RooAbsData::isWeighted
virtual Bool_t isWeighted() const
Definition: RooAbsData.h:126
RooStats::HypoTestResult::NullPValue
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
Definition: HypoTestResult.h:45
RooAbsArg::Print
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:320
NumEventsTestStat.h
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
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:3997
ToyMCSampler.h
TGeant4Unit::pc
static constexpr double pc
Definition: TGeant4SystemOfUnits.h:130
RooStats::SimpleLikelihoodRatioTestStat::SetNullParameters
void SetNullParameters(const RooArgSet &nullParameters)
Definition: SimpleLikelihoodRatioTestStat.h:109
RooStats::ModelConfig::Clone
virtual ModelConfig * Clone(const char *name="") const override
clone
Definition: ModelConfig.h:54
HypoTestInverterPlot.h
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:2333
TCanvas.h
ROOT::Math::normal_cdf
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
Definition: ProbFuncMathCore.cxx:234
TString
Basic string class.
Definition: TString.h:136
RooStats::SimpleLikelihoodRatioTestStat::SetAltParameters
void SetAltParameters(const RooArgSet &altParameters)
Definition: SimpleLikelihoodRatioTestStat.h:117
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
RooDataSet.h
test
Definition: test.py:1
RooStats::ProfileLikelihoodTestStat::EnableDetailedOutput
virtual void EnableDetailedOutput(bool e=true, bool withErrorsAndPulls=false)
Definition: ProfileLikelihoodTestStat.h:113
TFile.h
RooStats::SamplingDistribution::GetSamplingDistribution
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
Definition: SamplingDistribution.h:65
RooStats::SetAllConstant
bool SetAllConstant(const RooAbsCollection &coll, bool constant=true)
Definition: RooStatsUtils.h:89
RooAbsData::sumEntries
virtual Double_t sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
RooStats::HypoTestPlot
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
Definition: HypoTestPlot.h:22
RooStats::HypoTestResult::SetPValueIsRightTail
void SetPValueIsRightTail(Bool_t pr)
Definition: HypoTestResult.cxx:248
q
float * q
Definition: THbookFile.cxx:89
HybridCalculator.h
TROOT.h
RooStats::AsymptoticCalculator
Hypothesis Test Calculator based on the asymptotic formulae for the profile likelihood ratio.
Definition: AsymptoticCalculator.h:27
RooFit::Asimov
RooCmdArg Asimov(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:239
HypoTestInverter.h
RooStats::FrequentistCalculator
Does a frequentist hypothesis test.
Definition: FrequentistCalculator.h:25
RooStats::ToyMCSampler::SetNEventsPerToy
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
Definition: ToyMCSampler.h:151
RooStats::RatioOfProfiledLikelihoodsTestStat
TestStatistic that returns the ratio of profiled likelihoods.
Definition: RatioOfProfiledLikelihoodsTestStat.h:26
TSystem.h
FrequentistCalculator.h
RooStats::ModelConfig::GetPriorPdf
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return NULL if not existing)
Definition: ModelConfig.h:246
RooWorkspace::Print
void Print(Option_t *opts=0) const
Print contents of the workspace.
Definition: RooWorkspace.cxx:2194
ModelConfig.h
RooFit
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition: RooCFunction1Binding.h:29
SimpleLikelihoodRatioTestStat.h
RooAbsPdf.h
RatioOfProfiledLikelihoodsTestStat.h
AsymptoticCalculator.h
RooStats::ProofConfig
Holds configuration options for proof and proof-lite.
Definition: ProofConfig.h:46
RooStats::ProfileLikelihoodTestStat::SetPrintLevel
void SetPrintLevel(Int_t printlevel)
Definition: ProfileLikelihoodTestStat.h:103
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:96
RooStats::HypoTestResult::SetBackgroundAsAlt
void SetBackgroundAsAlt(Bool_t l=kTRUE)
Definition: HypoTestResult.h:90
RooRandom.h
RooArgSet::add
Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Add element to non-owning set.
Definition: RooArgSet.cxx:261
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
RooWorkspace::pdf
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
Definition: RooWorkspace.cxx:1277
RooStats::SimpleLikelihoodRatioTestStat
TestStatistic class that returns -log(L[null] / L[alt]) where L is the likelihood.
Definition: SimpleLikelihoodRatioTestStat.h:24
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:304
RooStats::HypoTestCalculatorGeneric::GetHypoTest
virtual HypoTestResult * GetHypoTest() const
inherited methods from HypoTestCalculator interface
Definition: HypoTestCalculatorGeneric.cxx:106
RooStats::HypoTestCalculatorGeneric
Common base class for the Hypothesis Test Calculators.
Definition: HypoTestCalculatorGeneric.h:34
RooRealVar.h
RooStats::HypoTestResult::AlternatePValue
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
Definition: HypoTestResult.h:48
RooStats::HypoTestResult::GetAltDistribution
SamplingDistribution * GetAltDistribution(void) const
Definition: HypoTestResult.h:71
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
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:234
TMath::Freq
Double_t Freq(Double_t x)
Computation of the normal frequency function freq(x).
Definition: TMath.cxx:265
RooStats::ToyMCSampler::SetTestStatistic
virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i)
Definition: ToyMCSampler.h:186
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:237
RooWorkspace
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
RooStats::ToyMCSampler::SetProofConfig
void SetProofConfig(ProofConfig *pc=NULL)
Definition: ToyMCSampler.h:235
RooStats::HypoTestCalculatorGeneric::GetTestStatSampler
TestStatSampler * GetTestStatSampler(void) const
Returns instance of TestStatSampler.
Definition: HypoTestCalculatorGeneric.h:66
RooAbsArg::getObservables
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition: RooAbsArg.h:294
RooStats::SamplingDistPlot::Draw
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
Definition: SamplingDistPlot.cxx:299
RooStats::ProfileLikelihoodTestStat::SetOneSidedDiscovery
void SetOneSidedDiscovery(Bool_t flag=true)
Definition: ProfileLikelihoodTestStat.h:88
Info
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition: TError.cxx:220
RooStats
Namespace for the RooStats classes.
Definition: Asimov.h:19
file
Definition: file.py:1
RooFit::Hybrid
@ Hybrid
Definition: RooGlobalFunc.h:70
RooStats::ToyMCSampler
ToyMCSampler is an implementation of the TestStatSampler interface.
Definition: ToyMCSampler.h:73
TObject::Write
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:798
TGraphAsymmErrors.h
TFile::Close
void Close(Option_t *option="") override
Close a file.
Definition: TFile.cxx:879
RooStats::HypoTestResult
HypoTestResult is a base class for results from hypothesis tests.
Definition: HypoTestResult.h:22
name
char name[80]
Definition: TGX11.cxx:110
ROOT::Math::normal_quantile_c
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
Definition: QuantFuncMathCore.cxx:126
MaxLikelihoodEstimateTestStat.h
RooStats::RatioOfProfiledLikelihoodsTestStat::EnableDetailedOutput
virtual void EnableDetailedOutput(bool e=true)
Definition: RatioOfProfiledLikelihoodsTestStat.h:74
RooAbsPdf
Definition: RooAbsPdf.h:43
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
RooStats::ProfileLikelihoodTestStat
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
Definition: ProfileLikelihoodTestStat.h:32
RooStats::ModelConfig::GetNuisanceParameters
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
Definition: ModelConfig.h:240
type
int type
Definition: TGX11.cxx:121
TGraphErrors.h
RooStats::ModelConfig::GetSnapshot
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return NULL if not existing)
Definition: ModelConfig.cxx:225
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
TMath::Quantiles
void Quantiles(Int_t n, Int_t nprob, Double_t *x, Double_t *quantiles, Double_t *prob, Bool_t isSorted=kTRUE, Int_t *index=0, Int_t type=7)
Computes sample quantiles, corresponding to the given probabilities.
Definition: TMath.cxx:1183
HypoTestInverterResult.h
RooStats::ModelConfig
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition: ModelConfig.h:30
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
If true, PDF can provide extended likelihood term.
Definition: RooAbsPdf.h:238
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:182
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:29
RooStats::MakeNuisancePdf
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
Definition: RooStatsUtils.cxx:141
gROOT
#define gROOT
Definition: TROOT.h:406
RooStats::ToyMCSampler::SetGenerateBinned
void SetGenerateBinned(bool binned=true)
Definition: ToyMCSampler.h:205
ProfileLikelihoodTestStat.h
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
RooStats::SamplingDistPlot::SetLogYaxis
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
Definition: SamplingDistPlot.h:93