Notebook

StandardHypoTestDemo¶

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)

Author: Lorenzo Moneta
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Thursday, May 15, 2025 at 06:41 AM.

In [1]:

%%cpp -d

#include "TFile.h"
#include "RooWorkspace.h"
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooStats/ModelConfig.h"
#include "RooRandom.h"
#include "TGraphErrors.h"
#include "TGraphAsymmErrors.h"
#include "TCanvas.h"
#include "TLine.h"
#include "TSystem.h"
#include "TROOT.h"

#include "RooStats/AsymptoticCalculator.h"
#include "RooStats/HybridCalculator.h"
#include "RooStats/FrequentistCalculator.h"
#include "RooStats/ToyMCSampler.h"
#include "RooStats/HypoTestPlot.h"

#include "RooStats/NumEventsTestStat.h"
#include "RooStats/ProfileLikelihoodTestStat.h"
#include "RooStats/SimpleLikelihoodRatioTestStat.h"
#include "RooStats/RatioOfProfiledLikelihoodsTestStat.h"
#include "RooStats/MaxLikelihoodEstimateTestStat.h"

#include "RooStats/HypoTestInverter.h"
#include "RooStats/HypoTestInverterResult.h"
#include "RooStats/HypoTestInverterPlot.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 enableDetailedOutput = false; // for detailed output
};

HypoTestOptions optHT;

Arguments are defined.

In [2]:

const char *infile = "";
const char *workspaceName = "combined";
const char *modelSBName = "ModelConfig";
const char *modelBName = "";
const char *dataName = "obsData";
int calcType = 0;
int testStatType = 3;
int ntoys = 5000;
bool useNC = false;
const char *nuisPriorName = 0;

In [3]:

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 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);

In [4]:

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

In [5]:

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) {
      // 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

In [6]:

TFile *file = TFile::Open(filename);

if input file was specified but not found, quit

In [7]:

if (!file) {
   cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
   return;
}

Tutorial starts here¶

get the workspace out of the file

In [8]:

RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);
if (!w) {
   cout << "workspace not found" << endl;
   return;
}
w->Print();

RooWorkspace(combined) combined contents

variables
---------
(Lumi,SigXsecOverSM,alpha_syst1,alpha_syst2,alpha_syst3,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)

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[ signal_channel1_scaleFactors * signal_channel1_shapes + background1_channel1_scaleFactors * background1_channel1_shapes + background2_channel1_scaleFactors * background2_channel1_shapes ] = 240
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.190787
RooSimultaneous::simPdf[ indexCat=channelCat channel1=model_channel1 ] = 0.190787

functions
--------
RooHistFunc::background1_channel1_Hist_alphanominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 100
RooStats::HistFactory::FlexibleInterpVar::background1_channel1_epsilon[ paramList=(alpha_syst2) ] = 1
RooProduct::background1_channel1_scaleFactors[ background1_channel1_epsilon * Lumi ] = 1
RooProduct::background1_channel1_shapes[ background1_channel1_Hist_alphanominal * mc_stat_channel1 * channel1_model_binWidth ] = 200
RooHistFunc::background2_channel1_Hist_alphanominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 0
RooStats::HistFactory::FlexibleInterpVar::background2_channel1_epsilon[ paramList=(alpha_syst3) ] = 1
RooProduct::background2_channel1_scaleFactors[ background2_channel1_epsilon * Lumi ] = 1
RooProduct::background2_channel1_shapes[ background2_channel1_Hist_alphanominal * mc_stat_channel1 * channel1_model_binWidth ] = 0
RooBinWidthFunction::channel1_model_binWidth[ HistFuncForBinWidth=signal_channel1_Hist_alphanominal HistFuncForBinWidth=signal_channel1_Hist_alphanominal ] = 2
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
RooHistFunc::signal_channel1_Hist_alphanominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 20
RooStats::HistFactory::FlexibleInterpVar::signal_channel1_epsilon[ paramList=(alpha_syst1) ] = 1
RooProduct::signal_channel1_scaleFactors[ signal_channel1_epsilon * SigXsecOverSM * Lumi ] = 1
RooProduct::signal_channel1_shapes[ signal_channel1_Hist_alphanominal * channel1_model_binWidth ] = 40

datasets
--------
RooDataSet::obsData(obs_x_channel1,channelCat)
RooDataSet::asimovData(obs_x_channel1,channelCat)

embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::signal_channel1_Hist_alphanominalDHist(obs_x_channel1)
RooDataHist::background1_channel1_Hist_alphanominalDHist(obs_x_channel1)
RooDataHist::background2_channel1_Hist_alphanominalDHist(obs_x_channel1)

parameter snapshots
-------------------
NominalParamValues = (nominalLumi=1[C],nom_alpha_syst1=0[C],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],obs_x_channel1=1.25,Lumi=1 +/- 0.1[C],alpha_syst1=0 +/- 1[C],alpha_syst2=0 +/- 1,alpha_syst3=0 +/- 1,gamma_stat_channel1_bin_0=1 +/- 0.05,gamma_stat_channel1_bin_1=1 +/- 0.1,SigXsecOverSM=1 +/- 0)

named sets
----------
ModelConfig_GlobalObservables:(nominalLumi,nom_alpha_syst1,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,channelCat)
ModelConfig_POI:(SigXsecOverSM)
globalObservables:(nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
observables:(obs_x_channel1,channelCat)

generic objects
---------------
RooStats::ModelConfig::ModelConfig

get the modelConfig out of the file

In [9]:

ModelConfig *sbModel = (ModelConfig *)w->obj(modelSBName);

get the modelConfig out of the file

In [10]:

RooAbsData *data = w->data(dataName);

input_line_85:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration
 RooAbsData *data = w->data(dataName);
 ^

make sure ingredients are found

In [11]:

if (!data || !sbModel) {
   w->Print();
   cout << "data or ModelConfig was not found" << endl;
   return;
}

input_line_86:2:7: error: reference to 'data' is ambiguous
 if (!data || !sbModel) {
      ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^

make b model

In [12]:

ModelConfig *bModel = (ModelConfig *)w->obj(modelBName);

case of no systematics remove nuisance parameters from model

In [13]:

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;
      RooStats::SetAllConstant(*nuisPar);
   }
   if (bModel) {
      const RooArgSet *bnuisPar = bModel->GetNuisanceParameters();
      if (bnuisPar)
         RooStats::SetAllConstant(*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() );
}

Info in <StandardHypoTestInvDemo>: The background model  does not exist
Info in <StandardHypoTestInvDemo>: Copy it from ModelConfig ModelConfig and set POI to zero
Info in <StandardHypoTestDemo>: Model ModelConfig has no snapshot  - make one using model poi

part 1, hypothesis testing

In [14]:

SimpleLikelihoodRatioTestStat *slrts = new SimpleLikelihoodRatioTestStat(*bModel->GetPdf(), *sbModel->GetPdf());

null parameters must include snapshot of poi plus the nuisance values

In [15]:

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

In [16]:

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");
}

input_line_91:3:42: error: reference to 'data' is ambiguous
   hypoCalc = new FrequentistCalculator(*data, *sbModel, *bModel);
                                         ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_91:5:37: error: reference to 'data' is ambiguous
   hypoCalc = new HybridCalculator(*data, *sbModel, *bModel);
                                    ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_91:7:41: error: reference to 'data' is ambiguous
   hypoCalc = new AsymptoticCalculator(*data, *sbModel, *bModel);
                                        ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^

check for nuisance prior pdf in case of nuisance parameters

In [17]:

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();
   std::unique_ptr<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");
   }

   ((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);

   // 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;
}

input_line_92:57:24: error: reference to 'data' is ambiguous
         int nEvents = data->numEntries();
                       ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_92:67:8: error: reference to 'data' is ambiguous
   if (data->isWeighted() && !generateBinned) {
       ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_92:67:8: error: unknown type name 'data'
   if (data->isWeighted() && !generateBinned) {
       ^
input_line_92:67:12: error: cannot use arrow operator on a type
   if (data->isWeighted() && !generateBinned) {
           ^
input_line_92:70:12: error: reference to 'data' is ambiguous
           data->numEntries(), data->sumEntries());
           ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^
input_line_92:70:32: error: reference to 'data' is ambiguous
           data->numEntries(), data->sumEntries());
                               ^
input_line_85:2:14: note: candidate found by name lookup is 'data'
 RooAbsData *data = w->data(dataName);
             ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:319:5: note: candidate found by name lookup is 'std::data'
    data(initializer_list<_Tp> __il) noexcept
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:290:5: note: candidate found by name lookup is 'std::data'
    data(_Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:300:5: note: candidate found by name lookup is 'std::data'
    data(const _Container& __cont) noexcept(noexcept(__cont.data()))
    ^
/usr/lib/gcc/x86_64-redhat-linux/11/../../../../include/c++/11/bits/range_access.h:310:5: note: candidate found by name lookup is 'std::data'
    data(_Tp (&__array)[_Nm]) noexcept
    ^

look at expected significances found median of S+B distribution

In [18]:

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;
   }
}

input_line_93:28:21: error: cannot initialize an array element of type 'void *' with an rvalue of type 'double (*)(double, double, double, bool, bool)'
      double pval = AsymptoticCalculator::GetExpectedPValues(htr->NullPValue(), htr->AlternatePValue(), -sig, false);
                    ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

write result in a file in case of toys

In [19]:

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();
}

input_line_95:2:3: error: use of undeclared identifier 'htr'
 (htr != __null && ((*(_Bool*)0x7f11c3335000)))
  ^
input_line_95:2:24: error: use of undeclared identifier '_Bool'
 (htr != __null && ((*(_Bool*)0x7f11c3335000)))
                       ^
input_line_95:2:30: error: expected expression
 (htr != __null && ((*(_Bool*)0x7f11c3335000)))
                             ^
Error in <HandleInterpreterException>: Error evaluating expression (htr != __null && ((*(_Bool*)0x7f11c3335000)))
Execution of your code was aborted.

Draw all canvases

In [20]:

gROOT->GetListOfCanvases()->Draw()