A C++ program for estimating detector sensitivities to long-lived particles: Displaced Decay Counter

A series of far-detector programs have been proposed for operation at various interaction points of the Large Hadron Collider during the upcoming runs. Investigating the potential and complementarity of these experiments for new-physics searches goes through the estimation of their sensitivity to specific long-lived particle models. Here, we present an integrated numerical tool written in the C++ language and called Displaced Decay Counter , which we have created to this end and which can be used in association with MadGraph5 , Pythia8 , or any other state-of-the-art Monte-Carlo collider simulation tool. Several far-detector models have been implemented within the program, accounting for the geometry and integrated luminosity of projected detectors. Additional or more accurate designs can be easily constructed through a dedicated interface. The functionality of this tool is exemplified through the discussion of three benchmark scenarios, which we considered for the validation of the implemented detector models.


Introduction
A decade of operation at the Large Hadron Collider (LHC) has left a picture of particle physics at energies below the TeV scale which seems well compatible with a Standard Model (SM) interpretation.In particular, no conclusive evidence has emerged to corroborate the existence of heavy color-or electroweakly-charged resonances leaving massive missing-energy signatures, as expected in various models of new-physics where a Z 2 symmetry distinguishes between the standard and exotic sectors -such as R-parity in supersymmetry (SUSY)-inspired models or T -parity in constructions involving a new strongly-coupled interaction.While the ensuing constraints leave substantial parameter space (possibly inaccessible to the LHC) in the scenarios of the previous type, it is desirable to exploit as much information as TeVscale colliders can provide us with for the searches of physics signatures beyond the SM.Long-lived particles (LLPs) taking mass at the electroweak scale form a well-motivated alternative to the prompt-decay picture, more challenging to test, and may thus have been overlooked in traditional collider search strategies.A long lifetime in an exotic particle may appear as the result of the protection by an approximate symmetry (e.g. of the 'matter-parity' type), compressed phase space, or a portal to the SM involving very massive mediators, for instance.Consequently, LLPs are widely predicted in various new-physics models [1].We refer the reader to Refs.[2][3][4][5][6] for some recent reviews on LLP searches.
Several initiatives have emerged in the last few years, both from the experimental and phenomenological communities, to promote the search for LLPs in LHC-based experiments.First, a number of LLP signatures can be tested at existing detectors, such as ATLAS and CMS (see e.g.Refs.[7][8][9][10] for some recent searches in this direction), and allow to constrain displaced leptons [11], displaced vertices with missing transverse momentum [12], leptons [13], or jets [14], or heavy charged tracks [15], to restrict to a few examples.It is then possible to recast such limits and make them available for testing arbitrary models [16][17][18][19][20].However, such searches usually focus on macroscopic, but short lifetimes, so as to account for the proximity of the detector with the interaction point.Far detectors offer a complementary search strategy, as they would be sensitive to exotic particles decaying further away.In addition, the SM background is anticipated to be very restricted, or even negligible.A collection of designs for such far detectors has flourished in recent years, some of which have been approved and are scheduled for operation in the current LHC Run3.While the detail of the efficiencies is in general unavailable, it is still enlightening to examine the exclusion potential of the projected detectors and their complementarity in constraining various types of new physics models.Simulations of this type have been carried out in numerous analyses, using a more or less detailed modelization of the detectors to estimate their acceptance in purely geometrical terms (see e.g.Refs.[21][22][23][24][25] for a few existing works).Specific studies on optimization of the geometry of such far detectors can also be found in e.g.Refs.[26,27].Moreover, recasting works on the sensitivities of these experiments can be found in e.g.Refs.[28][29][30].The purpose of this paper is to present an integrated tool allowing to assess the sensitivity of multiple far detectors in an efficient and systematic fashion: Displaced Decay Counter (DDC).
Several tools have been presented recently that accomplish somewhat comparable functions to those of DDC.First, MadDump [31] is a plugin for MadGraph5 aMC@NLO [32,33] and dedicated to the computation of signal-event rates of LLPs including their decays, scattering off SM particles, and their far detection, at beam-dump experiments.It takes into account the LLP kinematics as well as the geometry of the experiments.It can also be employed for the simulation of non-beamdump processes, in principle.Next, FORESEE [34] is a Python-based code that calculates the sensitivity reach of far detectors at the LHC and FCC-hh.This package computes the kinematical distributions of LLPs with the help of tables of kinematical distributions of various SM particles, allowing to derive the signal-event rates at the LHC far detectors.The event records are generated in the HepMC format.Although currently only a limited number of models have been implemented, the code can be extended to further LLP environments.A cut functionality is available but the efficiencies of decay products are not taken into account.Then, the authors of Ref. [35] implemented the tool ALPINIST, which estimates the signal-event rates at extracted-beam experiments, for axion-like particles (ALPs) with generic couplings.Making use of Mathematica, ROOT, and Python, the code reconstructs the final-state particles of the ALP inside a detector and performs the computation with fully reconstructed final states.Furthermore, the new package SensCalc was presented in Ref. [36].This code employs a pipeline of Mathematica notebooks in order to calculate the detector acceptance, sequentially deriving the LLP distributions, computing the signal-event rates, and plotting.The LLP distributions are either pre-defined or can be computed by running Monte-Carlo (MC) simulation.Portalphysics models -heavy neutral leptons, dark scalars, dark photons, and axion-like particleshave already been considered in SensCalc, and it could be extended to further LLP models.Again SensCalc does not reconstruct the final states.Finally, the authors of Ref. [37] very recently published a Python module FastSim that allows efficient signal-event studies of LLP decays in various far detectors.The tool includes idealized reconstruction criteria for displaced vertices, such as the number of observable charged particles from the LLP decays that should pass a certain momentum threshold and also intersect at least a certain number of detector planes.
The main functionality of the tool DDC under discussion is to compute the decay probability of an LLP within a pre-defined detector volume.A collection of detector models based on LHC designs is already pre-encoded and an integrated editor allows for further construction, possibly applying to other types of colliders.MCsimulation events in the LHEF [38] or HepMC [39] 1 format can be used as input, as well as CMND input files for internal generation through a Pythia8 [40] interface.LLPs are identified for testing via their PDG code.Compared to the existing codes discussed above, our package DDC completely relies on linked tools such as MadGraph5 and Pythia8 for performing MC simulation and obtaining the LLP distributions, without any attempt at exploiting pre-generated tabulated densities.Similarly to SensCalc, DDC is not restricted to a certain type of experiments, but should be easily applicable to a wide range of experimental facilities.Furthermore, DDC achieves a high degree of versatility: it is indeed compatible with several MC-simulation tools and can be used for a wide range of theoretical models or interaction types (lepton or hadron collisions, beam-dump experiments, neutrino-nucleus scattering, heavy ion collision, etc.).In addition, an editor eases the implementation of personalizable detector models by the user.A current disadvantage of DDC is the absence of implemented final-state-particle efficiencies, accounting for the detector sensitivity to various LLP decay products.Nevertheless the code should allow enough flexibility for a later inclusion of this feature.
In the following section, we describe the essential features of the far-detector simulator DDC2 .We then apply this framework to several LLP scenarios in section 3, in order to illustrate the functionality of this tool (and validate our implementation of detector models).A brief summary is provided in section 4.

The tool 2.1 General features
The C++ program DDC provides a framework based on HepMC and Pythia that is able to simulate detector response (up to a certain extent) for collider experiments looking for LLPs.Its basic input consists in: • a PDG code [41] representing a (long-lived) particle suspected to decay in a detector far from the interaction point, together with its mass, lifetime, and 'visible' branching ratio; multiple LLPs are supported as input, as long as they are produced independently (and not in a chain); • a set of Monte-Carlo events describing the production cross-section of this particle at a collider; LHEF and HepMC formats are accepted; alternatively, events can be generated internally via a call to Pythia: a Pythia input file in CMND format is then needed; finally, a total cross-section normalizing the set of events is also expected; • a list of known detector models corresponding to the detectors that the user wants to simulate, together with the considered integrated luminosities (default values are associated to each detector but can be overwritten).
The event generation is entirely left to the care of the user, as such a feature strongly depends on the model.Dedicated tools (MadGraph, Pythia, Herwig [42,43]) are available on the market.Event generation through an internal call to Pythia is nevertheless possible via the CMND input.Initial and final-state QCD radiation (showering), as well as hadronization, can also be run in-code via Pythia.Yet, such functionalities are not those targeted by our tool design and we will not comment further about them here.We should emphasize that in order to obtain reliable predicted results, it is imperative to simulate sufficiently many events such that not too few simulated LLPs travel inside the considered detector's direction (up to the azimuthal angle in some cases).
Given the input described above, DDC computes the decay probability within the volume of the tested detector models of all tracked particles (i.e.LLPs identified by their PDG code) appearing in the simulated events.The cross-section normalization, 'visible' branching ratios, and integrated luminosities then allow to derive a number of expected events.Detector efficiencies are not taken into account in the current version -the corresponding information is generally not available for far-detector experiments -even though an average efficiency can be included within the 'visible' branching ratio input and within the weighing factors associated to (local) detector volumes (see the description in section 2.3).A full MC simulation of the detector response to LLP decays within its volume (from which local detection efficiencies would be extracted) is a priori beyond the scope of DDC at present.In addition, the program allows for the implementation of event cuts (e.g.restricting the energy of the LLPs) for each detector model.

Getting started
The source code for DDC can be downloaded from https://github.com/wzeren/Displaced-Decay-Counter .Its operation requires installed versions of Pythia8 [40] and HepMC [39] 3 as prerequisites: we refer to the corresponding manuals for installation.A C++17 compiler is also needed.After unpacking DDC, the paths in DDC/Makefile should be adjusted to the local setup.The make command then produces the central executable main in DDC/bin.

Operation and output
After compilation and once a valid input is available, the DDC/bin/main executable can be run in command line as (we use the names of template input files as example) ./maininputEvents.datinputLLPs.datresults.txtResults are then saved in DDC/bin/results.txt,and some logs on the availability of the detectors and a summary of the input information are printed in the directory DDC/bin/Logs/, providing the simulated acceptance and observed number of events for each detector and each LLP (the latter might be distinguishable through their decay modes).

Code structure
The code content of DDC is organized around three main tasks: 1. read and interpret the input files; 2. encode the properties of detector models; 3. analyze the MC events against the former.
The LLP and event information is processed from the main routine via the class inputInterface.A class object is constructed from the paths to the input files.The input is then extracted and stored using the json interface: the characteristics of the LLPs are in particular stored in objects of the sub-class CLLP.
The characteristics of the detector models are stored in objects of the Detector class.The latter and its sub-classes are described in detail in the following subsection.The central property of these detector models rests with their class function computing a probability of decay within the detector based on geometric considerations, once an LLP trajectory (represented by its polar angle) and a boosted decay length have been identified.
The analysis of MC events is processed within the analysis class.Such an object is created in the main routine from the stored input, giving access to various class subroutines, which allow testing the events against the detector models.The event file is read (or produced with Pythia).For each event, LLPs are searched for by their PDG code within the event chain.The LLPs are checked against a list of detector models, which is generated from comparison of the requested input with the stored detector models.The outcome of this analysis is then copied in the output file after application of the relevant normalization factors.

Modeling the decay probability
The differential probability distribution of decay for an unstable particle with lifetime τ is described in its rest-frame by an exponential distribution as a function of the proper time t * of the particle: For a free particle flying in the laboratory frame with velocity ⃗ v away from its production point at the origin O, we may translate this decay law in terms of the distance ℓ flown from the origin: where γ ≡ (1 − ||⃗ v|| 2 ) −1/2 is the relativistic factor.Consequently, the probability that the particle decays in a detector that comes across its trajectory between distances ℓ 1 and ℓ 2 (ℓ 1 < ℓ 2 ) simply reads: Observing this simple fact, the estimation of the detector acceptance reduces to the simple geometrical exercise of determining the distances at which the LLP trajectory intersects the outer layers of the detector.Additivity forms a remarkable property of the probability law of Eq. ( 2): given a volume V disjointly covered by sub-volumes V i , V = ⊔ i V i , the probability that the particle decays in V amounts to the sum of probabilities that it decays in the sub-volumes V i .We will constantly exploit this property by covering the detector volume by a collection of sub-volumes.Finally, we observe that Eq. ( 2) only involves the outer surfaces of the detector volume.Each surface intersected by the particle trajectory contributes an exponential term depending only on the distance of the intersection from the interaction point; the sign of this contribution is determined by the orientation of the surface with respect to that of the trajectory.
We stress three underlying hypotheses in this picture of the interplay between LLPs and the detector: 1. the LLP is produced at the interaction point, so that displaced vertices (resulting from the decays of primary LLPs into secondary LLPs) are not taken into account; 2. the LLP lifetime is not modified by the interactions of this particle with the matter constituting the detector; 3. the LLP flies in a straight line from the interaction point, so that the impact of the electromagnetic fields, in particular surrounding the interaction point, on charged LLPs is neglected.

Defining detectors in cylindrical geometry
In order to minimize the time cost of event simulation, we design the (virtual) detectors in a fashion exploiting the cylindrical symmetry (of LLP production) around the beam axis.Detectors (which need not observe the cylindrical symmetry in their design) are thus constructed from 'shapes' in the cylindrical half-plane -the latter being defined by the beam axis (Oz) and the radial coordinate (Oh) -, to which an angular aperture δφ is associated.In this way, LLPs emitted at any azimuthal angle are used in the calculation of the acceptance.Should control over the cylindrical angle prove necessary in a later stage of development of DDC, it can be restored by associating angular cuts to the detector volumes.For a detector covering an azimuthal aperture of ∆φ, this type of implementation saves a factor ≈ 2π/∆φ in terms of event generation, as compared to a strict three-dimensional modelization.
A set of basic objects is defined in the program in order to facilitate this implementation of detectors in cylindrical geometry: Oriented cylindrical segments (class CylSeg) define our most basic geometrical structure.It is defined by two coordinate sets in the cylindrical plane (std::array<double,2> AA={zA,hA};) representing its extremities, together with a sign (stored as an integer), representing its orientation (incoming or outgoing).A decay probability function is attached to this class, taking as input the polar angle θ of the particle trajectory and the effective decay length γ∥⃗ v∥τ , and returning an exponential term weighed by the orientation ε O = ±1 of the segment if the trajectory intersects the segment at a distance ℓ from the interaction point (returning 0 otherwise): ε O exp − ℓ γ∥⃗ v∥τ .Given that the segment is 'straight' in the cylindrical plane, it means that we renounce implementing curved surfaces in these directions, and that the latter need to be approximated by polygonal shapes.

Cylindrical detector layers (class CylDetLayer)
correspond to a list of oriented cylindrical segments (std::vector<CylSeg> SegList) with a weighing factor W (double).The corresponding decay probability function returns the sum of the decay probability terms associated with the individual segments, weighed by W: W The list of segments is meant to delimit a shape in the cylindrical plane, while the weight W primarily represents the azimuthal aperture δφ associated with the detector layer W φ = δφ/2π.In a more sophisticated description of detectors, a sensitivity contribution of the detector layer to specific decay products of the decaying particle could be incorporated within this weight, along the line: represents the branching fraction of the LLP into the final state X .However, in the current implementation, we only consider the detector acceptance, so that W S reduces to BR [visible] and can be extracted as a global factor.Cuts (discussed below) offer another handle on the detector efficiency.
The detector layer can be directly constructed from a list of preliminarily defined segments and a weight: CylDetLayer::CylDetLayer(std::vector<CylSeg> Seglist, double wgh) .As such a definition could however accommodate comparatively exotic constructions that would not amount to a realistic detector layer, we also provide a constructor taking a list of coordinates and the weight as input: CylDetLayer::CylDetLayer(std::vector<std:: array<double,2>> ptlist, double wgh) .The constructor then builds the convex polygonal envelope of the considered list of points (ptlist), and automatically determines the orientation of each derived segment.
The cylindrical detector layers are meant to be the building blocks of our detector modelization.They should be used to fill the volume of a detector.A simple (though rarely optimal) possibility would consist in defining detector layers of elementary dimensions δz, δh, δφ and fill the detector volume with such elementary 'bricks'.The function CylDetLayer CylBrick(std::array<double,2> coord, double length, double height, double apphi, double wgh) can be used in such an approach: it defines a detector layer of rectangular shape in the cylindrical plane, of dimensions length×height, centered on one set of coordinates coord.Other elementary shapes fulfilling this purpose are left to the user's imagination.

Detectors (class Detector)
combine an identifier (std::string), an integrated luminosity (double; this is a default value which can be overwritten in the input), and a list of cylindrical detector layers (std::vector<CylDetLayer>).The associated function Detector::DetAcc(double th,double leff) simply sums the weighed decay probabilities of individual layers for a particle emitted at a polar angle θ =th and characterized by the decay length leff= γ∥⃗ v∥τ .A cut function is also associated with each detector.This object represents the higher level of description of the detector in cylindrical coordinates.Contrary to the case of smaller objects, we do not provide automatic building functions (from coordinates) for detectors.The user may refer to the built-in detectors that we provide for templates.

Built-in detectors
A list of virtual detectors is provided by default within our tool.These objects are meant to model existing designs of far detectors at the LHC 4 .The mismatch in volumes between our implementation and the actual project (i.e. the sum of all subvolumes belonging to only the implementation, or only to the physical detector), normalized to the total volume of the detector offers an estimate of the accuracy of the description.In practice, compensation between fictitious and missing detector sub-volumes leads to a better precision at the level of the acceptance.Nevertheless, the provided detector models in cylindrical coordinates are often perfectible.

MATHUSLA
[45] is a box-shaped far-detector project of massive 100 m×100 m×25 m dimensions for the CMS interaction point.Its depth is parallel with the beam axis, with near and far collinear sides at respectively z min = 68 m and z max = 168 m.The base and cover of the box stand at respectively y min = 60 m and y max = 85 m above the beam.The left-and right-hand side boundaries are slightly offset: x min = −42.41m and x max = 57.59m.A nominal LHC integrated luminosity of 3000 fb −1 is associated per default to this detector.We offer two models for this project: • MATHUSLA0 is a basic but popular approximation consisting of a single detector layer delimited by the coordinates z min = 68 m, z max = 168 m, h min = 60 m, h max = 85 m, and the angular aperture ∆φ ≈ 1.39 rad.The mismatch in volume with the actual design is sizable, leading to an uncertainty of order 30% at the level of the acceptance.• MATHUSLA1 is constructed out of 21 cylindrical layers, aiming at a mismatch in volume below 4% as compared to the actual design, and a precision on the acceptance at the percent level.
The projections of these models in the polar plane (orthogonal to the beam axis) are shown in the upper plot of Fig. 1.

FASER and FASER2
[21] are projected cylindrical detectors aligned with the beam axis in the far-forward region of the ATLAS interaction point (480 m).The dimensions of FASER (FASER2) are given by a length of 1.5 m (5 m) and a radius of 0.1 m (1 m).The assumed integrated luminosities are 150 fb −1 and 3000 fb −1 , respectively.Given that this geometry respects the cylindrical symmetry, there is no difficulty in modeling these detectors with a single detector layer.The corresponding virtual detectors that we provide are straightforwardly named FASER and FASER2.

CODEX-b
[46] is a projected cubic detector of (10 m) 3  of the compact dimensions of CODEX-b, this approximation works rather well: the mismatch in volume with the actual layout is of order 10%, and the precision of the description should amount to a few percent.• CODEXB1 is constructed out of ten detector layers and aims at a geometrical accuracy better than 2%.
The projections of these models in the polar plane (orthogonal to the beam axis) are shown in the lower plot of Fig. 1.

AL3X
[47] is a proposed cylindrical detector with inner radius 0.85 m, outer radius 5 m and a length of 12 m, extending from z min = 5.25 down the beam axis of the interaction point in the ALICE cavern.Again, the geometry does not raise difficulties and the model AL3X consists of a single detector layer.The default integrated luminosity is set to 250 fb −1 .
ANUBIS [48] is a proposal for a detector of cylindrical shape, with diameter 18 m, length 56 m, placed in one of the service shafts about 24 m above the ATLAS beam line 5 , which is parallel with the base of the cylinder.The central axis is at z C = 14 m.The targeted integrated luminosity is 3000 fb −1 .

MAPP1 and MAPP2
[51] are projected sub-detectors of the MoEDAL experiment at LHCb.Both have trapezoidal shapes with no side aligned with the beam axis.
The planned integrated luminosities are respectively 30 fb −1 and 300 fb −1 .As the cylindrical geometry is not optimal to describe such volumes, we simply employ the 'brick' strategy, filling the volumes with cylindrical objects of elementary size, with a targeted accuracy of about 4%.The resulting models are simply called MAPP1 and MAPP2.
FACET [52] is a cylindrical detector project placed at 101 m of the CMS interaction point.Its radius is 0.5 m and its length 18 m.This straightforward geometry is implemented as a single detector layer in the FACET model.It covers the polar angle between 1 and 4 mrad.The expected integrated luminosity is 3000 fb −1 .

Implementing new detectors
The user can implement new detectors, either in order to investigate the potential of new geometries or simply to improve on the features of the built-in detectors.To this end, an assistant can be invoked.The corresponding executable needs first be created by running make DetEditor in the main folder.Then, DetEditor is available in the DDC/src folder, and can be called from the command line as ./DetEditor.Its function is to prepare a skeleton code and link it to the rest of the program.The only needed input for the assistant is the name of the new detector.
The user can then access the barebone script in DDC/src/Detectors and edit the properties of the considered detector.See sec.2.3.2 for details on a detector's geometrical definition in the code.

Implementing cuts
In order to better model the detector response, it is possible to implement cuts on events containing LLPs.A cut function is attached to each detector in the DDC/src/Detectors folder.It is by default trivial but can be edited by the user.The input of this function is a collider event in HepMC format, leaving a maximal flexibility on the possible types of conditions (possibly beyond mimicking the detector response).Nevertheless, we note that, while the implementation of cuts is critical in the description of detectors located close to the main interaction points, it remains largely incidental in modelling the response of far detectors at the current stage of their development. 6 Physics benchmarks In this section, we illustrate the functionality of the tool DDC via three simple LLP models that we test against the detector projects described in the previous section.As we consider only detector models planned for the LHC, events are produced in proton-proton collisions with 13 or 14 TeV center-of-mass energy.Nevertheless, we stress that other types of colliders and detectors can be encoded for simulation within DDC in a straightforward manner.We use Pythia8.245 for the numerical studies in this section.The default settings are used except that we turn off multi-parton interactions in order to speed up the simulation process without affecting the results.

Light neutral fermion produced in bottom decays
The first scenario concerns a light neutral exotic fermion (denoted as N here) with mass under 5 GeV.Such a new-physics state can indeed escape collider constraints as long as it is essentially singlet under the electroweak gauge group (so that it is not produced in over-abundant proportions in Z-boson decays and is not accompanied by a comparatively light charged partner).Astrophysical and cosmological limits do not apply as long as the considered particle is not stable (or very long-lived).The decays of the light fermion often proceed through lepton-and/or baryonnumber violating processes.Examples of such a light decaying fermion include a massive sterile neutrino (see e.g.Ref. [53] for a study of long-lived sterile neutrinos produced in meson decays at the LHC) or a light bino in R-parity violating SUSY (RPV-SUSY) models (existing studies include for instance Refs.[23,24,54]).The production of such a particle at colliders may proceed through numerous channels, lepton-and baryon-number conserving or violating.
Here, we focus on a lepton-number violating, baryon-number conserving production in bottom hadron decays (under the understanding that N carries no charge under lepton number).For simplicity, we focus on production modes mediated by partonic operators of the form (N Γ l e)(bΓ q u) or (N Γ l ν e )(bΓ q d), involving leptons (e, ν) and quarks (u, d) of the first generation and where Γ l and Γ q represent (model-dependent) Dirac structures.Bottom hadrons H b may then decay into N and a lepton, with possible additional light hadrons H q .The most relevant H b 's for N -production are the comparatively long-lived B-mesons and Λ 0 b baryons (other H b 's have a much larger SMlike width, making an exotic disintegration likely uncompetitive) and conservation of the angular momentum forbids the two-body Λ 0 b → N ν.Thus, under the further hypothesis of dominant twobody decays, N -production is largely captured by the disintegration of B-mesons B + → N e + , B 0 → N ν, likely comparable in magnitude as long as the new-physics controlling the short-distance effect is SU (2) L -conserving (here, we specialize in the case where u and d are left-handed fields).We assume that the decay of the exotic fermion is not correlated with the production mode, so that the lifetime can be set independently from the production cross-section (this is typically the case in a R-parity-violating SUSY model, for instance, where different couplings mediate these different processes).Possible N decays include baryonic, semi-leptonic, and (radiative) leptonic channels [55], and need not be specified here, as long as we ignore the sensitivity of the detectors to specific decay modes.Finally, we specialize in the case of a mass m N = 1 GeV for the exotic fermion, and consider various values of the proper lifetime cτ N , with an assumed visible branching ratio of 100 %.
The generation of the parton-level pp → b b events and the ensuing hadronization processes at the center-of-mass energy √ s = 14 TeV are performed with Pythia8 (here used internally within DDC).The B ± meson production cross sections over the whole 4π solid angle at the 14-TeV LHC is estimated to be about 4.87 × 10 11 fb [53].The branching ratio of the exotic B-meson decay (into e ± N or ν e N ) is set to the fictitious value of 1 in the event generation (for an efficient use of computer resources) and the production cross-section must be a posteriori re-scaled to account for the actual branching ratio.The events are then analysed with DDC for multiple values of cτ N7 .We thus obtain the predicted numbers of signal events for each of the considered detectors at chosen integrated luminosities.Here as well as in the following subsections, we model the MATHUSLA, CODEX-b, and ANUBIS experiments with the implementations MATHUSLA1, CODEXB1, and ANUBIS1 described above.The results are presented in the upper plot of Fig. 2: the boundaries of the region with more than 3 signal events (95% confidence level in the absence of background) are plotted in the plane spanned by Br(B + /B 0 → e + /ν e N ) and cτ N for a fixed N mass of 1 GeV, considering the listed integrated luminosities.The respective relevance of each experiment in terms of parameter-space coverage in this benchmark can be appreciated directly by reading the plot.In the lower plot of Fig. 2 we analyze in which proportions the predictions depend on the modelization of the detectors.We define the variable ∆N S = |N 1 S − N 0 S |/N 0 S , where N (0/1) S denotes the corresponding signal-event numbers, in order to quantify the relative difference between the simplistic and more sophisticated modelizations, namely ANUBIS0 and ANUBIS1, CODEXB0 and CODEXB1, and MATHUSLA0 and MATHUSLA1.The wiggles in the small cτ N limit are due to insufficient simulation statistics in this regime.Expectedly, the compact dimension of CODEXB ensure a good percent-level performance of CODEXB0, while the discrepancies are larger for the more massive MATHUSLA (∼ 10%) and ANUBIS (∼ 30%).Still, we observe that the order of magnitude is properly captured by the simple approximations, which legitimates their common use in the literature.
Models of long-lived heavy neutral leptons (HNLs) or light RPV neutralinos produced in B-decays have been previously studied in e.g.Refs.[22][23][24][25]54] from the perspective of far-detector proposals.The results of DDC were compared to those obtained in these papers and showed good agreement, up to traceable discrepancies in the choice of B-meson parameters or the inclusion of three-body decays, for instance.In the case of CODEX-b, ANUBIS, and MATHUSLA, simple approximations were routinely employed, which we thus tested against the CODEXB0, ANUBIS0, and MATHUSLA0 designs.The compatibility of our results with those of earlier works argues in favor of the validity of the implemented detector models within DDC.
Below, we focus on the specific case of HNLs in order to compare our results with those obtained with FastSim.

Heavy neutral lepton in the minimal scenario
In the minimal scenario for Majorana neutrinos, the SM is extended with three HNLs which mix with the active neutrinos.The interaction Lagrangian for the HNLs and the SM electroweak gauge bosons is: where i, j = 1, 2, 3, and ℓ α , with α = e, µ, τ , denoting the charged leptons of the SM.U αj is the mixing-matrix element relating the (dominantly singlet) HNL N j to its subdominant SU (2) L component ν α .V L labels the PMNS matrix.For phenomenological purposes, a "3+1" scenario is commonly assumed where the sterile neutrino mixes with only one generation of the active neutrinos.Here, we restrict ourselves to the case of a single HNL N that mixes with the electron neutrino only, and consider its production in the rare decays of both charm mesons (D 0 , D ± , and D ± s ) and bottom mesons (B 0 , B ± , and B 0 s ).We follow Ref. [56] for the evaluation of the total number of such mesons produced in LHC collisions, as well as that of their decay rates into HNL's and the subsequent HNL decay rates into SM particles.The mixing angle plays a central role in enabling the latter two.All the two-and threebody semi-leptonic and leptonic decay channels of the HNL are considered (and computed at leading order in chiral perturbation theory) for the calculation of its total decay width, and among all the decay channels of the HNL, only the tri-neutrino ones are counted as the invisible final states.
In Fig. 3, we present the expected sensitivity reach of far-detector experiments in the (m N , |U e4 | 2 ) plane, corresponding to the boundary of the region where 3 signal events can be detected within the exposure time of the experiment; these are identified with the exclusion bounds at 95% confidence level in the expected absence of background events.The upper plot contains the results obtained with DDC for MATH-USLA (with the MATHUSLA1 implementation), FASER, FASER2, CODEX-b (with the CODEXB1 implementation), ANUBIS (with the ANUBIS1 implementation), MAPP1, and MAPP2.Furthermore, in the background, we indicate the current upper bounds extracted from Refs.[29,[57][58][59][60][61][62] (the light gray area), the lower bounds computed from considerations relative to the Big Bang Synthesis (BBN) [63] (the gray dashed curve), as well as the target region where a type-I seesaw is compatible with the active neutrino masses between 0.05 eV and 0.12 eV [64,65].In the lower plot, we compare the sensitivity reach of MATHUSLA USLA with corresponding results obtained in Ref. [37] with FastSim.The gray area in the background represents the regions currently excluded experimentally [29,[57][58][59][60][61][62].The lower bounds from Big Bang Nucleosynthesis (BBN) [63] is shown as a gray dashed curve, and the type-I seesaw target region for mν e between 0.05 eV and 0.12 eV [64,65] is displayed in the brown band in the lower left part of the plots.
computed with DDC with its counterpart employing FastSim [37].The latter is read from Fig. 6 of Ref. [37] and corresponds to the 4-signal-event curve.We observe that both tools lead to largely compatible sensitivity limits.The higher reach achieved by FastSim at large HNL mass can be essentially interpreted in terms of the more diverse production channels considered in Ref. [37], in particular through the B ± c mesons (in addition to W and Z-bosons as well as the τ leptons channels).In the mass-range covered by both DDC and FastSim, the slightly weaker sensitivity predicted by FastSim can be attributed to a number of effects including e.g. the precise number of mesons produced in the considered regime, as well as the finer treatment of the final-state particles and the LLP DV reconstruction within FastSim (see the discussion in Sec. 1).A detailed comparison goes beyond our purpose here, the general qualitative agreement being obvious.

Long-lived complex scalar mixed
with the SM Higgs boson and produced from B-meson decays In this example, we consider a complex neutral (electroweak singlet) scalar S, with mass m S , that mixes with the SM Higgs according to an angle θ.Such a complex scalar consequently acquires a coupling to SM particles through θ and may then be produced in e.g.rare decays of B-mesons such as B → K + S, which we consider here.Decays of the light scalar into SM final states are similarly expected.For the calculation of the production and decay rates of the light scalar S, we follow the module provided in the FORE-SEE GitHub repository [66].For simplicity, we shall assume that the trilinear coupling between the complex scalar S and two SM Higgs bosons h is vanishing.The following visible final states are considered in the decays of the light scalar: S → e − e + , µ − µ + , ππ, KK, gg, and ss.
Given the current bounds on θ applying in the relevant range of m S , such light exotic scalars of the studied type are expected to be long-lived.We perform the simulation and computation of the sensitivity prospects at LHC-based far-detector experiments using DDC.The corresponding results are presented in Fig. 4 in the (m S , θ) plane.Similarly to Fig. 3, we show the sensitivity reach of all the considered far detectors obtained with DDC in the upper plot, while the lower plot is dedicated to a comparison of our results with those obtained with FORESEE (3 signal events) and FastSim (4 signal events) in the case of FASER2 and MATH-USLA, respectively.The MATHUSLA sensitivity curve at 4 signal events corresponding to FastSim is borrowed from Fig. 5 of Ref. [37].The gray area corresponds to the present bounds in the modelparameter space (see Ref. [3] and the references therein).Again, we observe a qualitatively good agreement among the sensitivity results derived with DDC, FastSim, and FORESEE, which demonstrates the efficiency of the approach based on detector acceptance that we use in our tool.The minor differences on the exact position of the boundaries for the sensitivity regions may again be attributed to the detailed features implemented in each tool.In particular, FORESEE's requirement on the energy of the LLP, E LLP > 100 GeV (to ensure a negligible background [34,67,68]), results in a slightly more conservative coverage of the parameter space.As explained before, similar features can be implemented within DDC using e.g. its cut function, but we regard such a degree of refinement largely unnecessary at the current stage.

Long-lived scalar in an extended
Higgs sector As a third example, we consider an extended Higgs sector, containing new doublet-dominated states H, A, H ± at about 0.5 TeV and a lighter, mostly singlet CP-odd scalar A S , with mass in the 100 GeV range.In the limit where this pseudoscalar is purely singlet, its couplings to SM particles are suppressed, leading to a long lifetime, as well as a vanishing production rate in direct channels.However, some couplings of A S to BSM electroweakly-charged particles may remain sizable.For instance, the trilinear couplings [70] and mediate the heavy-Higgs decays A → A S h SM , H → A S Z, and H ± → A S W ± .These channels may then dominate the total width of A, H and H ± , which are easily produced in gluon fusion (ggf ) at hadron colliders, or in association with top and bottom quarks (but may be difficult to observe if decaying in the indicated fashion).These channels determine the production channels of the singlet.We refer the reader to Ref. [69] and references therein for a discussion of such a configuration in the context of the Next-to-Minimal SUSY SM (NMSSM), or more generally, in models with two Higgs doublets and a singlet.The lifetime of A S can be almost arbitrarily large in such a setup.In the NMSSM, the width may be dominated by the diphoton channel, mediated by chargino loops.When this mechanism is absent, or in the non-SUSY case, the A S decays, typically into SM fermion pairs, proceed through the vanishingly small doublet component of this state, or exploit larger Higgs-to-Higgs couplings in loop-induced channels.
Here, we assume that A S is long-lived and focus on the production mode gg → A → A S h SM .In practice, alternative production channels may be competitive, in particular gg → H → A S Z,

σ(pp
Fig. 5 Sensitivity plots for the extended Higgs models in the plane defined by the doublet production crosssection pp → A → A S h SM and the singlet proper lifetime cτ A S .The Higgs masses are chosen as (M A , M A S ) = (410, 70) GeV in the upper plot and (500, 200) GeV in the lower plot.The corresponding cross-section values in Ref. [69] are highlighted with a horizontal black dashed line.FASER has too weak a sensitivity to this scenario to appear in these plots.
or associated production with third-generation quarks, but we restrict ourselves to the single pseudoscalar channel as a simple handle on this scenario.Realistic cross-section values can be inferred from Table 2 of Ref. [69], considering the pp → A → (h SM → τ + τ − )(A S → γγ) search channel, after re-scaling the corresponding numbers by the branching ratios of the final Higgs states, Br(h SM → τ + τ − ) and Br(A S → γγ).Our scenario indeed differs from that of Ref. [69] in that A S no longer promptly decays (and we do not pay attention to the decays of h SM ).In practice, we will use the two benchmark points (m A , m A S ) = (410, 70) GeV and (500, 200) GeV (listed among other values in Table 2 of Ref. [69]), vary the pp → A → A S h SM cross section freely and assume a 100% visibility of the decay products of A S .
Again, we invoke Pythia8 for event generation at √ s = 13 TeV (in accordance with Ref. [69]) and let DDC compute the expected signal-event numbers at the various LHC far detectors.The results are displayed in Fig. 5. Again, the boundaries of the 3-signal events are shown in the plane spanned by σ(pp → A → A S h SM ) and cτ A S , highlighting the discovery potential of the far detectors for long-lived scalars produced by a single scalar resonance.The black horizontal dashed lines mark the cross-section values corresponding to Table 2 of Ref. [69], which may be deemed realistic for a Higgs-inspired model.We stress that the results of this benchmark can be easily generalized to other models where a single heavy resonance is produced via gluon fusion and decays to an LLP.

Long-lived fermion from a pair of heavy resonances
In the final example, we consider a RPV-SUSY-inspired scenario where the production of SUSY resonances is still dominated by R-parityconserving effects, while small R-parity-violating couplings mediate the decays of a long-lived LSP (lightest SUSY particle) bino χ0 1 .The latter may be arbitrarily long-lived, depending on the magnitude of the RPV couplings.The binos can be directly produced in pairs via s-channel contributions involving their mixings with higgsinos, or via t-channel diagrams involving squarks.However, if higgsinos and squarks are very heavy, which we assume for simplicity, the direct production crosssection is correspondingly suppressed.Instead, we consider an indirect production mechanism where right-handed sleptons, ẽR , exist in the ∼ 0.5 TeV range, are produced in pairs in pp-collisions (via their electroweak interactions) and decay into a lepton and a χ0 1 (via R-parity conserving gaugino couplings).We further assume a 100% visibility for the bino decays and the mass of this particle is set to 1 GeV.We do not pay particular attention to limits from prompt searches at the LHC here, as the scenario is simply meant for illustration, but mono-and multi-lepton searches can be relevant in this context, as the long-lived neutralino would leave a clear missing energy signature in the ATLAS or CMS detectors.
Once again, we generate events with Pythia8 for the slepton pair production at √ s = 13 TeV in pp collisions, with subsequent decays into binos and leptons.We normalize the cross-section to the next-to-leading and next-to-leading-logarithmic augmented order, as provided in Ref. [71].The sensitivity reach of the far-detectors is computed with DDC and presented in Fig. 6.It shows the limited detection capabilities of far detectors in such a scenario, as the 3-signal event threshold is barely reached in the cτ χ0 1 ∼ 10 −2 − 1 m range.This benchmark exemplifies the case of LLP production from a pair of heavy resonances, and we stress that it can be easily generalized to scenarios beyond the RPV-SUSY model that we considered here.
These simple scenarios have provided us with the opportunity of demonstrating straightforward functionalities of DDC.The input files for all these sample benchmark studies, as well as the plotting scripts, the final sensitivity plots, and the data points, can be found in the folder DDC/examples.

Conclusions
In this paper, we have described Displaced Decay Counter (DDC), a versatile C++ program for the calculation of detector acceptances, relying on Pythia8 and HepMC.DDC is primarily meant as a simulation tool for long-lived-particle models and the currently booming far-detector proposals.It accepts MC events as input in LHEF and HepMC format (a Pythia8 run card can also be used for internal event generation).The current version of DDC comes with a collection of far-detector designs, approximating several detectors currently planned or in construction at the LHC, and provides the user with an editor for the implementation of further far-detector models.
We have presented several benchmark scenarios to illustrate the usage of our tool.The first one, the already well-studied case of light long-lived neutral particles produced in B-meson decays, gave us the opportunity to confront the performance of DDC to older results and validate the implemented detector models.We then focused on the specific case of HNL's produced in the decays of charm and bottom mesons and noted the general agreement of our predictions with those derived for MATHUSLA with the tool FastSim.Another scenario, with a long-lived scalar singlet mixing with the SM Higgs, was considered, mainly for comparison purposes with FORESEE (in the case of the FASER2 experiment) and FastSim (in the case of MATHUSLA).Again, we stressed the close agreement of the sensitivity prospects derived with DDC, with those obtained with these other tools, despite the currently less refined treatment of final states in our implemented detector models (strictly focusing on geometrical acceptance).The further two models exemplify the cases of LLP production via a single or a pair of resonances in pp collisions.On a concluding note, we stress that DDC is highly versatile, offering modelling capabilities beyond those demonstrated here, and thus comes with a high potential for studying LLP models at colliders.

Fig. 1
Fig.1Comparison between the actual detector layouts (in red), the simple designs (in blue) and the more exact coverage (in green) that we provide in the polar plane.Upper: MATHUSLA (red) vs. MATHUSLA0 (blue) and MATHUSLA1 (green); Lower: CODEX-b (red) vs. CODEXB0 (blue) and CODEXB1 (green).The axes x and y (in meters) span the polar plane (i.e. are orthogonal to the collision axis (Oz)).

Fig. 2
Fig.2Upper : Ideal exclusion bounds of far detectors at 95% confidence level for a long-lived neutral fermion produced in B + and B 0 meson decays at the LHC far detectors.The 3-signal-event boundaries are shown in the plane Br(B + /B 0 → e + /νeN ) vs. cτ N , when the neutral fermion mass is fixed at 1 GeV.(They represent upper limits on the branching ratios.)Lower: Comparison between the implementations of simple approximations and more sophisticated implementations, in terms of the defined variable ∆N S , for ANUBIS0 and ANUBIS1, CODEXB0 and CODEXB1, and MATHUSLA0 and MATHUSLA1.

Fig. 3
Fig.3Upper : Sensitivity reach of all the far detectors considered in this work at 95% confidence level derived with DDC.The results are displayed in the parameter space of the HNL model spanned by the mass m N and mixingmatrix element squared |U e4 | 2 .Lower : Comparison of the sensitivity reach achieved by DDC in the case of MATH-USLA with corresponding results obtained in Ref.[37] with FastSim.The gray area in the background represents the regions currently excluded experimentally[29,[57][58][59][60][61][62].The lower bounds from Big Bang Nucleosynthesis (BBN)[63] is shown as a gray dashed curve, and the type-I seesaw target region for mν e between 0.05 eV and 0.12 eV[64,65] is displayed in the brown band in the lower left part of the plots.

Fig. 4
Fig.4The projected sensitivity of the far-detected experiments are shown in the parameter space of the complex Higgs scalar model with mixing.The upper plot summarizes the results of DDC, while the lower plot performs a comparison with FastSim and FORESEE in the case of MATHUSLA and FASER2, respectively.The background gray area represents the current experimentally excluded region (see Ref.[3] and the references therein).

Fig. 6 1 = 1
Fig. 6 Number of events for a light bino m χ0 1 = 1 GeV, produced in the decays of selectrons pair produced, with m ẽL = 500 GeV.The number of events are shown as a function of the lifetime of the bino.The black horizontal dashed line indicates the 3 signal event threshold.FASER and FASER2 have too weak sensitivity to this scenario to be displayed on this plot.