{
"cells": [
{
"cell_type": "markdown",
"id": "d88736e6",
"metadata": {},
"source": [
"# StandardHypoTestInvDemo\n",
"Standard tutorial macro for performing an inverted hypothesis test for computing an interval\n",
"\n",
"This macro will perform a scan of the p-values for computing the interval or limit\n",
"\n",
"Usage:\n",
"\n",
"```cpp\n",
"root>.L StandardHypoTestInvDemo.C\n",
"root> StandardHypoTestInvDemo(\"fileName\",\"workspace name\",\"S+B modelconfig name\",\"B model name\",\"data set\n",
"name\",calculator type, test statistic type, use CLS,\n",
" number of points, xmin, xmax, number of toys, use number counting)\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 two sided\n",
" = 3 Profile Likelihood one sided (i.e. = 0 if mu < mu_hat)\n",
" = 4 Profile Likelihood signed ( pll = -pll if mu < mu_hat)\n",
" = 5 Max Likelihood Estimate as test statistic\n",
" = 6 Number of observed event as test statistic\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 Sunday, December 21, 2025 at 02:55 PM."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "4c6fba34",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.545446Z",
"iopub.status.busy": "2025-12-21T13:55:19.545275Z",
"iopub.status.idle": "2025-12-21T13:55:19.566448Z",
"shell.execute_reply": "2025-12-21T13:55:19.565948Z"
}
},
"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 \"TROOT.h\"\n",
"#include \"TSystem.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",
"#include \"RooStats/NumEventsTestStat.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",
"using std::cout, std::endl;"
]
},
{
"cell_type": "markdown",
"id": "23205035",
"metadata": {},
"source": [
"structure defining the options"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "d0d6edff",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.568324Z",
"iopub.status.busy": "2025-12-21T13:55:19.568186Z",
"iopub.status.idle": "2025-12-21T13:55:19.881288Z",
"shell.execute_reply": "2025-12-21T13:55:19.880787Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_49:40:2: error: expected expression\n",
" %%cpp -d\n",
" ^\n",
"input_line_49:40:3: error: expected expression\n",
" %%cpp -d\n",
" ^\n",
"input_line_49:40:4: error: use of undeclared identifier 'cpp'\n",
" %%cpp -d\n",
" ^\n",
"input_line_49:40:9: error: use of undeclared identifier 'd'\n",
" %%cpp -d\n",
" ^\n",
"input_line_49:46:20: error: namespaces can only be defined in global or namespace scope\n",
"namespace RooStats {\n",
" ^\n"
]
}
],
"source": [
"struct HypoTestInvOptions {\n",
"\n",
" bool plotHypoTestResult = true; // plot test statistic result at each point\n",
" bool writeResult = true; // write HypoTestInverterResult in a file\n",
" TString resultFileName; // file with results (by default is built automatically using the workspace input file name)\n",
" bool optimize = true; // optimize evaluation of test statistic\n",
" bool useVectorStore = true; // convert data to use new roofit data store\n",
" bool generateBinned = false; // generate binned data sets\n",
" bool noSystematics = false; // force all systematics to be off (i.e. set all nuisance parameters as constat\n",
" // to their nominal values)\n",
" double nToysRatio = 2; // ratio Ntoys S+b/ntoysB\n",
" double maxPOI = -1; // max value used of POI (in case of auto scan)\n",
" bool enableDetailedOutput =\n",
" false; // enable detailed output with all fit information for each toys (output will be written in result file)\n",
" bool rebuild = false; // re-do extra toys for computing expected limits and rebuild test stat\n",
" // distributions (N.B this requires much more CPU (factor is equivalent to nToyToRebuild)\n",
" int nToyToRebuild = 100; // number of toys used to rebuild\n",
" int rebuildParamValues = 0; // = 0 do a profile of all the parameters on the B (alt snapshot) before performing a\n",
" // rebuild operation (default)\n",
" // = 1 use initial workspace parameters with B snapshot values\n",
" // = 2 use all initial workspace parameters with B\n",
" // Otherwise the rebuild will be performed using\n",
" int initialFit = -1; // do a first fit to the model (-1 : default, 0 skip fit, 1 do always fit)\n",
" int randomSeed = -1; // random seed (if = -1: use default value, if = 0 always random )\n",
"\n",
" int nAsimovBins = 0; // number of bins in observables used for Asimov data sets (0 is the default and it is given by\n",
" // workspace, typically is 100)\n",
"\n",
" bool reuseAltToys = false; // reuse same toys for alternate hypothesis (if set one gets more stable bands)\n",
" double confLevel = 0.95; // confidence level value\n",
"\n",
" std::string minimizerType =\n",
" \"\"; // minimizer type (default is what is in ROOT::Math::MinimizerOptions::DefaultMinimizerType()\n",
" std::string massValue = \"\"; // extra string to tag output file of result\n",
" int printLevel = 0; // print level for debugging PL test statistics and calculators\n",
"\n",
" bool useNLLOffset = false; // use NLL offset when fitting (this increase stability of fits)\n",
"};\n",
"%%cpp -d\n",
"\n",
"HypoTestInvOptions optHTInv;\n",
"\n",
"// internal class to run the inverter and more\n",
"\n",
"namespace RooStats {\n",
"\n",
"class HypoTestInvTool {\n",
"\n",
"public:\n",
" HypoTestInvTool();\n",
" ~HypoTestInvTool(){};\n",
"\n",
" HypoTestInverterResult *RunInverter(RooWorkspace *w, const char *modelSBName, const char *modelBName,\n",
" const char *dataName, int type, int testStatType, bool useCLs, int npoints,\n",
" double poimin, double poimax, int ntoys, bool useNumberCounting = false,\n",
" const char *nuisPriorName = nullptr);\n",
"\n",
" void AnalyzeResult(HypoTestInverterResult *r, int calculatorType, int testStatType, bool useCLs, int npoints,\n",
" const char *fileNameBase = nullptr);\n",
"\n",
" void SetParameter(const char *name, const char *value);\n",
" void SetParameter(const char *name, bool value);\n",
" void SetParameter(const char *name, int value);\n",
" void SetParameter(const char *name, double value);\n",
"\n",
"private:\n",
" bool mPlotHypoTestResult;\n",
" bool mWriteResult;\n",
" bool mOptimize;\n",
" bool mUseVectorStore;\n",
" bool mGenerateBinned;\n",
" bool mRebuild;\n",
" bool mReuseAltToys;\n",
" bool mEnableDetOutput;\n",
" int mNToyToRebuild;\n",
" int mRebuildParamValues;\n",
" int mPrintLevel;\n",
" int mInitialFit;\n",
" int mRandomSeed;\n",
" double mNToysRatio;\n",
" double mMaxPoi;\n",
" int mAsimovBins;\n",
" std::string mMassValue;\n",
" std::string\n",
" mMinimizerType; // minimizer type (default is what is in ROOT::Math::MinimizerOptions::DefaultMinimizerType()\n",
" TString mResultFileName;\n",
"};\n",
"\n",
"} "
]
},
{
"cell_type": "markdown",
"id": "dd896f56",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "7b2a62f6",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.885821Z",
"iopub.status.busy": "2025-12-21T13:55:19.885690Z",
"iopub.status.idle": "2025-12-21T13:55:19.902510Z",
"shell.execute_reply": "2025-12-21T13:55:19.902158Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_54:1:1: error: expected expression\n",
"/ end namespace RooStats\n",
"^\n"
]
}
],
"source": [
"%%cpp -d\n",
"/ end namespace RooStats\n",
"\n",
"RooStats::HypoTestInvTool::HypoTestInvTool()\n",
" : mPlotHypoTestResult(true), mWriteResult(false), mOptimize(true), mUseVectorStore(true), mGenerateBinned(false),\n",
" mEnableDetOutput(false), mRebuild(false), mReuseAltToys(false),\n",
" mNToyToRebuild(100), mRebuildParamValues(0), mPrintLevel(0), mInitialFit(-1), mRandomSeed(-1), mNToysRatio(2),\n",
" mMaxPoi(-1), mAsimovBins(0), mMassValue(\"\"), mMinimizerType(\"\"), mResultFileName()\n",
"{\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "f0a4b645",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "399c6e0b",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.908636Z",
"iopub.status.busy": "2025-12-21T13:55:19.908499Z",
"iopub.status.idle": "2025-12-21T13:55:19.910941Z",
"shell.execute_reply": "2025-12-21T13:55:19.910632Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_55:1:16: error: no member named 'HypoTestInvTool' in namespace 'RooStats'\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, bool value)\n",
" ~~~~~~~~~~^\n"
]
}
],
"source": [
"%%cpp -d\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, bool value)\n",
"{\n",
" //\n",
" // set boolean parameters\n",
" //\n",
"\n",
" std::string s_name(name);\n",
"\n",
" if (s_name.find(\"PlotHypoTestResult\") != std::string::npos)\n",
" mPlotHypoTestResult = value;\n",
" if (s_name.find(\"WriteResult\") != std::string::npos)\n",
" mWriteResult = value;\n",
" if (s_name.find(\"Optimize\") != std::string::npos)\n",
" mOptimize = value;\n",
" if (s_name.find(\"UseVectorStore\") != std::string::npos)\n",
" mUseVectorStore = value;\n",
" if (s_name.find(\"GenerateBinned\") != std::string::npos)\n",
" mGenerateBinned = value;\n",
" if (s_name.find(\"EnableDetailedOutput\") != std::string::npos)\n",
" mEnableDetOutput = value;\n",
" if (s_name.find(\"Rebuild\") != std::string::npos)\n",
" mRebuild = value;\n",
" if (s_name.find(\"ReuseAltToys\") != std::string::npos)\n",
" mReuseAltToys = value;\n",
"\n",
" return;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "2ccf9bb0",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "0b90011c",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.921584Z",
"iopub.status.busy": "2025-12-21T13:55:19.921407Z",
"iopub.status.idle": "2025-12-21T13:55:19.923870Z",
"shell.execute_reply": "2025-12-21T13:55:19.923559Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_56:1:16: error: no member named 'HypoTestInvTool' in namespace 'RooStats'\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, int value)\n",
" ~~~~~~~~~~^\n"
]
}
],
"source": [
"%%cpp -d\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, int value)\n",
"{\n",
" //\n",
" // set integer parameters\n",
" //\n",
"\n",
" std::string s_name(name);\n",
"\n",
" if (s_name.find(\"NToyToRebuild\") != std::string::npos)\n",
" mNToyToRebuild = value;\n",
" if (s_name.find(\"RebuildParamValues\") != std::string::npos)\n",
" mRebuildParamValues = value;\n",
" if (s_name.find(\"PrintLevel\") != std::string::npos)\n",
" mPrintLevel = value;\n",
" if (s_name.find(\"InitialFit\") != std::string::npos)\n",
" mInitialFit = value;\n",
" if (s_name.find(\"RandomSeed\") != std::string::npos)\n",
" mRandomSeed = value;\n",
" if (s_name.find(\"AsimovBins\") != std::string::npos)\n",
" mAsimovBins = value;\n",
"\n",
" return;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "cf2170a7",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "a5ca78f5",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.931805Z",
"iopub.status.busy": "2025-12-21T13:55:19.931654Z",
"iopub.status.idle": "2025-12-21T13:55:19.938762Z",
"shell.execute_reply": "2025-12-21T13:55:19.937833Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_57:1:16: error: no member named 'HypoTestInvTool' in namespace 'RooStats'\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, double value)\n",
" ~~~~~~~~~~^\n"
]
}
],
"source": [
"%%cpp -d\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, double value)\n",
"{\n",
" //\n",
" // set double precision parameters\n",
" //\n",
"\n",
" std::string s_name(name);\n",
"\n",
" if (s_name.find(\"NToysRatio\") != std::string::npos)\n",
" mNToysRatio = value;\n",
" if (s_name.find(\"MaxPOI\") != std::string::npos)\n",
" mMaxPoi = value;\n",
"\n",
" return;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "3b8db7ae",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "a92d4b5f",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.949741Z",
"iopub.status.busy": "2025-12-21T13:55:19.949555Z",
"iopub.status.idle": "2025-12-21T13:55:19.958550Z",
"shell.execute_reply": "2025-12-21T13:55:19.958218Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_58:1:16: error: no member named 'HypoTestInvTool' in namespace 'RooStats'\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, const char *value)\n",
" ~~~~~~~~~~^\n"
]
}
],
"source": [
"%%cpp -d\n",
"void RooStats::HypoTestInvTool::SetParameter(const char *name, const char *value)\n",
"{\n",
" //\n",
" // set string parameters\n",
" //\n",
"\n",
" std::string s_name(name);\n",
"\n",
" if (s_name.find(\"MassValue\") != std::string::npos)\n",
" mMassValue.assign(value);\n",
" if (s_name.find(\"MinimizerType\") != std::string::npos)\n",
" mMinimizerType.assign(value);\n",
" if (s_name.find(\"ResultFileName\") != std::string::npos)\n",
" mResultFileName = value;\n",
"\n",
" return;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "8270e02d",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "ac4293ee",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.964065Z",
"iopub.status.busy": "2025-12-21T13:55:19.963916Z",
"iopub.status.idle": "2025-12-21T13:55:19.966686Z",
"shell.execute_reply": "2025-12-21T13:55:19.966348Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_59:1:16: error: no member named 'HypoTestInvTool' in namespace 'RooStats'\n",
"void RooStats::HypoTestInvTool::AnalyzeResult(HypoTestInverterResult *r, int calculatorType, int testStatType,\n",
" ~~~~~~~~~~^\n"
]
}
],
"source": [
"%%cpp -d\n",
"void RooStats::HypoTestInvTool::AnalyzeResult(HypoTestInverterResult *r, int calculatorType, int testStatType,\n",
" bool useCLs, int npoints, const char *fileNameBase)\n",
"{\n",
"\n",
" // analyze result produced by the inverter, optionally save it in a file\n",
"\n",
" double lowerLimit = 0;\n",
" double llError = 0;\n",
"#if defined ROOT_SVN_VERSION && ROOT_SVN_VERSION >= 44126\n",
" if (r->IsTwoSided()) {\n",
" lowerLimit = r->LowerLimit();\n",
" llError = r->LowerLimitEstimatedError();\n",
" }\n",
"#else\n",
" lowerLimit = r->LowerLimit();\n",
" llError = r->LowerLimitEstimatedError();\n",
"#endif\n",
"\n",
" double upperLimit = r->UpperLimit();\n",
" double ulError = r->UpperLimitEstimatedError();\n",
"\n",
" // std::cout << \"DEBUG : [ \" << lowerLimit << \" , \" << upperLimit << \" ] \" << std::endl;\n",
"\n",
" if (lowerLimit < upperLimit * (1. - 1.E-4) && lowerLimit != 0)\n",
" std::cout << \"The computed lower limit is: \" << lowerLimit << \" +/- \" << llError << std::endl;\n",
" std::cout << \"The computed upper limit is: \" << upperLimit << \" +/- \" << ulError << std::endl;\n",
"\n",
" // compute expected limit\n",
" std::cout << \"Expected upper limits, using the B (alternate) model : \" << std::endl;\n",
" std::cout << \" expected limit (median) \" << r->GetExpectedUpperLimit(0) << std::endl;\n",
" std::cout << \" expected limit (-1 sig) \" << r->GetExpectedUpperLimit(-1) << std::endl;\n",
" std::cout << \" expected limit (+1 sig) \" << r->GetExpectedUpperLimit(1) << std::endl;\n",
" std::cout << \" expected limit (-2 sig) \" << r->GetExpectedUpperLimit(-2) << std::endl;\n",
" std::cout << \" expected limit (+2 sig) \" << r->GetExpectedUpperLimit(2) << std::endl;\n",
"\n",
" // detailed output\n",
" if (mEnableDetOutput) {\n",
" mWriteResult = true;\n",
" Info(\"StandardHypoTestInvDemo\", \"detailed output will be written in output result file\");\n",
" }\n",
"\n",
" // write result in a file\n",
" if (r != nullptr && mWriteResult) {\n",
"\n",
" // write to a file the results\n",
" const char *calcType = (calculatorType == 0) ? \"Freq\" : (calculatorType == 1) ? \"Hybr\" : \"Asym\";\n",
" const char *limitType = (useCLs) ? \"CLs\" : \"Cls+b\";\n",
" const char *scanType = (npoints < 0) ? \"auto\" : \"grid\";\n",
" if (mResultFileName.IsNull()) {\n",
" mResultFileName = TString::Format(\"%s_%s_%s_ts%d_\", calcType, limitType, scanType, testStatType);\n",
" // strip the / from the filename\n",
" if (!mMassValue.empty()) {\n",
" mResultFileName += mMassValue.c_str();\n",
" mResultFileName += \"_\";\n",
" }\n",
"\n",
" TString name = fileNameBase;\n",
" name.Replace(0, name.Last('/') + 1, \"\");\n",
" mResultFileName += name;\n",
" }\n",
"\n",
" // get (if existing) rebuilt UL distribution\n",
" TString uldistFile = \"RULDist.root\";\n",
" TObject *ulDist = nullptr;\n",
" bool existULDist = !gSystem->AccessPathName(uldistFile);\n",
" if (existULDist) {\n",
" TFile *fileULDist = TFile::Open(uldistFile);\n",
" if (fileULDist)\n",
" ulDist = fileULDist->Get(\"RULDist\");\n",
" }\n",
"\n",
" TFile *fileOut = new TFile(mResultFileName, \"RECREATE\");\n",
" r->Write();\n",
" if (ulDist)\n",
" ulDist->Write();\n",
" Info(\"StandardHypoTestInvDemo\", \"HypoTestInverterResult has been written in the file %s\", mResultFileName.Data());\n",
"\n",
" fileOut->Close();\n",
" }\n",
"\n",
" // plot the result ( p values vs scan points)\n",
" std::string typeName = \"\";\n",
" if (calculatorType == 0)\n",
" typeName = \"Frequentist\";\n",
" if (calculatorType == 1)\n",
" typeName = \"Hybrid\";\n",
" else if (calculatorType == 2 || calculatorType == 3) {\n",
" typeName = \"Asymptotic\";\n",
" mPlotHypoTestResult = false;\n",
" }\n",
"\n",
" const char *resultName = r->GetName();\n",
" TString plotTitle = TString::Format(\"%s CL Scan for workspace %s\", typeName.c_str(), resultName);\n",
" HypoTestInverterPlot *plot = new HypoTestInverterPlot(\"HTI_Result_Plot\", plotTitle, r);\n",
"\n",
" // plot in a new canvas with style\n",
" TString c1Name = TString::Format(\"%s_Scan\", typeName.c_str());\n",
" TCanvas *c1 = new TCanvas(c1Name);\n",
" c1->SetLogy(false);\n",
"\n",
" plot->Draw(\"CLb 2CL\"); // plot all and Clb\n",
"\n",
" // if (useCLs)\n",
" // plot->Draw(\"CLb 2CL\"); // plot all and Clb\n",
" // else\n",
" // plot->Draw(\"\"); // plot all and Clb\n",
"\n",
" const int nEntries = r->ArraySize();\n",
"\n",
" // plot test statistics distributions for the two hypothesis\n",
" if (mPlotHypoTestResult) {\n",
" TCanvas *c2 = new TCanvas(\"c2\");\n",
" if (nEntries > 1) {\n",
" int ny = TMath::CeilNint(TMath::Sqrt(nEntries));\n",
" int nx = TMath::CeilNint(double(nEntries) / ny);\n",
" c2->Divide(nx, ny);\n",
" }\n",
" for (int i = 0; i < nEntries; i++) {\n",
" if (nEntries > 1)\n",
" c2->cd(i + 1);\n",
" SamplingDistPlot *pl = plot->MakeTestStatPlot(i);\n",
" pl->SetLogYaxis(true);\n",
" pl->Draw();\n",
" }\n",
" }\n",
" gPad = c1;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "ae147cb8",
"metadata": {},
"source": [
" internal routine to run the inverter\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "f23a003f",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.975297Z",
"iopub.status.busy": "2025-12-21T13:55:19.975108Z",
"iopub.status.idle": "2025-12-21T13:55:19.994278Z",
"shell.execute_reply": "2025-12-21T13:55:19.986088Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_60:1:35: error: no member named 'HypoTestInvTool' in namespace 'RooStats'\n",
"HypoTestInverterResult *RooStats::HypoTestInvTool::RunInverter(RooWorkspace *w, const char *modelSBName,\n",
" ~~~~~~~~~~^\n"
]
}
],
"source": [
"%%cpp -d\n",
"HypoTestInverterResult *RooStats::HypoTestInvTool::RunInverter(RooWorkspace *w, const char *modelSBName,\n",
" const char *modelBName, const char *dataName, int type,\n",
" int testStatType, bool useCLs, int npoints,\n",
" double poimin, double poimax, int ntoys,\n",
" bool useNumberCounting, const char *nuisPriorName)\n",
"{\n",
"\n",
" std::cout << \"Running HypoTestInverter on the workspace \" << w->GetName() << std::endl;\n",
"\n",
" w->Print();\n",
"\n",
" RooAbsData *data = w->data(dataName);\n",
" if (!data) {\n",
" Error(\"StandardHypoTestDemo\", \"Not existing data %s\", dataName);\n",
" return nullptr;\n",
" } else\n",
" std::cout << \"Using data set \" << dataName << std::endl;\n",
"\n",
" if (mUseVectorStore) {\n",
" RooAbsData::setDefaultStorageType(RooAbsData::Vector);\n",
" data->convertToVectorStore();\n",
" }\n",
"\n",
" // get models from WS\n",
" // get the modelConfig out of the file\n",
" ModelConfig *bModel = (ModelConfig *)w->obj(modelBName);\n",
" ModelConfig *sbModel = (ModelConfig *)w->obj(modelSBName);\n",
"\n",
" if (!sbModel) {\n",
" Error(\"StandardHypoTestDemo\", \"Not existing ModelConfig %s\", modelSBName);\n",
" return nullptr;\n",
" }\n",
" // check the model\n",
" if (!sbModel->GetPdf()) {\n",
" Error(\"StandardHypoTestDemo\", \"Model %s has no pdf \", modelSBName);\n",
" return nullptr;\n",
" }\n",
" if (!sbModel->GetParametersOfInterest()) {\n",
" Error(\"StandardHypoTestDemo\", \"Model %s has no poi \", modelSBName);\n",
" return nullptr;\n",
" }\n",
" if (!sbModel->GetObservables()) {\n",
" Error(\"StandardHypoTestInvDemo\", \"Model %s has no observables \", modelSBName);\n",
" return nullptr;\n",
" }\n",
" if (!sbModel->GetSnapshot()) {\n",
" Info(\"StandardHypoTestInvDemo\", \"Model %s has no snapshot - make one using model poi\", modelSBName);\n",
" sbModel->SetSnapshot(*sbModel->GetParametersOfInterest());\n",
" }\n",
"\n",
" // case of no systematics\n",
" // remove nuisance parameters from model\n",
" if (optHTInv.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 || bModel == sbModel) {\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(\"_with_poi_0\"));\n",
" RooRealVar *var = dynamic_cast(bModel->GetParametersOfInterest()->first());\n",
" if (!var)\n",
" return nullptr;\n",
" double oldval = var->getVal();\n",
" var->setVal(0);\n",
" bModel->SetSnapshot(RooArgSet(*var));\n",
" var->setVal(oldval);\n",
" } else {\n",
" if (!bModel->GetSnapshot()) {\n",
" Info(\"StandardHypoTestInvDemo\", \"Model %s has no snapshot - make one using model poi and 0 values \",\n",
" modelBName);\n",
" RooRealVar *var = dynamic_cast(bModel->GetParametersOfInterest()->first());\n",
" if (var) {\n",
" double oldval = var->getVal();\n",
" var->setVal(0);\n",
" bModel->SetSnapshot(RooArgSet(*var));\n",
" var->setVal(oldval);\n",
" } else {\n",
" Error(\"StandardHypoTestInvDemo\", \"Model %s has no valid poi\", modelBName);\n",
" return nullptr;\n",
" }\n",
" }\n",
" }\n",
"\n",
" // check model has global observables when there are nuisance pdf\n",
" // for the hybrid case the globals are not needed\n",
" if (type != 1) {\n",
" bool hasNuisParam = (sbModel->GetNuisanceParameters() && sbModel->GetNuisanceParameters()->getSize() > 0);\n",
" bool hasGlobalObs = (sbModel->GetGlobalObservables() && sbModel->GetGlobalObservables()->getSize() > 0);\n",
" if (hasNuisParam && !hasGlobalObs) {\n",
" // try to see if model has nuisance parameters first\n",
" RooAbsPdf *constrPdf = RooStats::MakeNuisancePdf(*sbModel, \"nuisanceConstraintPdf_sbmodel\");\n",
" if (constrPdf) {\n",
" Warning(\"StandardHypoTestInvDemo\", \"Model %s has nuisance parameters but no global observables associated\",\n",
" sbModel->GetName());\n",
" Warning(\"StandardHypoTestInvDemo\",\n",
" \"\\tThe effect of the nuisance parameters will not be treated correctly \");\n",
" }\n",
" }\n",
" }\n",
"\n",
" // save all initial parameters of the model including the global observables\n",
" RooArgSet initialParameters;\n",
" std::unique_ptr allParams{sbModel->GetPdf()->getParameters(*data)};\n",
" allParams->snapshot(initialParameters);\n",
"\n",
" // run first a data fit\n",
"\n",
" const RooArgSet *poiSet = sbModel->GetParametersOfInterest();\n",
" RooRealVar *poi = (RooRealVar *)poiSet->first();\n",
"\n",
" std::cout << \"StandardHypoTestInvDemo : POI initial value: \" << poi->GetName() << \" = \" << poi->getVal()\n",
" << std::endl;\n",
"\n",
" // fit the data first (need to use constraint )\n",
" TStopwatch tw;\n",
"\n",
" bool doFit = mInitialFit;\n",
" if (testStatType == 0 && mInitialFit == -1)\n",
" doFit = false; // case of LEP test statistic\n",
" if (type == 3 && mInitialFit == -1)\n",
" doFit = false; // case of Asymptoticcalculator with nominal Asimov\n",
" double poihat = 0;\n",
"\n",
" if (mMinimizerType.empty())\n",
" mMinimizerType = ROOT::Math::MinimizerOptions::DefaultMinimizerType();\n",
" else\n",
" ROOT::Math::MinimizerOptions::SetDefaultMinimizer(mMinimizerType.c_str());\n",
"\n",
" Info(\"StandardHypoTestInvDemo\", \"Using %s as minimizer for computing the test statistic\",\n",
" ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str());\n",
"\n",
" if (doFit) {\n",
"\n",
" // do the fit : By doing a fit the POI snapshot (for S+B) is set to the fit value\n",
" // and the nuisance parameters nominal values will be set to the fit value.\n",
" // This is relevant when using LEP test statistics\n",
"\n",
" Info(\"StandardHypoTestInvDemo\", \" Doing a first fit to the observed data \");\n",
" RooArgSet constrainParams;\n",
" if (sbModel->GetNuisanceParameters())\n",
" constrainParams.add(*sbModel->GetNuisanceParameters());\n",
" RooStats::RemoveConstantParameters(&constrainParams);\n",
" tw.Start();\n",
" std::unique_ptr fitres{sbModel->GetPdf()->fitTo(\n",
" *data, InitialHesse(false), Hesse(false), Minimizer(mMinimizerType.c_str(), \"Migrad\"), Strategy(0),\n",
" PrintLevel(mPrintLevel), Constrain(constrainParams), Save(true), Offset(RooStats::NLLOffsetMode()))};\n",
" if (fitres->status() != 0) {\n",
" Warning(\"StandardHypoTestInvDemo\",\n",
" \"Fit to the model failed - try with strategy 1 and perform first an Hesse computation\");\n",
" fitres = std::unique_ptr{sbModel->GetPdf()->fitTo(\n",
" *data, InitialHesse(true), Hesse(false), Minimizer(mMinimizerType.c_str(), \"Migrad\"), Strategy(1),\n",
" PrintLevel(mPrintLevel + 1), Constrain(constrainParams), Save(true), Offset(RooStats::NLLOffsetMode()))};\n",
" }\n",
" if (fitres->status() != 0)\n",
" Warning(\"StandardHypoTestInvDemo\", \" Fit still failed - continue anyway.....\");\n",
"\n",
" poihat = poi->getVal();\n",
" std::cout << \"StandardHypoTestInvDemo - Best Fit value : \" << poi->GetName() << \" = \" << poihat << \" +/- \"\n",
" << poi->getError() << std::endl;\n",
" std::cout << \"Time for fitting : \";\n",
" tw.Print();\n",
"\n",
" // save best fit value in the poi snapshot\n",
" sbModel->SetSnapshot(*sbModel->GetParametersOfInterest());\n",
" std::cout << \"StandardHypoTestInvo: snapshot of S+B Model \" << sbModel->GetName()\n",
" << \" is set to the best fit value\" << std::endl;\n",
" }\n",
"\n",
" // print a message in case of LEP test statistics because it affects result by doing or not doing a fit\n",
" if (testStatType == 0) {\n",
" if (!doFit)\n",
" Info(\"StandardHypoTestInvDemo\", \"Using LEP test statistic - an initial fit is not done and the TS will use \"\n",
" \"the nuisances at the model value\");\n",
" else\n",
" Info(\"StandardHypoTestInvDemo\", \"Using LEP test statistic - an initial fit has been done and the TS will use \"\n",
" \"the nuisances at the best fit value\");\n",
" }\n",
"\n",
" // build test statistics and hypotest calculators for running the inverter\n",
"\n",
" SimpleLikelihoodRatioTestStat slrts(*sbModel->GetPdf(), *bModel->GetPdf());\n",
"\n",
" // null parameters must includes snapshot of poi plus the nuisance values\n",
" RooArgSet nullParams(*sbModel->GetSnapshot());\n",
" if (sbModel->GetNuisanceParameters())\n",
" nullParams.add(*sbModel->GetNuisanceParameters());\n",
" if (sbModel->GetSnapshot())\n",
" slrts.SetNullParameters(nullParams);\n",
" RooArgSet altParams(*bModel->GetSnapshot());\n",
" if (bModel->GetNuisanceParameters())\n",
" altParams.add(*bModel->GetNuisanceParameters());\n",
" if (bModel->GetSnapshot())\n",
" slrts.SetAltParameters(altParams);\n",
" if (mEnableDetOutput)\n",
" slrts.EnableDetailedOutput();\n",
"\n",
" // ratio of profile likelihood - need to pass snapshot for the alt\n",
" RatioOfProfiledLikelihoodsTestStat ropl(*sbModel->GetPdf(), *bModel->GetPdf(), bModel->GetSnapshot());\n",
" ropl.SetSubtractMLE(false);\n",
" if (testStatType == 11)\n",
" ropl.SetSubtractMLE(true);\n",
" ropl.SetPrintLevel(mPrintLevel);\n",
" ropl.SetMinimizer(mMinimizerType.c_str());\n",
" if (mEnableDetOutput)\n",
" ropl.EnableDetailedOutput();\n",
"\n",
" ProfileLikelihoodTestStat profll(*sbModel->GetPdf());\n",
" if (testStatType == 3)\n",
" profll.SetOneSided(true);\n",
" if (testStatType == 4)\n",
" profll.SetSigned(true);\n",
" profll.SetMinimizer(mMinimizerType.c_str());\n",
" profll.SetPrintLevel(mPrintLevel);\n",
" if (mEnableDetOutput)\n",
" profll.EnableDetailedOutput();\n",
"\n",
" profll.SetReuseNLL(mOptimize);\n",
" slrts.SetReuseNLL(mOptimize);\n",
" ropl.SetReuseNLL(mOptimize);\n",
"\n",
" if (mOptimize) {\n",
" profll.SetStrategy(0);\n",
" ropl.SetStrategy(0);\n",
" ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);\n",
" }\n",
"\n",
" if (mMaxPoi > 0)\n",
" poi->setMax(mMaxPoi); // increase limit\n",
"\n",
" MaxLikelihoodEstimateTestStat maxll(*sbModel->GetPdf(), *poi);\n",
" NumEventsTestStat nevtts;\n",
"\n",
" AsymptoticCalculator::SetPrintLevel(mPrintLevel);\n",
"\n",
" // create the HypoTest calculator class\n",
" HypoTestCalculatorGeneric *hc = nullptr;\n",
" if (type == 0)\n",
" hc = new FrequentistCalculator(*data, *bModel, *sbModel);\n",
" else if (type == 1)\n",
" hc = new HybridCalculator(*data, *bModel, *sbModel);\n",
" // else if (type == 2 ) hc = new AsymptoticCalculator(*data, *bModel, *sbModel, false, mAsimovBins);\n",
" // else if (type == 3 ) hc = new AsymptoticCalculator(*data, *bModel, *sbModel, true, mAsimovBins); // for using\n",
" // Asimov data generated with nominal values\n",
" else if (type == 2)\n",
" hc = new AsymptoticCalculator(*data, *bModel, *sbModel, false);\n",
" else if (type == 3)\n",
" hc = new AsymptoticCalculator(*data, *bModel, *sbModel,\n",
" true); // for using Asimov data generated with nominal values\n",
" else {\n",
" Error(\"StandardHypoTestInvDemo\", \"Invalid - calculator type = %d supported values are only :\\n\\t\\t\\t 0 \"\n",
" \"(Frequentist) , 1 (Hybrid) , 2 (Asymptotic) \",\n",
" type);\n",
" return nullptr;\n",
" }\n",
"\n",
" // set the test statistic\n",
" TestStatistic *testStat = nullptr;\n",
" if (testStatType == 0)\n",
" testStat = &slrts;\n",
" if (testStatType == 1 || testStatType == 11)\n",
" testStat = &ropl;\n",
" if (testStatType == 2 || testStatType == 3 || testStatType == 4)\n",
" testStat = &profll;\n",
" if (testStatType == 5)\n",
" testStat = &maxll;\n",
" if (testStatType == 6)\n",
" testStat = &nevtts;\n",
"\n",
" if (testStat == nullptr) {\n",
" Error(\"StandardHypoTestInvDemo\", \"Invalid - test statistic type = %d supported values are only :\\n\\t\\t\\t 0 (SLR) \"\n",
" \", 1 (Tevatron) , 2 (PLR), 3 (PLR1), 4(MLE)\",\n",
" testStatType);\n",
" return nullptr;\n",
" }\n",
"\n",
" ToyMCSampler *toymcs = (ToyMCSampler *)hc->GetTestStatSampler();\n",
" if (toymcs && (type == 0 || type == 1)) {\n",
" // look if pdf is number counting or extended\n",
" if (sbModel->GetPdf()->canBeExtended()) {\n",
" if (useNumberCounting)\n",
" Warning(\"StandardHypoTestInvDemo\", \"Pdf is extended: but number counting flag is set: ignore it \");\n",
" } else {\n",
" // for not extended pdf\n",
" if (!useNumberCounting) {\n",
" int nEvents = data->numEntries();\n",
" Info(\"StandardHypoTestInvDemo\",\n",
" \"Pdf is not extended: number of events to generate taken from observed data set is %d\", nEvents);\n",
" toymcs->SetNEventsPerToy(nEvents);\n",
" } else {\n",
" Info(\"StandardHypoTestInvDemo\", \"using a number counting pdf\");\n",
" toymcs->SetNEventsPerToy(1);\n",
" }\n",
" }\n",
"\n",
" toymcs->SetTestStatistic(testStat);\n",
"\n",
" if (data->isWeighted() && !mGenerateBinned) {\n",
" Info(\"StandardHypoTestInvDemo\", \"Data set is weighted, nentries = %d and sum of weights = %8.1f but toy \"\n",
" \"generation is unbinned - it would be faster to set mGenerateBinned to true\\n\",\n",
" data->numEntries(), data->sumEntries());\n",
" }\n",
" toymcs->SetGenerateBinned(mGenerateBinned);\n",
"\n",
" toymcs->SetUseMultiGen(mOptimize);\n",
"\n",
" if (mGenerateBinned && sbModel->GetObservables()->getSize() > 2) {\n",
" Warning(\"StandardHypoTestInvDemo\", \"generate binned is activated but the number of observable is %d. Too much \"\n",
" \"memory could be needed for allocating all the bins\",\n",
" sbModel->GetObservables()->getSize());\n",
" }\n",
"\n",
" // set the random seed if needed\n",
" if (mRandomSeed >= 0)\n",
" RooRandom::randomGenerator()->SetSeed(mRandomSeed);\n",
" }\n",
"\n",
" // specify if need to re-use same toys\n",
" if (mReuseAltToys) {\n",
" hc->UseSameAltToys();\n",
" }\n",
"\n",
" if (type == 1) {\n",
" HybridCalculator *hhc = dynamic_cast(hc);\n",
" assert(hhc);\n",
"\n",
" hhc->SetToys(ntoys, ntoys / mNToysRatio); // can use less ntoys for b hypothesis\n",
"\n",
" // remove global observables from ModelConfig (this is probably not needed anymore in 5.32)\n",
" bModel->SetGlobalObservables(RooArgSet());\n",
" sbModel->SetGlobalObservables(RooArgSet());\n",
"\n",
" // check for nuisance prior pdf in case of nuisance parameters\n",
" if (bModel->GetNuisanceParameters() || sbModel->GetNuisanceParameters()) {\n",
"\n",
" // fix for using multigen (does not work in this case)\n",
" toymcs->SetUseMultiGen(false);\n",
" ToyMCSampler::SetAlwaysUseMultiGen(false);\n",
"\n",
" RooAbsPdf *nuisPdf = nullptr;\n",
" if (nuisPriorName)\n",
" nuisPdf = w->pdf(nuisPriorName);\n",
" // use prior defined first in bModel (then in SbModel)\n",
" if (!nuisPdf) {\n",
" Info(\"StandardHypoTestInvDemo\",\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(\"StandardHypoTestInvDemo\",\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(\"StandardHypoTestInvDemo\", \"Cannot run Hybrid calculator because no prior on the nuisance \"\n",
" \"parameter is specified or can be derived\");\n",
" return nullptr;\n",
" }\n",
" }\n",
" assert(nuisPdf);\n",
" Info(\"StandardHypoTestInvDemo\", \"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(\"StandardHypoTestInvDemo\",\n",
" \"Prior nuisance does not depend on nuisance parameters. They will be smeared in their full range\");\n",
" }\n",
"\n",
" hhc->ForcePriorNuisanceAlt(*nuisPdf);\n",
" hhc->ForcePriorNuisanceNull(*nuisPdf);\n",
" }\n",
" } else if (type == 2 || type == 3) {\n",
" if (testStatType == 3)\n",
" ((AsymptoticCalculator *)hc)->SetOneSided(true);\n",
" if (testStatType != 2 && testStatType != 3)\n",
" Warning(\"StandardHypoTestInvDemo\",\n",
" \"Only the PL test statistic can be used with AsymptoticCalculator - use by default a two-sided PL\");\n",
" } else if (type == 0) {\n",
" ((FrequentistCalculator *)hc)->SetToys(ntoys, ntoys / mNToysRatio);\n",
" // store also the fit information for each poi point used by calculator based on toys\n",
" if (mEnableDetOutput)\n",
" ((FrequentistCalculator *)hc)->StoreFitInfo(true);\n",
" } else if (type == 1) {\n",
" ((HybridCalculator *)hc)->SetToys(ntoys, ntoys / mNToysRatio);\n",
" // store also the fit information for each poi point used by calculator based on toys\n",
" // if (mEnableDetOutput) ((HybridCalculator*) hc)->StoreFitInfo(true);\n",
" }\n",
"\n",
" // Get the result\n",
" RooMsgService::instance().getStream(1).removeTopic(RooFit::NumericIntegration);\n",
"\n",
" HypoTestInverter calc(*hc);\n",
" calc.SetConfidenceLevel(optHTInv.confLevel);\n",
"\n",
" calc.UseCLs(useCLs);\n",
" calc.SetVerbose(true);\n",
"\n",
" if (npoints > 0) {\n",
" if (poimin > poimax) {\n",
" // if no min/max given scan between MLE and +4 sigma\n",
" poimin = int(poihat);\n",
" poimax = int(poihat + 4 * poi->getError());\n",
" }\n",
" std::cout << \"Doing a fixed scan in interval : \" << poimin << \" , \" << poimax << std::endl;\n",
" calc.SetFixedScan(npoints, poimin, poimax);\n",
" } else {\n",
" // poi->setMax(10*int( (poihat+ 10 *poi->getError() )/10 ) );\n",
" std::cout << \"Doing an automatic scan in interval : \" << poi->getMin() << \" , \" << poi->getMax() << std::endl;\n",
" }\n",
"\n",
" tw.Start();\n",
" HypoTestInverterResult *r = calc.GetInterval();\n",
" std::cout << \"Time to perform limit scan \\n\";\n",
" tw.Print();\n",
"\n",
" if (mRebuild) {\n",
"\n",
" std::cout << \"\\n***************************************************************\\n\";\n",
" std::cout << \"Rebuild the upper limit distribution by re-generating new set of pseudo-experiment and re-compute \"\n",
" \"for each of them a new upper limit\\n\\n\";\n",
"\n",
" allParams = std::unique_ptr{sbModel->GetPdf()->getParameters(*data)};\n",
"\n",
" // define on which value of nuisance parameters to do the rebuild\n",
" // default is best fit value for bmodel snapshot\n",
"\n",
" if (mRebuildParamValues != 0) {\n",
" // set all parameters to their initial workspace values\n",
" allParams->assign(initialParameters);\n",
" }\n",
" if (mRebuildParamValues == 0 || mRebuildParamValues == 1) {\n",
" RooArgSet constrainParams;\n",
" if (sbModel->GetNuisanceParameters())\n",
" constrainParams.add(*sbModel->GetNuisanceParameters());\n",
" RooStats::RemoveConstantParameters(&constrainParams);\n",
"\n",
" const RooArgSet *poiModel = sbModel->GetParametersOfInterest();\n",
" bModel->LoadSnapshot();\n",
"\n",
" // do a profile using the B model snapshot\n",
" if (mRebuildParamValues == 0) {\n",
"\n",
" RooStats::SetAllConstant(*poiModel, true);\n",
"\n",
" sbModel->GetPdf()->fitTo(*data, InitialHesse(false), Hesse(false),\n",
" Minimizer(mMinimizerType.c_str(), \"Migrad\"), Strategy(0), PrintLevel(mPrintLevel),\n",
" Constrain(constrainParams), Offset(RooStats::NLLOffsetMode()));\n",
"\n",
" std::cout << \"rebuild using fitted parameter value for B-model snapshot\" << std::endl;\n",
" constrainParams.Print(\"v\");\n",
"\n",
" RooStats::SetAllConstant(*poiModel, false);\n",
" }\n",
" }\n",
" std::cout << \"StandardHypoTestInvDemo: Initial parameters used for rebuilding: \";\n",
" RooStats::PrintListContent(*allParams, std::cout);\n",
"\n",
" tw.Start();\n",
" SamplingDistribution *limDist = calc.GetUpperLimitDistribution(true, mNToyToRebuild);\n",
" std::cout << \"Time to rebuild distributions \" << std::endl;\n",
" tw.Print();\n",
"\n",
" if (limDist) {\n",
" std::cout << \"Expected limits after rebuild distribution \" << std::endl;\n",
" std::cout << \"expected upper limit (median of limit distribution) \" << limDist->InverseCDF(0.5) << std::endl;\n",
" std::cout << \"expected -1 sig limit (0.16% quantile of limit dist) \"\n",
" << limDist->InverseCDF(ROOT::Math::normal_cdf(-1)) << std::endl;\n",
" std::cout << \"expected +1 sig limit (0.84% quantile of limit dist) \"\n",
" << limDist->InverseCDF(ROOT::Math::normal_cdf(1)) << std::endl;\n",
" std::cout << \"expected -2 sig limit (.025% quantile of limit dist) \"\n",
" << limDist->InverseCDF(ROOT::Math::normal_cdf(-2)) << std::endl;\n",
" std::cout << \"expected +2 sig limit (.975% quantile of limit dist) \"\n",
" << limDist->InverseCDF(ROOT::Math::normal_cdf(2)) << std::endl;\n",
"\n",
" // Plot the upper limit distribution\n",
" SamplingDistPlot limPlot((mNToyToRebuild < 200) ? 50 : 100);\n",
" limPlot.AddSamplingDistribution(limDist);\n",
" limPlot.GetTH1F()->SetStats(true); // display statistics\n",
" limPlot.SetLineColor(kBlue);\n",
" new TCanvas(\"limPlot\", \"Upper Limit Distribution\");\n",
" limPlot.Draw();\n",
"\n",
" /// save result in a file\n",
" limDist->SetName(\"RULDist\");\n",
" TFile *fileOut = new TFile(\"RULDist.root\", \"RECREATE\");\n",
" limDist->Write();\n",
" fileOut->Close();\n",
"\n",
" // update r to a new updated result object containing the rebuilt expected p-values distributions\n",
" // (it will not recompute the expected limit)\n",
" if (r)\n",
" delete r; // need to delete previous object since GetInterval will return a cloned copy\n",
" r = calc.GetInterval();\n",
"\n",
" } else\n",
" std::cout << \"ERROR : failed to re-build distributions \" << std::endl;\n",
" }\n",
"\n",
" return r;\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "681c0101",
"metadata": {},
"source": [
" Definition of a helper function: "
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "a19132a3",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:19.996136Z",
"iopub.status.busy": "2025-12-21T13:55:19.995930Z",
"iopub.status.idle": "2025-12-21T13:55:19.999714Z",
"shell.execute_reply": "2025-12-21T13:55:19.999309Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_61:5:4: error: use of undeclared identifier 'StandardHypoTestInvDemo'\n",
" StandardHypoTestInvDemo(fileName, resultName, \"\", \"\", \"\", 0, 0, useCLs);\n",
" ^\n"
]
}
],
"source": [
"%%cpp -d\n",
"void ReadResult(const char *fileName, const char *resultName = \"\", bool useCLs = true)\n",
"{\n",
" // read a previous stored result from a file given the result name\n",
"\n",
" StandardHypoTestInvDemo(fileName, resultName, \"\", \"\", \"\", 0, 0, useCLs);\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "124e0827",
"metadata": {},
"source": [
" Arguments are defined. "
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "9f30fe48",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:20.001424Z",
"iopub.status.busy": "2025-12-21T13:55:20.001291Z",
"iopub.status.idle": "2025-12-21T13:55:20.212169Z",
"shell.execute_reply": "2025-12-21T13:55:20.211606Z"
}
},
"outputs": [],
"source": [
"const char *infile = nullptr;\n",
"const char *wsName = \"combined\";\n",
"const char *modelSBName = \"ModelConfig\";\n",
"const char *modelBName = \"\";\n",
"const char *dataName = \"obsData\";\n",
"int calculatorType = 0;\n",
"int testStatType = 0;\n",
"bool useCLs = true;\n",
"int npoints = 6;\n",
"double poimin = 0;\n",
"double poimax = 5;\n",
"int ntoys = 1000;\n",
"bool useNumberCounting = false;\n",
"const char *nuisPriorName = nullptr;"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "0cadb747",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:20.214187Z",
"iopub.status.busy": "2025-12-21T13:55:20.214058Z",
"iopub.status.idle": "2025-12-21T13:55:20.426245Z",
"shell.execute_reply": "2025-12-21T13:55:20.425674Z"
}
},
"outputs": [],
"source": [
"/*\n",
"\n",
" Other Parameter to pass in tutorial\n",
" apart from standard for filename, ws, modelconfig and data\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 < mu_hat)\n",
" = 4 Profiel Likelihood signed ( pll = -pll if mu < mu_hat)\n",
" = 5 Max Likelihood Estimate as test statistic\n",
" = 6 Number of observed event as test statistic\n",
"\n",
" useCLs scan for CLs (otherwise for CLs+b)\n",
"\n",
" npoints: number of points to scan , for autoscan set npoints = -1\n",
"\n",
" poimin,poimax: min/max value to scan in case of fixed scans\n",
" (if min > max, try to find automatically)\n",
"\n",
" ntoys: number of toys to use\n",
"\n",
" useNumberCounting: set to true when using number counting events\n",
"\n",
" 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.\n",
"\n",
" extra options are available as global parameters of the macro. They major ones are:\n",
"\n",
" plotHypoTestResult plot result of tests at each point (TS distributions) (default is true)\n",
" writeResult write result of scan (default is true)\n",
" rebuild rebuild scan for expected limits (require extra toys) (default is false)\n",
" 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",
"\n",
"\n",
"*/\n",
"\n",
"TString filename(infile);\n",
"if (filename.IsNull()) {\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": "1d3cf7cb",
"metadata": {},
"source": [
"Try to open the file"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "dcec550f",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:20.433811Z",
"iopub.status.busy": "2025-12-21T13:55:20.433657Z",
"iopub.status.idle": "2025-12-21T13:55:20.755788Z",
"shell.execute_reply": "2025-12-21T13:55:20.755365Z"
}
},
"outputs": [],
"source": [
"TFile *file = TFile::Open(filename);"
]
},
{
"cell_type": "markdown",
"id": "e4285f82",
"metadata": {},
"source": [
"if input file was specified but not found, quit"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "2efb7a52",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:20.757775Z",
"iopub.status.busy": "2025-12-21T13:55:20.757642Z",
"iopub.status.idle": "2025-12-21T13:55:20.969592Z",
"shell.execute_reply": "2025-12-21T13:55:20.969115Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_81:7:1: error: unknown type name 'HypoTestInvTool'\n",
"HypoTestInvTool calc;\n",
"^\n"
]
}
],
"source": [
"if (!file) {\n",
" cout << \"StandardRooStatsDemoMacro: Input file \" << filename << \" is not found\" << endl;\n",
" return;\n",
"}\n",
"\n",
"HypoTestInvTool calc;"
]
},
{
"cell_type": "markdown",
"id": "8ef20381",
"metadata": {},
"source": [
"set parameters"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "056b4dd0",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:20.971987Z",
"iopub.status.busy": "2025-12-21T13:55:20.971788Z",
"iopub.status.idle": "2025-12-21T13:55:21.200567Z",
"shell.execute_reply": "2025-12-21T13:55:21.194277Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_83:2:3: error: use of undeclared identifier 'calc'\n",
" (calc.SetParameter(\"PlotHypoTestResult\", optHTInv.plotHypoTestResult))\n",
" ^\n",
"input_line_83:2:43: error: use of undeclared identifier 'optHTInv'\n",
" (calc.SetParameter(\"PlotHypoTestResult\", optHTInv.plotHypoTestResult))\n",
" ^\n",
"Error in : Error evaluating expression (calc.SetParameter(\"PlotHypoTestResult\", optHTInv.plotHypoTestResult))\n",
"Execution of your code was aborted.\n"
]
}
],
"source": [
"calc.SetParameter(\"PlotHypoTestResult\", optHTInv.plotHypoTestResult);\n",
"calc.SetParameter(\"WriteResult\", optHTInv.writeResult);\n",
"calc.SetParameter(\"Optimize\", optHTInv.optimize);\n",
"calc.SetParameter(\"UseVectorStore\", optHTInv.useVectorStore);\n",
"calc.SetParameter(\"GenerateBinned\", optHTInv.generateBinned);\n",
"calc.SetParameter(\"NToysRatio\", optHTInv.nToysRatio);\n",
"calc.SetParameter(\"MaxPOI\", optHTInv.maxPOI);\n",
"calc.SetParameter(\"EnableDetailedOutput\", optHTInv.enableDetailedOutput);\n",
"calc.SetParameter(\"Rebuild\", optHTInv.rebuild);\n",
"calc.SetParameter(\"ReuseAltToys\", optHTInv.reuseAltToys);\n",
"calc.SetParameter(\"NToyToRebuild\", optHTInv.nToyToRebuild);\n",
"calc.SetParameter(\"RebuildParamValues\", optHTInv.rebuildParamValues);\n",
"calc.SetParameter(\"MassValue\", optHTInv.massValue.c_str());\n",
"calc.SetParameter(\"MinimizerType\", optHTInv.minimizerType.c_str());\n",
"calc.SetParameter(\"PrintLevel\", optHTInv.printLevel);\n",
"calc.SetParameter(\"InitialFit\", optHTInv.initialFit);\n",
"calc.SetParameter(\"ResultFileName\", optHTInv.resultFileName);\n",
"calc.SetParameter(\"RandomSeed\", optHTInv.randomSeed);\n",
"calc.SetParameter(\"AsimovBins\", optHTInv.nAsimovBins);"
]
},
{
"cell_type": "markdown",
"id": "d547ea30",
"metadata": {},
"source": [
"enable offset for all roostats"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "7b3abb08",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:21.211026Z",
"iopub.status.busy": "2025-12-21T13:55:21.210800Z",
"iopub.status.idle": "2025-12-21T13:55:21.453687Z",
"shell.execute_reply": "2025-12-21T13:55:21.453238Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"input_line_85:2:3: error: use of undeclared identifier 'optHTInv'\n",
" (optHTInv)\n",
" ^\n",
"Error in : Error evaluating expression (optHTInv)\n",
"Execution of your code was aborted.\n"
]
}
],
"source": [
"if (optHTInv.useNLLOffset)\n",
" RooStats::SetNLLOffsetMode(\"initial\");\n",
"\n",
"RooWorkspace *w = dynamic_cast(file->Get(wsName));\n",
"HypoTestInverterResult *r = nullptr;\n",
"std::cout << w << \"\\t\" << filename << std::endl;\n",
"if (w != nullptr) {\n",
" r = calc.RunInverter(w, modelSBName, modelBName, dataName, calculatorType, testStatType, useCLs, npoints, poimin,\n",
" poimax, ntoys, useNumberCounting, nuisPriorName);\n",
" if (!r) {\n",
" std::cerr << \"Error running the HypoTestInverter - Exit \" << std::endl;\n",
" return;\n",
" }\n",
"} else {\n",
" // case workspace is not present look for the inverter result\n",
" std::cout << \"Reading an HypoTestInverterResult with name \" << wsName << \" from file \" << filename << std::endl;\n",
" r = dynamic_cast(file->Get(wsName)); //\n",
" if (!r) {\n",
" std::cerr << \"File \" << filename << \" does not contain a workspace or an HypoTestInverterResult - Exit \"\n",
" << std::endl;\n",
" file->ls();\n",
" return;\n",
" }\n",
"}\n",
"\n",
"calc.AnalyzeResult(r, calculatorType, testStatType, useCLs, npoints, infile);\n",
"\n",
"return;"
]
},
{
"cell_type": "markdown",
"id": "9bf1a274",
"metadata": {},
"source": [
"Draw all canvases "
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "32ae256f",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2025-12-21T13:55:21.461607Z",
"iopub.status.busy": "2025-12-21T13:55:21.461474Z",
"iopub.status.idle": "2025-12-21T13:55:21.671445Z",
"shell.execute_reply": "2025-12-21T13:55:21.670995Z"
}
},
"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
}