Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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)
int type
Definition TGX11.cxx:121
Definition test.py:1
␛[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby␛[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please read http://roofit.sourceforge.net/license.txt
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) 0x5611717c59e0 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) 0x5611717c6290 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();
}
}
const Bool_t kFALSE
Definition RtypesCore.h:92
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition TError.cxx:220
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:187
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:231
char name[80]
Definition TGX11.cxx:110
float * q
#define gROOT
Definition TROOT.h:406
R__EXTERN TSystem * gSystem
Definition TSystem.h:559
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:313
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition RooAbsArg.h:340
Int_t getSize() const
RooAbsArg * first() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:49
virtual Double_t sumEntries() const =0
Return effective number of entries in dataset, i.e., sum all weights.
virtual Bool_t isWeighted() const
Definition RooAbsData.h:126
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Bool_t canBeExtended() const
If true, PDF can provide extended likelihood term.
Definition RooAbsPdf.h:238
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:91
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:29
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:39
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Hypothesis Test Calculator based on the asymptotic formulae for the profile likelihood ratio.
Does a frequentist hypothesis test.
Same purpose as HybridCalculatorOriginal, but different implementation.
Common base class for the Hypothesis Test Calculators.
TestStatSampler * GetTestStatSampler(void) const
Returns instance of TestStatSampler.
virtual HypoTestResult * GetHypoTest() const
inherited methods from HypoTestCalculator interface
This class provides the plots for the result of a study performed with any of the HypoTestCalculatorG...
HypoTestResult is a base class for results from hypothesis tests.
void SetBackgroundAsAlt(Bool_t l=kTRUE)
void SetPValueIsRightTail(Bool_t pr)
virtual Double_t AlternatePValue() const
Return p-value for alternate hypothesis.
virtual Double_t NullPValue() const
Return p-value for null hypothesis.
void Print(const Option_t *="") const
Print out some information about the results Note: use Alt/Null labels for the hypotheses here as the...
SamplingDistribution * GetAltDistribution(void) const
ModelConfig is a simple class that holds configuration information specifying how a model should be u...
Definition ModelConfig.h:30
virtual void SetSnapshot(const RooArgSet &set)
Set parameter values for a particular hypothesis if using a common PDF by saving a snapshot in the wo...
virtual ModelConfig * Clone(const char *name="") const override
clone
Definition ModelConfig.h:54
const RooArgSet * GetParametersOfInterest() const
get RooArgSet containing the parameter of interest (return NULL if not existing)
const RooArgSet * GetNuisanceParameters() const
get RooArgSet containing the nuisance parameters (return NULL if not existing)
const RooArgSet * GetObservables() const
get RooArgSet for observables (return NULL if not existing)
const RooArgSet * GetSnapshot() const
get RooArgSet for parameters for a particular hypothesis (return NULL if not existing)
RooAbsPdf * GetPdf() const
get model PDF (return NULL if pdf has not been specified or does not exist)
RooAbsPdf * GetPriorPdf() const
get parameters prior pdf (return NULL if not existing)
ProfileLikelihoodTestStat is an implementation of the TestStatistic interface that calculates the pro...
virtual void EnableDetailedOutput(bool e=true, bool withErrorsAndPulls=false)
Holds configuration options for proof and proof-lite.
Definition ProofConfig.h:46
TestStatistic that returns the ratio of profiled likelihoods.
void SetLogYaxis(Bool_t ly)
changes plot to log scale on y axis
void Draw(Option_t *options=0)
Draw this plot and all of the elements it contains.
This class simply holds a sampling distribution of some test statistic.
const std::vector< Double_t > & GetSamplingDistribution() const
Get test statistics values.
TestStatistic class that returns -log(L[null] / L[alt]) where L is the likelihood.
void SetNullParameters(const RooArgSet &nullParameters)
void SetAltParameters(const RooArgSet &altParameters)
ToyMCSampler is an implementation of the TestStatSampler interface.
void SetProofConfig(ProofConfig *pc=NULL)
virtual void SetTestStatistic(TestStatistic *testStatistic, unsigned int i)
void SetGenerateBinned(bool binned=true)
virtual void SetNEventsPerToy(const Int_t nevents)
Forces the generation of exactly n events even for extended PDFs.
The RooWorkspace is a persistable container for RooFit projects.
RooAbsData * data(const char *name) const
Retrieve dataset (binned or unbinned) with given name. A null pointer is returned if not found.
void Print(Option_t *opts=0) const
Print contents of the workspace.
TObject * obj(const char *name) const
Return any type of object (RooAbsArg, RooAbsData or generic object) with given name)
RooAbsPdf * pdf(const char *name) const
Retrieve p.d.f (RooAbsPdf) with given name. A null pointer is returned if not found.
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition TFile.h:54
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:3997
void Close(Option_t *option="") override
Close a file.
Definition TFile.cxx:879
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition TNamed.cxx:140
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
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
Basic string class.
Definition TString.h:136
TString & Replace(Ssiz_t pos, Ssiz_t n, const char *s)
Definition TString.h:682
const char * Data() const
Definition TString.h:369
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2331
virtual Bool_t AccessPathName(const char *path, EAccessMode mode=kFileExists)
Returns FALSE if one can access a file using the specified access mode.
Definition TSystem.cxx:1294
double normal_cdf(double x, double sigma=1, double x0=0)
Cumulative distribution function of the normal (Gaussian) distribution (lower tail).
double normal_quantile_c(double z, double sigma)
Inverse ( ) of the cumulative distribution function of the upper tail of the normal (Gaussian) distri...
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
Definition Asimov.h:19
bool SetAllConstant(const RooAbsCollection &coll, bool constant=true)
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
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
Definition file.py:1
Author
Lorenzo Moneta

Definition in file StandardHypoTestDemo.C.