Logo ROOT  
Reference Guide
StandardHypoTestDemo.C File Reference

Detailed Description

Standard tutorial macro for hypothesis test (for computing the discovery significance) using all RooStats hypotheiss tests calculators and test statistics.

View in nbviewer Open in SWAN 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) 0x5649b009ef40 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) 0x5649b009f820 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 constat
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 byt 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 includes 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:261
RooWorkspace.h
RooStats::SimpleLikelihoodRatioTestStat::EnableDetailedOutput
virtual void EnableDetailedOutput(bool e=true)
Definition: SimpleLikelihoodRatioTestStat.h:167
RooRealVar::setVal
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:216
RooStats::RatioOfProfiledLikelihoodsTestStat::SetSubtractMLE
void SetSubtractMLE(bool subtract)
Definition: RatioOfProfiledLikelihoodsTestStat.h:145
TLine.h
RooStats::HybridCalculator
Definition: HybridCalculator.h:28
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:176
TNamed::SetName
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
RooAbsData
Definition: RooAbsData.h:46
RooStats::SamplingDistribution
Definition: SamplingDistribution.h:39
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:115
RooStats::HypoTestResult::NullPValue
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
Definition: HypoTestResult.h:57
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:3946
ToyMCSampler.h
TGeant4Unit::pc
static constexpr double pc
Definition: TGeant4SystemOfUnits.h:136
RooStats::SimpleLikelihoodRatioTestStat::SetNullParameters
void SetNullParameters(const RooArgSet &nullParameters)
Definition: SimpleLikelihoodRatioTestStat.h:121
RooStats::ModelConfig::Clone
virtual ModelConfig * Clone(const char *name="") const override
clone
Definition: ModelConfig.h:66
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
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:2311
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
Definition: TString.h:136
RooStats::SimpleLikelihoodRatioTestStat::SetAltParameters
void SetAltParameters(const RooArgSet &altParameters)
Definition: SimpleLikelihoodRatioTestStat.h:129
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:125
TFile.h
RooStats::SamplingDistribution::GetSamplingDistribution
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
Definition: SamplingDistribution.h:87
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
Definition: HypoTestPlot.h:28
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
Definition: AsymptoticCalculator.h:27
RooFit::Asimov
RooCmdArg Asimov(Bool_t flag=kTRUE)
Definition: RooGlobalFunc.cxx:233
HypoTestInverter.h
RooStats::FrequentistCalculator
Definition: FrequentistCalculator.h:31
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
RooStats::RatioOfProfiledLikelihoodsTestStat
Definition: RatioOfProfiledLikelihoodsTestStat.h:32
TSystem.h
FrequentistCalculator.h
RooStats::ModelConfig::GetPriorPdf
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return NULL if not existing)
Definition: ModelConfig.h:258
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
RooAbsPdf.h
RatioOfProfiledLikelihoodsTestStat.h
AsymptoticCalculator.h
RooStats::ProofConfig
Definition: ProofConfig.h:52
RooStats::ProfileLikelihoodTestStat::SetPrintLevel
void SetPrintLevel(Int_t printlevel)
Definition: ProfileLikelihoodTestStat.h:115
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
RooStats::HypoTestResult::SetBackgroundAsAlt
void SetBackgroundAsAlt(Bool_t l=kTRUE)
Definition: HypoTestResult.h:102
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
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
Definition: SimpleLikelihoodRatioTestStat.h:30
RooAbsData::numEntries
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Definition: RooAbsData.cxx:307
RooStats::HypoTestCalculatorGeneric::GetHypoTest
virtual HypoTestResult * GetHypoTest() const
inherited methods from HypoTestCalculator interface
Definition: HypoTestCalculatorGeneric.cxx:106
RooStats::HypoTestCalculatorGeneric
Definition: HypoTestCalculatorGeneric.h:40
RooRealVar.h
RooStats::HypoTestResult::AlternatePValue
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
Definition: HypoTestResult.h:60
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
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:192
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
RooStats::ToyMCSampler::SetProofConfig
void SetProofConfig(ProofConfig *pc=NULL)
Definition: ToyMCSampler.h:241
RooStats::HypoTestCalculatorGeneric::GetTestStatSampler
TestStatSampler * GetTestStatSampler(void) const
Returns instance of TestStatSampler.
Definition: HypoTestCalculatorGeneric.h:78
RooAbsArg::getObservables
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
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:297
RooStats::ProfileLikelihoodTestStat::SetOneSidedDiscovery
void SetOneSidedDiscovery(Bool_t flag=true)
Definition: ProfileLikelihoodTestStat.h:100
Info
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition: TError.cxx:220
RooStats
Definition: Asimov.h:19
file
Definition: file.py:1
RooFit::Hybrid
@ Hybrid
Definition: RooGlobalFunc.h:70
RooStats::ToyMCSampler
Definition: ToyMCSampler.h:79
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:795
TGraphAsymmErrors.h
TFile::Close
void Close(Option_t *option="") override
Close a file.
Definition: TFile.cxx:876
RooStats::HypoTestResult
Definition: HypoTestResult.h:28
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:86
RooAbsPdf
Definition: RooAbsPdf.h:40
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
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
Definition: RooRealVar.h:36
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
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
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:171
RooArgSet
Definition: RooArgSet.h:28
RooStats::MakeNuisancePdf
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
Definition: RooStatsUtils.cxx:118
gROOT
#define gROOT
Definition: TROOT.h:406
RooStats::ToyMCSampler::SetGenerateBinned
void SetGenerateBinned(bool binned=true)
Definition: ToyMCSampler.h:211
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:105