Reduce, Reuse, Reinterpret: an end-to-end pipeline for recycling particle physics results

Abstract

Searches for new physics at the Large Hadron Collider have constrained many models of physics beyond the Standard Model. Many searches also provide resources that allow them to be reinterpreted in the context of other models. We describe a reinterpretation pipeline that examines previously untested models of new physics using supplementary information from ATLAS Supersymmetry (SUSY) searches in a way that provides accurate constraints even for models that differ meaningfully from the benchmark models of the original analysis. The public analysis information, such as public analysis routines and serialized probability models, is combined with common event generation and simulation toolkits MadGraph, Pythia8, and Delphes into workflows steered by toml configuration files, and bundled into the mapyde Python package. The use of mapyde is demonstrated by constraining previously untested SUSY models with compressed sleptons and electroweakinos using ATLAS results.

1 Introduction

Direct searches for new phenomena at the Large Hadron Collider (LHC) have constrained many models of beyond-the-Standard-Model (BSM) physics. A typical LHC search focuses on an experimental signature chosen for its ability to discriminate between SM and BSM sources, and benchmarks the results of the search using one or more specific models. In some cases the chosen benchmarks are representative models that solve a particular problem in particle physics, while in other cases the benchmark models are simplified models [LHCNewPhysicsWorkingGroup:2011mji] that represent a broader (but still limited) parameter space within a theoretical framework such as supersymmetry (SUSY). In either case, the experimental results as published by LHC collaborations are only strictly applicable to their benchmark models, leaving the vast majority of BSM parameter space unexplored. Leveraging the full power of LHC data in the search for new physics requires tools that facilitate the re-use of those experimental results to test models that were not considered in the original search. The challenge to experiments is to publish or provide enough information to enable such efforts, and the challenge to the rest of the community is to use that information to extend our understanding of what BSM theories are still viable. These challenges form the reinterpretation problem.

Existing toolkits solve the reinterpretation problem in various ways, as described in Section III of Ref. [LHCReinterpretationForum:2020xtr] and in Refs. [Cranmer:2021urp, Bailey:2022tdz]. Some toolkits implement existing analysis workflows in independent software frameworks which are more simulation-based: checkmate [Dercks:2016npn], MadAnalysis5 [Conte:2018vmg, Araz:2020lnp, Dumont:2014tja], GAMBIT’s ColliderBit [GAMBIT:2017yxo, Kvellestad:2019vxm, GAMBIT:2018gjo, zenodo:gambit], rivet [Bierlich:2019rhm, Bierlich:2020wms], and Contur [Buckley:2021neu] (which interprets rivet outputs). Others match simplified versions of full models to experimental results using efficiency maps, relying more on experimental data: SModelS [Alguero:2021dig] and recast-based approaches [Cranmer:2010hk] such as in Refs. [zenodo:LHCreinterpretation, llpRepo, RECAST1, RECAST2, RECAST3]. The CMS collaboration provides Simplified Likelihood [CMS-NOTE-2017-001] correlation/covariance matrices for some analyses, while the ATLAS collaboration has started to provide full probability models [ATL-PHYS-PUB-2019-029] in which additional data is encoded, such as correlations and background estimates. The recast framework [Cranmer:2010hk] facilitates full production and processing of simulated events using the same software tools used by ATLAS in the physics results of interest. This has the advantage of high accuracy and precision, but requires significant computational resources to produce fully simulated and reconstructed samples.

Recent progress in releasing full public probability models has further expanded the possibilities for reinterpreting LHC analyses. It is no longer a technical challenge to assess the sensitivity of an LHC analysis with a public serialized probability model, including all uncertainties and correlations, to an arbitrary model of new physics. The only challenge that remains is to evaluate the acceptance and efficiency of the analysis to the new model. This can often be done using public event generation and detector simulation tools, informed by the experimental details of the reference analysis.

In this paper we present a new pipeline for calculating the constraints on new physics from existing analyses with public probability models. This pipeline is built using mapyde [mapyde], a pure-Python package that chains public tools in HEP with a single configuration file. Additionally, mapyde provides a user-friendly interface to configure, run, and extend the toolchain with other tools as needed. The mapyde toolkit is described in Section 2, including the software employed, the configuration, and deployment in containerized environments. Two example uses of the mapyde toolkit are provided in Section 3, in which we reproduce existing ATLAS results, test a new simplified model of slepton-wino-bino production, and run a pMSSM-like scan of electroweak SUSY model parameters to test SUSY models with mixed wino-bino-higgsino states.

2 The mapyde toolkit

Figure 1: Overview of the mapyde toolchain and the role of each component.

This paper describes the implementation of, and results obtained with, version v0.5.0 of mapyde. This package can be installed via pip [pypi] or conda [anaconda], as shown in Figure 2. Figure 3 provides some examples of how the authors intend this software to be executed. The mapyde package provides support for using the following tools:

• MadGraph5_aMC@NLO (MadGraph) [Alwall:2014hca, Frederix:2018nkq] (event generation)

• Pythia8 [Bierlich:2022pfr] (parton shower, hadronization, decays)

• Delphes [deFavereau:2013fsa, Selvaggi:2014mya, Mertens:2015kba] (detector simulation)

• SimpleAnalysis [simpleanalysis, atlas_simpleanalysis] (analysis description)

• pyhf [pyhf, pyhf_joss] (probability model fitting)

and was named after the first three tools in a typical mapyde simulation pipeline: MadGraph, Pythia8, and Delphes (MaPyDe). Additional tools, such as those listed in Section 1, can be supported in the analysis pipeline of this flexible framework with a little extra custom configuration, and they can be more natively supported in future versions of mapyde upon request or by pull requests from external contributors. A flow chart illustrating the role of different tools in the mapyde pipeline is shown in Figure 1.

⬇

python -m pip install mapyde # via pypi

pipx install mapyde # via pypi

conda install -c conda-forge mapyde # via conda-forge

mamba install mapyde # via conda-forge

Figure 2: Snippet illustrating various ways to install mapyde.

A mapyde pipeline is constructed as a series of containerized transforms, with the mapyde package providing configuration and steering of the individual stages. The use of containers allows the mapyde package to remain lightweight and easy to install, while still enabling event generation/simulation and subsequent data analysis steps to run on any system that provides Docker or Singularity/Apptainer. Default containers are provided for MadGraph+Pythia8, for Delphes, and for probability model fitting with pyhf. The ATLAS collaboration provides containers that implement the SimpleAnalysis routines for many SUSY searches, including those reinterpreted in Section 3. Containers for a variety of different MadGraph+Pythia8 releases are provided in the mapyde GitHub container registry [MapydeRegistry]. Users can also provide their own containers for any stage of the analysis.

⬇

$ mapyde config parse user.toml # parse and print configuration to screen

$ mapyde –prefix cards # print the on-disk path for configuration cards

$ mapyde run all user.toml # run the default workflow for user.toml

$ mapyde run madgraph user.toml # run the madgraph step for user.toml

Figure 3: Some example commands for running mapyde v0.5.0 from the command-line. From top to bottom, (a) parse the entire configuration and print a JSON-serialized representation to the screen, (b) print the prefix for where to find configuration cards on the installed machine, (c) run the default workflow steps & tools provided by mapyde for the given configuration, and (d) run only MadGraph step of the workflow using the provided configuration.

The mapyde workflow is controlled by a user-generated dictionary that encodes the configuration for all stages of the analysis. Templates for this dictionary are provided, and custom templates can be used to provide a consistent reference for subsequent analysis runs that modify a small number of configuration parameters. The configuration dictionary offers direct access to selected parameters in the MadGraph run card, in addition to allowing the user to point towards custom process and model parameter cards. Analyses performed with mapyde can use either a command-line interface or take advantage of direct access to Python functions. The command-line interface is implemented with the click Python package [click], and it takes a toml [toml] configuration file as input, which is translated internally into the configuration dictionary. As shown in Figure 4, mapyde also provides a Text User Interface (TUI) using Textual [textual] to support users less familiar with the command-line interface. The Python interface to mapyde can either parse a toml input file or accept the configuration dictionary directly. An example toml configuration file is provided in Appendix C.

Figure 4: Screenshot of an iTerm2 terminal showing the Textual User Interface option of mapyde called from the command-line via mapyde --tui.

In the containers provided by mapyde, MadGraph and Pythia8 are bundled together and are both launched through a MadGraph control file that is generated based on the mapyde configuration dictionary. The MadGraph outputs are stored as Les Houches Event (LHE) records [Alwall:2006yp], and are passed directly into Pythia8 for parton showering and hadronization. The Pythia8 outputs in hepmc format are then passed into Delphes, or directly into SimpleAnalysis when evaluating the acceptance of an analysis selection without detector simulation. Outputs from Delphes in root format are passed to an intermediate stage that either analyzes the Delphes output or transforms it into a format able to be processed by the next step in the pipeline. In the results presented below, we use a script that converts the Delphes output into a root file that can serve as input to the SimpleAnalysis phase, called Delphes2SA.py. The output of SimpleAnalysis is then processed into json format in another custom script called SA2JSON.py. That script encodes a mapping from the branch names in the SimpleAnalysis output file to the signal region names in the public probability model provided with the published analysis. The result of that transform can then be used to patch the serialized probability model and perform hypothesis tests with the generated signal. In the cases presented below, we perform the hypothesis test with the muscan.py script, whose output is a json file containing both the 95% confidence level upper limits on the cross section (expressed as a ratio of the cross section to the model prediction, called the signal strength and denoted $\mu_{\mathrm{sig}}$), as well as the full configuration dictionary used for the job. These outputs can then be used to determine if a given set of parameters used to generate events for the hypothesis test are compatible with existing limits or not. All outputs from intermediate stages of the analysis chain are kept to facilitate re-running the pipeline from an arbitrary starting point.

3 Reinterpreting Compressed SUSY Searches from ATLAS

We demonstrate the utility of the mapyde analysis chain by reproducing and reinterpreting the results of ATLAS searches for supersymmetry. The searches described in Ref. [ATLAS:2019lng] are optimized for SUSY models with “compressed” mass spectra, where the next-to-lightest SUSY partner (NLSP) and lightest SUSY partner (LSP) are separated by $\mathcal{O}(1-10)$ GeV in mass. The small mass splittings imply low-momentum (soft) Standard Model decay products, since most of the momentum of the NLSP is given to the LSP. The ATLAS searches considered here focus on final states including two low-momentum (“soft”) charged leptons, substantial missing transverse momentum ($E_{\mathrm{T}}^{\mathrm{miss}}$) from the invisible LSP’s (which in this case are the lightest neutralinos, $\widetilde{\chi}_{1}^{0}$), and one or more energetic jets that boost the SUSY system. We focus specifically on two searches from that paper: a search optimized for light sleptons (SUSY partners of the SM leptons), and a search optimized for light electroweakinos (SUSY partners of the SM electroweak bosons).

3.1 Implementation

We use mapyde to generate, shower, and simulate SUSY events, to analyze them using SimpleAnalysis, and to interpret the results using pyhf. The implementation corresponds to mapyde version 0.5.0, with configuration cards and scripts provided in a public GitHub repository [mapyde-user]. The configuration uses MadGraph version 2.9.3, Pythia8 version 8.306, Delphes version 3.5.0, and SimpleAnalysis version 1.1.0. The SimpleAnalysis code runs on dedicated containers provided by ATLAS [SAGitLabRegistry], where we use the “EwkCompressed2018” selection corresponding to the analyses from Ref. [ATLAS:2019lng]. We use pyhf version 0.7.2 to patch the public probability models provided in the HEPData [HepData] repository [hepdata.91374] for Ref. [ATLAS:2019lng] with the signal yields from mapyde, and to compute upper limits on $\mu_{\mathrm{sig}}$. When reproducing the results for the same benchmark models from Ref. [ATLAS:2019lng] the pyhf output is compared with limit contours provided in HEPData.

The signal samples produced in our mapyde workflow differ from those used in Ref. [ATLAS:2019lng] in three significant ways. First, events in the ATLAS samples contained up to two jets in the matrix element in addition to a pair of SUSY particles, with different jet multiplicity processes merged in the parton shower using the CKKW-L algorithm [Lonnblad:2011xx]. In our samples, we produce only one-jet events in MadGraph, and allow Pythia8 to model the contributions from additional QCD emissions. Second, the ATLAS samples use MadSpin [Artoisenet:2012st] to perform the three-body decays of the electroweakinos to SM leptons and a $\widetilde{\chi}_{1}^{0}$, while we perform decays using Pythia8. (The mapyde pipeline does support the use of MadSpin, but it is not used here.) Third, and most importantly, the ATLAS samples use ATLAS simulation and reconstruction software to transform the Pythia8 output into ROOT files containing physics objects, while we use Delphes.

The impact of the first two differences is evaluated by comparing the acceptance of the event selection applied to events at particle level. ATLAS provides the acceptance of the event selection by processing particle-level events with SimpleAnalysis as part of the public event record of the search. We perform a similar calculation by passing the hepmc event record from Pythia8 directly into SimpleAnalysis using mapyde. The resulting signal region yields are then compared to the product of the ATLAS acceptance, cross section, branching ratio [Fiaschi:2023tkq, Fuks:2013vua, Debove:2010kf, Fuks:2012qx, Fiaschi:2018hgm, Beenakker:1999xh], and integrated luminosity. On average we find that the mapyde and ATLAS acceptances are very similar, except for very compressed points where on average the ATLAS acceptance is approximately 10% higher.

The final difference between the ATLAS simulation framework and mapyde is evaluated using signal yields and model constraints after detector simulation and tuned by modifying the efficiencies in the Delphes configuration card. We focus in particular on the treatment of electrons and muons in Delphes, since these also required special handling in the ATLAS analysis. Further details of the lepton efficiency tuning in Delphes are described below.

3.2 Compressed Sleptons

We first reproduce the results of the ATLAS search for compressed sleptons. The ATLAS search is optimized using a slepton-bino model, in which the slepton NLSP decays to a bino-like LSP and a charged SM lepton. We generate events with pairs of charged sleptons (including scalar superpartners of both left- and right-handed SM leptons) and compare the leading-order (LO) cross sections calculated by MadGraph in an inclusive sample (without additional emissions) with the next-to-leading-order (NLO) cross sections reported by ATLAS, which do not include electroweak corrections. We find a mass-independent NLO:LO $k$-factor of 1.18, which is used to scale the one-jet MadGraph samples produced with mapyde.

After the acceptance corrections described above, the lepton efficiencies in Delphes are tuned to reproduce the ATLAS signal yields as documented in Figures 14 and 15. The efficiencies of electrons and muons after object selection are provided as part of the public record for Ref. [ATLAS:2019lng] and are the starting point for the modified Delphes configuration used to reproduce the ATLAS results. We find that setting all electron and muon efficiencies in Delphes to the upper range of values reported in Fig. 3 of [ATLAS:2019lng] is sufficient to reproduce the results of the ATLAS slepton and electroweakino searches to adequate precision, as shown in Figures 5 and 10. The only exception is the lowest-$p_{\mathrm{T}}$ bin for both electrons and muons, where the mapyde efficiencies needed to reproduce the ATLAS limits are roughly 15% higher than the values reported by ATLAS. Figure 5 also shows the model constraints for the slepton-bino search when using the default Delphes configuration card, which does a reasonable job of describing the high-splitting regions, but fails to describe the most compressed mass points, which are most sensitive to the low-$p_{\mathrm{T}}$ lepton efficiencies provided by ATLAS.

Figure 5: Constraints on the slepton-bino model from: the ATLAS paper [ATLAS:2019lng] (blue dots); mapyde before tuning the Delphes lepton efficiencies (gray line); and mapyde after tuning (blue line). The color map shows the relative difference in the limits on the SUSY signal strength, $\mu_{\mathrm{SUSY}}$, between mapyde and ATLAS results.

With the framework tuned to the ATLAS response for compressed SUSY events, we now use mapyde to assess the sensitivity of the ATLAS search to a new model. We consider a “slepton-wino-bino” simplified model, where a light slepton decays to a wino and a SM charged lepton, as illustrated in Figure 6. We set the slepton branching ratio $\tilde{\ell}\to\ell\widetilde{\chi}_{2}^{0}$ to 100%, with the $\widetilde{\chi}_{2}^{0}$ decaying to a $\widetilde{\chi}_{1}^{0}$ and to two fermions through an off-shell $Z$ boson. In such events, the $p_{\mathrm{T}}$ values of the SM leptons will be largely determined by the $\widetilde{\ell}$–$\widetilde{\chi}_{2}^{0}$ mass gap, in contrast to the slepton-bino model where the lepton $p_{\mathrm{T}}$ is driven by the $\widetilde{\ell}$–$\widetilde{\chi}_{1}^{0}$ mass gap. In the limit of vanishing $\widetilde{\chi}_{2}^{0}$–$\widetilde{\chi}_{1}^{0}$ mass differences, the slepton-wino-bino model described here is phenomenologically identical to the slepton-bino model.

Figure 6: Feynman diagram of slepton pair production, with slepton ($\tilde{\ell}$) decays to wino-like NLSP’s ($\tilde{\chi}_{2}^{0}$), which then decay to bino-like LSP’s ($\tilde{\chi}_{1}^{0}$).

(a)

(b)

(c)

Figure 7: Mass hierarchies for large $\Delta m(\widetilde{\ell},\widetilde{\chi}_{1}^{0})$ splittings for different splittings of $\Delta m(\widetilde{\chi}_{2}^{0},\widetilde{\chi}_{1}^{0})$ being (a) large, (b) intermediate, and (c) small. In all figures, $m(\widetilde{\ell})=130$ GeV.

We parameterize the model in terms of the slepton mass, the $\widetilde{\ell}$-$\widetilde{\chi}_{1}^{0}$ mass gap, and the $\widetilde{\chi}_{2}^{0}$-$\widetilde{\chi}_{1}^{0}$ mass gap, to allow us to display the model constraints in the same plane as the slepton-bino results from ATLAS. The resulting constraints (both expected and observed) are shown in Figures 8 and 9 for a range of $\widetilde{\chi}_{2}^{0}$-$\widetilde{\chi}_{1}^{0}$ splittings. At small $\Delta m(\widetilde{\chi}_{2}^{0},\widetilde{\chi}_{1}^{0})$ the results resemble those of the ATLAS slepton-bino model, since the $\widetilde{\chi}_{2}^{0}$ is nearly degenerate with the $\widetilde{\chi}_{1}^{0}$ leading to virtually no kinematic difference arising from the additional soft decay products. By contrast, for large $\widetilde{\chi}_{2}^{0}$-$\widetilde{\chi}_{1}^{0}$ splittings (dark blue contour lines), the requirement that the slepton must decay through a $\widetilde{\chi}_{2}^{0}$ results in stronger constraints at large $\widetilde{\ell}$-$\widetilde{\chi}_{1}^{0}$ differences. For each model considered, the contour is prevented from covering arbitrarily small values of $\Delta m(\widetilde{\ell},\widetilde{\chi}_{2}^{0})$ by the fact that such models produce softer SM leptons and lower lepton efficiences. Similarly, the contour is bounded from above by the reduced cross-section due to larger slepton masses, and by ATLAS analysis selections, which were optimized for compressed decays, and which are sensitive to the momenta of the fermions from the $\widetilde{\chi}_{2}^{0}$ decay. For models with large $\Delta m(\widetilde{\ell},\widetilde{\chi}_{1}^{0})$, the second-lightest neutralino $\widetilde{\chi}_{2}^{0}$ has more allowed phase-space illustrated in Figure 7.

(a)

(b)

Figure 8: Expected (a) and observed (b) constraints on the slepton-wino-bino model. The model is parameterized by the slepton mass, the slepton-$\widetilde{\chi}_{1}^{0}$ mass splitting, and the $\widetilde{\chi}_{2}^{0}$-$\widetilde{\chi}_{1}^{0}$ splitting, and compared against the slepton-bino results from Ref. [ATLAS:2019lng].

(a)

(b)

Figure 9: Expected (a) and observed (b) constraints on the slepton-wino-bino model, shown in logarithmic scale. The model is parameterized by the slepton mass, the slepton-$\widetilde{\chi}_{1}^{0}$ mass splitting, and the $\widetilde{\chi}_{2}^{0}$-$\widetilde{\chi}_{1}^{0}$ splitting, and compared against the slepton-bino results from Ref. [ATLAS:2019lng].

This class of models, while being less “simplified” than the slepton-bino model traditionally studied by LHC experiments, includes models that nominally escape constraint by both existing slepton-bino searches as well as direct searches for electroweakinos. In particular, the wino-bino constraints reported in Ref. [ATLAS:2019lng] constrain electroweakinos with $\widetilde{\chi}_{2}^{0}$-$\widetilde{\chi}_{1}^{0}$ splittings as low as 1 GeV only for electroweakino masses near 100 GeV. In the slepton-wino-bino results shown in Fig. 8 we find that the ATLAS search excludes models with 1 GeV $\widetilde{\chi}_{2}^{0}$-$\widetilde{\chi}_{1}^{0}$ splittings up to slepton masses of well over 200 GeV, which implies limits on $\widetilde{\chi}_{1}^{0}$ masses also exceed 200 GeV for small slepton-$\widetilde{\chi}_{1}^{0}$ splittings.

3.3 Compressed Electroweakinos

The ATLAS search for compressed electroweakinos considered two different simplified models: one with “pure” Higgsino-like $\widetilde{\chi}_{2}^{0}$, $\widetilde{\chi}_{1}^{\pm}$, and $\widetilde{\chi}_{1}^{0}$ states, and another with wino-like $\widetilde{\chi}_{2}^{0}$/$\widetilde{\chi}_{1}^{\pm}$ and bino-like $\widetilde{\chi}_{1}^{0}$. We focus on the Higgsino model to validate the mapyde output. Following a procedure similar to the validation of the ATLAS slepton-bino results described above, we generate a grid of Higgsino model points to reproduce the ATLAS results. The mapyde electroweakino samples use the same $k$-factor and lepton efficiencies as the slepton search. The branching ratio of the Higgsino-like $\widetilde{\chi}_{2}^{0}$ to leptons is taken from LHC SUSY Cross Section Working Group [susyxs]. The comparison between the mapyde and ATLAS results is shown in Figure 10, where the mapyde exclusion contour follows the corresponding ATLAS contour.

Figure 10: Constraints on the “pure Higgsino” simplified model from ATLAS (blue dots) and mapyde (blue line). The color map shows the relative difference in the limits on the SUSY signal strength, $\mu_{\mathrm{SUSY}}$, between mapyde and ATLAS results.

(a)

(b)

Figure 11: Example feynman diagrams of possible electroweakino decays with (a) $\widetilde{\chi}_{2}^{\pm}\to\widetilde{\chi}_{1}^{0}$ and (b) $\widetilde{\chi}_{2}^{\pm}\to\widetilde{\chi}_{2}^{0}$.

We next use mapyde to assess the ATLAS sensitivity to MSSM SUSY models in which the bino, wino, and Higgsino mass terms ($M_{1}$, $M_{2}$, and $\mu$, respectively) are all relatively low, leading to “wino-bino-Higgsino” models that have potentially rich phenomenologies. Some example models are illustrated in Figure 11. In these models, the $W$-boson is treated as off-shell for the “compressed” phase-space under study in this paper. We use a pMSSM scanning tool (EasyScan_HEP [Shang:2023gfy, Han:2016gvr]) to generate particle spectra for models with $M_{1}$, $M_{2}$, and $\mu$ ranging from $-500$ GeV to 500 GeV, sampled with flat priors. Particle masses are calculated using Spheno [Porod:2003um, Porod:2011nf] and stored in the SUSY Les Houches Accord (SLHA) format [Allanach:2008qq]. Additional tools are run as described in Appendix B and are used to calculate values for the selection criteria shown in Figure 12. Since our goal is to investigate models with compressed electroweakino mass spectra accessible to ATLAS Run-2 searches, we select models that satisfy $m(\widetilde{\chi}_{1}^{0})>100$ GeV, $m(\widetilde{\chi}_{3}^{0})<300$ GeV, and $(m(\widetilde{\chi}_{3}^{0})-m(\widetilde{\chi}_{1}^{0}))<50$ GeV for further study. We further require that any selected models have valid output from the spectrum generator Spheno [Porod:2003um, Porod:2011nf], have a valid Higgs mass as computed by FeynHiggs [Bahl:2018qog, Bahl:2017aev, Bahl:2016brp, Hahn:2013ria, Frank:2006yh, Degrassi:2002fi, Heinemeyer:1998np, Heinemeyer:1998yj], and to satisfy dark matter constraints ($\Omega h^{2}<0.12$) implemented in MicrOMEGAS [Belanger:2020gnr], flavor physics constraints implemented in SuperIso [Arbey:2018msw], and optionally muon $g-2$ constraints implemented in GM2Calc [Athron:2015rva, Athron:2021evk]. The SLHA records for the selected 81 models are used as inputs for event generation with mapyde. Branching ratios for electroweakino decays are taken directly from the Spheno output.

⬇

(SP_m_h!=-1) & (SPfh_m_h!=-1) # spheno, feynhiggs: ok

& (SI_BR_Bs_to_mumu!=-1) & (GM2_gmuon!=-1) # superiso, gm2calc: ok

& (MO_Omega!=-1) & (MO_Omega < 0.12) # micromegas: ok, DM relic density

& (SP_m_chi_10>100) & (abs(SP_m_chi_30)<300) # N1 > 100 GeV, N3 < 300 GeV

& ((abs(SP_m_chi_30)-abs(SP_m_chi_10))<50) # m(N3, N1) < 50 GeV

Figure 12: The mask used to define the selection of models for assessing the ATLAS sensitivity to Higgsino-Win We next use mapyde to assess the ATLAS sensitivity to MSSM SUSY models in which the bino, wino, and Higgsino o models.

The ATLAS constraints on the wino-bino-higgsino model scan are parameterized in terms of the mass of the $\widetilde{\chi}_{2}^{0}$ and the mass difference $\Delta m=m(\widetilde{\chi}_{2}^{0})-m(\widetilde{\chi}_{1}^{0})$, to facilitate comparisons with the pure-Higgsino constraints from ATLAS. The models under study are binned in a coarse grid in the $\Delta m$ vs $m(\widetilde{\chi}_{2}^{0})$ plane, and the fraction of excluded models is calculated for each bin. The expected and observed constraints on these models are shown in Figure 13, with the ATLAS Higgsino results overlaid for comparison. In general the selected pMSSM points are more constrained than those of a pure-Higgsino model at larger mass splittings, likely due to the presence of the additional $\widetilde{\chi}_{3}^{0}$ and its own decays to final states similar to that of the $\widetilde{\chi}_{2}^{0}$.

(a)

(b)

Figure 13: Expected (a) and observed (b) constraints on wino-bino-higgsino models containing at least three light neutralinos. Black dots represent the specific model points used to calculate the fraction of excluded models within a bin, shown in color. White bins contain no models from the current scan. Results from an ATLAS search [ATLAS:2019lng] are overlaid for comparison.

4 Conclusion

The combination of rich physics results from the LHC with public supplementary material, including analysis routines and probability models, has opened the door to new ways of quantifying constraints on unexplored models of BSM physics with already-published results. We demonstrated one such pipeline, implemented in the mapyde Python package, that facilitates the re-use of LHC analyses. We illustrated the utility of a user-friendly Python package, mapyde, by probing previously-untested models of supersymmetry, including a simplified model of sleptons that undergo cascade decays to wino-like and bino-like electroweakinos, and a parameter scan of highly-mixed higgsino, wino, and bino states with a rich set of possible decay chains. In both cases we find that existing searches are able to constrain non-trivial regions of parameter space that would have been difficult or impossible to predict by simple extrapolations of existing results.

5 Acknowledgments

We thank Sam English, Bryn Lonsbrough, Len Morales Zaragoza, Zach Dethloff, and Ryota Johnson for suffering with early versions of scripts that would later develop into mapyde. We thank Zach Marshall, Jeff Shahinian, and Jason Nielsen for reading early drafts of this work and providing useful comments. Hance and Stark are supported by the Department of Energy Office of Science grant DE-SC0010107.

Appendix A Lepton Efficiencies in Delphes

Figure 14: A snippet from the mapyde-provided Delphes configuration for electron efficiencies. The portion that is commented out is what comes from the default ATLAS configuration. These numbers come from Ref. [ATLAS:2019lng].

Figure 15: A snippet from the mapyde-provided Delphes configuration for muon efficiencies. The portion that is commented out is what comes from the default ATLAS configuration. These numbers come from Ref. [ATLAS:2019lng].

Appendix B EasyScan_HEP Configuration

We used an EasyScan_HEPini configuration file that defined the scan ranges in Figure 16 as well as additional programs to run, described in the text above. The specific versions used are:

• EasyScan_HEP (v1.0.0): pMSSM scanning and program control [Shang:2023gfy, Han:2016gvr]

• Spheno (v4.0.4): spectrum generator [Porod:2003um, Porod:2011nf]

• FeynHiggs (v2.16.0): Higgs mass calculation [Bahl:2018qog, Bahl:2017aev, Bahl:2016brp, Hahn:2013ria, Frank:2006yh, Degrassi:2002fi, Heinemeyer:1998np, Heinemeyer:1998yj]

• MicrOMEGAS (v5.2.1): Dark Matter calculations (e.g., relic density) [Belanger:2020gnr]

• SuperIso (v4.0): Flavor Physics observables [Arbey:2018msw]

• GM2Calc (v2.0.0): $g-2$ calculation [Athron:2015rva, Athron:2021evk]

⬇

[scan]

Scan method: random

# ID Prior Min MAX

Input parametes: tanb, Flat, 1, 60

M_1, Flat, -500, 500

M_2, Flat, -500, 500

M_3, Flat, 2000, 2000

AT, Flat, 2000, 2000

Ab, Flat, 2000, 2000

Atau, Flat, 2000, 2000

MU, Flat, -500, 500

mA, Flat, 2000, 2000

meL, Flat, 2000, 2000

mtauL, Flat, 2000, 2000

meR, Flat, 2000, 2000

%**** main.tex Line 350 ****mtauR, Flat, 2000, 2000

mqL1, Flat, 2000, 2000

mqL3, Flat, 2000, 2000

muR, Flat, 2000, 2000

mtR, Flat, 2000, 2000

mdR, Flat, 2000, 2000

mbR, Flat, 2000, 2000

Figure 16: A portion of the easyscan.ini configuration defining the random sampling for the electroweakinos scan.

Appendix C Input toml Configuration for a Slepton Sample

The configuration dictionary used within mapyde can be created as a Python dictionary and passed to mapyde through the Python interface, or can be generated from a toml configuration file. We provide an example toml configuration file below. The code in this section represents a single file, but is described in blocks for easier interpretation.

C.1 The base block

The base block provides mapyde with paths for inputs and outputs of the analysis pipeline. Users can take advantage of pre-generated configuration cards Pythia8 and Delphes, but will usually at least need to define their own process cards for generating events with MadGraph.

⬇

[base]

%**** main.tex Line 375 ****path = "{{PWD}}"

output = "output"

logs = "logs"

data_path = "{{MAPYDE_DATA}}"

cards_path = "{{MAPYDE_CARDS}}"

scripts_path = "{{MAPYDE_SCRIPTS}}"

process_path = "{{MAPYDE_CARDS}}/process/"

param_path = "{{MAPYDE_CARDS}}/param/"

run_path = "{{MAPYDE_CARDS}}/run/"

pythia_path = "{{MAPYDE_CARDS}}/pythia/"

delphes_path = "{{MAPYDE_CARDS}}/delphes/"

madspin_path = "{{MAPYDE_CARDS}}/madspin/"

likelihoods_path = "{{MAPYDE_LIKELIHOODS}}"

Figure 17: The base block of an example toml configuration file for generating slepton events.

Note that in Figure 17, the template is left undefined, which means that the base template shipped with mapyde will be loaded in. This behavior can be disabled by setting template = false. The (nested) template inheriting behavior for toml configuration parsing enables defining grid scans more cleanly with less duplicated code.

C.2 The madgraph block

The madgraph block defines some top-level options for running madgraph, such as:

• whether or not to skip the MadGraph job, sometimes useful when re-running an analysis chain and re-using LHE inputs from a previous job

• the name of the param card to use, in this example making use of a toml feature that allows references to other toml parameters

• the number of jobs to run in parallel when generating MadGraph events and when processing those events in Pythia8

• the name of the container that provides the version of MadGraph to use, in this case version 2.9.3 (taken from the GitHub container registry)

These options, and options from subsequent MadGraph-related blocks, are processed by mapyde to produce a MadGraph run script. The madgraph.generator block only specifies the name of that run script.

⬇

[madgraph]

skip = false

params = "SleptonBino"

paramcard = "{{madgraph[’params’]}}.slha"

cores = 1

batch = false

version = "madgraph:2.9.3"

\par[madgraph.generator]

output = "run.mg5"

Figure 18: The madgraph block of an example toml configuration file for generating slepton events.

C.3 The madgraph.masses block

The madgraph.masses block enables the use of user-defined substitutions in the SLHA param card.

⬇

[madgraph.masses]

MSLEP = 250

MN1 = 240

Figure 19: The madgraph.masses block of an example toml configuration file for generating slepton events.

In the example above, the slepton masses in the SLHA file have been replaced with the string ’MSLEP’, which mapyde then substitutes with the value of 250 GeV at run-time. This allows a single param-card template to be used for jobs that want to keep most parameters fixed, but vary one or more parameters as part of a parameter scan.

C.4 The madgraph.run blocks

The madgraph.run blocks exposes the MadGraph run card to the mapyde configuration dictionary. The madgraph.run block defines the name of the MadGraph run card, which will be accessed via the run_path defined in the base block as well as some commonly-changed run parameters: nevents, iseed, and ecms. The madgraph.run.options block defines modifications to be made to the MadGraph run card. In the example below, default values for mmjj, ptj, and ptj1min are overwritten by the values provided in the toml configuration file.

⬇

[madgraph.run]

card = "default_LO.dat"

ecms = 13000

nevents = 50000

seed = 0

\par[madgraph.run.options]

mmjj = 500

ptj = 20

ptj1min = 50

Figure 20: The madgraph.run blocks of an example toml configuration file for generating slepton events.

C.5 The madgraph.proc block

The madgraph.proc block defines the name of the process card that contains the hard process information for MadGraph. Figure 21 shows an example using an existing process card, while  Figure 22 creates a new one from the specified contents.

⬇

[madgraph.proc]

name = "isrslep"

card = "{{madgraph[’proc’][’name’]}}"

contents = false

Figure 21: The madgraph.proc block of an example full toml configuration file for generating slepton events.

⬇

[madgraph.proc]

name = "isrslep"

card = false

contents = """\ set …

define susystrong …

define …

generate p p > chsleptons chsleptons j / susystrong @1

output -f

"""

Figure 22: The madgraph.proc block demonstrating how to create a process card on-the-fly using contents to generate slepton events.

C.6 The madspin block

The madspin block defines the name of the MadSpin card, in jobs where MadSpin is run. The running of MadSpin can be skipped if it is not needed.

⬇

[madspin]

skip = true

card = ”

Figure 23: The madgraph.run block of an example toml configuration file for generating slepton events.

C.7 The pythia block

The pythia block holds all configuration information for Pythia8. It defines the name of the Pythia8 configuration file, controls whether Pythia8 is run or not, and allows the passing of some additional options that will be appended to the Pythia8 configuration card.

⬇

[pythia]

%**** main.tex Line 525 ****skip = false

card = "pythia8_card.dat"

additional_opts = ""

Figure 24: The pythia block of an example toml configuration file for generating slepton events.

C.8 The delphes block

The delphes block holds all configuration information for the Delphes stage. It defines the name of the tcl configuration card, the name of the container that provides Delphes(version), and the name of the output file. The input is assumed to be the hepmc output from a previous Pythia8 job.

⬇

[delphes]

skip = false

card = "delphes_card_ATLAS_lowptleptons_sleptons_notrackineffic.tcl"

version = "delphes"

output = "delphes/delphes.root"

Figure 25: The delphes block of an example toml configuration file for generating slepton events.

C.9 The analysis block

The analysis block defines any transform that comes after the Delphes stage. In the example code below, the script Delphes2SA.py will transform the root output from Delphes into a root file that can be processed by SimpleAnalysis. This stage also allows the implementation of scale factors (kfactor) that are applied to the event weights from the Delphes output. Scaling of the Delphes outputs to a total integrated luminosity is also facilitated by the lumi flag. The cross section used to calculate event weights can be over-ridden with the XSoverride flag, which is multiplied by the kfactor to define the final cross section. Users can also provide their own scripts to run at this stage, which can choose to implement or ignore the options in this example configuration. Running multiple scripts can be accomplished in at least two different ways: by wrapping multiple scripts and having mapyde call the wrapper, or by calling mapyde multiple times with separate configuration files that differ only in the script name.

The code defined with the script flag is run within the same container as the Delphes stage, to allow the use of Delphes libraries for analyzing the output root file.

⬇

[analysis]

script = "Delphes2SA.py"

XSoverride = -1

kfactor = 1.18

output = "analysis/Delphes2SA.root"

lumi = 139000

Figure 26: The analysis block of an example toml configuration file for generating slepton events.

C.10 The simpleanalysis block

The simpleanalysis block defines the inputs, outputs, and analysis code to be run in the SimpleAnalysis stage. The output filenames will be given the name of the SimpleAnalysis algorithm (here EwkCompressed2018, also see Figure 27), with the optional outputtag appended. This is particularly useful in cases where the input is specified as hepmc, indicating the output from Pythia8 should be analyzed, rather than the output from Delphes. In such cases, specifying an outputtag like _hepmc can help to distinguish between those results and the results of running SimpleAnalysis on Delphes output. Currently mapyde only supports the definition of a single SimpleAnalysis stage in the pipeline, so multiple SimpleAnalysis transforms, e.g., one for Delphes inputs and one for hepmc inputs, need to be performed using separate configurations.

⬇

#include "SimpleAnalysisFramework/AnalysisClass.h"

DefineAnalysis(EwkOneLeptonTwoBjets2018)

// Wh->1l+bb+met analysis (Run2 data)

// https://arxiv.org/pdf/1909.09226.pdf

void EwkOneLeptonTwoBjets2018::Init()

{

        // Define signal/control regions

        addRegions({"SR_h_Low", "SR_h_Med", "SR_h_High"} );

        addRegions({"SR_h_Low_Incl", "SR_h_Med_Incl", "SR_h_High_Incl"} );

        addRegions({"SR_h_Low_bin1", "SR_h_Low_bin2", "SR_h_Low_bin3"});

        addRegions({"SR_h_Med_bin1", "SR_h_Med_bin2", "SR_h_Med_bin3"});

        addRegions({"SR_h_High_bin1", "SR_h_High_bin2", "SR_h_High_bin3"});

        addRegions({"CR_tt_Low", "CR_tt_Med", "CR_tt_High"});

        addRegions({"WCR", "STCR"});

        addRegions({"VR_off_Low", "VR_off_Med", "VR_off_High"});

        addRegions({"VR_on_Low", "VR_on_Med", "VR_on_High"});

        // Preselection for debugging

        addRegions({"presel_1L", "presel_2J", "presel_bb", "presel_met", "presel_mbb"});

}

void EwkOneLeptonTwoBjets2018::ProcessEvent(AnalysisEvent *event)

{

        //

        // Soft Objects

        //

        auto baseEle = event->getElectrons(7,2.47, ELooseBLLH);

        auto baseMuon = event->getMuons(6,2.70, MuMedium | MuNotCosmic | MuZ05mm | MuQoPSignificance);

        auto baseJets = event->getJets(20.,4.5);

        auto met_Vect  = event->getMET();

        float met       = met_Vect.Pt();

        auto weights        = event->getMCWeights();

        //

        // Baseline objects

        //

        auto radiusCalcLep = [] (const AnalysisObject& lep,const AnalysisObject&) {

                                                         return (0.04 + 10/lep.Pt()) > 0.4 ? 0.4 : 0.04 + 10/lep.Pt();

                                                 };

        // ele-mu OR

        baseEle = overlapRemoval(baseEle, baseMuon, 0.01);

        // ele-jet OR

        baseJets = overlapRemoval(baseJets, baseEle, 0.2);

        baseEle = overlapRemoval(baseEle, baseJets, radiusCalcLep);

        // mu-jet OR

        baseJets = overlapRemoval(baseJets, baseMuon, 0.2, LessThan3Tracks);

        baseMuon = overlapRemoval(baseMuon, baseJets, radiusCalcLep);

        // Get signal electrons with FCLoose for pT < 200 GeV and FCHighPtCaloOnly for pT > 200

        // filterObjects() only allows lower pT cut, so work around that

        auto signalEle = filterObjects(baseEle,7., 2.47, ETightLH | ED0Sigma5  | EZ05mm  | EIsoFCLoose);

        AnalysisObjects signalEleLowPt;

    for (const auto& ele : signalEle) {

      if ((ele.Pt() < 200.)) signalEleLowPt.push_back(ele);

        }

        auto signalEleHighPt = filterObjects(baseEle,200., 2.47, ETightLH | ED0Sigma5  | EZ05mm  | EIsoFCHighPtCaloOnly);

        signalEle = signalEleLowPt + signalEleHighPt;

        auto signalMuon = filterObjects(baseMuon,6., 2.7, MuD0Sigma3 | MuZ05mm | MuIsoFCLoose);

        auto signalLept =  signalEle + signalMuon;

        auto signalJets = filterObjects(baseJets, 30., 2.80, JVT120Jet);

        auto signalBJets = filterObjects(signalJets, 30., 2.8, BTag77MV2c10);

    unsigned int N_baseEle = baseEle.size();

        unsigned int N_baseMuon = baseMuon.size();

        unsigned int N_baseLept = N_baseEle + N_baseMuon;

        unsigned int N_signalEle = signalEle.size();

        unsigned int N_signalMuon = signalMuon.size();

        unsigned int N_signalLept = N_signalEle + N_signalMuon;

        unsigned int N_signalJets = signalJets.size();

        unsigned int N_signalBJets = signalBJets.size();

        float mt=0, m_CT=0,  mbb=0, mlb1=0;

        if (signalLept.size()==1   && signalBJets.size()==2 ) {

                mt   = calcMT(signalLept[0], met_Vect);

                m_CT = calcMCT(signalBJets[0],signalBJets[1],met_Vect);

                mbb= (signalBJets[0]+signalBJets[1]).M();

                mlb1= (signalBJets[0]+signalLept[0]).M();

        }

        if(N_signalLept == 1 && N_baseLept == 1) accept("presel_1L");

        if(N_signalJets<=3 && N_signalJets >= 2) accept("presel_2J");

        if(N_signalBJets == 2) accept("presel_bb");

        if(met> 220. && mt> 50.) accept("presel_met");

        if(mbb>100 && mbb<=140) accept("presel_mbb");

        // Preselection

        if(N_baseLept != 1) return;

        if(N_signalLept != 1) return;

        if(N_signalJets>3 || N_signalJets < 2) return;

        if(N_signalBJets != 2) return;

        if(mt< 50.) return;

        if(met< 220.) return;

        ntupVar("AnalysisType", (baseMuon.size() == 1) ? 2 : 1);

        ntupVar("met",met);

        ntupVar("mct2",m_CT);

        ntupVar("mbb",mbb);

        ntupVar("mlb1",mlb1);

        ntupVar("mt",mt);

        ntupVar("nBJet",met);

        ntupVar("nJet30",(int)signalJets.size());

    ntupVar("nBJet30_MV2c10",(int)signalBJets.size());

        ntupVar("lep1Pt", signalLept[0].Pt());

    ntupVar("lep1Eta", signalLept[0].Eta());

    ntupVar("lep1Phi", signalLept[0].Phi());

        ntupVar("jet1Pt", signalJets[0].Pt());

    ntupVar("jet1Eta", signalJets[0].Eta());

    ntupVar("jet1Phi", signalJets[0].Phi());

        ntupVar("nLep_signal", (int)signalLept.size());

    ntupVar("nLep_base", (int)N_baseLept);

        // Signal regions

        if(met>240 && mbb>100 && mbb<=140 && m_CT>180) {

                if(mt>100 && mt<160) {

                        accept("SR_h_Low");

                        if(m_CT > 180 && m_CT <= 230) {

                                accept("SR_h_Low_bin1");

                        }else if(m_CT > 230 && m_CT <= 280) {

                                accept("SR_h_Low_bin2");

                        }else{

                                accept("SR_h_Low_bin3");

                        }

                }

                else if(mt>160 && mt<240) {

                        accept("SR_h_Med");

                        if(m_CT > 180 && m_CT <= 230) {

                                accept("SR_h_Med_bin1");

                        }else if(m_CT > 230 && m_CT <= 280) {

                                accept("SR_h_Med_bin2");

                        }else{

                                accept("SR_h_Med_bin3");

                        }

                }

                else if(mt>240 && mlb1 > 120) {

                        accept("SR_h_High");

                        if(m_CT > 180 && m_CT <= 230) {

                                accept("SR_h_High_bin1");

                        }else if(m_CT > 230 && m_CT <= 280) {

                                accept("SR_h_High_bin2");

                        }else{

                                accept("SR_h_High_bin3");

                        }

                }

        }

        if(met>240 && mbb>100 && mbb<=140 && m_CT>180) {

                if(mt > 100) {

                        accept("SR_h_Low_Incl");

                }

                if(mt >= 160) {

                        accept("SR_h_Med_Incl");

                }

                if(mt > 240 && mlb1 > 120) {

                        accept("SR_h_High_Incl");

                }

        }

        // Control  regions W

        if(met>240 && mbb<80 && mbb > 50 && mt<100 && mt>50 && m_CT>180) {

                accept("WCR");

        }

        // Control  regions tt

        if(met>240 && mt>100 && m_CT<180 && ((mbb < 100 && mbb > 50)|| mbb > 140)) {

                if(mt < 160) {

                        accept("CR_tt_Low");

                }else if(mt < 240) {

                        accept("CR_tt_Med");

                }else if(mlb1 > 120) {

                        accept("CR_tt_High");

                }

        }

        // Control  regions ST

        if(met>240 && m_CT>180 && mt > 100 && mbb > 195) {

                accept("STCR");

        }

        // validation regions on

        if(mbb>100 && mbb<140 && m_CT<180 && mt > 100) {

                if(mt < 160 && met > 240) {

                        accept("VR_on_Low");

                }else if(met > 240 && mt < 240 && mt > 160) {

                        accept("VR_on_Med");

                }else if(mt > 240) {

                        accept("VR_on_High");

                }

        }

        // validation regions off

        if(((mbb > 50 && mbb<80) ||( mbb>160  && mbb<195))&& m_CT>180 && mt > 100) {

                if(mt<160) {

                        accept("VR_off_Low");

                }else if(met > 240 && mt < 240 && mt > 160) {

                        accept("VR_off_Med");

                }else if(mt > 240 && ((mbb > 50 && mbb < 75)||(mbb > 165 && mbb < 195))) {

                        accept("VR_off_High");

                }

        }

        return;

}

Figure 27: A snippet of SimpleAnalysisCodes/src/ANA-SUSY-2019-08.cxx [SAGitLab] showing how the analysis name is defined.

⬇

[simpleanalysis]

name = "EwkCompressed2018"

input = ""

outputtag = ""

Figure 28: The simpleanalysis block of an example toml configuration file for generating slepton events.

C.11 The sa2json block

The sa2json block controls the translation of the SimpleAnalysis outputs in root file format to a json patch file that can be used to update the probability model for statistical inference. This stage requires an understanding of the SimpleAnalysis output format, and how to relate that output to the names of signal regions defined in the public probability model from HEPData. An example is provided by the SA2JSON.py script in mapyde, but each serialized probability model will need a custom translation script. The mapyde developers welcome contributions of such scripts to the mapyde code base.

The image tag controls the name of the container to be used when running the transform, while the input and output options point to the SimpleAnalysis outputs and output json filename, respectively. The options entry collects command-line options for the SA2JSON.py script. In this case, -c tells SA2JSON.py to perform a special selection for the compressed electroweak SUSY searches provided by the EwkCompressed2018 SimpleAnalysis selection.

⬇

[sa2json]

inputs = "{{simpleanalysis[’name’]}}{{simpleanalysis[’outputtag’]}}.root"

image = "pyplotting:latest"

output = "{{simpleanalysis[’name’]}}{{simpleanalysis[’outputtag’]}}_patch.json"

options = "-c"

Figure 29: The sa2json block of an example toml configuration file for generating slepton events.

C.12 The pyhf block

The pyhf block controls the options for running pyhf on an input probability model. The image flag controls the container to use for running pyhf. The serialized model, provided in json format, is specified with the likelihood flag. The gpu-options and other-options flags both pass command-line arguments to the script muscan.py, provided as part of the mapyde package, which performs a scan of the signal strength, $\mu_{\mathrm{sig}}$. If the container includes support for a GPU (such as the CUDA libraries provided in pyplotting:latest-cuda11 container), and a GPU is present, then pyhf will use the jax backend and perform the calculations using the GPU. In the example below, the -c option tells muscan.py to not use a GPU even if it is available (and to use the CPU instead), while the -b jax option tells muscan.py to use the jax backend for calculations.

⬇

[pyhf]

script = "muscan.py"

skip = false

likelihood = "Slepton_bkgonly.json"

image = "pyplotting:latest"

gpu-options = "-c -B jax"

other-options = ""

%**** main.tex Line 625 ****

Figure 30: The pyhf block of an example toml configuration file for generating slepton events.