Event generation for beam dump experiments

A wealth of new physics models which are motivated by questions such as the nature of dark matter, the origin of the neutrino masses and the baryon asymmetry in the universe, predict the existence of hidden sectors featuring new particles. Among the possibilities are heavy neutral leptons, vectors and scalars, that feebly interact with the Standard Model (SM) sector and are typically light and long lived. Such new states could be produced in high-intensity facilities, the so-called beam dump experiments, either directly in the hard interaction or as a decay product of heavier mesons. They could then decay back to the SM or to hidden sector particles, giving rise to peculiar decay or interaction signatures in a far-placed detector. Simulating such kind of events presents a challenge, as not only short-distance new physics (hard production, hadron decays, and interaction with the detector) and usual SM phenomena need to be described but also the travel has to be accounted for as determined by the geometry of the detector. In this work, we describe a new plugin to the {\sc MadGraph5\_aMC@NLO} platform, which allows the complete simulation of new physics processes relevant for beam dump experiments, including the various mechanisms for the production of hidden particles, namely their decays or scattering off SM particles, as well as their far detection, keeping into account spatial correlations and the geometry of the experiment.


Introduction
Experiments at the LHC have yet not reported any sign of physics Beyond the Standard Model (BSM). Nevertheless, the problem of reconciling our description of the fundamental interactions and particles with long-standing problems, such as the matter-antimatter asymmetry in the Universe, the evidence for dark matter from many astrophysical and cosmological observations and the origin of the neutrino masses, becomes ever more pressing. Many ideas have been proposed, some of which addressing one problem at the time, others, more ambitious, providing solutions to two or more open questions at the same time. In this context, a recurrent theme is the hypothesis that a hidden sector involving new light particles, might be coupled to the Standard Model via, for instance, one of the three portals (scalar, fermion and vector) in a feeble way. Such scenarios can provide not only dark matter candidates, but also other states, such as heavy neutral leptons, vectors, scalars, which could be long-lived and also possibly decay back to SM particles.
To prove the existence and measure the properties of such elusive particles is extremely difficult. The situation is in fact similar to neutrino production and detection 1 : as the energy does not pose a hindrance, one is lead to consider high-intensity facilities and design experimental setups that maximise the rates. In short, one needs very intense beams and then let such beams cross a heavy and instrumented target to detect their scattering or, if necessary, to create a decay tunnel as long as possible to observe their decay products. A first example is NA62 [3,4] which can run in a beam dump mode and is expected to collect data in this configuration soon. The DUNE [5] experiment will be operational in 2026 to study neutrino oscillations. As a by-product it could also search for hidden sector particles. The SHiP experiment [6,7] has been designed on purpose to search for such light and feebly interacting particles originated in interactions of 400 GeV/c protons produced by the CERN SPS [8]. More recently, other proposals have been put forward to also exploit proton collisions at the LHC experiments with detectors placed not very far from the collision points, namely, the CODEX-b [9], MATHUSLA [10,11] and FASER [12,13] experiments.
In the present paper we address the issue of how to efficiently simulate the production of a flux of particles belonging to the hidden sector and their subsequent interactions and/or decays. In the following, we will generically call a Beam Dump Facility (BDF) every experiment where a known flux of primary SM probes strikes on a fixed target and a detector is placed in an optimal position with respect to the target, with the aim of detecting either neutrinos or a new kind of feebly interacting particles produced off the primary beam interaction. The case of a detector placed close to a collider experiment such as those proposed in [9][10][11][12][13] can be equally treated within our framework without any modifications. In practice, sensitivity studies of such experiments to new physics phenomena, rely on the simulation of two distinct processes, one where the new particles are produced and the other where the new particles (or their decay products) interact with a detector placed at some macroscopic distance, from tens of meters to thousands of kilometers. The production of a feebly interacting particle in a beam dump can proceed through at least three phenomenologically different phenomena: i) its prompt production in the high energy scattering of the primary beam particle with a nucleus in the target; ii) as the result of the decay of SM particles produced in the primary collision or in the cascade process in the target; iii) through the bremsstrahlung process of primary or secondary particles in the target. The detection, on the other hand, will proceed either through the decay in flight of new particle back to visible SM final states or directly through the scattering with the matter in the detector. The aim of this paper is to provide an implementation that allows the simulation of the complete chain of subprocesses, from the production to the final detection at a BDF in one go. Our starting point are FeynRules [14][15][16] for the implementation of the new physics model lagrangian and MadGraph5 aMC@NLO [17,18], MG5aMC for short, for providing the necessary short-distance physics elements, the automatic production of particle-level unweighted events and the framework. To achieve maximal flexibility we provide the implementation as a MG5aMC plugin, in line with other recently developed applications [19][20][21]. Figure 1 shows a sketch of the elements of the simulation which are automatically combined in our implementation. These functionalities are available so that samples of unweighted events in a standard format can be generated in a single step and eventually passed to the simulation of the detector response. For the rest of the paper we dub the MG5aMC plugin for the simulation of hidden particle effects at beam dump facilities with the short-hand name MadDump.
An important aspect of our implementation is that it provides the elements of the simulation that are related to BSM physics in a single framework. This entails a number of advantages. First, it eliminates the possibility of making mistakes in the generation or in the combination of event samples for the production and the detection stages. This is particularly relevant when scanning over the parameters of a BSM model, where, although every step is simple in principle, the combinatorics and the bookkeeping would make the whole construction cumbersome. Second, by using functionalities already present in MG5aMC it allows to fully automate the scanning over the BSM parameters. Third, once implemented in FeynRules and available in the UFO format, the same BSM model and parameter points can be constrained in different contexts within the same framework, for instance, in collider physics using MadAnalysis5 [22][23][24] recasting capabilities or in MadDM [21,25,26].
In MadDump the primary flux of Standard Model probes that can generate the hidden particles has to be provided by the user. As shown in Figure 1, it can be either the original beam hitting the target, or the flux of hadrons following the hard interaction interaction that can produce the hidden particle through their decays. In the former case, the user has just to  Figure 1: Schematic view of relevant processes that might happen at a BDF and can be computed/simulated within MG5aMC .
provide the specific particle code of the probe and its energy in the laboratory frame, while for the latter case the flux can be given as an event file featuring the decaying particle momenta and specifying the particle identifier. Note that, for the case of hidden particles generated in the target from meson decays, our approach is more flexible than just directly linking event generators like Pythia8 [27,28] or HERWIG7 [29] as it allows, in principle, to later include other effects, such as the cascade production of secondary particles which in some cases could be relevant. MadDump is able to handle event files in all formats of the most used event generators. Another important ingredient is the geometry of the BDF, which is provided by the user in a dedicated file. The paper is organised as follows: in section 2 we introduce the algorithms at the core of MadDump. In section 3 we present illustrative examples of possible applications, considering physics cases relevant for the SHiP experiment. Conclusions and the perspectives of the present work are given in section 4. In three appendices A, B and C we provide many details on the numerical techniques employed and the associated uncertainties, while appendix D documents the scripts that produce the results presented in section 3.

Approach
The first important aspect of our implementation is the idea of considering the beam dump experiment as a two-step process: • Production: Hidden particle flux generation upon interaction of the beam with the target; • Detection: Interaction of the hidden particles (or their decay products) in the (possibly far-placed) detector.
While both steps depend on the details of the new physics model and therefore they have to be considered together, it is possible to factorise the simulation into two independent steps: the results of the Production phase simulation are used to build a (two-dimensional) parametrisation of the incoming hidden particle flux hitting the detector and leading to different signatures in the detector. By disentangling the Production from the Detection phase and the corresponding event generation into two subsequent steps, the possibility of following the full history, from production to the final signature in the detector, of each event is lost. However, the gain in efficiency in the event generation is enormous, an element which is a key aspect in the simulation of a high-intensity experiment.
The second important aspect of MadDump is that it has been designed as a plugin of MG5aMC . In other words, it heavily relies on already existing modules which are at the core of MG5aMC , such as the phase space integration provided by MadEvent and the decay package MadSpin, integrating them with functionalities that are specifically required for BDF's, so the various steps of the simulation can be undertaken to obtain the final result in one go. Among the key new functionalities, we stress • the determination of doubly differential scatter data in the Production phase of nonstandard particle beams and their support in the Detection phase; • the support of HepMC as input format with the aim of making easier the interplay with other Monte Carlo generators like Pythia8 or HERWIG7.
The third aspect is the underlying idea of factorising SM physics from the BSM one, whenever possible. The former, while accessible via standard MC tools, is in general quite involved and needs the modeling of many effects. Being strongly dependent on the particular experimental setup, a dedicated simulation of the target and/or detector effects is almost always needed. However, while cumbersome, this part of the simulation does not have any dependence on the new physics model considered and can be taken care once for all. On the other hand, the new physics short-distance part by definition depends on the details of the model and therefore has to be generated/considered for each different data interpretation. Fortunately, it can be described quite easily from first principles and dealt with by usual or especially developed MG5aMC modules.

Production
In a typical beam dump experiment, a collimated and mono-energetic beam of protons or electrons impinges on a thick target, at rest in the laboratory frame. A copious number of SM particles is generated both in primary and subsequent secondary interactions inside the target, which is designed to maximise the particles yields. The production of hidden sector particles may proceed according to different mechanisms. In the following, we focus on two cases, i.e., • prompt production in primary or secondary beam interactions • rare meson decays.
Without loss of generality, we consider the case of proton beam dump experiment, keeping in mind that other situations can be dealt with by MadDump in a completely analogous way. Depending on the specific BSM physical case (model and parameter point of interest), the prompt production may be described in perturbation theory and it can be treated directly in MadDump. In this case, the main input is the BSM Lagrangian (in the UFO format) which fixes the hidden sector model and its interactions with the SM particles. For example, in a model where a new massive vector mediator couples to quarks and DM fermionic particles, the main production mechanism resembles Drell-Yan production and decay at fixed target experiments. This description must, however, be consistent with the typical scales characterising the model. For instance, in the above example the mass of the mediator must be larger then the QCD scale for the computation to be reliable. As stated above, the main goal of MadDump is to handle BSM interactions and then embed them consistently into a complete and modular simulation chain which can take into account the rest of the SM interactions, possibly also using other inputs. In particular, an accurate simulation of the cascade production of hadronic particles, mesons and baryons, is expected to be handled by other MC tools or dedicated simulations, which can fully include parton showers, hadronisation, nuclear effects, meson decays and so on. This part can be very important if hidden particles are produced in the decays of mesons. In that case, the meson production is assumed to be simulated independently. MadDump, on the other hand, by parsing the event files 2 describing the beam-target event takes care of the decay of mesons into hidden particles, employing for example, an effective field theory approach, which can also be implemented at the level of the UFO. In this way, mesons are considered on the same footing of the elementary particles in the model and their decays occur through interaction vertices that can be handled by MG5aMC in the usual way. Examples of both prompt production and meson decay studies are given in section 3.
Either way, by the hard-interaction or via the decay of mesons, hidden particles are created, which fly out of the target close to the forward direction. The hidden particles produced during the beam dump, however, do not form a standard collimated and mono-energetic particle beam. On the contrary, they have a spatial distribution, they are produced in different points inside the volume of the target, and a phase space, i.e., a non trivial four-momenta spread, distribution. Assuming that the hidden particles travel freely until they eventually enter in the active region of the detector, after macroscopic distances that can go from meters to hundreds or thousands of kilometers, we can describe the beam of hidden particles by means of a multi-differential flux function where the vector x denotes the collection of all the other relevant kinematical variables (angles, spatial distribution of the hidden particles production point within target, etc), but the energy. The flux function in general not known a-priori and/or in an analytical form since it represents the result of scattering/decay processes in the Production phase. This distribution is implicitly determined by the simulation of the Production phase and in turn it can be extracted from a sufficiently large sample of production events. In practice, since the flux depends on the particular BSM model and the specific production mechanism, it cannot be fit it once and for all (as for the proton pdf). An on-the-fly fitting procedure is needed that is fast, robust and flexible.

Detection
The final detection of the hidden particles might occur according to the two distinct physical processes • the hidden particle interacts with the active volume of the detector, resulting in a neutrino-like signature; • the hidden particle decays to SM particles inside a dedicated decay tunnel (included of what we dub "detector"), resulting in a "displaced vertex" signature.
The interaction of the hidden particle with the detector turns out to be the most difficult part to simulate. One can exploit some approximations and different Monte Carlo techniques to obtain a generator with a satisfactory level of accuracy. On the contrary, as we will discuss later, in the displaced decay case, the same complications do not arise and the situation is much easier to handle. Let us discuss first the interaction case.

Interactions of hidden particles in the detector
The outcoming flux of hidden particles from the Production phase, eq.(2.1), corresponds to the incoming hidden particles flux of the Detection phase. Our strategy is to parametrise the flux by using Production event samples and use it as a generalised partonic distribution function (pdf) for the needed computations in the Detection phase. In doing so, we will also able to parametrise not only the acceptance of the detector but also some of the efficiencies/features that are model dependent. The total interaction cross section with the fiducial volume of the detector is obtained by convoluting of the flux function φ(E, x) with the elementary cross sectionσ I for the "partonic" sub-process HP + X → HP + X , where X represents the SM matter particle in the detector and "HP" the hidden particle. Our implementation is able to handle: • elastic electron scattering, X = e − • deep inelastic scattering of nucleons (DIS), X = u, d, s, c.
The geometrical detector acceptance sets the integration limits in the convolution integral. This is equivalent to introduce a weight function W (E, x) which is 0 if the point does not pass the acceptance cut and 1 otherwise. In this way, it is possible to restore the integration limits to their full ranges. This simple idea is the basis of more sophisticated re-weighting strategies in Monte Carlo integration. We can exploit these techniques with the aim of modeling in a realistic way the detector efficiency. For example, due to its shape and composition, the particles entering the detector may travel a longer or shorter path inside its volume. Correspondingly, the probability that the particles interact inside the detector will be greater or lesser resulting in an efficiency function depending on kinematical variables of the incoming particles. We can effectively describe this effect giving a suitable weight W (E, x) to each incoming hidden particle particle penalising those which will travel a shorter path. To this aim, we have introduced in our framework the possibility of introducing a re-weighting procedure and we have extensively tested/used it to describe some common detector effects. The user can easily customise the weight function in order to refine the simulation at will. Our master formula for the Detection cross section σ D is given by: The crucial point here is that the partonic cross sectionσ D does not depend on the other variables but the energy, as it follows directly from the Lorentz invariance of the interaction among point-like particles. Due to this property we are allowed to formally perform the integral over the x vector before performing the convolution with the partonic cross section leading to the introduction of an effective one-dimensional energy pdf: Before proceeding further, we discuss some practical implications of the above formula. The 1D-functionφ(E) can be obtained through a simple 1D fit of the energy histogram of the input DM production events, after having re-weighted them by the weight function W (E, x). According to eq.(2.4), this is the only ingredient needed to compute the total cross section, which in turn is crucial to extract the hidden particles yields in the detector. Up to this point, the simulation of a collimated (along the beam axis) but not mono-energetic beam of hidden particles particle impinging on the detector is achieved. However, this is not sufficient to develop a full event generator of unweighted "interaction events". A complete event reconstruction that gives access to correlations (between energies, positions, angles), may be of great importance. It can allow to study possible kinematical cuts with the aim of maximising the signal yield with respect to the backgrounds, and to accurately model detector effects. For example, due to different energy-angular correlations of the DM particles wrt the neutrino ones, it is possible to have a signal-enriched sample of events for off-axis detector configurations, as pointed out in ref. [31]. Hence, in principle this information might be useful to design optimised experimental configurations. Moreover, the output events can be further post-processed exploiting for example parton shower programs or other dedicated tools in order to have a better estimates of the detector efficiency. In principle, given the "trivial" dependence of the partonic cross section on the parameters x, we can assign them a posteriori on a event-by-event base according to the distribution where E is the energy of the current event. In practice, the procedure underlying the above formula is limited by the computational issue of performing a robust and reliable multi-dimensional fit, since the incoming particle flux is not known analytically. As mentioned above, our fitting procedure relies on the point-like approximation of the target in the primary interaction. In a typical beam dump experiment, the distance between the production target and the near detector is greater than the characteristic size of the target, so that the point-like target approximation is a reasonable first approximation. Under this assumption, the complexity of the problem reduces considerably. As depicted in Figure 2, just three kinematical variables are needed to describe the incoming flux impinging the detector: the energy E, the polar angle θ and the azimuthal angle φ. Furthermore, the physics occurring at the production point is invariant under a rotation around the beam axis resulting in flat distributions for the azimuthal correlations. Hence, the only relevant correlations are the E − θ ones. In Figure 3 we show a typical plot of the production scatter data, which enter into the neutrino detector, in the E − θ plane for a SHiP-like configuration [6,7]. Figure 3: Input 2D-Scatter plot in the DM E − ϑ plan. As example, we consider DM particles produced by a Dark Photon-like mediator, which enter a SHiP-like detector In order to generate E, θ values distributed in the same way, we have developed a numerical 2D-fitting algorithm which is fast, robust and automated. The main design concepts are based on the adaptive algorithms exploited in Monte Carlo integrators like VEGAS [32] and FOAM [33]. Our algorithm produces a 2D-mesh of bins for the 2D-histogram of the input points in the (E, θ)-plane in such a way that the histogram heights are flat. It is based on a deterministic procedure: we apply a sequence of alternate splittings, one along the x-axis and one along the y-axis, according to a democratic principle of equal weights. Further details on this technique can be found in Appendix A. As an example, in Figure 4, we show the 2D-mesh associated to the scatter data in Figure 3. Starting from this mesh, we can generate new E, θ points with the same distribution as the original ones. The result is plotted in Figure 5, where the goodness of the procedure can be appreciated by inspection.
We pass now to describe the generation of the azimuthal angle φ. As already stated above, at fixed θ, the φ distribution of the flux is flat due to the symmetry with respect to beam axis at the production point. The only kind of angular correlations which can be introduced are of geometrical nature and depend on the shape of the detector. An off-axis detector or even simply a detector on-axis with a rectangular surface exposed to the hidden particle flux, produces θ−dependent boundaries on φ. To be definite let us consider the on-axis detector as in the SHiP experiment. In this case, at fixed θ value, the projection of the points at the  detector surface obtained by varying the φ angle is simply a circle. If the circle is entirely contained in the detector surface, one can generate the azimuthal angle flat in the full range [0, 2π). Otherwise, one must take into account the geometrical intersections between the circle and the surface, and generate a φ value flat only in the 'inside' region. In a simplified notation, our prescription states that we pick a φ value uniformly in the range In practice, it may happen that there are more than one interval for the azimuthal angle, as it is indeed the case for the rectangular detector we considered above. In that situation, for some values of θ there may be (it depends on the customary sizes of the sides) up to four intervals in which φ can lie, which correspond to the four corners of the rectangle. In Figure 6, we compare the angular correlations as obtained by the original scatter points which enter into the detector and the ones we reconstructed following our strategy. Again, the agreement is fairly good. We stress that the above procedure is fully consistent with the re-weighting strategy outlined at the beginning, which we summarise here: 1. we first build the one-dimensional energy pdf on top of the original (unweighted) scatter points re-weighted by the effective function W , which takes into account, for example, that for some θ values, there are φ values not allowed; 2. we then reconstruct the missing variables in their actual ranges.
We exploited the same re-weighting strategy to modeled the full 3D-geometry of the detector. Indeed, up to this point, the same weight has been assigned to each particle direction. More technical details on this procedure are reported in Appendix C. Depending on the specific new physics model, the interaction between hidden particles and the SM matter in the detector can be based on different types of processes: • elastic scattering off electron; • DIS-like scattering off nuclei; • elastic and coherent scattering off nuclei.
In principle, the processes included in the above list can be easily simulated in our framework if a suitable model file is supplied, at least for the first two cases.

Displaced decays
In the case of displaced decays, the algorithm simplifies considerably. The decay process does not require the regeneration of the events and then the issue of their full reconstruction does not arise. Indeed, the decay can be generated event-by-event on top of the incoming flux of the unweighted events. The probability that a given particle decays in a specified decay channel i after having traveled a distance l from the production point is given by where Br i is the branching ratio for the i-channel, γ and β are, respectively, the Lorentz factor and the velocity in the laboratory frame of the decaying particle, and Γ is its width. The displacement from the production to the decay point can be determined starting from the partial decay widths. The latter, for the given BSM model, are computed on-the-fly according to the actual parameters of the simulation ( the so-called "auto-width" option).
parameter value # proton-on-target 2 · 10 20 detector configuration on-axis distance target-detector 5650 cm detector density 3.72 g/cm 3 detector shape parallelepiped x-side: 187 cm y-side: 69 cm z-side: 200 cm detector efficiency 1 Table 1: Main input values used for the simulation of the SHiP detector geometry, as reported in ref. [7].
This feature in combination with the "auto-scan" mode, both provided by MG5aMC, allows for a complete simulation scan over the relevant parameter space. In the present version of MadDump, the displaced decay events are not forced to be generated inside the actual decay vessel. It is only as the last step of the simulation, that a rejection of the events that occur outside the detection volume happens. The algorithm could be improved by reweighting each event according to the distance that the decaying particle could actually travel inside the decay vessel.

Illustrative examples
In this section we provide some illustrative applications of MadDump, considering three new physics models and making the corresponding predictions for the SHiP experimental setup [7]. We stress that the SHiP facility configuration used in the following is based on the one developed for the Technical Proposal in 2015 [7]. Since then, the SHiP Collaboration has continuously improved its setup aiming at higher sensitivity in the different channels. The newest setup as well as the corresponding background estimates are not yet available. The analyses reported in this paper will have to be redone once the updated information becomes available. Table 1 summaries the relevant input parameters which specify also the geometry of the apparatus. In particular, we determine the number of detectable events for each model in a typical run and then using the expected rates for the background processes reported in [7] we estimate the experimental sensitivities at 99.73% confidence level, which corresponds to the 3σ contour, and compare them with existing limits.

Quark-DM scattering: leptophobic portals
Let us first consider the case where an hidden particle that could be a dark matter candidate interacts with the visible sector via a new leptophobic force. This is a good benchmark model to study quark-DM scattering, in particular we will focus on signatures of deep inelastic scattering.
3.1.1 Vector portal: Baryonic U (1) B The simplest possibility for a leptophobic force mediated by a spin 1 particle is provided by models where the baryon number U (1) B is gauged such as where the actual mass generation mechanism is not relevant here. The quarks are the only SM fermions charged under this new gauge symmetry thus: while for the DM particle χ where the only important requirement on ψ χ is that it is long-lived enough to reach the detector. The U (1) B is anomalous and the cancellation of anomalies could lead to additional strong constraints as discussed in [34,35] . However, these constraints depend on whether the anomalies are canceled by fermions chiral or not under SM gauge symmetries, hence they are UV-dependent and we will not include them while comparing sensitivity of various low energy probes. The existing bounds on the Z coupling in the 1-10 GeV mass range come from the Z exchange induced invisible decays of quarkonia such as Υ → χχ and J/ψ → χχ) see [36]. Monojet searches at hadron colliders set a bound on g Z , the strongest one coming from a CDF search [37,38] at Tevatron, Moreover, existing and previous neutrino facilities like MiniBooNE could have sensitivity to few GeV leptophobic Z as discussed in [31,39,40] where it is shown that a reanalysis of existing data could set the strongest bounds on an ample region of the parameter space. In this study, we exploited the prompt production mode of MadDump for the generation of dark matter particles in the s-channel for an almost on-shell Z . We use our own version for the UFO file for this model, which is available in the MadDump directory. Considering the relevant parameter, the mass of the Z , in the range [2,10] GeV with equal steps of 1 GeV, we have generated 100 k production events (we remark that in each event the multiplicity of dark Figure 7: Sensitivity plot to a GeV leptophobic force at a SHiP like experiment obtained with MadDump assuming 9×10 5 neutrino deep inelastic background events [7]. matter particle is 2), while we generate 10k DIS interaction events between the dark matter particles and the detector. The fraction of events which passes the detector acceptance ranges from ∼ 14% for m Z = 2 GeV to ∼ 1% for m Z = 7 GeV. Exploiting a small workstation with sixteen cores, we have got the following timings for the complete simulation of each benchmark point: • production : ∼ 1 • fit : ∼ 1 • interaction : ∼ 3 .
Based on [7] we assume 9 × 10 5 background events. In Fig. 7 we present the potential sensitivity of SHiP for 2 × 10 20 POT to this scenario and we compare it with the above mentioned existing constraints.

Leptophobic scalar and pseudo-scalar portal
Another interesting possibility to consider is a leptophobic force mediated by a scalar or pseudo-scalar particle.
We consider the following simplified model with i = u, d, s, c, b, t; as before, χ is a Dirac fermion stable or long-lived enough to cross the SHiP detector. In the case of a real scalar S the interaction with quarks could arise via the renormalisable interaction which induces a singlet-Higgs mixing sin θ such as the singlet S inherits couplings to the SM fermions: where y i is the SM Yukawa coupling of the fermion i. A different flavor structure from the SM for the singlet-fermion coupling could be arranged via the dimension five operators, that is: and where Λ is the cutoff above which either extra Higgs bosons or vector-like leptons are expected.
Depending on the origin of the interaction among S and the quarks, bounds from Higgs invisible decay and/or electroweak precision measurements could be relevant. However, we will not discuss them due to their dependence on the UV-completion. Moreover, we assume that CP is a good symmetry of the Lagrangian (see for instance [41] for a discussion of possible constraints). We consider the benchmark scenario where the scalar S couples to up and down quarks (see [41]). We have used the same setup of the previous case, 100 k prompt production events and 10 k DIS interaction events. As for the UFO file for the scalar model, we exploited the general model in ref. [42]. The program scanned over the relevant parameter, the mass of the scalar mediator S, in the range [2,10] GeV with equal steps of 1 GeV. The fraction of events which passes the detector acceptance ranges from ∼ 16% for m Z = 2 GeV to ∼ 4% for m Z = 7 GeV. The timings are analogous to those of the previous case.
In Fig. 8 we show the sensitivity at SHiP considering as before DIS as signal events. In this case the only existing bounds come from the CDF monojets bounds [7,37] and we notice that SHiP could constrain new regions of the parameter space with this proposed analysis.

Electron-DM scattering: the dark photon
As benchmark model to study DM-electron scattering we consider a new gauge boson associated to an abelian gauge symmetry U (1) , A , kinetically mixed with the photon [43], namely Figure 8: Sensitivity plot to a GeV scalar mediator at a SHiP like experiment obtained with MadDump assuming 9 × 10 5 neutrino DIS events compared to existing bounds from CDF [7,37].
a dark photon (DP). The relevant Lagrangian corresponds to: where is the DP-photon kinetic mixing. We further assume the existence of a particle χ either a scalar or a fermion charged under the new gauge symmetry U (1) and stable at least compared to the scale of SHiP, hence we also add the following Lagrangian: We choose as benchmark point m χ = m A /3 and α D = g 2 D 4π = 0.5 as in [44]. For this case study, we have used the decay-interaction mode of MadDump. The incoming meson fluxes has been generated with Pythia8, having care to store only the final state mesons which decayed directly in photons. This means that in a decay chain η → 3π 0 → 6γ only the 3π 0 are stored, while a η meson is stored in the list if it decayed directly into 2γ. With this caveat, we can limit ourselves to consider only one decay channel for each of the mesons included in our analysis: • π 0 → 2γ; • η → 2γ; • ω → γπ 0 .
We exploited the general UFO model for spin-1 as reference model for the DP [42,45]. The meson decays has been modeled applying a standard effective field theory (EFT) approach. Indeed, for the most interesting case in which the DP is almost on shell, the decay process can be approximated by the tree-level vertex depicted in Fig. 9 at first order in the EFT expansion. We added the meson particles and the minimum set of new interactions required to deal with Figure 9: Effective field theory approach to production of DP from rare meson decay.
their decays directly in the UFO file, on top of the reference model. We simulated 100 k proton on target (POT) events, which, in terms of meson yields/POT, resulted in: ∼ 6/POT for pions, ∼ 0.3/POT for η's and ∼ 0.07/POT and for ω's. We considered one meson species at time, and scanned over the relevant parameter space, for masses of the DP below the corresponding meson mass. For the case of the pions, which are the most numerous particles, the time per scan point has been ∼ 20 on a 4-cores CPU. The most time consuming tasks are the I/O operations related to the meson decay process, which took ∼ 14 of the whole time per benchmark point. In Fig. 10 we present in the (m A , ) plane the SHiP sensitivity compared to existing bounds, described in details in [44]. In the region of interest, the strong experimental constraints come from the monophoton BaBar search [46] and NA64 [47] via a missing energy analysis. Assuming α D = 0.5, experiments looking at electron-DM scattering such as MiniBooNE [48], LSND [49], and E137 [50] achieve a better sensitivity than NA64 so their reach is also presented here. An even stronger reach for m A 300 MeV could be reached by NOνA experiment at Fermilab by recasting existing data as shown in [51]. For electron scattering events, according to the SHiP technical proposal [7], we consider the following selection cuts: where θ χ−e is the angle between the incoming DM particle and the outgoing electron. We assume 284 electron-neutrino scattering events as background following [7]. We simulated  Figure 10: SHiP sensitivity for a Sub-GeV dark photon emitted by rare meson decays of π 0 , η and ω obtained with MadDump. For all production modes both electron-DM scattering and DIS are considered considering respectively 284 and 9000 background events. The gray region is excluded by current probes, the strongest are BaBar monophoton search [46] and MiniBooNE [48]. The dark gauge couplings is chosen to be α D = 0.5 and the ration between DP and DM mass is 3. events both for electron-φ χ scattering and DIS events comparing their sensitivity in Fig. 10. As expected, the electron sensitivity is significantly better than the one achievable with DIS. Our prediction here is conservative because we do not include potentially important contributions to the production stage like the decays from mesons produced in the cascade process and the prompt production.

Conclusion
In this paper we have presented a new MG5aMC plugin called MadDump that allows the generation of events where the production of a particle and its detection are separated by a long distance. In order to install it, it is enough to type "install maddump" within MG5aMC 3 . The main input provided by the user are the geometry of the experiment and the physics model under investigation. With these ingredients event generation corresponding to one or more benchmark scenarios can be performed automatically. We have shown illustrative examples based on different BSM scenarios and production/detection mechanisms and computed the corresponding number of events that would be produced at the SHiP experiment. The framework is fully general and can be applied to any BSM model and experiment at a beam dump facility that aims to test it. Our tool could be employed for a number of studies, from the search of new feebly interacting particles to the study of elusive SM processes like tau neutrino cross section at present and future beam dump experiments.

Acknowledgments
We are thankful to Eduardo Cortina, Giovanni De Lellis, Antonia Di Crescenzo, Jan Hajer and Roberta Volpe for useful comments on the manuscript.

A Techniques for event generation
Consider the problem of generating unweighted points in a 2D-space according to the distribution P (x, y)dxdy = f (x, y)g(x, y)dxdy , (A.1) where g(x, y) is a modulation function whose expression is supposed to be known analytically. More precisely, we are mainly interested in the problem of generating unweighted y values given a fixed x = x according to the profile function φ(y)dy = P (x, y) dy .

(A.2)
When the function f (x, y) is given in closed form the problem reduces to generating points according to a given function P (x, y) and it can be accomplished by standard Monte Carlo techniques, using for example the classic hit-or-miss algorithm. Here we consider the more interesting situation in which the function f (x, y) is only available numerically indirectly from sample of events. One can re-interpret it as a fitting problem. The function f (x, y) is given in an approximated way, with a level of precision which can be in principle reduced at will (by generating more points) but at each step it is finite, as a 2D histogram built out of events. In its standard formulation, the task of obtaining a fit from an histogram consists basically of two different parts: • the choice of the model, i.e. an n-parameter family of functions together with the cost function of the fit; • minimisation of the cost function.
In so doing the result is given by a function supplemented with extra information on the accuracy of the fit (covariance matrices, goodness of the fit, etc). Fitting a function, even in the "simple" 2D case, is however not always trivial. In particular, aside the technical aspects underlying the minimisation procedure, a certain amount of knowledge of the function to be fitted is required (in order to choose a reasonable class of models). For our purposes, no apriori assumptions can be formulated on the behavior of the function f (x, y), as, in general, it can result from very different classes of physical processes. For this reason, we look for a procedure that allows to automatize the process. Though an analytical fit has advantages (including also the possibility of smoothing a discrete data sample in to continuous distribution), such a level of accuracy is not strictly required in order to perform the generation of the unweighted events and we do not adopt it. Our approach is based on importance sampling and variance reduction methods implemented in Monte Carlo integrator algorithms. The strategy is based on the flattening of the integrand function via a numerical adaptation of the integration grid. Moreover, once the grids are available, they can be used to regenerate unweighted points according to the integral function. In our implementation, we have devised a complete deterministic procedure to construct a grid, very closely the above concept of adaptive grid. In what follows, we will give a detailed description of our algorithm together with some validation examples.

A.1 Grid construction
As a first example, let us consider the case in which the modulation function g(x, y) reduces to the identity map. As simple as it may appear (in this case P (x, y)dxdy = f (x, y)dxdy and an unweighted generator for that distribution is assumed to be known), it allows us to clarify a few useful points. First, it may happen that generating events with the grid is more efficient (for what concern both time and space resources), or more usable in some sense, than exploiting the original generator. This is in fact the case in our applications, in which the unweighted generator has a very complex structure and the 2D events we interested in are a tiny part of the whole result. For this reason, it is not only an illustrative case, leading to a clearer illustration of the basic concepts, but it is relevant per se. We assume to have at our disposal a sample of N unweighted points (x, y) distributed according to the function f (x, y). Our aim is to generate unweighted events distributed according to the same distribution. Starting from the available points, the profile of the function f (x, y) is given by the heights of a 2D histogram with bins of equal size. The idea is that of resizing the bins in such a way to flatten the histogram, or, equivalently, to have the same number of points lying in each bin. In this way, the distribution of the bins will follow the behavior of the function: they will be denser and smaller near the region where f (x, y) is peaked and sparse and bigger where it is flat. The resulting 2D map will retain almost the full information of the 3D plot, and it is very similar to the idea of a contour plot. In order to obtain such a parametrisation, we employ a decision-tree-like algorithm, which is very simple and efficient. Before describing it, a technical remark is needed. Adaptive grids are often constructed with lines parallel to the axis coordinates. This is efficient in all the situations where the function can be expressed in the factorised form A great improvement is given by an approach in which irregular grids, made of cells of different sizes, are allowed, as in the case of the FOAM algorithm [33]. The cells adapt better to the behavior of the function reproducing it in a more faithful way, for example near circular peaks. The cell represents the basic object of our algorithm. A cell can be split in two cells along the x-direction (horizontally) or the y-direction (vertically). Given these basic ingredients, the algorithm proceeds as follows: 1. start from a cell containing all the available points; 2. alternate an horizontal split and a vertical split, in such a way that, in each of the two splits, half of the point fall in a subcell and half in the other one; 3. repeat step 2 for each subcell until the number of point for cell is lesser/equal than a prefixed value (exit condition parameter ).
It is clear that the above procedure gives the grid we are looking for. The exit condition parameter n min controls the grain of the mesh. The choice of its value is based on the compromise between having it small for a finer grain and having a sufficient number of points per bin to be statistical significant. We now restore the proper role of the modulation function g(x, y), which as mentioned above, can be arbitrary yet to be expressed in an analytical form. It can be viewed as a reweighting of the original sample of points: where we have conventionally set to 1 the common weight of the unweighted sample. We have P uniform = g(x, y)f (x, y)dxdy = g(ξ, η)dξdη = g f where the notation · pdf denotes the average wrt the pdf in the subscript. Under the hypothesis f (x, y) is a distribution function, a well-defined change of variables is implicitly given by the relation f (x, y)dxdy = dξdη where the function f is the Jacobian of the transformation. Furthermore, if also g has a definite sign it is possible to perform an extra change of variables This relation proves the equivalence between unweighted generation of the product distribution f · g and the generation reweighted by g starting from a sample of unweighted points generated according to f . We are now ready to generalise the previous case. We require that the rebinning procedure leads to a grid with (almost) the same weight w for each bin, with the definition of the weight w(b) of the bin b given by The generalisation of the algorithm is straightforward 1. start from a cell containing all the available points; 2. alternate an horizontal split and a vertical split, in such a way that, in each of the two splits, each subcell have half the weight; 3. repeat step 2 for each subcell until its weight is greater than the prefixed value w (exit condition parameter ).
The exit condition parameter can be chosen of the form where n min has the same meaning as before and α is a dimensional factor which can be adjusted in the direction of refining the grain or increasing the number of points per bin. By construction, the above procedure cannot handle distributions which vanish on some regions inside the fitted domain. This limitation is particularly severe in the case the distribution presents a falling-down tail and vanishes inside the fitted region. Indeed, even if the cells become larger and larger when approaching the tail, there is always a non null probability to generate points inside them, also in the empty regions. In this way, unphysical points are generated. In order to milder this limitation, we implemented a further refinement step after the mesh has been constructed. The peripheral cells, i.e. the cells which share a side with the frame of the fitted regions, are reshaped in a such a way to limit the cell to the actual region populated by the input points. We refer again to Fig. 4 in section 2 to appreciate the reliability of this improvement.

A.2 Example
As a validation example, we consider the situation in which both the function f (x, y) and the modulation g(x, y) are given analytically in order to show that the algorithm works correctly. Furthermore, we test its robustness considering a highly non-trivial case in which the modulation affects and distorts in a severe manner the original function. We take a simple Gaussian function, (see Fig. 11a): f (x, y) = e −(x 2 +y 2 ) with the following modulation g(x, y) = 2 x 2 + y 2 cos 2 (2y) + 1.
The product function is shown in the 3D plot in Fig. 11b. In Fig. 12a and Fig. 12b we report the corresponding 2D meshes obtained with our algorithm for 100 k and 1 M input points respectively. We stress that the starting point has been the generation of a sample of unweighted points distributed according to the Gaussian function. We have reweighted the points according to the modulation function and then we have applied our algorithm for weighted events. The algorithm reproduces faithfully the behaviour of the function with a level of accuracy which, as expected, improves with the number of input points.
(a) Simple gaussian distribution as f function. (b) Full distribution after modulation is applied. Figure 11: 3D-plot of the analytical distributions occuring in the validation example.

A.3 Generation of unweighted events
Generating a 2D sample starting from the 2D mesh with the same distribution as f (x, y) is trivial. By construction -25 -  • the probability of generating a point in a given cell is proportional to the inverse of its area, • inside a cell, the probability of generating a point is uniform, and therefore it is enough to generate an equal number of points uniformly in each cell. Let us now turn to the issue of generating an y value according to the profile function Eq. (A.2) at a given x = x point. We introduce a small resolution parameter related to the x variable such that the x = x value is fixed within the interval [x − , x + ]. Then, the thin stripe centered in x = x with width 2 and parallel to the y-axis will intercept the mesh in a subset S of cells. We associate a normalised weight to each cell in S proportional to the ratio of the overlapping area between the stripe and the cell over the total area of the cell. Then, we pick a cell according to the value of these weights by generating a uniform random number in the interval [0, 1]. Finally, we generate a uniform y value within the cell. This construction solves our problem, i.e. the y values are distributed according to the profile function φ. The procedure is independent of in the limit → 0. In practice, this means that magnitude of should be chosen as a fraction of the minimum x-width of the cells of the mesh.

A.4 Example
Since an example of the generation of the entire 2D sample has been already shown in section 2, here we focus on the constrained one-dimensional generation. Let us consider again the previous example and fix a x value, for instance x = 0. We generate y values distributed as the profile function (A.2) according to the above procedure. In Figure 13, we plot the comparison between the generated points and the analytical profile function φ using our meshes with 100 k points ( 13a) and 1 M points ( 13b). The generated histograms are in good agreement with the analytical curve and they reproduce well also the sharp deep in y = 0. The result improves by exploiting the mesh with a greater number of points giving a solid indication that the procedure is asymptotically converging to the true distribution.
(a) Using the 100 k mesh.
(b) Using the 1 M mesh Figure 13: Generation of an unweighted sample according to the 1D profile function at x = 0. The distribution is normalised to 1.

B Systematic uncertainties
In summary, our fitting procedure is a way to approximate the probability density function associated to a given 2D scatter data with a piecewise function, i.e a histogram, with an automatic choice of the bins. In the example of the previous section, we have provided a non-trivial numerical proof of concept of it. Moreover, since in that case, the analytical distribution is known a-priori, we have a full control on the systematics uncertainties and the convergence of the method. This is not the case in the practical applications, where the probability distribution is available only in the form of scatter data. As a consequence, estimating the systematics of the approximation becomes more difficult. We follow a pragmatic approach which should provide a guideline for the user to tame the systematics according to his own scopes. Despite the fact that this systematics can be made smaller and smaller by providing more and more statistics (initial input events), in practice a compromise between the accuracy goal and the actual computational resources needs to be found. Let us start from some basic and general considerations. First, one can always separate the prediction of the total rates (including the geometrical acceptance of the detector) from those of more exclusive observables, as the angular distributions. Since the physical interaction cross section depends only on the energy, the inclusive total rates depend only on the effective dark sector particle energy distribution introduced in section 2.2.1 eq.(2.4). In our approach, we do not rely on the 2D fit to obtain this quantity. Instead, we perform a dedicated 1D fit exploiting a smoother class of functions. In this way we have a better control on the result and also on its uncertainties. Indeed, a 1D fit is a simpler operation and, since we are integrating over angles, we have access to a higher level of statistics. The 1D fit works as follows. Starting from the input weighted data, we first build a 1D histogram and we assign to each energy bin the usual Poisson uncertainty. We then fit the histogram using a weighted cubic splines fitting. In order to assess the error on the fit, we add the possibility to vary the values within the histograms uncertainty bands. This can be done setting the parameter rescale fac in the fit2D card card file. which takes values in the interval [−1, 1], where 0 stands for the central value, 1 for the upper limit, -1 for the lower one. A realistic study case is given in the following subsection.
We pass now to discussing the case of the 2D fit. At a fixed number of input events, the algorithm depends mainly on a unique parameters which gives the exit condition of the splitting loop. Naively, it represents the number of points which lie in each bin. Hence, there is a competition between choosing it small, to have a better description of the shape of the distribution (more bins with a smaller size), or choosing it large, to avoid to be overwhelmed by the statistical fluctuations (less bins with a larger size). The user can change this parameter by setting the value of the npoints cell variable in the fit2D card card file (the default value is 50). We have introduced in the code a consistency test (that can be enabled by setting to True the flag fit syst, again in the fit2D card card file) which compares the mesh obtained by varying the central value of this parameter by a factor of 1/2 and a factor of 2. To this aim, we consider the classifier where x represents a generic event, and P i , i = 1, 2, are the probability densities we are comparing. The values of D ranges over the interval [0, 1]. An average value D ∼ 1 means that P 2 (x) underestimates P 1 (x), while for D ∼ 0 we have the opposite. For P 1 (x) = P 2 (x), D(x) = 1/2. Hence, in the case the average value D ∼ 1/2 and its standard deviation σ D is small, we cannot distinguish between the two probability functions. The study of the D classifier put on a quantitative foot the qualitative results given by the visual inspection of the mesh grids. Its mean and standard deviation give us a measure of the global goodness of the fit. Furthermore, we can perform also a more local comparison of the angular shape at fixed value of the energy variable. We postpone a more detailed discussion to the following subsection in which we apply the above analyses to a realistic case study.

B.1 Study case: leptophobic model
As a study case, we consider the example of the leptophobic model presented in the main part of this work. Since the generation of the input DM events can be simulated directly internally in MadDump, we have a direct access to input samples of different number of events. Our setup is outlined in the script reported in the listing appendix. We select the value m = 2 GeV for the mass of the DM mediator. We analyse first the uncertainties on the total rates. We reported the predictions for the DM yields in Tab. 2 for input samples of increasing statistics. The uncertainties on the predictions correspond to the 1D fit variation around the histogram error bars, as stated above. The three results are consistent within their uncertainties and, as expected, the accuracy improves increasing the statistics. We observe, that in this case, for a sample of 100 k input events, which corresponds to 100 k ×0.28 = 28 k DM particles passing the geometrical acceptance, we already get a result accurate at the few percent level. Note that here we absorb the multiplicity factor of 2 within the definition of the geometrical acceptance The default value of the npoints cell parameter is 50. We study what happens by varying it by a factor of 1/2 and a factor of 2. We briefly explain the strategy we followed for the comparison. We produce the three meshes corresponding to the three values npoints cell = 25, 50, 100. We use the result corresponding to the central value as our reference point and we denoted by P (x) the corresponding probability density. This means that we evaluate P starting from the 2D mesh as follows  Table 2: Estimates of the the total DM yields for increasing numbers of input events. The uncertainties refer to the variation around the 1D fit error bars as explained in the main text. The number of input here does not include the geometrical acceptance , which is reported in the second column.
where n cells is the number of cells of the mesh and A(x) is the area of the cells in which the point x lies. We build the two classifiers where P l and P h are the probabilities densities corresponding respectively to the lower and the higher values of npoints cell. We generate random (uniformly distributed) points x = (E, ϑ) and compute the mean and standard deviation of the two classifiers. We consider the simple unbiased estimator given by the uniform average We reported our results in the third and fifth columns of Tab. 3, respectively for D l and D h . All the mean values are fairly consistent with 0.5, which is the indication that the different meshes are consistent among themselves. The standard deviation is lower for D h , which is what is indeed expected since, with a lower npoints cell, the fit is more sensitive to the statistical fluctuations of the original scatter data. Furthermore, we observe that the standard deviation decreases increasing the statistics from 10 k to 100 k events but, then, there are not any improvements from 100 k to 1 M. The explanation for this behavior is related to the vanishing tail of the distributions. Indeed, since we are fitting using piecewise functions, the accuracy of the method is worse in the long vanishing tails, where we decided to exploit a cut prescription instead of spreading uniformly the weights on a very big cell. Then, the error is dominated by the fluctuations near the boundary regions, where, however, the probability densities is approaching zero. To test quantitatively this argument, we re-weight the events accordingly to our reference probability P (x), and we consider the weighted estimators 10 k 0.27 0.55 ± 0.20 0.54 ± 0.03 0.48 ± 0.12 0.549 ± 0.013 100 k 0.282 0.54 ± 0.12 0.525 ± 0.015 0.50 ± 0.08 0.518 ± 0.005 1 M 0.282 0.55 ± 0.13 0.511 ± 0.005 0.47 ± 0.09 0.510 ± 0.004 Table 3: Comparison of the classifiers D f , with f = l, h calculated as weighted and unweighted averages.
where the weights are given by The results are reported in the fourth and sixth columns of Tab. 3. The errors drop significantly and for the largest input sample we see that it is not possible from the practical point of view to distinguish the probability densities given by the three meshes, so that one, for instance, might choose to use the mesh corresponding to npoints cell=25 since it is the finest (it has more cells wrt the other two).
Another important aspect concerns the convergence of our method to reproduce the original data set. With this, we mean the minimum number of regenerated events N gen needed to have a distribution which is consistent with that of the input data. Indeed, for N gen small, we expect the regenerated distribution is dominated by the statistical fluctuations. On the other hand, we expect that after having reached the desired level of agreement, further generations of events will not spoil the convergence. The naive expectation would be to have N min gen N input . For a quantitative analysis, we rely again on the classifiers introduced above. In the following, we outline our strategy. We fix the mesh associated to the central value npoints cell=50 as our reference for the 2D-dimensional data distribution. Starting from this mesh, we regenerate samples of events with increasing statistics. We perform a second fit on top of each regenerated sample obtaining new meshes. We assume that these meshes represents the approximate bi-dimensional distributions of the samples, according to eq. (B.2). We finally compute the two averages for each of the classifier D f (x), where now f = N gen labels the regenerated samples. We report our result for the input sample of 100 k events. In this case, the effective number of input data events is 30 k, i.e. the events passing the geometrical acceptance cuts of the detector. The results are reported in Tab. 4. They confirm on the quantitative ground our naive expectations N min gen N input , leading to the prescription N min gen ∼ 2 − 3N input . We conclude this discussion on the fit systematics showing some lesser inclusive results at the level of differential distributions wrt to the polar angle. We adopted the following strategy. We select an energy value E. Then, we pick the cells in the 2D mesh fit in which E is included. We consider as energy resolution parameter some multiple (the default is 2) of the minimum energy width of the selected cells. This allows us to consider an energy bin centered in E N gen  with width given by the resolution parameter. Then, we consider the input data and the regenerated ones which lie in this bin and we compared their angular distributions given by standard 1D-histogram. Our choice of the energy bin size guarantees that the statistics is comparable for any starting E values. We have analyzed and compared the results obtained with the two input samples 100 k and 1 M. They are shown in Fig. 14 and Fig. ??, respectively. In both cases, we have used regenerated samples with N gen = 300 k. We found a good agreement both between data distributions for different statistics and between data and our fits which is of the order of ∼ 25 − 30% for the 100 k case and ∼ 10% for the 1 M one. As expected, in the region corresponding to the bulk of the events, the corresponding energy bin size are smaller. For example, in our study case, we observe that in the central energy range 40 GeV < E < 150 GeV the energy bin size is of order ∼ 1 GeV for the 100 k input sample and the situation fairly improves for the 1 M one. On the contrary, for the value E = 20 GeV which lies on the tail of the distribution, we need a bigger bin size, ∼ 16 − 18 GeV. The analyses performed in this section are encoded in MadDump and the user can reproduce the same studies for his particular situation. As a rule of thumb, to take cum grano salis, 500 k events entering the detector can be consider a reasonable amount of input statistics. Together with the default settings of the internal MadDump parameters and the choice of N min gen ∼ 2 − 3N input this should lead to an uncertainty of 1% for the total rates and 5 − 10% on the angular distributions, which is usually lesser than the other systematics of the simulation.
C A method to take into account the depth of the detector Consider the scattering of a flux of incident particles impinging on a thick target. Let us fix the geometry of the problem and consider for instance a parallelepiped shape for the fiducial volume of the target. In general, the incident flux is neither collimated nor mono-energetic. For the sake of simplicity, we consider the flux to be originating from a point-like source placed on the target axis and mono-energetic. Consider a cartesian 3D reference frame in which the z-axis is along the "depth" of the target and x and y are parallel to the other dimensions.
Let us subdivide the target in thin sheets along the z-axis. The number of scattered events in each sheet is given by where S(z) is the surface of the shell at z, F is the flux of incident particles impinging on the sheet, ρ is the number density of the target particles and σ is the interaction cross section. In the formula above, we assumed that the cross section is constant all over the surface of the sheet. Furthermore, we assume that the cross section is constant over the whole fiducial volume of the target and we consider a uniform target. Then, integrating also over the depth of the target, we get the master formula for the number of scattered events We are interested in the situation in which the cross section is very small, i.e. we can neglect (at least in first approximation) the variation of the flux F due to the scattered particles. This means that the dependence of ∂ 2 F /∂x∂y on z is purely geometrical: from a given configuration at a point z * it is possible to construct the flux at a new z point by prolongating the flying direction of the particles in the flux. For this reason, in the above example, even though the surface of the sheet is constant, the number of particles impinging on the different sheets along the z-axis is different: S dxdy ∂ 2 F ∂x∂y (x, y, z) = n(z) .

(C.3)
In Monte Carlo integration/generation this translates in employing different weights for the bunch of events describing the incident flux. It is still possible to use unweighted events if one considers a variant of the hit-or-miss rejection method. Here, by "unweighted events" we mean that the points have been generated according to the dxdy ∂ 2 F ∂x∂y distribution. Since the integrand is positive definite, enlarging the integration region we obtain the inequality S(z) ≥ S =⇒ S dxdy ∂ 2 F ∂x∂y (x, y, z) < S(z) dxdy ∂ 2 F ∂x∂y (x, y, z) ≡ n * (z) .

(C.4)
In particular, we can choose S(z) such that the above integral n * is constant. This means that we enlarge the surface according to the radial projection starting from the point-source of the flux. Then, for the new integral, we can employ unweighted events for the flux: A full event is given once a z variable or equivalently a travel distance along the flying direction of the event is generated uniformly between the minimum and the maximum value inside the largest volume. Then, we accept or reject the event whether it lies or not in the true fiducial volume. While correct, this method is not efficient as generated events can be rejected. An alternative approach entails applying a simple reweighting procedure. Intuitively, we just have to penalize the events that would be produced by particles crossing the fiducial volume of the detector over smaller paths. Indeed, a given event may contribute or not to n(z) in Eq. (C.3) depending whether at z it is inside or not the integration region: ∂x∂y (x(z), y(z), z 0 ), z < d(x, y) 0, z > d(x, y) where we denoted with d(x, y) the z-distance after which the event goes out of the integration region. However, the weights retain a dependence on z due to the presence in the argument of the function of (x(z),y(z)), which represent the coordinates of the particle when it crossed the sheet at z 0 . If we replace x, y-coordinates by angular ones (θ, φ) which gives the flying direction of the events, the weight will not have any residual dependence on z but the theta function Θ(z − d(θ, φ)). This can be simply taken into account by reweighting the events as In terms of the travel distance inside the fiducial volume of the detector r(θ, φ), we have d(θ, φ) = r(θ, φ) cos(θ). (C.8) Then, using this reweighting strategy, we can reconstruct the full event by generating uniformly the z or, equivalently, the traveled distance variable according to the actual minimummaximum allowed by the geometry of the target. Notice that the number of scattered events in the two cases is given by: N = ρσ(z 1 − z 0 )(n * − n rejected ) ρσ(z 1 − z 0 ) n(z) , (C.9) which implies the integral condition: (n * − n rejected ) = n(z) . (C.10)

D Listings
In this appendix we report the input script files used for the examples presented in the main text.