Introduction

This resource contains the full statistical model related to the paper Search for higgsinos decaying to two Higgs bosons and missing transverse momentum in proton-proton collisions at √s= 13 TeV, by the CMS collaboration. The results of the analysis are interpreted in three Simplified Model Space (SMS) signal models, each a function of one or two mass parameters. The instructions below include the commands needed to reproduce the exclusion limits for the models named TChiHH-G, TChiHH, and T5HH after the respective SMS models, and corresponding to Figures 11, 13, and 14 in the paper. For additional technical documentation please consult the Combine documentation.

Datacards

While a full description of the analysis can be found in the paper, a short description of the datacards contained in this repository follows.

Datacards are held in the datacards directory. Three sets of datacards are provided, one for each of the three SNS models studied. Each set is contained in a subdirectory within the datacards directory, with a name matching the relevant model.

Within each datacards/[model] directory, a datacard is provided for each signal mass hypothesis. The naming conventions for datacards in the T5HH and TChiHH-G directories is:

1D[model][massNLSP]_LSP1.txt

where [model] is one of T5HH or TChiHH-G and [massNLSP] is the mass in GeV of the next-to-lightest supersymmetric particle (NLSP).

The naming convention for datacards in the TChiHH directory is:

2DTChiHH[massNLSP]_LSP[massLSP].txt

where [massNLSP] and [massLSP] are the masses of the NLSP and of the LSP, respectively, in GeV.

Each datacard is the combination of one counting experiment per analysis category.

The naming convention for categories is:

• for signal regions, referring to Fig. 10 in the paper: SR[bin_number]_[subcat]_[metcut], where:
  • [bin_number] runs from 1 to 22;
  • [subcat] can be one of dr(loose|tight)_(3|4)b or (1|2)H;
  • [metcut] is coded as met[1-4] and refers to the \(p_T^{miss}\) cut.
• for regions involved in the ABCD backround estimation, referring to Fig. 4 and to the text, the naming convention is: [region]_[subcat]_[metcut], where:
  • [region] is (CSR|CSB|SB)
  • [subcat] can be one of dr(loose|tight)_(2|3|4)b or (0|1|2)H
  • [metcut] is coded as met[1-4] and refers to the \(p_T^{miss}\) cut.
• the regions CR0Hb_[metcut] are control regions for the SR bins 17 to 22.

The signal process is always identified as sig. The sum of all background processes is always identified as bkg.

Nuisance parameters describing systematic uncertainties are described doc/sysNote.pdf and in input/systematics_SUS-20-004.html.

The parameter of interest (POI) is a signal strength modifier r of the expected cross section for the relevant model.

Software instructions

General installation instructions for Combine can be found in the Combine documentation.

A container image is provided to ensure reproducible results. The results in this README are obtained using v9.2.1:

docker run --name combine -it gitlab-registry.cern.ch/cms-cloud/combine-standalone:v9.2.1

A slim version of the container image is also available at gitlab-registry.cern.ch/cms-cloud/combine-standalone:v9.2.1-slim. Versions of packages in the slim container image do not match exactly with the ones in the default container, so small differences in the output of commands with respect to the ones shown below are to be expected.

You can copy files (such as the datacards and other inputs for combine) using docker cp as documented here.

For the commands below, you may require running ulimit -s unlimited; ulimit -u unlimited to avoid memory issues.

Scripts

The scripts subdirectory contains modules to facilitate the extraction of the asymptotic limits, plots of the profile likelihood, and the likelihood for any given value of the signal strength. Instructions for runnng these follow.

runAsymptoticLimits.py

$ python3 scripts/runAsymptoticLimits.py -h
usage: runAsymptoticLimits [-h] [--datacard [DATACARD ...]] [--processes PROCESSES] [--dryrun] sourceDir

Run combine AsymptoticLimits method on datacards, saving the limit
TTrees in root files.

positional arguments:
  sourceDir             Directory containing the datacards

optional arguments:
  -h, --help            show this help message and exit
  --datacard [DATACARD ...]
                        Datacard(s) to process; if omitted, process all .txt files in sourceDir
  --processes PROCESSES
                        Number of processes; by default the available cores.
  --dryrun              Just show the command to be run.

For example, running python3 scripts/runAsymptoticLimits.py datacards/TChiHH will result in several ROOT file, one per datacard, containing the expected and observed limit.

One can browse the expected/observed limits on r by doing:

$ root -l higgsCombine1DTChiHH1100_LSP1.AsymptoticLimits.mH120.root
root [0] 
Attaching file higgsCombine1DTChiHH1100_LSP1.AsymptoticLimits.mH120.root as _file0...
(TFile *) 0x4c9d980
root [1] limit->Scan("limit:quantileExpected")
************************************
*    Row   *     limit * quantileE *
************************************
*        0 * 0.7407836 * 0.0250000 *
*        1 * 1.1093192 * 0.1599999 *
*        2 * 1.7890625 *       0.5 *
*        3 * 2.9585924 * 0.8399999 *
*        4 * 4.7229156 * 0.9750000 *
*        5 * 1.3800972 *        -1 *
************************************

The expected limit is the one with quantileExpected value of 0.5. The different quantiles correspond to -2 \(\sigma\) (2.5% quantile), -1 \(\sigma\) (16% quantile), +1 \(\sigma\) (84% quantile), +2 \(\sigma\) (97.5% quantile). The observed limis is identified by quantileExpected=-1.

likeliTools.py

$ python3 scripts/likeliTools.py -h
usage: likeliTools [-h] [--nPoints NPOINTS] [--rMin RMIN] [--rMax RMAX] [--title TITLE] datacard

Examine the likelihood function for a Combine datacard.  The
class 'likelihood' can be instantiated from the user's program,
or this module can be invoked directly from the comand line to
produce a profile likelihood plot.

The likelihood class constructor reads the datacard workspace
.root file, or creates that file if given a .txt datacard file.
It then performs the maximum likelihood fit and adds the result
to the workspace. 

The deltaLnL method can be called to obtain the value of the
negative log(likelihood) difference with respect to the
maximimum likelihood solution for a given value of the signal
strength.

The profile method produces a profile likelihood plot.

The plotProfile method could be called if the MultiDimFit scan
root file is already available.

Finally, when called stand-alone, this program makes the profile
likelihood plot pdf file.

positional arguments:
  datacard           Path to datacard .txt or .root file

optional arguments:
  -h, --help         show this help message and exit
  --nPoints NPOINTS  Number of points for profile plot
  --rMin RMIN        Minimum r for profile plot
  --rMax RMAX        Maximum r for profile plot
  --title TITLE      Title for profile plot

Example usage:

python3 scripts/likeliTools.py datacards/TChiHH-G/1DTChiHH400_LSP1.txt --rMin 0 --rMax 0.5`

The script will produce a pdf file with the requested scan in the parameter of interest r. The values of \(-\Delta{}\mathrm{log}(\mathcal{L})\) are also saved in the limit TTree in the resulting ROOT file and can be browsed as:

root -l higgsCombine.1DTChiHH400_LSP1_scan.MultiDimFit.mH120.root
root [0] 
Attaching file higgsCombine.1DTChiHH400_LSP1_scan.MultiDimFit.mH120.root as _file0...
(TFile *) 0x4919a70
root [1] limit->Scan("r:deltaNLL")
************************************
*    Row   *         r *  deltaNLL *
************************************
*        0 * 0.1289093 *         0 *
*        1 * 0.0041666 * 1.6181876 *
*        2 * 0.0125000 * 1.4027500 *
*        3 * 0.0208333 *   1.20927 *
...

altSignal.py

Alternative signal models can be implemented by substituting their yields in the 22 signal and 39 control regions for those appearing in the "rate" line of the datacard; the value to be entered for the background process in all bins is 1. The altSignal.py script facilitates this process.
A datacard with parameters corresponding to similar kinematic conditions to those of the alternate signal should be chosen, since there is some dependence of the nuisance parameters on the signal model (see the Note on systematics in file doc/sysNote.pdf) for specifics.

$ python3 scripts/altSignal.py -h
usage: altSignal [-h] datacard replacement newcard

Replace the rate line of a datacard .txt file with a provided
alternate, or with the rate line of an alternate datacard.      

positional arguments:
  datacard     Path to input datacard .txt file
  replacement  Replacement rate line or datacard .txt file
  newcard      Output datacard .txt file

optional arguments:
  -h, --help   show this help message and exit

Instructions to reproduce Figs 11, 13, and 14

Fig. 11

To reproduce Fig 11, run:

python3 scripts/runAsymptoticLimits.py datacards/TChiHH-G

This will result in a set of ROOT files, each corresponding to a signal mass point, holding the expected and observed limits on r in the limit TTree object. For example, for NLSP mass point 1000 GeV, the produced ROOT tile is higgsCombine1DTChiHH1000_LSP1.AsymptoticLimits.mH120.root, and its content looks like:

$ root -l higgsCombine1DTChiHH1000_LSP1.AsymptoticLimits.mH120.root
root [0] 
Attaching file higgsCombine1DTChiHH1000_LSP1.AsymptoticLimits.mH120.root as _file0...
(TFile *) 0x501ae30
root [1] limit->Scan("limit:quantileExpected")
************************************
*    Row   *     limit * quantileE *
************************************
*        0 * 0.5126953 * 0.0250000 *
*        1 * 0.7495880 * 0.1599999 *
*        2 *  1.171875 *       0.5 *
*        3 * 1.8912311 * 0.8399999 *
*        4 * 2.9233350 * 0.9750000 *
*        5 * 0.9210153 *        -1 *
************************************

To translate the limits on r to those on the cross section, shown in Fig. 11, we need to multiply the limit on r by the expected cross section of the TChiHH-G model tabulated in HEPDATA. For example, for the 1000 GeV mass point, the expected cross section is 0.96885 fb, thus the expected limit is 1.171875*0.96885=1.1354 fb.

Fig. 13

To reproduce Fig 13, run:

python3 scripts/runAsymptoticLimits.py datacards/TChiHH

This will result in a set of ROOT files, each corresponding to a signal NLSP and LSP mass pair, holding the expected and observed limits on r in the limit TTree object. For example, for NLSP value of 800 GeV and LSP value of 675 GeV, the produced ROOT tile is higgsCombine2DTChiHH800_LSP675.AsymptoticLimits.mH120.root, and its content looks like:

$ root -l higgsCombine2DTChiHH800_LSP675.AsymptoticLimits.mH120.root
root [0] 
Attaching file higgsCombine2DTChiHH800_LSP675.AsymptoticLimits.mH120.root as _file0...
(TFile *) 0x458e4f0
root [1] limit->Scan("limit:quantileExpected")
************************************
*    Row   *     limit * quantileE *
************************************
*        0 *   693.375 * 0.0250000 *
*        1 * 944.02636 * 0.1599999 *
*        2 *      1376 *       0.5 *
*        3 * 2012.2362 * 0.8399999 *
*        4 * 2801.2753 * 0.9750000 *
*        5 * 1655.5354 *        -1 *
************************************

To translate the limits on r to those on the cross section, shown in Fig. 13, we need to multiply the limit on r by the expected cross section of the TChiHH model tabulated in HEPDATA, column \(\sigma_{thy}\). For example, for the NLSP value of 800 GeV, LSP value of 675 GeV, the expected cross section is 0.58621 fb, thus the expected limit is 1376*0.58621=806.63 fb.

Fig. 14

To reproduce Fig 14, run:

python3 scripts/runAsymptoticLimits.py datacards/T5HH

This will result in a set of ROOT files, each corresponding to a signal mass point, holding the expected and observed limits on r in the limit TTree object. For example, for NLSP mass point 1000 GeV, the produced ROOT tile is higgsCombine1DT5HH1000_LSP1.AsymptoticLimits.mH120.root, and its content looks like:

$ root -l higgsCombine1DT5HH1000_LSP1.AsymptoticLimits.mH120.root
root [0] 
Attaching file higgsCombine1DT5HH1000_LSP1.AsymptoticLimits.mH120.root as _file0...
(TFile *) 0x2382760
root [1] limit->Scan("limit:quantileExpected")
************************************
*    Row   *     limit * quantileE *
************************************
*        0 * 0.0021362 * 0.0250000 *
*        1 * 0.0024566 * 0.1599999 *
*        2 * 0.0034179 *       0.5 *
*        3 * 0.0055978 * 0.8399999 *
*        4 * 0.0088948 * 0.9750000 *
*        5 * 0.0026889 *        -1 *
************************************

To translate the limits on r to those on the cross section, shown in Fig. 14, we need to multiply the limit on r by the expected cross section of the T5HH model tabulated in HEPDATA, column \(\sigma_{Theory}\). For example, for the 1000 GeV mass point, the expected cross section is 325.39 fb, thus the expected limit is 0.0034179*325.39=1.1122 fb.