{
"cells": [
{
"cell_type": "markdown",
"id": "8b6790ec",
"metadata": {},
"source": [
"# StandardHypoTestDemo\n",
"Standard tutorial macro for hypothesis test (for computing the discovery significance) using all\n",
"RooStats hypothesis tests calculators and test statistics.\n",
"\n",
"Usage:\n",
"\n",
"```cpp\n",
"root>.L StandardHypoTestDemo.C\n",
"root> StandardHypoTestDemo(\"fileName\",\"workspace name\",\"S+B modelconfig name\",\"B model name\",\"data set\n",
"name\",calculator type, test statistic type, number of toys)\n",
"\n",
" type = 0 Freq calculator\n",
" type = 1 Hybrid calculator\n",
" type = 2 Asymptotic calculator\n",
" type = 3 Asymptotic calculator using nominal Asimov data sets (not using fitted parameter values but nominal ones)\n",
"\n",
"testStatType = 0 LEP\n",
" = 1 Tevatron\n",
" = 2 Profile Likelihood\n",
" = 3 Profile Likelihood one sided (i.e. = 0 if mu_hat < 0)\n",
"```\n",
"\n",
"\n",
"\n",
"\n",
"**Author:** Lorenzo Moneta \n",
"This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:18 PM."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "f454562e",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:48.933518Z",
"iopub.status.busy": "2024-03-19T19:18:48.933173Z",
"iopub.status.idle": "2024-03-19T19:18:48.980840Z",
"shell.execute_reply": "2024-03-19T19:18:48.979743Z"
}
},
"outputs": [],
"source": [
"%%cpp -d\n",
"\n",
"#include \"TFile.h\"\n",
"#include \"RooWorkspace.h\"\n",
"#include \"RooAbsPdf.h\"\n",
"#include \"RooRealVar.h\"\n",
"#include \"RooDataSet.h\"\n",
"#include \"RooStats/ModelConfig.h\"\n",
"#include \"RooRandom.h\"\n",
"#include \"TGraphErrors.h\"\n",
"#include \"TGraphAsymmErrors.h\"\n",
"#include \"TCanvas.h\"\n",
"#include \"TLine.h\"\n",
"#include \"TSystem.h\"\n",
"#include \"TROOT.h\"\n",
"\n",
"#include \"RooStats/AsymptoticCalculator.h\"\n",
"#include \"RooStats/HybridCalculator.h\"\n",
"#include \"RooStats/FrequentistCalculator.h\"\n",
"#include \"RooStats/ToyMCSampler.h\"\n",
"#include \"RooStats/HypoTestPlot.h\"\n",
"\n",
"#include \"RooStats/NumEventsTestStat.h\"\n",
"#include \"RooStats/ProfileLikelihoodTestStat.h\"\n",
"#include \"RooStats/SimpleLikelihoodRatioTestStat.h\"\n",
"#include \"RooStats/RatioOfProfiledLikelihoodsTestStat.h\"\n",
"#include \"RooStats/MaxLikelihoodEstimateTestStat.h\"\n",
"\n",
"#include \"RooStats/HypoTestInverter.h\"\n",
"#include \"RooStats/HypoTestInverterResult.h\"\n",
"#include \"RooStats/HypoTestInverterPlot.h\"\n",
"\n",
"#include \n",
"\n",
"using namespace RooFit;\n",
"using namespace RooStats;\n",
"\n",
"struct HypoTestOptions {\n",
"\n",
" bool noSystematics = false; // force all systematics to be off (i.e. set all nuisance parameters as constant)\n",
" double nToysRatio = 4; // ratio Ntoys Null/ntoys ALT\n",
" double poiValue = -1; // change poi snapshot value for S+B model (needed for expected p0 values)\n",
" int printLevel = 0;\n",
" bool generateBinned = false; // for binned generation\n",
" bool useProof = false; // use Proof\n",
" bool enableDetailedOutput = false; // for detailed output\n",
"};\n",
"\n",
"HypoTestOptions optHT;"
]
},
{
"cell_type": "markdown",
"id": "c32fb9d1",
"metadata": {},
"source": [
" Arguments are defined. "
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "04131225",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:48.992349Z",
"iopub.status.busy": "2024-03-19T19:18:48.991939Z",
"iopub.status.idle": "2024-03-19T19:18:49.367983Z",
"shell.execute_reply": "2024-03-19T19:18:49.365460Z"
}
},
"outputs": [],
"source": [
"const char *infile = \"\";\n",
"const char *workspaceName = \"combined\";\n",
"const char *modelSBName = \"ModelConfig\";\n",
"const char *modelBName = \"\";\n",
"const char *dataName = \"obsData\";\n",
"int calcType = 0;\n",
"int testStatType = 3;\n",
"int ntoys = 5000;\n",
"bool useNC = false;\n",
"const char *nuisPriorName = 0;"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "c06d62c5",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:49.373225Z",
"iopub.status.busy": "2024-03-19T19:18:49.372838Z",
"iopub.status.idle": "2024-03-19T19:18:49.590644Z",
"shell.execute_reply": "2024-03-19T19:18:49.589515Z"
}
},
"outputs": [],
"source": [
"bool noSystematics = optHT.noSystematics;\n",
"double nToysRatio = optHT.nToysRatio; // ratio Ntoys Null/ntoys ALT\n",
"double poiValue = optHT.poiValue; // change poi snapshot value for S+B model (needed for expected p0 values)\n",
"int printLevel = optHT.printLevel;\n",
"bool generateBinned = optHT.generateBinned; // for binned generation\n",
"bool useProof = optHT.useProof; // use Proof\n",
"bool enableDetOutput = optHT.enableDetailedOutput;"
]
},
{
"cell_type": "markdown",
"id": "caef90a9",
"metadata": {},
"source": [
"Other Parameter to pass in tutorial\n",
"apart from standard for filename, ws, modelconfig and data"
]
},
{
"cell_type": "markdown",
"id": "773bd456",
"metadata": {},
"source": [
"type = 0 Freq calculator\n",
"type = 1 Hybrid calculator\n",
"type = 2 Asymptotic calculator"
]
},
{
"cell_type": "markdown",
"id": "abd4076a",
"metadata": {},
"source": [
"testStatType = 0 LEP\n",
"= 1 Tevatron\n",
"= 2 Profile Likelihood\n",
"= 3 Profile Likelihood one sided (i.e. = 0 if mu < mu_hat)"
]
},
{
"cell_type": "markdown",
"id": "1a4ddd83",
"metadata": {},
"source": [
"ntoys: number of toys to use"
]
},
{
"cell_type": "markdown",
"id": "db23b701",
"metadata": {},
"source": [
"useNumberCounting: set to true when using number counting events"
]
},
{
"cell_type": "markdown",
"id": "b270efd8",
"metadata": {},
"source": [
"nuisPriorName: name of prior for the nuisance. This is often expressed as constraint term in the global model\n",
"It is needed only when using the HybridCalculator (type=1)\n",
"If not given by default the prior pdf from ModelConfig is used."
]
},
{
"cell_type": "markdown",
"id": "9c868c9e",
"metadata": {},
"source": [
"extra options are available as global parameters of the macro. They major ones are:"
]
},
{
"cell_type": "markdown",
"id": "de531c1a",
"metadata": {},
"source": [
"generateBinned generate binned data sets for toys (default is false) - be careful not to activate with\n",
"a too large (>=3) number of observables\n",
"nToyRatio ratio of S+B/B toys (default is 2)\n",
"printLevel"
]
},
{
"cell_type": "markdown",
"id": "62cf7864",
"metadata": {},
"source": [
"disable - can cause some problems\n",
"ToyMCSampler::SetAlwaysUseMultiGen(true);"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "5196a25f",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:49.605679Z",
"iopub.status.busy": "2024-03-19T19:18:49.605278Z",
"iopub.status.idle": "2024-03-19T19:18:49.824108Z",
"shell.execute_reply": "2024-03-19T19:18:49.822697Z"
}
},
"outputs": [],
"source": [
"SimpleLikelihoodRatioTestStat::SetAlwaysReuseNLL(true);\n",
"ProfileLikelihoodTestStat::SetAlwaysReuseNLL(true);\n",
"RatioOfProfiledLikelihoodsTestStat::SetAlwaysReuseNLL(true);"
]
},
{
"cell_type": "markdown",
"id": "edb7245f",
"metadata": {},
"source": [
"RooRandom::randomGenerator()->SetSeed(0);"
]
},
{
"cell_type": "markdown",
"id": "b219a7a0",
"metadata": {},
"source": [
"to change minimizers\n",
"~~~{.bash}\n",
"ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);\n",
"ROOT::Math::MinimizerOptions::SetDefaultMinimizer(\"Minuit2\");\n",
"ROOT::Math::MinimizerOptions::SetDefaultTolerance(1);\n",
"~~~"
]
},
{
"cell_type": "markdown",
"id": "a92432aa",
"metadata": {},
"source": [
"-------------------------------------------------------\n",
"First part is just to access a user-defined file\n",
"or create the standard example file if it doesn't exist"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "b565eabe",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:49.830982Z",
"iopub.status.busy": "2024-03-19T19:18:49.830575Z",
"iopub.status.idle": "2024-03-19T19:18:50.049110Z",
"shell.execute_reply": "2024-03-19T19:18:50.047613Z"
}
},
"outputs": [],
"source": [
"const char *filename = \"\";\n",
"if (!strcmp(infile, \"\")) {\n",
" filename = \"results/example_combined_GaussExample_model.root\";\n",
" bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code\n",
" // if file does not exists generate with histfactory\n",
" if (!fileExist) {\n",
" // Normally this would be run on the command line\n",
" cout << \"will run standard hist2workspace example\" << endl;\n",
" gROOT->ProcessLine(\".! prepareHistFactory .\");\n",
" gROOT->ProcessLine(\".! hist2workspace config/example.xml\");\n",
" cout << \"\\n\\n---------------------\" << endl;\n",
" cout << \"Done creating example input\" << endl;\n",
" cout << \"---------------------\\n\\n\" << endl;\n",
" }\n",
"\n",
"} else\n",
" filename = infile;"
]
},
{
"cell_type": "markdown",
"id": "78d00530",
"metadata": {},
"source": [
"Try to open the file"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "8effe78b",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:50.068928Z",
"iopub.status.busy": "2024-03-19T19:18:50.068512Z",
"iopub.status.idle": "2024-03-19T19:18:50.594727Z",
"shell.execute_reply": "2024-03-19T19:18:50.593651Z"
}
},
"outputs": [],
"source": [
"TFile *file = TFile::Open(filename);"
]
},
{
"cell_type": "markdown",
"id": "8e77683c",
"metadata": {},
"source": [
"if input file was specified but not found, quit"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "ee7c0dcb",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:50.603458Z",
"iopub.status.busy": "2024-03-19T19:18:50.603054Z",
"iopub.status.idle": "2024-03-19T19:18:50.813707Z",
"shell.execute_reply": "2024-03-19T19:18:50.812473Z"
}
},
"outputs": [],
"source": [
"if (!file) {\n",
" cout << \"StandardRooStatsDemoMacro: Input file \" << filename << \" is not found\" << endl;\n",
" return;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "08816b5c",
"metadata": {},
"source": [
"-------------------------------------------------------\n",
"Tutorial starts here\n",
"-------------------------------------------------------"
]
},
{
"cell_type": "markdown",
"id": "2a4fad32",
"metadata": {},
"source": [
"get the workspace out of the file"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "f120a706",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:50.818196Z",
"iopub.status.busy": "2024-03-19T19:18:50.817857Z",
"iopub.status.idle": "2024-03-19T19:18:51.905307Z",
"shell.execute_reply": "2024-03-19T19:18:51.897374Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"RooWorkspace(combined) combined contents\n",
"\n",
"variables\n",
"---------\n",
"(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)\n",
"\n",
"p.d.f.s\n",
"-------\n",
"RooGaussian::alpha_syst1Constraint[ x=alpha_syst1 mean=nom_alpha_syst1 sigma=1 ] = 1\n",
"RooGaussian::alpha_syst2Constraint[ x=alpha_syst2 mean=nom_alpha_syst2 sigma=1 ] = 1\n",
"RooGaussian::alpha_syst3Constraint[ x=alpha_syst3 mean=nom_alpha_syst3 sigma=1 ] = 1\n",
"RooRealSumPdf::channel1_model[ signal_channel1_scaleFactors * signal_channel1_shapes + background1_channel1_scaleFactors * background1_channel1_shapes + background2_channel1_scaleFactors * background2_channel1_shapes ] = 240\n",
"RooPoisson::gamma_stat_channel1_bin_0_constraint[ x=nom_gamma_stat_channel1_bin_0 mean=gamma_stat_channel1_bin_0_poisMean ] = 0.019943\n",
"RooPoisson::gamma_stat_channel1_bin_1_constraint[ x=nom_gamma_stat_channel1_bin_1 mean=gamma_stat_channel1_bin_1_poisMean ] = 0.039861\n",
"RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.1 ] = 1\n",
"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\n",
"RooSimultaneous::simPdf[ indexCat=channelCat channel1=model_channel1 ] = 0.190787\n",
"\n",
"functions\n",
"--------\n",
"RooHistFunc::background1_channel1_Hist_alphanominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 100\n",
"RooStats::HistFactory::FlexibleInterpVar::background1_channel1_epsilon[ paramList=(alpha_syst2) ] = 1\n",
"RooProduct::background1_channel1_scaleFactors[ background1_channel1_epsilon * Lumi ] = 1\n",
"RooProduct::background1_channel1_shapes[ background1_channel1_Hist_alphanominal * mc_stat_channel1 * channel1_model_binWidth ] = 200\n",
"RooHistFunc::background2_channel1_Hist_alphanominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 0\n",
"RooStats::HistFactory::FlexibleInterpVar::background2_channel1_epsilon[ paramList=(alpha_syst3) ] = 1\n",
"RooProduct::background2_channel1_scaleFactors[ background2_channel1_epsilon * Lumi ] = 1\n",
"RooProduct::background2_channel1_shapes[ background2_channel1_Hist_alphanominal * mc_stat_channel1 * channel1_model_binWidth ] = 0\n",
"RooBinWidthFunction::channel1_model_binWidth[ HistFuncForBinWidth=signal_channel1_Hist_alphanominal HistFuncForBinWidth=signal_channel1_Hist_alphanominal ] = 2\n",
"RooProduct::gamma_stat_channel1_bin_0_poisMean[ gamma_stat_channel1_bin_0 * gamma_stat_channel1_bin_0_tau ] = 400\n",
"RooProduct::gamma_stat_channel1_bin_1_poisMean[ gamma_stat_channel1_bin_1 * gamma_stat_channel1_bin_1_tau ] = 100\n",
"ParamHistFunc::mc_stat_channel1[ ] = 1\n",
"RooHistFunc::signal_channel1_Hist_alphanominal[ depList=(obs_x_channel1) depList=(obs_x_channel1) ] = 20\n",
"RooStats::HistFactory::FlexibleInterpVar::signal_channel1_epsilon[ paramList=(alpha_syst1) ] = 1\n",
"RooProduct::signal_channel1_scaleFactors[ signal_channel1_epsilon * SigXsecOverSM * Lumi ] = 1\n",
"RooProduct::signal_channel1_shapes[ signal_channel1_Hist_alphanominal * channel1_model_binWidth ] = 40\n",
"\n",
"datasets\n",
"--------\n",
"RooDataSet::obsData(obs_x_channel1,channelCat)\n",
"RooDataSet::asimovData(obs_x_channel1,channelCat)\n",
"\n",
"embedded datasets (in pdfs and functions)\n",
"-----------------------------------------\n",
"RooDataHist::signal_channel1_Hist_alphanominalDHist(obs_x_channel1)\n",
"RooDataHist::background1_channel1_Hist_alphanominalDHist(obs_x_channel1)\n",
"RooDataHist::background2_channel1_Hist_alphanominalDHist(obs_x_channel1)\n",
"\n",
"parameter snapshots\n",
"-------------------\n",
"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)\n",
"\n",
"named sets\n",
"----------\n",
"ModelConfig_GlobalObservables:(nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)\n",
"ModelConfig_NuisParams:(alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)\n",
"ModelConfig_Observables:(obs_x_channel1,channelCat)\n",
"ModelConfig_POI:(SigXsecOverSM)\n",
"globalObservables:(nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)\n",
"observables:(obs_x_channel1,channelCat)\n",
"\n",
"generic objects\n",
"---------------\n",
"RooStats::ModelConfig::ModelConfig\n",
"\n"
]
}
],
"source": [
"RooWorkspace *w = (RooWorkspace *)file->Get(workspaceName);\n",
"if (!w) {\n",
" cout << \"workspace not found\" << endl;\n",
" return;\n",
"}\n",
"w->Print();"
]
},
{
"cell_type": "markdown",
"id": "5b541870",
"metadata": {},
"source": [
"get the modelConfig out of the file"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "89bf409b",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:51.913508Z",
"iopub.status.busy": "2024-03-19T19:18:51.913096Z",
"iopub.status.idle": "2024-03-19T19:18:52.246649Z",
"shell.execute_reply": "2024-03-19T19:18:52.245496Z"
}
},
"outputs": [],
"source": [
"ModelConfig *sbModel = (ModelConfig *)w->obj(modelSBName);"
]
},
{
"cell_type": "markdown",
"id": "f93c3747",
"metadata": {},
"source": [
"get the modelConfig out of the file"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "a8e49546",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:52.251654Z",
"iopub.status.busy": "2024-03-19T19:18:52.251322Z",
"iopub.status.idle": "2024-03-19T19:18:52.510206Z",
"shell.execute_reply": "2024-03-19T19:18:52.494861Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_85:2:2: warning: 'data' shadows a declaration with the same name in the 'std' namespace; use '::data' to reference this declaration\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n"
]
}
],
"source": [
"RooAbsData *data = w->data(dataName);"
]
},
{
"cell_type": "markdown",
"id": "825c45a6",
"metadata": {},
"source": [
"make sure ingredients are found"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "1bd24c49",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:52.528685Z",
"iopub.status.busy": "2024-03-19T19:18:52.528305Z",
"iopub.status.idle": "2024-03-19T19:18:52.758231Z",
"shell.execute_reply": "2024-03-19T19:18:52.752672Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_86:2:7: error: reference to 'data' is ambiguous\n",
" if (!data || !sbModel) {\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n"
]
}
],
"source": [
"if (!data || !sbModel) {\n",
" w->Print();\n",
" cout << \"data or ModelConfig was not found\" << endl;\n",
" return;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "c07ef32d",
"metadata": {},
"source": [
"make b model"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "67f389f3",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:52.770350Z",
"iopub.status.busy": "2024-03-19T19:18:52.769948Z",
"iopub.status.idle": "2024-03-19T19:18:53.014318Z",
"shell.execute_reply": "2024-03-19T19:18:53.012694Z"
}
},
"outputs": [],
"source": [
"ModelConfig *bModel = (ModelConfig *)w->obj(modelBName);"
]
},
{
"cell_type": "markdown",
"id": "846f34e5",
"metadata": {},
"source": [
"case of no systematics\n",
"remove nuisance parameters from model"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "977d388e",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:53.025464Z",
"iopub.status.busy": "2024-03-19T19:18:53.025066Z",
"iopub.status.idle": "2024-03-19T19:18:53.355763Z",
"shell.execute_reply": "2024-03-19T19:18:53.354613Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Info in : The background model does not exist\n",
"Info in : Copy it from ModelConfig ModelConfig and set POI to zero\n",
"Info in : Model ModelConfig has no snapshot - make one using model poi\n"
]
}
],
"source": [
"if (noSystematics) {\n",
" const RooArgSet *nuisPar = sbModel->GetNuisanceParameters();\n",
" if (nuisPar && nuisPar->getSize() > 0) {\n",
" std::cout << \"StandardHypoTestInvDemo\"\n",
" << \" - Switch off all systematics by setting them constant to their initial values\" << std::endl;\n",
" RooStats::SetAllConstant(*nuisPar);\n",
" }\n",
" if (bModel) {\n",
" const RooArgSet *bnuisPar = bModel->GetNuisanceParameters();\n",
" if (bnuisPar)\n",
" RooStats::SetAllConstant(*bnuisPar);\n",
" }\n",
"}\n",
"\n",
"if (!bModel) {\n",
" Info(\"StandardHypoTestInvDemo\", \"The background model %s does not exist\", modelBName);\n",
" Info(\"StandardHypoTestInvDemo\", \"Copy it from ModelConfig %s and set POI to zero\", modelSBName);\n",
" bModel = (ModelConfig *)sbModel->Clone();\n",
" bModel->SetName(TString(modelSBName) + TString(\"B_only\"));\n",
" RooRealVar *var = dynamic_cast(bModel->GetParametersOfInterest()->first());\n",
" if (!var)\n",
" return;\n",
" double oldval = var->getVal();\n",
" var->setVal(0);\n",
" // bModel->SetSnapshot( RooArgSet(*var, *w->var(\"lumi\")) );\n",
" bModel->SetSnapshot(RooArgSet(*var));\n",
" var->setVal(oldval);\n",
"}\n",
"\n",
"if (!sbModel->GetSnapshot() || poiValue > 0) {\n",
" Info(\"StandardHypoTestDemo\", \"Model %s has no snapshot - make one using model poi\", modelSBName);\n",
" RooRealVar *var = dynamic_cast(sbModel->GetParametersOfInterest()->first());\n",
" if (!var)\n",
" return;\n",
" double oldval = var->getVal();\n",
" if (poiValue > 0)\n",
" var->setVal(poiValue);\n",
" // sbModel->SetSnapshot( RooArgSet(*var, *w->var(\"lumi\") ) );\n",
" sbModel->SetSnapshot(RooArgSet(*var));\n",
" if (poiValue > 0)\n",
" var->setVal(oldval);\n",
" // sbModel->SetSnapshot( *sbModel->GetParametersOfInterest() );\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "651c6525",
"metadata": {},
"source": [
"part 1, hypothesis testing"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "59d4e54e",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:53.365021Z",
"iopub.status.busy": "2024-03-19T19:18:53.364622Z",
"iopub.status.idle": "2024-03-19T19:18:53.581807Z",
"shell.execute_reply": "2024-03-19T19:18:53.580382Z"
}
},
"outputs": [],
"source": [
"SimpleLikelihoodRatioTestStat *slrts = new SimpleLikelihoodRatioTestStat(*bModel->GetPdf(), *sbModel->GetPdf());"
]
},
{
"cell_type": "markdown",
"id": "5119f394",
"metadata": {},
"source": [
"null parameters must include snapshot of poi plus the nuisance values"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "a3b037aa",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:53.586683Z",
"iopub.status.busy": "2024-03-19T19:18:53.586275Z",
"iopub.status.idle": "2024-03-19T19:18:53.804299Z",
"shell.execute_reply": "2024-03-19T19:18:53.803043Z"
}
},
"outputs": [],
"source": [
"RooArgSet nullParams(*bModel->GetSnapshot());\n",
"if (bModel->GetNuisanceParameters())\n",
" nullParams.add(*bModel->GetNuisanceParameters());\n",
"\n",
"slrts->SetNullParameters(nullParams);\n",
"RooArgSet altParams(*sbModel->GetSnapshot());\n",
"if (sbModel->GetNuisanceParameters())\n",
" altParams.add(*sbModel->GetNuisanceParameters());\n",
"slrts->SetAltParameters(altParams);\n",
"\n",
"ProfileLikelihoodTestStat *profll = new ProfileLikelihoodTestStat(*bModel->GetPdf());\n",
"\n",
"RatioOfProfiledLikelihoodsTestStat *ropl =\n",
" new RatioOfProfiledLikelihoodsTestStat(*bModel->GetPdf(), *sbModel->GetPdf(), sbModel->GetSnapshot());\n",
"ropl->SetSubtractMLE(false);\n",
"\n",
"if (testStatType == 3)\n",
" profll->SetOneSidedDiscovery(1);\n",
"profll->SetPrintLevel(printLevel);\n",
"\n",
"if (enableDetOutput) {\n",
" slrts->EnableDetailedOutput();\n",
" profll->EnableDetailedOutput();\n",
" ropl->EnableDetailedOutput();\n",
"}\n",
"\n",
"/* profll.SetReuseNLL(mOptimize);*/\n",
"/* slrts.SetReuseNLL(mOptimize);*/\n",
"/* ropl.SetReuseNLL(mOptimize);*/\n",
"\n",
"AsymptoticCalculator::SetPrintLevel(printLevel);\n",
"\n",
"HypoTestCalculatorGeneric *hypoCalc = 0;"
]
},
{
"cell_type": "markdown",
"id": "a312fc93",
"metadata": {},
"source": [
"note here Null is B and Alt is S+B"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "1a83970e",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:53.809732Z",
"iopub.status.busy": "2024-03-19T19:18:53.809247Z",
"iopub.status.idle": "2024-03-19T19:18:54.019657Z",
"shell.execute_reply": "2024-03-19T19:18:54.018505Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_91:3:42: error: reference to 'data' is ambiguous\n",
" hypoCalc = new FrequentistCalculator(*data, *sbModel, *bModel);\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n",
"input_line_91:5:37: error: reference to 'data' is ambiguous\n",
" hypoCalc = new HybridCalculator(*data, *sbModel, *bModel);\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n",
"input_line_91:7:41: error: reference to 'data' is ambiguous\n",
" hypoCalc = new AsymptoticCalculator(*data, *sbModel, *bModel);\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n"
]
}
],
"source": [
"if (calcType == 0)\n",
" hypoCalc = new FrequentistCalculator(*data, *sbModel, *bModel);\n",
"else if (calcType == 1)\n",
" hypoCalc = new HybridCalculator(*data, *sbModel, *bModel);\n",
"else if (calcType == 2)\n",
" hypoCalc = new AsymptoticCalculator(*data, *sbModel, *bModel);\n",
"\n",
"if (calcType == 0) {\n",
" ((FrequentistCalculator *)hypoCalc)->SetToys(ntoys, ntoys / nToysRatio);\n",
" if (enableDetOutput)\n",
" ((FrequentistCalculator *)hypoCalc)->StoreFitInfo(true);\n",
"}\n",
"if (calcType == 1) {\n",
" ((HybridCalculator *)hypoCalc)->SetToys(ntoys, ntoys / nToysRatio);\n",
" // n. a. yetif (enableDetOutput) ((HybridCalculator*) hypoCalc)->StoreFitInfo(true);\n",
"}\n",
"if (calcType == 2) {\n",
" if (testStatType == 3)\n",
" ((AsymptoticCalculator *)hypoCalc)->SetOneSidedDiscovery(true);\n",
" if (testStatType != 2 && testStatType != 3)\n",
" Warning(\"StandardHypoTestDemo\",\n",
" \"Only the PL test statistic can be used with AsymptoticCalculator - use by default a two-sided PL\");\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "61c04f06",
"metadata": {},
"source": [
"check for nuisance prior pdf in case of nuisance parameters"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "4ed5159a",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:54.034434Z",
"iopub.status.busy": "2024-03-19T19:18:54.034043Z",
"iopub.status.idle": "2024-03-19T19:18:54.241858Z",
"shell.execute_reply": "2024-03-19T19:18:54.240888Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_92:57:24: error: reference to 'data' is ambiguous\n",
" int nEvents = data->numEntries();\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n",
"input_line_92:67:8: error: reference to 'data' is ambiguous\n",
" if (data->isWeighted() && !generateBinned) {\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n",
"input_line_92:67:8: error: unknown type name 'data'\n",
" if (data->isWeighted() && !generateBinned) {\n",
" ^\n",
"input_line_92:67:12: error: cannot use arrow operator on a type\n",
" if (data->isWeighted() && !generateBinned) {\n",
" ^\n",
"input_line_92:70:12: error: reference to 'data' is ambiguous\n",
" data->numEntries(), data->sumEntries());\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n",
"input_line_92:70:32: error: reference to 'data' is ambiguous\n",
" data->numEntries(), data->sumEntries());\n",
" ^\n",
"input_line_85:2:14: note: candidate found by name lookup is 'data'\n",
" RooAbsData *data = w->data(dataName);\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:318:5: note: candidate found by name lookup is 'std::data'\n",
" data(initializer_list<_Tp> __il) noexcept\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:289:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:299:5: note: candidate found by name lookup is 'std::data'\n",
" data(const _Container& __cont) noexcept(noexcept(__cont.data()))\n",
" ^\n",
"/usr/include/c++/9/bits/range_access.h:309:5: note: candidate found by name lookup is 'std::data'\n",
" data(_Tp (&__array)[_Nm]) noexcept\n",
" ^\n"
]
}
],
"source": [
"if (calcType == 1 && (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters())) {\n",
" RooAbsPdf *nuisPdf = 0;\n",
" if (nuisPriorName)\n",
" nuisPdf = w->pdf(nuisPriorName);\n",
" // use prior defined first in bModel (then in SbModel)\n",
" if (!nuisPdf) {\n",
" Info(\"StandardHypoTestDemo\",\n",
" \"No nuisance pdf given for the HybridCalculator - try to deduce pdf from the model\");\n",
" if (bModel->GetPdf() && bModel->GetObservables())\n",
" nuisPdf = RooStats::MakeNuisancePdf(*bModel, \"nuisancePdf_bmodel\");\n",
" else\n",
" nuisPdf = RooStats::MakeNuisancePdf(*sbModel, \"nuisancePdf_sbmodel\");\n",
" }\n",
" if (!nuisPdf) {\n",
" if (bModel->GetPriorPdf()) {\n",
" nuisPdf = bModel->GetPriorPdf();\n",
" Info(\"StandardHypoTestDemo\",\n",
" \"No nuisance pdf given - try to use %s that is defined as a prior pdf in the B model\",\n",
" nuisPdf->GetName());\n",
" } else {\n",
" Error(\"StandardHypoTestDemo\", \"Cannot run Hybrid calculator because no prior on the nuisance parameter is \"\n",
" \"specified or can be derived\");\n",
" return;\n",
" }\n",
" }\n",
" assert(nuisPdf);\n",
" Info(\"StandardHypoTestDemo\", \"Using as nuisance Pdf ... \");\n",
" nuisPdf->Print();\n",
"\n",
" const RooArgSet *nuisParams =\n",
" (bModel->GetNuisanceParameters()) ? bModel->GetNuisanceParameters() : sbModel->GetNuisanceParameters();\n",
" std::unique_ptr np{nuisPdf->getObservables(*nuisParams)};\n",
" if (np->getSize() == 0) {\n",
" Warning(\"StandardHypoTestDemo\",\n",
" \"Prior nuisance does not depend on nuisance parameters. They will be smeared in their full range\");\n",
" }\n",
"\n",
" ((HybridCalculator *)hypoCalc)->ForcePriorNuisanceAlt(*nuisPdf);\n",
" ((HybridCalculator *)hypoCalc)->ForcePriorNuisanceNull(*nuisPdf);\n",
"}\n",
"\n",
"/* hypoCalc->ForcePriorNuisanceAlt(*sbModel->GetPriorPdf());*/\n",
"/* hypoCalc->ForcePriorNuisanceNull(*bModel->GetPriorPdf());*/\n",
"\n",
"ToyMCSampler *sampler = (ToyMCSampler *)hypoCalc->GetTestStatSampler();\n",
"\n",
"if (sampler && (calcType == 0 || calcType == 1)) {\n",
"\n",
" // look if pdf is number counting or extended\n",
" if (sbModel->GetPdf()->canBeExtended()) {\n",
" if (useNC)\n",
" Warning(\"StandardHypoTestDemo\", \"Pdf is extended: but number counting flag is set: ignore it \");\n",
" } else {\n",
" // for not extended pdf\n",
" if (!useNC) {\n",
" int nEvents = data->numEntries();\n",
" Info(\"StandardHypoTestDemo\",\n",
" \"Pdf is not extended: number of events to generate taken from observed data set is %d\", nEvents);\n",
" sampler->SetNEventsPerToy(nEvents);\n",
" } else {\n",
" Info(\"StandardHypoTestDemo\", \"using a number counting pdf\");\n",
" sampler->SetNEventsPerToy(1);\n",
" }\n",
" }\n",
"\n",
" if (data->isWeighted() && !generateBinned) {\n",
" Info(\"StandardHypoTestDemo\", \"Data set is weighted, nentries = %d and sum of weights = %8.1f but toy \"\n",
" \"generation is unbinned - it would be faster to set generateBinned to true\\n\",\n",
" data->numEntries(), data->sumEntries());\n",
" }\n",
" if (generateBinned)\n",
" sampler->SetGenerateBinned(generateBinned);\n",
"\n",
" // use PROOF\n",
" if (useProof) {\n",
" ProofConfig pc(*w, 0, \"\", kFALSE);\n",
" sampler->SetProofConfig(&pc); // enable proof\n",
" }\n",
"\n",
" // set the test statistic\n",
" if (testStatType == 0)\n",
" sampler->SetTestStatistic(slrts);\n",
" if (testStatType == 1)\n",
" sampler->SetTestStatistic(ropl);\n",
" if (testStatType == 2 || testStatType == 3)\n",
" sampler->SetTestStatistic(profll);\n",
"}\n",
"\n",
"HypoTestResult *htr = hypoCalc->GetHypoTest();\n",
"htr->SetPValueIsRightTail(true);\n",
"htr->SetBackgroundAsAlt(false);\n",
"htr->Print(); // how to get meaningful CLs at this point?\n",
"\n",
"delete sampler;\n",
"delete slrts;\n",
"delete ropl;\n",
"delete profll;\n",
"\n",
"if (calcType != 2) {\n",
" HypoTestPlot *plot = new HypoTestPlot(*htr, 100);\n",
" plot->SetLogYaxis(true);\n",
" plot->Draw();\n",
"} else {\n",
" std::cout << \"Asymptotic results \" << std::endl;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "65132b36",
"metadata": {},
"source": [
"look at expected significances\n",
"found median of S+B distribution"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "7679e14b",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:54.246949Z",
"iopub.status.busy": "2024-03-19T19:18:54.246533Z",
"iopub.status.idle": "2024-03-19T19:18:54.451137Z",
"shell.execute_reply": "2024-03-19T19:18:54.450124Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"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)'\n",
" double pval = AsymptoticCalculator::GetExpectedPValues(htr->NullPValue(), htr->AlternatePValue(), -sig, false);\n",
" ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"
]
}
],
"source": [
"if (calcType != 2) {\n",
"\n",
" SamplingDistribution *altDist = htr->GetAltDistribution();\n",
" HypoTestResult htExp(\"Expected Result\");\n",
" htExp.Append(htr);\n",
" // find quantiles in alt (S+B) distribution\n",
" double p[5];\n",
" double q[5];\n",
" for (int i = 0; i < 5; ++i) {\n",
" double sig = -2 + i;\n",
" p[i] = ROOT::Math::normal_cdf(sig, 1);\n",
" }\n",
" std::vector values = altDist->GetSamplingDistribution();\n",
" TMath::Quantiles(values.size(), 5, &values[0], q, p, false);\n",
"\n",
" for (int i = 0; i < 5; ++i) {\n",
" htExp.SetTestStatisticData(q[i]);\n",
" double sig = -2 + i;\n",
" std::cout << \" Expected p -value and significance at \" << sig << \" sigma = \" << htExp.NullPValue()\n",
" << \" significance \" << htExp.Significance() << \" sigma \" << std::endl;\n",
" }\n",
"} else {\n",
" // case of asymptotic calculator\n",
" for (int i = 0; i < 5; ++i) {\n",
" double sig = -2 + i;\n",
" // sigma is inverted here\n",
" double pval = AsymptoticCalculator::GetExpectedPValues(htr->NullPValue(), htr->AlternatePValue(), -sig, false);\n",
" std::cout << \" Expected p -value and significance at \" << sig << \" sigma = \" << pval << \" significance \"\n",
" << ROOT::Math::normal_quantile_c(pval, 1) << \" sigma \" << std::endl;\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "6b963c56",
"metadata": {},
"source": [
"write result in a file in case of toys"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "8af2f08c",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:54.464585Z",
"iopub.status.busy": "2024-03-19T19:18:54.464182Z",
"iopub.status.idle": "2024-03-19T19:18:54.668979Z",
"shell.execute_reply": "2024-03-19T19:18:54.667976Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_95:2:3: error: use of undeclared identifier 'htr'\n",
" (htr != __null && ((*(_Bool*)0x7fb7d8006000)))\n",
" ^\n",
"input_line_95:2:24: error: use of undeclared identifier '_Bool'\n",
" (htr != __null && ((*(_Bool*)0x7fb7d8006000)))\n",
" ^\n",
"input_line_95:2:30: error: expected expression\n",
" (htr != __null && ((*(_Bool*)0x7fb7d8006000)))\n",
" ^\n",
"Error in : Error evaluating expression (htr != __null && ((*(_Bool*)0x7fb7d8006000)))\n",
"Execution of your code was aborted.\n"
]
}
],
"source": [
"bool writeResult = (calcType != 2);\n",
"\n",
"if (enableDetOutput) {\n",
" writeResult = true;\n",
" Info(\"StandardHypoTestDemo\", \"Detailed output will be written in output result file\");\n",
"}\n",
"\n",
"if (htr != NULL && writeResult) {\n",
"\n",
" // write to a file the results\n",
" const char *calcTypeName = (calcType == 0) ? \"Freq\" : (calcType == 1) ? \"Hybr\" : \"Asym\";\n",
" TString resultFileName = TString::Format(\"%s_HypoTest_ts%d_\", calcTypeName, testStatType);\n",
" // strip the / from the filename\n",
"\n",
" TString name = infile;\n",
" name.Replace(0, name.Last('/') + 1, \"\");\n",
" resultFileName += name;\n",
"\n",
" TFile *fileOut = new TFile(resultFileName, \"RECREATE\");\n",
" htr->Write();\n",
" Info(\"StandardHypoTestDemo\", \"HypoTestResult has been written in the file %s\", resultFileName.Data());\n",
"\n",
" fileOut->Close();\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "4960a2f0",
"metadata": {},
"source": [
"Draw all canvases "
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "2b59ef4f",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:18:54.677747Z",
"iopub.status.busy": "2024-03-19T19:18:54.677363Z",
"iopub.status.idle": "2024-03-19T19:18:54.894845Z",
"shell.execute_reply": "2024-03-19T19:18:54.893556Z"
}
},
"outputs": [],
"source": [
"gROOT->GetListOfCanvases()->Draw()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "ROOT C++",
"language": "c++",
"name": "root"
},
"language_info": {
"codemirror_mode": "text/x-c++src",
"file_extension": ".C",
"mimetype": " text/x-c++src",
"name": "c++"
}
},
"nbformat": 4,
"nbformat_minor": 5
}