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 hypotheiss 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:120
RooCmdArg Asimov(Bool_t flag=kTRUE)
Double_t Freq(Double_t x)
Computation of the normal frequency function freq(x).
Definition: TMath.cxx:265
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) 0x5611cdb023c0 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) 0x5611cdb0b500 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
SimpleLikelihoodRatioTestStat *slrts = new SimpleLikelihoodRatioTestStat(*bModel->GetPdf(), *sbModel->GetPdf());
// 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);
ProfileLikelihoodTestStat *profll = new ProfileLikelihoodTestStat(*bModel->GetPdf());
RatioOfProfiledLikelihoodsTestStat *ropl =
new RatioOfProfiledLikelihoodsTestStat(*bModel->GetPdf(), *sbModel->GetPdf(), sbModel->GetSnapshot());
ropl->SetSubtractMLE(false);
if (testStatType == 3)
profll->SetOneSidedDiscovery(1);
profll->SetPrintLevel(printLevel);
if (enableDetOutput) {
slrts->EnableDetailedOutput();
profll->EnableDetailedOutput();
ropl->EnableDetailedOutput();
}
/* profll.SetReuseNLL(mOptimize);*/
/* slrts.SetReuseNLL(mOptimize);*/
/* ropl.SetReuseNLL(mOptimize);*/
AsymptoticCalculator::SetPrintLevel(printLevel);
HypoTestCalculatorGeneric *hypoCalc = 0;
// 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 =
(bModel->GetNuisanceParameters()) ? bModel->GetNuisanceParameters() : sbModel->GetNuisanceParameters();
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->SetPValueIsRightTail(true);
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) {
SamplingDistribution *altDist = htr->GetAltDistribution();
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:88
void Info(const char *location, const char *msgfmt,...)
void Error(const char *location, const char *msgfmt,...)
void Warning(const char *location, const char *msgfmt,...)
char name[80]
Definition: TGX11.cxx:109
float * q
Definition: THbookFile.cxx:87
#define gROOT
Definition: TROOT.h:415
R__EXTERN TSystem * gSystem
Definition: TSystem.h:560
RooArgSet * getObservables(const RooArgSet &set, Bool_t valueOnly=kTRUE) const
Return the observables of this pdf given a set of observables.
Definition: RooAbsArg.h:240
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:281
Int_t getSize() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:39
virtual Double_t sumEntries() const =0
virtual Bool_t isWeighted() const
Definition: RooAbsData.h:99
virtual Int_t numEntries() const
Definition: RooAbsData.cxx:307
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:28
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:35
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:278
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
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:48
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:3923
void Close(Option_t *option="") override
Close a file.
Definition: TFile.cxx:856
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131
const char * Data() const
Definition: TString.h:364
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
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:1287
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:20
bool SetAllConstant(const RooAbsCollection &coll, bool constant=true)
Definition: RooStatsUtils.h:88
RooAbsPdf * MakeNuisancePdf(RooAbsPdf &pdf, const RooArgSet &observables, const char *name)
static constexpr double pc
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:1190
Definition: file.py:1
Author
Lorenzo Moneta

Definition in file StandardHypoTestDemo.C.