Simulating low-energy neutrino interactions with MARLEY

Monte Carlo event generators are a critical tool for the interpretation of data obtained by neutrino experiments. Several modern event generators are available which are well-suited to the GeV energy scale used in studies of accelerator neutrinos. However, theoretical modeling differences make their immediate application to lower energies difficult. In this paper, I present a new event generator, MARLEY, which is designed to better address the simulation needs of the low-energy (tens of MeV and below) neutrino community. The code is written in C++14 with an optional interface to the popular ROOT data analysis framework. The current release of MARLEY (version 1.2.0) emphasizes simulations of the reaction $^{40}$Ar$(\nu_e, e^{-})^{40}$K$^{*}$ but is extensible to other channels with suitable user input. This paper provides detailed documentation of MARLEY's implementation and usage, including guidance on how generated events may be analyzed and how MARLEY may be interfaced with external codes such as Geant4. Further information about MARLEY is available on the official website at http://www.marleygen.org.


Introduction
Monte Carlo event generators are a widely-used tool in nuclear and particle physics. These computer programs implement probabilistic models of physics processes and produce corresponding sets of events: lists of particles (represented by their charges, 4-momenta, etc.) involved in simulated interactions.
While helpful as an aid to theoretical calculations, event generators are also often used by experiments for designing detectors, estimating efficiencies and backgrounds, assessing systematic uncertainties, and interpreting the results of measurements. Examples of event generators include PYTHIA [5] and Herwig [6] for high-energy particle collisions, HYDJET++ [7] for relativistic heavy ion collisions specifically, FREYA [8] for fission, and DECAY4 [9] for radioactive decays of unstable isotopes.
For studies of neutrinos, much community effort has been directed toward the development of event generators suitable for use by accelerator-based experiments at facilities like J-PARC [10] and Fermilab [11]. These experiments employ beams of primarily muon-flavor neutrinos which are produced over a broad energy range between the low hundreds of MeV to the tens of GeV. Neutrino scattering on atomic nuclei is the primary means of detection, and several modern event generators provide widely-used models of the relevant physics, including GENIE [12], GiBUU [13], NEUT [14], and NuWro [15].
Although important differences exist between each of these generators, all share a similar conceptual treatment of neutrino-nucleus interactions. Each scattering event on a complex nucleus is taken to involve a neutrino striking a single bound nucleon. 1 A removal energy and initial 3-momentum are assigned to the struck nucleon using a model of the nuclear ground state. Traditionally this is done using a variant of the Fermi gas model (e.g., that of ref. [17]), but implementations of more sophisticated treatments, such as the Correlated Basis Function approach [18,19], are beginning to become available [20].
With the initial state fully defined, the generator then simulates production of the particles that emerge from the neutrino-nucleon interaction vertex. A variety of models must be used at this stage due to competition between multiple nucleon-level processes (e.g, quasi-elastic and deep inelastic scattering) which may occur at accelerator neutrino energies.
Because they are subject to the strong force, outgoing hadrons from the primary neutrino interaction will often rescatter within the nuclear medium. These final-state interactions (FSIs) can have a pronounced effect on the kinematics and multiplicities of the hadrons that ultimately exit the nucleus. Intranuclear hadron transport and FSIs are handled in GENIE, NEUT, and NuWro using variations of the intranuclear cascade (INC) model [21][22][23]. This model assumes that propagation of a hadron h through the nucleus may be described in terms of a mean free path where E is the hadron energy, r is its radial position within the nucleus, and ρ is the number density of nucleons. The total cross section for h scattering on a free nucleon, σ hN , may be used directly but is typically modified with approximate corrections for nuclear effects. The GiBUU code simulates intranuclear hadron transport using a semi-classical model which considers the time evolution of the phase space density for each hadronic species. The equations describing the behavior of distinct kinds of hadrons are coupled through the common nuclear mean field and through a collision term which represents the influence of FSIs. The numerical implementation adopts a test-particle ansatz to solve the relativistic Boltzmann-Uehling-Uhlenbeck (BUU) equation. Extensive documentation of GiBUU's treatment of various nuclear reactions, including neutrino-nucleus scattering, is available in refs. [13,24]. The INC models used by other generators for FSIs may largely be regarded as approximate simplifications of the BUU approach [25, sec. 2].
When all produced particles either escape the nucleus or are re-absorbed by it, the hadron transport stage of the simulation is complete. This typically marks the end of the physics workflow needed to generate a single neutrino-nucleus scattering event. The finished events may be used both for standalone calculations of kinematic distributions and as input to later stages of an experiment's detector simulation chain. 2 In addition to accelerator-based neutrino oscillation experiments, there is also worldwide interest in pursuing detailed measurements of lower-energy (tens-of-MeV and below) neutrinos produced by supernovae [27][28][29], by the Sun [30,31], and by terrestrial facilities via pion and muon decays at rest [32][33][34][35][36]. Several distinct reaction modes will enable neutrino detection in these measurements. Elastic scattering on electrons and protons [37] is flavor-blind but produces signal events down to arbitrarily low neutrino energies. Coherent elastic neutrino-nucleus scattering (CEvNS), a neutralcurrent (NC) process in which a neutrino scatters off of a complex nucleus and leaves it in its ground state [38], shares these properties and was recently observed for the first time by the COHERENT experiment [39,40]. Apart from reactions involving complex nuclei, the only inelastic channel available at tens-of-MeV energies is charged-current (CC) absorption ofν e on free protons via the inverse beta decay (IBD) reaction 3ν e + p → n + e + . ( The IBD cross section is precisely known [41,42] and dominates the expected signal for supernova neutrinos in water Cherenkov and liquid scintillator detectors [28]. The low-energy neutrino interaction modes listed so far do not present any special difficulties from an event generator perspective. Relatively simple expressions are adequate to compute differential cross sections for all but the most precise calculations of these processes, and, with the recent addition of a CEvNS model [43], the GENIE event generator currently provides an implementation of all of them that may be suitable for use by low-energy neutrino experiments. In contrast to other low-energy processes, inelastic neutrino scattering on all but the lightest complex nuclei (e.g., deuterium) is theoretically cumbersome, requiring an elaborate treatment of nuclear physics in order to fully describe the interactions. Despite the significant challenges involved, however, obtaining realistic simulations of low-energy inelastic neutrino-nucleus scattering is highly desirable for a variety of scientific applications. Principal among these is detection of low-energy astrophysical neutrinos by the upcoming Deep Underground Neutrino Experiment (DUNE) [44]. Thanks to the experiment's planned use of four ten-kiloton liquid argon time projection chambers (LArTPCs), DUNE has the potential to perform detailed measurements of these neutrinos via the charged-current reaction ν e + 40 Ar → e − + 40 K * .
This process is anticipated to provide most of the signal in DUNE for the neutrino burst associated with a galactic core-collapse supernova, granting the experiment the potential for unique sensitivity among large detectors to the ν e component of the supernova neutrino flux [45,46]. DUNE also shows substantial promise as a detector for studying solar neutrinos [30]. Primary sensitivity to low-energy ν e is shared with DUNE by the 79-ton Helium and Lead Observatory (HALO) [47][48][49], but in this case the detection technique is indirect. Inelastic neutrino reactions on lead nuclei, such as will sometimes lead to the production of neutrino-induced neutrons (NINs) via de-excitations of the residual nucleus. HALO's lead neutrino target is instrumented with 3 He-based neutron counters, which will be used to search for NINs in the event of a nearby core-collapse supernova. Due to its widespread use as a radiation shielding material, lead is also a potential source of background NINs in precision CEvNS measurements. To better constrain theoretical modeling of this background, the COHERENT experiment is pursuing direct measurements of NIN production on lead, iron, and copper [35].
Despite the successes of standard neutrino event generators in describing accelerator neutrino data, their prevailing treatment of inelastic neutrino-nucleus scattering is likely to be inadequate when applied to interactions at energies of tens of MeV and below. This is due in part to approximations made in their modeling of nuclear structure. For few-MeV neutrinos, inelastic neutrino-nucleus cross sections are governed by the energetically-accessible transitions to low-lying discrete energy levels of the daughter nucleus. At somewhat higher energies, excitations of collective vibrational modes of the nucleus, known as giant resonances [79,80], also begin to play a major role. Both of these details are entirely missing from the conventional Fermi gas model of the nuclear ground state. While these deficiencies of the Fermi gas model become less problematic as the neutrino energy increases to several hundred MeV and beyond, it has been pointed out that low-energy nuclear excitations are still expected to exert a noticeable influence on ∼1 GeV neutrino-nucleus differential cross sections at very forward scattering angles [81,82].
A second modeling difference that limits the suitability of typical neutrino event generators for the low-energy regime is the description of hadronic final-state interactions. Rather than taking an INC-or BUU-like dynamical approach to FSIs, in which intranuclear scattering is explicitly modeled, typical low-energy calculations [60,[83][84][85][86][87][88][89][90] opt instead for a statistical treatment, in which only bulk properties of the nuclear system (such as its excitation energy and spin-parity) are employed to predict its behavior. The latter strategy is usually justified by assuming that the struck nucleon from the primary neutrino interaction will scatter repeatedly within the nuclear medium without being directly knocked out. These multiple intranuclear collisions lead to the transferred energy being widely shared among the constituent nucleons and thus to thermal equilibration: de-excitations of the resultant compound nucleus may be treated independently of the manner in which it was formed.
While theoretical predictions of neutrino-nucleus cross sections at high (low) energies tend to rely exclusively on a dynamical (statistical) model of FSIs, a more complete treatment may be achieved by combining the two approaches. The current release of GiBUU enables such calculations by providing an optional interface to the Statistical Multifragmentation Model (SMM) [91] code. Disintegration of the residual nucleus is simulated by SMM in a post-processing step which takes otherwise complete GiBUU events as input. Although the two codes have been used together to examine other processes [92][93][94][95], their joint application to neutrino interactions remains unstudied.
At present, native support for statistical nuclear de-excitation models in the other three neutrino generators mentioned above is limited to some simple approximations used in GENIE's treatment of nucleon emission [96]. However, official GENIE interfaces to INCL++ [97,98] and to the Bertini Cascade implementation [99] in Geant4 [100,101] are in the late stages of development [102,103]. These will provide enhancements to GENIE FSI modeling which are similar in scope to those obtained with SMM for GiBUU. Unofficial interfaces are also being explored by outside groups, with at least one attempt [104] having been made to apply the TALYS [105,106] nuclear de-excitation model to events generated using both GENIE and NuWro.
Although interfacing with these external tools provides an FSI treatment more compatible with standard low-energy approaches, a key omission remains problematic. With the exception of TALYS, the usual configurations of all of the remaining codes 5 lack a means of simulating discrete γ-ray transitions between low-lying nuclear energy levels. 6 Gamma-rays created in this way represent a major component of the final state for inelastic scattering of solar neutrinos on complex nuclei. For inelastic NC reactions, de-excitation γ-rays may often be the only final-state particles which are experimentally observable.
In this work, I present a new neutrino event generator, MARLEY, 7 which implements a model of inelastic neutrino-nucleus scattering designed specifically for the low-energy regime. The early version of MARLEY described herein is primarily focused on simulations of charged-current absorption of ν e on 40 Ar [108]. However, preparation of additional input data would allow other reaction channels and nuclear targets to be handled without difficulty within the existing code framework. Section 2 provides an overview of the theoretical treatment of neutrino scattering used in MARLEY. Section 3 presents the MARLEY approach to random sampling, which makes ample use of modern improvements to the C++ language and, as discussed in section 3.3, implements a new inverse transform sampling algorithm [109] in a physics event generator for the first time. Section 4 discusses the MARLEY event generation workflow and implementation details. Sections 5-7 outline how to install, configure, run, and interpret the output of the code. Section 8 explains how MARLEY can be interfaced with external software toolkits, using the popular Geant4 [100,101] particle transport package as an example. Finally, section 9 considers prospects for future improvements to MARLEY.

MARLEY treatment of neutrino scattering
This section provides a brief overview of the physics models currently implemented in MARLEY, with more details available in ref. [108]. Natural units ( = c = 1) are used throughout.
As has been done in many previous calculations of low-energy neutrino cross sections [60,[83][84][85][86][87][88][89][90], MARLEY treats neutrino-nucleus scattering events as proceeding via a two-step process. In the first step, a two-to-two scattering reaction involving the neutrino and the target nucleus is simulated, and the final nucleus is left in a state with a well-defined excitation energy, spin, and parity. In the second step, which is handled independently from the first, the final nucleus de-excites. At low excitation energies, where discrete level data are available, the de-excitations are modeled using tabulated γ-ray branching ratios. At higher excitation energies, where the final nucleus becomes unbound, the formation of a compound nucleus is assumed, and the decay widths for all open channels are calculated using the Hauser-Feshbach statistical model [110].

Nuclear 2 → 2 scattering model
To simulate two-to-two neutrino-nucleus scattering processes at low energies, MARLEY 1.2.0 evaluates the nuclear matrix elements in the long-wavelength limit (in which the four-momentum transfer q → 0) and the slow-nucleon limit (in which |p N |/m N → 0, where p N is the initial 3-momentum of the struck nucleon and m N is its mass). The combination of these two limits is sometimes referred to as the allowed approximation. Under this approach, the differential cross section in the center-of-momentum (CM) frame for a transition to a particular nuclear final state is given by the expression Here G F is the Fermi constant, Mandelstam s is the square of the total energy in the CM frame, and E i (E f ) is the total energy of the initial (final) nucleus. The final-state lepton has total energy E , 3-momentum p , speed β = |p |/E , and scattering angle θ , which is defined with respect to the incident neutrino direction. The first factor in square brackets, E i E f /s, arises due to nuclear recoil and is often neglected. The symbol F CC is used to introduce extra factors needed when computing the cross section for charged-current scattering. Discussion of F CC is deferred to section 2.1.2. The spin-reduced Fermi and Gamow-Teller nuclear matrix elements may be written in the form where g V (g A ) is the vector (axial-vector) weak coupling constant of the nucleon, and J i (J f ) is the initial (final) nuclear spin. The Fermi matrix element B(F) is subject to the spin-parity selection rule where Π i (Π f ) is the initial (final) nuclear parity. The Gamow-Teller matrix element B(GT) likewise obeys the selection rule The Fermi (O F ) and Gamow-Teller (O GT ) operators from eqs. (6) and (7) are defined for charged-current and neutral-current scattering processes via the relations where σ is the Pauli vector, A is the nucleon number, and the isospin lowering (raising) operator t − (t + ) should be chosen 8 for an incident neutrino (antineutrino). The symbol t 3 denotes the third component of isospin, and an operator suffixed by (n) is understood to act only on the nth nucleon. The weak nuclear charge Q W for a nucleus with neutron number N and proton number Z is given in terms of the weak mixing angle θ W by

Coherent elastic neutrino-nucleus scattering
For NC reactions on a spin-zero (J i = 0) target nucleus, the differential cross section from eq. (5) reduces to a particularly simple form when describing scattering that leaves the nucleus in its ground state. The Gamow-Teller selection rule (eq. (9)) ensures that B(GT) vanishes identically, while the Fermi matrix element for the ground-state-to-ground-state transition becomes For this special case, eq. (5) may be rewritten in terms of the lab-frame kinetic energy T f of the recoiling ground-state nucleus as where M is the mass of the nuclear target. The maximum value of T f allowed by the reaction kinematics is given in terms of the initial lab-frame neutrino energy E ν by This process is known as coherent elastic neutrino-nucleus scattering (CEvNS) [38,39]. The nuclear recoil correction factor is usually neglected 9 in the CEvNS literature (see, e.g., [111]). Due to its adoption of the allowed approximation, the MARLEY treatment of CEvNS currently does not include a q 2 -dependent nuclear form factor that accounts for imperfect coherence in the cross section [112].

Coulomb corrections
In eq. (5), the symbol F CC is used to include extra factors needed solely for CC scattering. It is defined by where V ud is the Cabibbo-Kobayashi-Maskawa matrix element connecting the up and down quarks. The Coulomb correction factor F C accounts for the electromagnetic interaction between the outgoing lepton and nucleus in an approximate way. Three prescriptions for computing F C are available in MARLEY: the Fermi function, the effective momentum approximation (EMA), and the modified effective momentum approximation (MEMA).
Fermi function. For very low energies of the outgoing lepton (such as those observed in beta decay), using the Fermi function [113,114] as the Coulomb correction factor is a standard approach. A minor complication emerges, however, because the derivation of the Fermi function 10 assumes that the final nucleus is at rest, while the differential cross section in eq. (5) is evaluated in the CM frame and accounts for nuclear recoil. To work around this discrepancy, MARLEY uses the relative speed β rel of the two final-state particles [116] to evaluate the Fermi function in the rest frame of the final nucleus using the Lorentz-invariant expression In eqs. (17) and (18), p , m f , and Z f denote the 4-momentum, mass, and proton number of the final nucleus. Likewise, k , m , and z represent the 4-momentum, mass, and electric charge (in units of the elementary charge) of the final-state lepton. The quantity S is defined in terms of the fine structure constant α by is the nuclear radius (in natural units), and the Sommerfeld parameter η is given by Effective momentum approximation. For higher energies of the outgoing lepton, the Fermi function is known to overestimate the magnitude of the Coulomb corrections, and an alternative approach called the effective momentum approximation (EMA) becomes more appropriate [117]. Let the symbol K (E) denote the momentum (total energy) of the outgoing lepton in the rest frame of the final nucleus: Then the effective values of these variables are those that exist in the presence of the nuclear Coulomb potential, which is taken to be that at the center of a uniformly-charged sphere: In the MARLEY implementation of the EMA, the Coulomb correction factor is given by the ratio 11 Modified effective momentum approximation. In ref. [117], an adjustment to the standard EMA prescription is proposed which improves the accuracy of the approximation in cases where the final lepton mass cannot be neglected. Under this modified effective momentum approximation (MEMA), the Coulomb correction factor defined in eq. (25) is replaced by Default behavior. The final-state Coulomb interaction increases the charged-current cross section for neutrinos and decreases it for antineutrinos. Because the (M)EMA is known to overestimate the size of this effect at low energies while the Fermi function does the same at high energies, previous calculations [63,118] have adopted a simple prescription for combining the two approaches: in any particular case, adopt the method that yields the smallest Coulomb correction. For MARLEY, this amounts to defining the Coulomb correction factor F C by Although eq. (27) represents the default MARLEY approach to Coulomb corrections, this behavior may be altered by the user in the job configuration file (see section 6.7.5).

Total cross section
With the definitions given above, integration of eq. (5) over cos θ becomes trivial, leading to the total cross section

Nuclear de-excitation model
De-excitations from high-lying nuclear levels are simulated in MARLEY using the Hauser-Feshbach statistical model (HFSM) [110]. This model assumes that the decaying nuclear state may be adequately described as a thermally-equilibrated compound nucleus with a definite excitation energy (E x ), spin (J), and parity (Π). In the MARLEY HFSM implementation, emissions of γ-rays and light nuclear fragments (1 ≤ A ≤ 4) are treated as a sequence of binary decays while fission, production of heavy nucleon clusters (A ≥ 5), and simultaneous multiparticle evaporation are neglected.

Nuclear fragment emission
According to the HFSM, the distribution of final excitation energies E x that may result from the emission of a fragment a (with parity π a and separation energy S a ) from the compound nucleus is described by the differential decay width where J is the final nuclear spin; s, , and j are the spin, orbital, and total angular momentum quantum numbers of the emitted fragment; is the value of the final-state nuclear parity needed to enforce parity conservation; and ε is the total kinetic energy of the decay products in the rest frame of the initial nucleus. The maximum accessible final-state excitation energy is related to the total kinetic energy ε via The functions ρ i and ρ f represent the density of nuclear levels in the vicinity of the initial and final states, while the transmission coefficient T j quantifies how readily the fragment may be emitted from the nucleus. Because the value of T j at fixed ε falls off rapidly with increasing , MARLEY truncates the infinite sum over orbital angular momenta using an upper limit max . By default, max = 5 is used. This value may be adjusted by the user as described in section 6.7.4. At low excitation energies, where individual nuclear levels can be resolved, MARLEY treats the level density ρ f as a sum of delta functions, with one term per level. For a specific nuclear level, the partial decay width for emission of fragment a may be written in the form where the symbol δ π , which enforces parity conservation, is equal to one if eq. (30) is satisfied and zero if it is not. At higher excitation energies, MARLEY computes ρ f according to the Back-shifted Fermi gas model (BFM) from version 3 of the Reference Input Parameter Library (RIPL-3) [119]. The "BFM effective" values of the level density parameters for this model are adopted from a global fit of nuclear level data for 289 nuclides reported in ref. [120]. A full description of the BFM as implemented in MARLEY is given in appendix B of ref. [108]. The initial level density ρ i is always evaluated according to the BFM regardless of excitation energy.
The partial decay width for a transition to the continuum of nuclear levels via emission of a fragment a may be computed by integrating eq. (29) over the final excitation energy interval . The lower bound of the continuum E min,c x is taken by MARLEY to be the excitation energy of the highest tabulated discrete nuclear level for the nuclide of interest. In cases where no such level data are available, the continuum is taken to start at the nuclear ground state (E min,c x = 0).

Fragment transmission coefficients
The fragment transmission coefficients T j that appear in eqs. (29) and (33) are computed by numerically solving the radial Schrödinger equation where u j is the fragment's radial wavefunction, is its kinetic energy in the rest frame of the final nucleus, 12 and is the magnitude of its 3-momentum in the rest frame of the initial nucleus. In eqs. (35) and (36), m a (M ) denotes the mass of the emitted fragment (final nucleus). The optical potential U used by MARLEY for nucleon emission is the global parameterization of Koning and Delaroche [121]. For complex nuclear fragments, a folding approach similar to that of Madland [122] is used to construct the optical potential by weighting the individual neutron and proton potentials. More details about the MARLEY nuclear optical potential are available in appendix C of ref. [108].
Far from the nucleus, the optical potential approaches the Coulomb potential, and the fragment radial wavefunction approaches the limiting form 12 The label lab for this quantity reflects it status as the laboratory-frame kinetic energy in the time-reversed process wherein the fragment is absorbed to form the compound nucleus. See appendix A of ref. [108].
where H ± are the Coulomb wavefunctions [123, ch. 33]. The Sommerfeld parameter is evaluated in terms of the proton number z (Z ) of the emitted fragment (final nucleus) and the relative speed The energy-averaged S-matrix element S j that appears in eq. (37) is related to the transmission coefficient To approximate S j , MARLEY first obtains a numerical solution u j (r) of eq. (34) using Numerov's method [124][125][126]. This method computes u j iteratively on a regular grid with fixed radial step size ∆. If one defines the function a j (r) to be equal to the non-derivative terms enclosed in square brackets in eq. (34), i.e., then the Numerov solution (accurate to order ∆ 4 ) u n ≈ u j (r n ) at the nth grid point is given by the recurrence relation and boundary conditions u 0 = 0 (44) Here I have defined the abbreviations h ≡ ∆ 2 /12 and a n ≡ a j (r n ) for n ≥ 1, with a 0 ≡ 0. The MARLEY calculation of S j proceeds through iterations of the Numerov method until the first grid point r A is encountered such that that is, the difference between the nuclear optical potential U and the Coulomb potential V C falls below a small threshold V thresh . Iterations continue further until a second grid point r B is encountered such that In MARLEY 1.2.0, the parameter values ∆ = 0.1 fm / c, V thresh = 1 keV and S = 1.2 are used.
Comparing the full solution obtained in this way to the asymptotic form from eq. (37) at r A and r B leads to the expression Numerical values of the Coulomb wavefunctions are obtained by interfacing with the GNU Scientific Library [1,2].

Gamma-ray emission
In the Hauser-Feshbach formalism, γ-ray emission is described by the differential decay width where λ ≥ 1 is the multipolarity, E γ ≈ E x − E x is the energy of the emitted γ-ray, 13 and labels the type of transition as either electric (E) or magnetic (M). The infinite sum over multipolarities in eq. (49) is truncated at λ max . The default cutoff value λ max = 5 may be configured by the user as described in section 6.7.4. The calculation of the level densities ρ i and ρ f is identical to the approach used for nuclear fragment emission (see section 2.2.1). In particular, the level density ρ f used for γ-ray transitions to discrete levels is once again treated as a sum of delta functions. The partial decay width for γ-ray emission to a particular nuclear level then becomes If J + J < 1, the width Γ γ vanishes.
Similarly to the fragment emission case, calculation of the partial decay width for γ-ray transitions to the continuum of nuclear levels is performed by integrating eq.

Gamma-ray transmission coefficients
The transmission coefficients T Xλ that appear in eqs. (49) and (51) are typically expressed in terms of a strength function f Xλ (E γ ) such that In MARLEY 1.2.0, the expression used to evaluate the strength function is taken from the RIPL-3 Standard Lorentzian model [119]. According to this model, γ-ray emissions of type Xλ are assumed to take place via de-excitation of the corresponding giant multipole resonance, which has centroid excitation energy E Xλ , width Γ Xλ , and peak cross section σ Xλ . A full listing of the values of these parameters is given in Table II of ref. [108]. Neutrino

Neutrino-electron elastic scattering
In addition to neutrino-nucleus scattering, MARLEY is also capable of simulating neutrino-electron elastic scattering. For a target atom with proton number Z, the differential cross section in the CM frame for this process is computed according to where m e is the electron mass and E ν (θ ν ) is the energy (scattering angle) of the neutrino. The coupling constants g 1 and g 2 depend on the neutrino species and are given in table 1. Electron binding energies are neglected.

Random sampling implementation
Like any other Monte Carlo event generator, MARLEY must make extensive use of pseudorandom numbers and produce samples from a variety of discrete and continuous probability distributions. With the advent of C++11, a suite of high-quality random number generation tools were adopted as part of the C++ standard library [127], and MARLEY relies heavily on these new features. All random numbers used by MARLEY are obtained using the C++ standard library object std::mt19937 64, which provides a 64-bit implementation of the Mersenne Twister algorithm developed by M. Matsumoto and T. Nishimura [128,129].

Discrete distributions
All sampling from discrete distributions is handled using instances of the C++ standard library object std::discrete distribution. In cases where the sampling weights for each possible outcome already exist in memory, the usual method for initializing this object is to supply iterators that point to the beginning and end of a collection of sampling weights. For example, line 5 of the code snippet 14 1 # include < random > 2 # include < vector > 3 4 std :: vector < double > weights = { 1. , 2. }; 5 std :: d iscr ete_d istr ibut ion dist1 ( weights . begin () , weights . end () ) ; initializes a std::discrete distribution object dist1 which will sample int values of 0 (1) with probability 1/3 (2/3).
When the sampling weights are stored as object data members, however, initializing the distribution becomes more complicated. In particular, a naïve attempt using the approach from the example above  triggers a compilation error on line 12. Because the iterators returned by As.begin() and As.end() refer to objects of type A instead of the weights (of type double), they cannot be used to initialize dist2.
MARLEY works around this difficulty by implementing the iterator to member interface proposed by T. Becker [130]. This interface converts an iterator that points to an object (A) into an iterator that points to one of that object's data members (weight). In this example, including the appropriate MARLEY header file (include/marley/IteratorToMember.hh) and replacing line 12 with compiles successfully and yields the same sampling behavior for dist2 as for dist. For cases in which the input iterators refer to object pointers instead of the objects themsleves, MARLEY provides a similar interface via the marley::IteratorToPointerMember class template, which is compatible with both bare (e.g., A*) and smart pointers (e.g., std::unique ptr<A>).
This approach is used in several places in the MARLEY source code to initialize discrete distributions using sampling weights stored as object data members, e.g., relative intensities owned by marley::Gamma objects representing distinct γ-ray de-excitations from a particular nuclear energy level.

Continuous 1D distributions: accept/reject approach
To sample from continuous one-dimensional distributions, MARLEY implements two general schemes, both of which take the bounds of a sampling interval [a, b] and an arbitrary probability density function f (x) (for which no particular normalization is assumed) as input. The first scheme uses a simple rejection method which relies on an accurate knowledge of the global maximum f max of f (x) within the sampling interval. If f max is known in advance, it may be supplied along with the other input parameters. Otherwise, it is estimated within a specified tolerance by minimizing −f (x) using Brent's method [131]. Once the value of f max has been obtained, pairs of uniformly-distributed variables x ∈ [a, b] and y ∈ [0, f max ] are repeatedly sampled. This continues until y ≤ f (x), at which point the sampled x value is accepted.

Continuous 1D distributions: inverse transform approach
In cases where the global maximum f max is unknown and its estimation via Brent's method is either unreliable (because a local maximum may be found instead of the global one) or inefficient (because f is computationally expensive to evaluate), MARLEY employs a second sampling scheme based on the "fast inverse transform sampling" algorithm originally proposed in ref. [109]. Although MARLEY uses numerical techniques similar to those in the original Matlab code [132], which was written as an extension of the Chebfun package [133,134], the C++ implementation described herein is original. To the author's knowledge, MARLEY represents the first application of this algorithm to physics event generation.

Algorithm
To obtain a random sample x from f (x) using inverse transform sampling, MARLEY first constructs a polynomial approximantf of f using a grid of N + 1 ordered pairs (x j , f j ) for j ∈ {0, 1, . . . , N }, where the x j are the Chebyshev points of the second kind and the f j are the function values At points other than the x j (where f j is used directly), the polynomial approximantf is given by the barycentric formula [135] f with weights The f j are related to the coefficients α k that appear in an N th order expansion of f in Chebyshev polynomials T k , i.e., by the type-I discrete cosine transform (DCT-I) To approximate the cumulative density function (CDF) MARLEY uses the formulas to integrate eq. (60) term-by-term, yielding an (N + 1)th order Chebyshev expansion where the coefficients β k are given by One may obtain a polynomial approximantF of the cumulative density function F by applying a second DCT-I (which is its own inverse) to the expansion coefficients β k : Using the N + 1 Chebyshev points one may approximate F (x) for any x ∈ [a, b] using eq. (58) with the substitutions f → F , j → , and N → N + 1.
With the approximate CDFF (x) constructed in this manner, MARLEY obtains a random sample x from f (x) by generating a uniform random number ξ ∈ [0, 1]. Bisection is then used to find the sampled value of x ∈ [a, b] which satisfies the relation Figure 1 shows an application of this sampling technique to the modeling of compound nuclear decays in MARLEY. In the left-hand plot, the black histogram shows the spectrum of proton kinetic energies T p obtained from a simulation of 200,000 decays 40 K → p + 39 Ar to the excitation energy continuum of the daughter 39 Ar nucleus. The initial 40 K state had excitation energy E x = 45 MeV and spin-parity J π = 4 + . Decay modes other than proton emission to the continuum were switched off for simplicity, and the kinetic energy T p is reported in the rest frame of the mother 40 K nucleus. The contents of each histogram bin are normalized to provide a calculation of the differential decay width dΓ p /dT p , where the average value of this quantity in the wth bin is approximated by the Monte Carlo estimator

Validation
Here N events = 200,000 is the total number of simulated decay events, n w is the number of these that fall within the wth bin, ∆T w p is the width of the wth bin, and Γ p is the total width for proton emission to the continuum.
The blue curve shows a direct calculation of the differential decay width via where M (M ) is the mass of the initial (final) nucleus and dΓ p /dE x is evaluated according to eq. (29) with the fragment species a = p. The agreement seen between the exact calculation and the simulated events is achieved by sampling the latter from a CDF constructed using a Chebyshev polynomial approximant to dΓ p /dE x with grid size N = 64. The polynomial approximant with N = 64 and the exact calculation are indistinguishable on the scale of the plot. For reference, a lower-order (N = 4) Chebyshev polynomial approximation to the distribution is also shown by the dotted purple line. The right-hand plot in fig. 1 shows simulation results obtained using an identical procedure, except that the decay process 40 K → γ + 40 K is considered using the differential width from eq. (49). Three lower-order Chebyshev approximants are shown which provide improved agreement with the exact calculation as the grid size N grows.

Implementation details
Although the barycentric interpolation scheme described here is used in MARLEY solely for inverse transform sampling, the C++ implementation is very general and may find useful applications elsewhere. The marley::ChebyshevInterpolatingFunction class constructs the polynomial interpolantf for an arbitrary input function f (x), represented by a std::function<double(double)>. If the grid size N is not specified in the constructor, then an adaptive technique is used to choose a grid size sufficiently large to represent f (x) at close to machine precision. Starting with N = 2, the value of N is doubled (and the Chebyshev expansion coefficients are recomputed) until the stopping criterion (with machine epsilon ε) is satisfied or N reaches a large maximum value. The evaluate(double x) member function returns the value off (x), and the cdf() member function returns a new marley ::ChebyshevInterpolatingFunction representingF (x). The DCT-I calculations in eqs. (61) and (67) are carried out using the fast Fourier transform C library FFTPACK4 [136], which is included in the MARLEY source code distribution. This library is based on the original FFTPACK Fortran code developed by P. Swarztrauber [137,138].
No simultaneous sampling of multiple variables from a joint probability distribution is needed to implement the physics models in the current version of MARLEY.

Event generation workflow
The flowchart in fig. 2 illustrates the procedure used by MARLEY to generate events. In the following paragraphs, each stage in the process will be described. Unless otherwise noted by providing an explicit namespace specifier (e.g., std::), all C++ classes referred to using typewriter font in this section are defined within the marley namespace.

Generator initialization
At the beginning of an event generation job, a JSONConfig object is used to parse and interpret the settings stored in a configuration file (whose format is described in section 6). The create generator member function is then used to construct a Generator object which handles the actual simulation of events. The steps below are followed to initialize the Generator.

Set random number seed
If the user has specified an integer value for the random number seed in the configuration file (see section 6.1), then this value is used to initialize a std::mt19937 64 object owned by the Generator. If not, then the system time since the Unix epoch is used as a seed.

Set neutrino direction
By default, MARLEY generates events in a reference frame in which the incident neutrino is traveling in the positive z direction. If the user has specified a different neutrino direction (see section 6.2), then this information is passed to a ProjectileDirectionRotator object owned by the Generator. A rotation matrix is precalculated which will allow the coordinate system to be appropriately transformed at the end of each iteration of the event loop.

Define incident neutrino spectrum
Based on the user's description of the neutrino energy spectrum (see section 6.5), one of several derived classes of the NeutrinoSource abstract base class is instantiated and stored as a member of the Generator object. In typical MARLEY use cases, the NeutrinoSource describes the energy distribution of incident neutrinos. This behavior is desirable so that the energy spectrum P (E ν ) of reacting neutrinos where φ(E ν ) is the NeutrinoSource spectrum and σ(E ν ) is the total cross section, is fully consistent with the MARLEY physics models. In unusual situations where cross section weighting is not wanted (e.g., event generation with a uniform distribution of reacting neutrino energies), it may be disabled in the job configuration file as described in section 6.7.3. This corresponds to the substitution σ(E ν ) → 1 in eq. (73).

Configure reactions
Each distinct 2 → 2 scattering mode that may be simulated during a MARLEY job is represented by an object that implements the abstract base class Reaction. This class includes member functions called total xs and diff xs, which respectively compute (in units of MeV −2 per target atom) the reaction total cross section σ r (E ν ) and differential cross section d σ r (E ν )/d cos θ , where θ is the lepton scattering angle in the CM frame. The create event member function simulates the scattering process using these quantities (see sections 4.2.3 and 4.2.4) and returns the results in a newly-created Event object.
Based on the set of one or more scattering modes enabled by the user in the configuration file (see section 6.4), a vector of pointers to Reaction objects is initialized during job startup and stored as the reactions member of the Generator. This vector is populated by calls to the factory method Reaction::load from file, which processes a reaction specification given in an input file (see section 6.4 and appendix B). After all reaction input files have been fully processed, the abundance-weighted sum σ(E ν ) of the total cross sections σ r (E ν ) for all of the enabled reactions is used together with the NeutrinoSource spectrum φ(E ν ) to construct the probability density function for reacting neutrino energies P (E ν ) shown in eq. (73). The weights needed to compute σ(E ν ) are the nuclide fractions in the neutrino target, as described in section 4.1.6.
In MARLEY 1.2.0, two concrete derived classes of Reaction are implemented. The Nuclear Reaction class implements the neutrino-nucleus scattering model described in section 2.1, while the ElectronReaction class does the same for the neutrino-electron elastic scattering model from section 2.3.
While simulating a neutrino-nucleus scattering event, MARLEY represents the nuclear state in terms of the following quantities: the proton number Z, nucleon number A, excitation energy above the ground state E x , total spin J, and parity Π. The values of these variables are tracked throughout the event loop, both during 2 → 2 scattering (section 4.2) and during simulation of nuclear de-excitations (section 4.3).
A distinction is made in MARLEY between bound versus unbound nuclear states. A nuclear state is considered to be unbound if the excitation energy E x exceeds the separation energy for at least one nuclear fragment with mass number A ≤ 4. Separation energies are computed using atomic and particle mass data from refs. [139,140] tabulated in the file data/mass table.js. The singleton class MassTable provides a C++ API for accessing and manipulating these data. If an atomic mass value is not tabulated for a particular nuclide, the MassTable class computes an estimate using a formula for the liquid-drop model mass excess developed by Myers and Swiatecki [120,141].
A set of transition matrix elements to the final nuclear states that may be populated by means of a 2 → 2 neutrino-nucleus scatter are represented in MARLEY by a vector of MatrixElement objects. Each of these objects is labeled as either a Fermi or Gamow-Teller matrix element (see section 2.1) using the type member variable, while the strength member stores the corresponding B(F) or B(GT) value. These quantities are given in the reaction input file together with the final nuclear excitation energy, which is stored as the MatrixElement member level energy . A Nuclear Reaction is constructed using a std::shared ptr to a vector of MatrixElement objects, allowing shared ownership between multiple reactions as appropriate.
While parsing an input file repesenting a neutrino-nucleus reaction, the static method Reaction ::load from file constructs one MatrixElement object for each listed nuclear transition. After the full list of transitions has been processed, a StructureDatabase object owned by the Generator is consulted to determine whether discrete level data are available for the final-state nucleus of interest. The first such query will trigger loading of these data in the form of DecayScheme objects as described in section 4.1.5. If a suitable DecayScheme is available, then every MatrixElement which represents a transition to a bound nuclear final state is matched to a Level object owned by the DecayScheme. This matching is performed by selecting the Level whose excitation energy most closely matches the value of the MatrixElement's member variable level energy . Duplicate matches occurring for two MatrixElement objects belonging to the same NuclearReaction will lead to an exception being thrown. A successful match will result in the MatrixElement's member variable final level being loaded with a pointer to the matched Level. If the spin-parity of the matched Level is not compatible with the relevant selection rules (from either eq. (8) or eq. (9)), then a warning message will be printed to the screen.
Matching of this kind is not performed for MatrixElement objects representing transitions to unbound nuclear states. In such cases (or when an appropriate DecayScheme is not available), a null pointer is assigned to the final level member.

Load nuclear structure data
To simulate nuclear transitions at low excitation energies, MARLEY relies on tabulated nuclear structure data files. These files contain listings of discrete nuclear energy levels and branching ratios for γ-ray transitions between them. The file format is documented in appendix A.
For MARLEY 1.2.0, the recommended structure data files are largely taken (with reformatting) from version 1.6 of the TALYS nuclear reaction code [105,106]. The data are organized by element in the folder data/structure/. This folder also includes a text file named nuclide index.txt which serves as an index for the entire dataset: each line of the text file has the format NucPDG DataFileName where NucPDG is an integer PDG code (see section 7.1) representing a nuclear species and DataFile Name is the name of the corresponding file (assumed to appear in data/structure/) in which its discrete level data are given.
A StructureDatabase object owned by the Generator provides an interface for accessing and manipulating nuclear structure data within MARLEY. The StructureDatabase manages a lookup table indexed by nuclear PDG code that stores pointers to DecayScheme objects. Each Decay Scheme represents a table of nuclear levels and γ-ray branching ratios for a particular nuclide. A fully-initialized DecayScheme owns one Level object for every discrete level listed for the nuclide of interest in the associated structure data file. Each Level itself owns one Gamma object for each listed γ-ray transition that may originate from it.
External access to the DecayScheme objects is provided by the StructureDatabase member function get decay scheme. This function takes a nuclear PDG code or a (Z, A) pair as input. If a corresponding DecayScheme already exists in the lookup table, a pointer to it is immediately returned. If a suitable DecayScheme has not been constructed yet, the StructureDatabase consults the nuclide index (which is automatically loaded from nuclide index.txt when needed and cached for repeated use) in an attempt to find a matching structure data file. If a match is found, Decay Scheme objects are constructed and added to the lookup table for every nuclide listed in the file. A pointer to the requested DecayScheme is then returned if it was successfully loaded. A null pointer is stored in the lookup table and returned when the attempt to load the requested DecayScheme fails.
Additional information indexed by nuclear PDG code is also stored in the StructureDatabase, including (1) ground-state nuclear spin-parities 15 loaded from the file data/structure/gs spin parity table.txt and (2) three varieties of objects (represented by the abstract base classes Level DensityModel, OpticalModel, and GammaStrengthFunctionModel) used to compute quantities of interest for the Hauser-Feshbach statistical model (see sections 2.2 and 4.3.1). In MARLEY 1.2.0, these three abstract base classes each have a single concrete implementation which is used for calculations. 16 The BackshiftedFermiGasModel class is derived from LevelDensityModel and computes nuclear level densities (e.g., ρ i and ρ f in eq. (29)) as described in section 2.2.1. The KoningDelarocheOpticalModel class is derived from OpticalModel and computes nuclear fragment transmission coefficients using the procedure outlined in section 2.2.2. Finally, the Standard LorentzianModel class is derived from GammaStrengthFunctionModel and computes the γ-ray strength function from eq. (54) as described in section 2.2.4.
Beyond the nuclide-specific items managed by the StructureDatabase, there are also two settings which are common to all nuclei: the cutoff values max and λ max used to truncate the sums in eq. (29) and eq. (49), respectively, (see sections 2.2.1, 2.2.3 and 6.7.4) and a list (elements of which are represented by the Fragment class) of the nuclear fragments to consider when simulating decays using the Hauser-Feshbach model.

Define target
In a MARLEY simulation, the isotopic composition of the material illuminated by the incident neutrinos is represented by a Target object owned by the Generator. In the absence of an explicit list of abundances specified by the user (see section 6.3), MARLEY assumes equal amounts of each nuclide that appears in the initial state of at least one configured Reaction. Using the information stored in the Target object, the Generator computes the abundance-weighted total cross section per target atom via where f r is the nuclide fraction for the initial-state atom involved in the rth Reaction.

Enable/disable nuclear de-excitations
By default, MARLEY simulates both the prompt 2 → 2 scattering reaction and the subsequent nuclear de-excitations for every event. In applications where only the prompt reaction is important, the user may disable de-excitations as described in section 6.7.2. The choice of whether or not to simulate de-excitations is represented as a boolean member variable owned by the Generator object.

Event loop: 2 → 2 scattering
After the Generator object has been fully initialized, the marley command-line executable enters an event loop. This loop iterates until a fixed number of events requested by the user (see section 6.6) has been reached or the loop is interrupted (by an error condition or by the user pressing ctrl+C). A single event loop iteration corresponds to a call to the Generator member function create event. This function returns an Event object constructed according to the following steps.

Neutrino energy selection
In the first step of the event loop, the Generator employs rejection sampling (see section 3.2) to select a reacting neutrino energy from the probability distribution P (E ν ) defined in eq. (73). By default, Brent's method is used during the first event loop iteration to obtain a numerical estimate P max of the maximum value of the probability density function P (E ν ). The value of P max is cached and reused during neutrino energy sampling for subsequent events.
In cases where Brent's method fails to converge to the global maximum of P (E ν ), the resulting distribution of E ν in the generated events will be biased, with too few events produced in energy regions where P (E ν ) is larger than P max . At the start of each pass through the event loop, MARLEY verifies that all values of P (E ν ) computed during rejection sampling never exceed P max . If the cached value of P max is ever found to be an underestimate of the true global maximum, an error message is printed (see listing 2 in section 5.4.1 for an example), and the value of P max is updated to match the largest value of P (E ν ) encountered.
Users can override the default use of Brent's method by supplying their own values of P max in the job configuration file. Instructions for doing so can be found in section 6.7.7. For convenience in troubleshooting, the error message printed by MARLEY in response to a rejection sampling problem contains a recommended value of P max to use in this way.

Reaction mode sampling
Once a reacting neutrino energy E ν has been chosen, the Generator samples a specific reaction mode r from the discrete probability distribution Control then passes to the create event member function of the rth Reaction object owned by the Generator. This function handles selection of 2 → 2 scattering kinematics and initializes an Event object.

Selection of a nuclear transition
If the sampled reaction mode involves scattering on a nuclear target, then simulation of the prompt 2 → 2 reaction is handled by the NuclearReaction class, and the outgoing nucleus may be left in one of potentially many final states. A specific final state is chosen for a nuclear reaction by selecting a MatrixElement according to the discrete distribution where σ r (E ν , Λ) is the partial reaction cross section computed for the Λth MatrixElement owned by the NuclearReaction of interest. That is, σ r (E ν , Λ) is given by the expression in eq. (28) with the quantity B(F) + B(GT) set equal to the strength member of the Λth MatrixElement. The total cross section σ r (E ν ) is obtained by summing over all of the owned MatrixElement objects which represent energetically-accessible nuclear transitions: If the sampled MatrixElement has been matched to a known discrete nuclear Level (i.e., if its member pointer final level is non-null, see section 4.1.4), then the excitation energy, spin, and parity of the final nucleus are determined directly from the Level object.
If a match is not available (e.g., because the associated nuclear state is unbound), then the final nuclear excitation energy E x is set equal to the sampled MatrixElement's member variable level energy . The final-state spin J and parity Π are determined according to the selection rules given in either eq. (8) or eq. (9) as appropriate for the MatrixElement type. In the case of a Gamow-Teller transition involving an initial nucleus with nonzero spin (as determined using the table of ground-state spin-parities managed by the StructureDatabase, see section 4.1.5), multiple J values will satisfy the spin selection rule in eq. (9). In such cases, MARLEY 1.2.0 makes a rough approximation: equipartition of spin is assumed, and a definite value of J is sampled from the discrete distribution where ρ (which is calculated using the same level density treatment as in section 2.2.1) is the density of final-state nuclear levels with spin J and parity Π in the vicinity of excitation energy E x . In eq. (78), the sum in the denominator runs over all final spin values K that satisfy the Gamow-Teller spin selection rule from eq. (9).

Outgoing lepton direction
Due to the simple kinematics of 2 → 2 scattering, the 4-momenta of the outgoing particles are fully determined by specifying their masses (including the final nuclear excitation energy) and the direction of the outgoing lepton. This direction is represented in MARLEY using the lepton's azimuthal scattering angle φ and polar scattering cosine cos θ , both measured with respect to the incident neutrino direction in the center-of-momentum (CM) frame.
For all reactions currently implemented in MARLEY 1.2.0, the differential cross section is independent of φ , and the value of this variable is therefore sampled uniformly on the half-open interval [0, 2π).
For neutrino-nucleus reactions, a value of cos θ is obtained via rejection sampling (see section 3.2) on the interval [−1, 1] from the probability density function where d σ r (E ν , Λ) / d cos θ is the CM-frame differential cross section for the reaction r proceeding via a nuclear transition described by the Λth MatrixElement. This differential cross section is computed as in eq. (5), except that only one of the two matrix element terms is used. For a Fermi transition, B(F) is set equal to the strength member variable of the MatrixElement, while B(GT) is set to zero. For Gamow-Teller transitions, the reverse is done. In the case of neutrino-electron elastic scattering, rejection sampling is still performed using the distribution from eq. (79), but the differential cross section is calculated as in eq. (55). For nuclear transitions proceeding purely via a Fermi or Gamow-Teller operator, the angular dependence in the CM frame reduces to (see eq. (5)) where β ∈ [0, 1] is the speed of the outgoing lepton. The linear expressions seen here for either case allow the maximum of the distribution (needed for rejection sampling) to be found analytically. Nuclear transitions that do not correspond to one of these two cases are neglected in the current version of MARLEY. The maximum of the cos θ distribution for neutrino-electron elastic scattering is also found analytically. After a CM frame direction is sampled for the outgoing lepton, the active Reaction initializes Particle objects with the initial and final 4-momenta for the 2 → 2 scattering event. The values of these in the laboratory frame are stored in a newly-constructed Event object together with the excitation energy and spin-parity of the final energy level of the target.
At this point, program control passes back to the Generator, which applies two additional operations to the Event object before the event loop is complete. Once an Event object has been created by a Reaction, further operations on it are handled by derived instances of the EventProcessor abstract base class. These define a member function called process event which takes references to the Event and to the Generator as arguments.

Event loop: Nuclear de-excitations
For reactions involving a transition to an excited state of a nuclear target, the subsequent deexcitations are simulated using an EventProcessor called NucleusDecayer. This class manages a de-excitation loop that is initialized using the proton number (Z), nucleon number (A), excitation energy (E x ), spin (J), and parity (Π) of the outgoing nuclear state stored in the Event object. The values of these variables are updated during each iteration of the de-excitation loop, which simulates a sequence of binary decays until the nuclear ground state (E x = 0) is reached.
Each binary decay step begins with a check to determine whether the current nuclear state is bound or unbound (see section 4.1.4).

Unbound nuclear states
In the case of an unbound nuclear state, the NucleusDecayer constructs a HauserFeshbachDecay object, which provides a Monte Carlo implementation of the Hauser-Feshbach statistical model (see section 2.2) for decays of a compound nucleus. Upon construction, the HauserFeshbachDecay object uses the properties of the initial nuclear state (Z, A, E x , J, Π) and information from the Generator's owned StructureDatabase to initialize a member vector of pointers to ExitChannel objects.
The abstract base class ExitChannel represents a nuclear de-excitation via emission of a γray or a particular nuclear fragment. Two pairs of abstract derived classes virtually inherit from ExitChannel. The classes in the first pair, DiscreteExitChannel and ContinuumExitChannel, specialize to handling transitions to a specific discrete nuclear level and to the continuum, respectively. The classes in the second pair, FragmentExitChannel and GammaExitChannel, respectively represent de-excitations involving fragment and γ-ray emission. Concrete ExitChannel objects belong to a class that inherits from exactly one member of each pair, with the four possibilities being FragmentDiscreteExitChannel, FragmentContinuumExitChannel, GammaDiscreteExitChannel, and GammaContinuumExitChannel.
Every ExitChannel object owns a member variable called width , which is initialized during construction with the partial decay width of interest (in MeV −1 ) via a call to the protected member function compute total width. The FragmentDiscreteExitChannel and GammaDiscreteExit Channel classes implement this function using the expressions from eqs. (33) and (51), respectively. The FragmentContinuumExitChannel and GammaContinuumExitChannel classes both share the ContinuumExitChannel::compute total width implementation, which integrates the differential decay width over the energetically-accessible continuum (see sections 2.2.1 and 2.2.3) to obtain the width value where the emitted particle species u ∈ {p, n, d, t, h, α, γ}. The numerical integration is performed using Clenshaw-Curtis quadrature [142] by an instance of the Integrator class. Evaluation of the differential decay width dΓ u /dE x is delegated to the pure virtual member function differential width, which is implemented in FragmentContinuumExitChannel and Gamma ContinuumExitChannel according to the expressions given in eq. (29) and eq. (49), respectively. When the HauserFeshbachDecay object has been fully initialized, its owned vector of Exit Channel pointers will include an individual object derived from DiscreteExitChannel for each accessible transition to a discrete nuclear level present in the StructureDatabase. It will also contain a pointer to a single object derived from ContinuumExitChannel for each particle species that may be emitted via a transition to the continuum.
To simulate a binary decay step, the NucleusDecayer calls the do decay method of the new HauserFeshbachDecay object. This causes a particular ExitChannel e to be sampled with probability Here Γ e is the partial width for the eth ExitChannel (given by the width member variable) and the sum in the denominator runs over all of the owned ExitChannel objects. Control then passes to the do decay member function of the sampled ExitChannel to determine the properties of the final nuclear state. If a DiscreteExitChannel has been sampled, the final nuclear excitation energy (E x ), spin (J ), and parity (Π ) are retrieved from the associated Level object. For a ContinuumExitChannel, a ChebyshevInterpolatingFunction approximation to the dif- via inverse transform sampling (see section 3.3). The differential decay width is then evaluated (via a call to the differential width member function) for the selected final excitation energy E x . Each individual term τ b in the differential decay width sum (see either eq. (29) or eq. (49) as appropriate for the emitted particle species) is stored in a SpinParityWidth object together with its corresponding J and Π values. The final nuclear spin-parity is determined by choosing the bth SpinParityWidth object with probability where the differential decay width in the denominator is evaluated at the sampled value of E x . Once values of the variables defining the final nuclear state (E x , J , and Π ) have been determined, the procedure required to complete an iteration of the event loop (see section 4.3.3) is similar to that used for decays of bound nuclear levels.

Bound nuclear states
If the nucleus is in a bound excited state, i.e., it has a nonzero excitation energy E x which is below all of the nuclear fragment emission thresholds, then the NucleusDecayer consults the StructureDatabase to determine whether a DecayScheme has been previously configured for the nuclide of interest (see section 4.1.5). If this is the case, then simulation of the remaining deexcitations is delegated to the DecayScheme member function do cascade. During each binary decay step handled by this function, a γ-ray transition γ j originating from the current Level is chosen with probability where I j is the relative intensity of the jth Gamma owned by the current Level and the sum in the denominator runs over all of the owned Gamma objects. The properties of the final nuclear state, known immediately from the Level pointed to by the sampled Gamma object (via its end level member variable) are used to complete the de-excitation loop iteration as described in section 4.3.3. If a suitable DecayScheme cannot be found in the StructureDatabase, then MARLEY follows the same procedure as for unbound nuclear levels (see section 4.3.1).

Finishing a de-excitation step
At the end of each iteration of the de-excitation loop, a direction for the emitted γ-ray or nuclear fragment is chosen by sampling a polar cosine cos θ and an azimuthal angle φ isotropically in the rest frame of the decaying nucleus. This information is used together with the daughter particle PDG codes and masses to initialize two new Particle objects via a call to the utility function marley kinematics::two body decay. This function handles the elementary kinematical calculations needed to obtain the outgoing particle 4-momenta in the laboratory frame. The Particle representing the emitted γ-ray or fragment is appended to the vector of final particles owned by the Event object being processed. The Particle representing the daughter nucleus, on the other hand, replaces the mother nucleus Particle in the Event. After these adjustments have been made to the Event object and the values of the variables tracked in the de-excitation loop (Z, A, E x , J, Π) have been updated, simulation of the current binary decay step is complete.
The de-excitation loop will continue to simulate additional binary decay steps until the nuclear ground state is reached, or, in cases where no discrete level data are available for a particular final-state nucleus, until the excitation energy falls below the cutoff E cut x = 1 keV.

Event loop: Rotation of coordinates
For convenience during internal calculations, all MARLEY events are initially generated in a frame in which the incident neutrino direction lies along the positive z axis. If the user has specified a different neutrino direction in the job configuration file (see section 6.2), then an instance of the ProjectileDirectionRotator class (which inherits from EventProcessor) is used to process the Event after the de-excitation loop terminates. All particle 3-momenta in the event are rotated into a new reference frame in which the neutrino is traveling along the direction specified by the user. After this rotation is applied, generation of the new MARLEY event is complete.

Installation and usage
The current release of MARLEY has been tested on both Linux and macOS platforms and is expected to work in any Unix-like environment in which the prerequisites are installed. Building and running MARLEY on Windows is not currently supported. The installation instructions presented in this section assume use of the Bash shell [143] and the availability of several standard commandline tools.

Prerequisites
The following prerequisites are required to build the MARLEY source code: • A C++-14-compliant compiler. Two compilers are officially supported: -The GNU Compiler Collection (GCC) [144], version 4.9.4 or greater -The Clang frontend for the LLVM compiler infrastructure [145], version 3.5.2 or greater • The GNU Scientific Library (GSL) [1,2] On Linux architectures, all of these prerequisites will likely be available for installation through the standard package manager. On macOS, the use of Homebrew [147] to install GSL is recommended. This may be done by executing the terminal command brew install gsl after Homebrew has been installed on the target system.
Although it is not required in order to build or use MARLEY, the popular ROOT data analysis framework [3,4] provides convenient tools for plotting and analyzing simulation results. Users who wish to use the optional interface between the two codes (see sections 7.2.4, 7.3.4 and 7.3.5) should ensure that ROOT is installed before building MARLEY.
Other than GSL, MARLEY's only required external dependencies are the C++ standard library and the symbols defined in the dirent.h and sys/stat.h headers of the C POSIX library [148].

Obtaining and building the code
The source code for MARLEY 1.2.0 may be downloaded as a compressed archive file from the releases webpage (https://github.com/MARLEY-MC/marley/releases). Both zip and tar.gz file formats are available. After downloading the source code, the user should unpack the archive file in the desired installation folder. For the tar.gz format, this may be done via the command tar xvfz marley -1.2.0. tar . gz After unpacking the source code, the user may build MARLEY by navigating to the build/ subdirectory and invoking GNU Make: cd marley -1.2.0/ build make At build time, MARLEY verifies that GSL is installed by checking for the presence of the gsl-config script on the system path. Similarly, the optional MARLEY interface to ROOT is automatically enabled or disabled based on whether root-config script is present. Users who have an existing installation of ROOT but desire to build MARLEY without ROOT support may manually disable the interface by setting the IGNORE ROOT variable while invoking make, e.g., via make IGNORE_ROOT = yes Users wishing to contribute improvements to MARLEY may prefer to clone the official source code repository instead of downloading a release archive file. Instructions for doing so are available in the "developer documentation" section of the official website [149].

Configuring the runtime environment
At runtime, the marley command-line executable relies on the system environment variable MARLEY to store the path to the root folder of the source code. If generation of events is attempted without this variable being set, then the program will halt after printing the error message Although the user may manually set the value of the MARLEY variable, use of the Bash shell script setup marley.sh to configure the system environment is recommended. In addition to storing the path to the source code, this script makes several other changes to environment variables for user convenience, including adding the build/ folder to the system search paths for executables 17 and dynamic libraries. 18 The setup marley.sh script appears in the root source code folder and should be executed ("sourced") using the source command. From within the build/ folder, for example, one should source the script via source ../ setup_marley . sh Sourcing the setup script does not produce any console output. The instructions given in the remainder of this section assume that the setup marley.sh script has already been run in the current terminal session. 17 The PATH environment variable 18 Either the LD LIBRARY PATH (Linux) or the DYLD LIBRARY PATH (macOS) environment variable

Running a simulation
The typical procedure for running a MARLEY simulation is to invoke the executable via a command of the form 19

marley CONFIG_FILE
where CONFIG FILE is the name (with any needed path specification) of a job configuration file. Section 6 gives a full description of the file format and configuration parameters. Three example job configuration files are included with MARLEY in the examples/config/ folder: annotated.js A heavily-commented example that provides documentation of the configuration file format similar to the contents of section 6 COPY ME.js An example intended to serve as the basis for new job configuration files written by users minimal.js An example of the simplest possible job configuration: default settings are used for all parameters except those that must be explicitly specified As it executes, the marley program will print a variety of logging messages to the screen. Several rows at the bottom of the screen are reserved for a status display that tracks the progress of the simulation. Listing 1 shows the format of the status display. Following the first two rows, which report the current event count and the elapsed time since the simulation began, zero or more rows record the total amount of data written to each output file. The final row displays an estimate of the time at which the simulation job will be completed. If it becomes necessary to end the simulation before all requested events have been generated, the user may interrupt program execution by pressing ctrl+C. In response, the marley executable will terminate gracefully after writing the current event to any open output files. As discussed in section 6.6, two of the available output file formats allow for a simulation job interrupted in this way to be resumed from where it left off.

Troubleshooting rejection sampling problems
To select the energy E ν of the reacting neutrino in each event, MARLEY employs the rejection sampling technique discussed in section 3.2 and the probability density function P (E ν ) given in section 4.1.3. The validity of the technique depends on obtaining an accurate estimate of the global maximum P max of P (E ν ) within the sampling region of interest. Underestimates of P max will lead to bias in the neutrino energy distribution of the simulated events, while significant overestimates will adversely impact sampling efficiency.
Although numerical estimation of P max in MARLEY is reasonably robust for a variety of realistic neutrino spectra, it is not foolproof: there exist pathological energy distributions (e.g., those with multiple sharp peaks) for which automatic detection of the global maximum often fails. To protect against this problem, MARLEY checks that each calculation of P (E ν ) performed during neutrino energy sampling yields a value that does not exceed P max . If a value larger than P max is ever encountered, a set of warning and error messages similar to those in listing 2 will be printed to the screen. A new estimate of P max , given by the problematic value of P (E ν ) increased by a small "safety factor," will be adopted in subsequent event generation. energy_pdf_max : 0.00135331 , 7 If this error message persists after raising energy_pdf_max to a relatively high value , please contact the MARLEY developers for troubleshooting help .
Listing 2: Example warning and error messages printed in response to a rejection sampling problem encountered when selecting reacting neutrino energies While this adaptive approach to improving estimation of P max may resolve rejection sampling problems in the later part of the simulation, it cannot correct for bias in the energy distribution of the events that have already been generated. A simple strategy for overcoming this difficulty is for the user to restart the simulation from the beginning. By adding a line to the job configuration file containing the energy pdf max key (see section 6.7.7) and the improved estimate of P max recommended by the rejection sampling error message (see line 6 of listing 2), automatic estimation of P max will be disabled in favor of adopting the user-specified value.
For severe underestimations of P max , several repetitions of this "interrupt and rerun" strategy (each with a larger value of energy pdf max than the one before) may be needed before correct behavior is achieved. If rejection sampling errors persist across multiple fix attempts, or if the corrected value of P max results in prohibitively slow performance, a different approach to running the simulation (e.g., splitting a problematic spectrum into energy regions which are handled by separate MARLEY jobs) should be pursued.

The marley-config utility
For a variety of applications, e.g., analyzing the standard MARLEY output files (section 7.3) and interfacing MARLEY with external codes (section 8.1), it may be desirable to make use of MARLEY classes within an external C++ program. To facilitate compilation of such programs, a Bash script named marley-config is placed in the build/ folder whenever MARLEY is successfully built. 20 This script may be queried to obtain various pieces of information about the build.

Directory contents
Relative path --bindir executables build/ --datadir input data data/ --incdir header files include --libdir libraries build/ --srcdir source files src/ --topdir top-level MARLEY directory ./ Each marley-config query is executed by providing the script with one or more command-line options, each beginning with a prefix of two hyphens (--). The command-line options may appear in any order and may be repeated. If the --help option is present on the command line, then all other arguments are ignored, and a multi-line usage message is printed. All other recognized options produce console output that will be combined into a single line and printed in the order that the corresponding options appeared on the command line.
Three categories of command-line options are recognized by the marley-config script. The compilation options are used to retrieve compiler flags (compatible with both GCC and Clang) and other helpful information for building external programs that link to the MARLEY shared libraries (see section 5.5.1). Brief descriptions of the four compilation options currently recognized by marley-config are given below.
--cflags Prints the C++ compiler flags needed to build code that includes MARLEY header files. If MARLEY has been built with ROOT support, then the compiler flags will include -DUSE ROOT, which defines the preprocessor macro USE ROOT. This macro may be used together with the preprocessor directives #ifdef or #ifndef to test for ROOT-dependent MARLEY features in compiled code. An example is given in section 5.5.1.
--cxx-std Prints the version of the C++ Standard used by the compiler when MARLEY was built. The format used is the same as for parameters passed via the -std flag to the compiler. For example, if MARLEY was built using the -std=c++14 flag, then marley-config will print c++14 for this option.
--libs Prints the compiler flags needed for linking to the MARLEY shared libraries --use-root Prints the string yes if MARLEY was built with ROOT support and no if it was not The directory options are used to print the full paths to various subfolders of the MARLEY installation. Table 2 lists the available directory options in the first column. The second and third columns list, respectively, a description of the directory's contents and its path relative to the root MARLEY folder. For an installation of MARLEY in which the top-level directory of the source code is /home/marley, executing marley -config --topdir --incdir --bindir will produce the output / home / marley / home / marley / include / home / marley / build A final marley-config option category, the version options, may be used to identify the installed version of MARLEY. Two options in this category are currently recognized by the marley-config script: --git-revision and --version.
The --git-revision option is used to print a unique hash value reported by Git at build time. If MARLEY was built using a cloned Git repository (see section 5.2), and the git executable is present on the system path when GNU Make is invoked, then the hash identifier for the latest commit will be retrieved 21 and saved in the marley-config script. If uncommitted changes exist in the Git index or working tree when MARLEY is built, then the string -dirty will be appended to the hash value. In cases where a hash value is unavailable (e.g., because MARLEY was installed without cloning a Git repository), the marley-config script will print the string unknown version in response to this option.
If a tagged release of MARLEY was built, 22 then the --version option may be used to print the corresponding version number (e.g., 1.2.0). Otherwise, the output produced by marley-config in response to --version is the same as for --git-revision.

Compiling programs using marley-config
To illustrate usage of the marley-config utility, three simple C++ programs that employ MARLEY classes are provided in the folder examples/executables/minimal/. Each one of these may be compiled using GCC by executing a command of the form g ++ -o EXECUTABLE_NAME $ ( marley -config --cflags --libs ) SOURCE_NAME from within examples/executables/minimal/. Here EXECUTABLE NAME is the desired file name for the compiled executable and SOURCE NAME is the name of one of the example C++ source files.

Generator configuration
Before entering the event loop, the MARLEY executable reads in various settings (see section 4.1) from a job configuration file. The format of the job configuration file is based on JavaScript Object Notation (JSON) [151], which represents data structures as a nested hierarchy of key-value pairs. A JSON object is an unordered set of key-value pairs surrounded by curly braces ({}). Each key is an arbitrary string delimited by double quotes ("") and separated from its corresponding value by a colon. Each value may be a JSON object itself, a double-quoted string, a number, an array, or one of the words true, false, or null (without surrounding quotes). A JSON array is an ordered collection of values surrounded by square brackets ([]). Elements of an array need not have the same data type, e.g., the values true and 1.05 may be members of the same array. Array elements are separated by commas, as are object key-value pairs. A valid JSON document is itself an object and is thus delimited by a pair of curly braces.
For the convenience of users, the job configuration file format used by MARLEY permits three minor deviations from the JSON standard: • Single-word keys (no whitespace) may be given without surrounding double quotes.
• A trailing comma is allowed after the final element in objects and arrays.
A filename extension of .js is recommended for MARLEY job configuration files because the default JavaScript syntax highlighting used by many text editors is well-suited to the JSON-like format. 23 To parse job configuration files, MARLEY relies on a modified version of the SimpleJSON library [152], which was incorporated into the code base as the class marley::JSON. An example MARLEY job configuration file is shown in listing 4. Explanations of each parameter needed to fully configure the generator are given in the following subsections.

Random number seed
The optional seed key may be used to specify an integer seed for the random number generator. Any value corresponding to a C++ long int (guaranteed by the standard to be at least 32 bits long) may be used as a seed. This includes negative integers, although users should note that these will be reinterpreted as unsigned integers by MARLEY. If the seed key is omitted from the job configuration file, then the current Unix time will be used as the random number seed.

Neutrino direction
For simplicity, MARLEY initially generates all events in a frame where the incident neutrino is traveling along the positive z direction (see section 4.4). However, it is often useful to simulate events for neutrinos traveling in an arbitrary direction in the lab frame. This functionality is enabled via use of the direction key to define a lab-frame direction 3-vector for the incident neutrinos. The Cartesian components of this vector are specified as shown on line 7 of listing 4. They need not have any particular normalization, but at least one component must be nonzero. If the direction key is omitted, then the positive z direction is assumed.
In certain special cases, e.g., studies of the Diffuse Supernova Neutrino Background (DSNB) [153], it may be desirable to generate MARLEY events with an isotropic distribution of initial neutrino directions. This behavior may be enabled by providing the string "isotropic" as the value for the direction key instead of a direction 3-vector.

Target composition
The nuclidic composition of the neutrino target material in a MARLEY simulation may be specified by defining a JSON object in the job configuration file labeled with the target key. The JSON object must define two arrays of equal size. The first of these, which is labeled with the nuclides key, includes one or more nuclear PDG codes (see section 7.1) with one code per species. The second array, which is labeled with the atom fractions key, contains the corresponding nuclide fractions in the target material. After being automatically normalized to sum to unity, these will be used as the abundance weights f r that appear in eqs. (74) and (75). If any of the elements of the atom fractions array is negative, MARLEY will halt and print an error message.
Lines 10-13 of listing 4 provide an example configuration for a neutrino target composed of pure 40 Ar. If the target key is omitted from the job configuration file, then the default behavior mentioned in section 4.1.6 will be used: every distinct nuclide that appears in the initial state of at least one configured reaction (see section 6.4) will be included in the neutrino target with equal abundance.

Reactions
To define the set of neutrino scattering processes that should be considered in a MARLEY simulation, the user must specify a list of reaction input files as a JSON string array labeled by the reactions key. This key is required to be present in all MARLEY job configuration files. For neutrinonucleus reactions, each reaction input file provides the values of the allowed nuclear matrix elements B(F) and B(GT) needed to compute the differential cross section in eq. (5) for a particular target nucleus and interaction mode. For neutrino-electron elastic scattering, the reaction input file provides a list of nuclides for which this process should be enabled. Details of the file format are discussed in appendix B.
In MARLEY 1.2.0, the code is capable of simulating reactions belonging to four distinct interaction modes: charged-current neutrino-nucleus scattering (ν CC), charged-current antineutrino-nucleus scattering (ν CC), neutral-current (anti)neutrino-nucleus scattering (NC), and (anti)neutrino-electron elastic scattering (ES). As shown in eq. (10), the nuclear scattering modes may be distinguished by the isospin operator that appears in the nuclear matrix elements (e.g., t − for ν CC).

Example reaction input files
Three example reaction input files for ν CC on 40 Ar, intended for use in simulations of the process 40 Ar(ν e , e − ) 40 K * , are included in the data/react folder as part of the standard MARLEY distribution. At high excitation energies, all three files include the same set of theoretical matrix elements taken from a QRPA 24 calculation by Cheoun, Ha, and Kajino [154]. At low excitation energies, the files contain different sets of nuclear matrix elements extracted from experimental data. The references used to obtain the data-driven matrix elements for each file are as follows: ve40ArCC Bhattacharya1998.react Measurement of 40 Ti β + decay reported in ref. [155] ve40ArCC Liu1998.react Independent 40 Ti β + decay measurement reported in ref. [156] ve40ArCC Bhattacharya2009.react Forward-angle (p,n) scattering data reported in ref. [157].
A description of the data evaluation procedure used to prepare these files is available in ref. [108,sec. III].
Two other reaction input files are included in the data/react folder. The file ES.react enables simulations of neutrino-electron scattering on an atomic 40 Ar target. Simulation of this process for other nuclides may be enabled in ES.react by appending additional nuclear PDG codes.
The file CEvNS40Ar.react provides a simple example of an NC configuration. Although MARLEY is capable of simulating inelastic NC reactions if provided with suitable input, CEvNS40Ar.react includes only one B(F) matrix element appropriate for simulating coherent elastic neutrino nucleus scattering on 40 Ar (see section 2.1.1). Altering the nuclear PDG code given in this file (see appendix B) will immediately enable this process to be simulated for any spin-zero target nucleus.

Configuration example
Line 16 of listing 4 provides an example reactions setting that enables ν CC and ES reactions on 40 Ar to be simulated simultaneously. As shown in the example, a full path specification is not needed for files that appear in the ${MARLEY}/data/react/ folder, where ${MARLEY} is the value of the MARLEY environment variable at runtime (see section 5.3). If multiple files representing the same reaction mode and the same target nuclide are listed for the reactions key, a warning message will be printed and only the configuration from the first such file will be used in the simulation.
Every element of the reactions array must be a simple string literal containing a single file name. The use of wildcard characters, regular expressions, and environment variables is not currently allowed by the MARLEY job configuration file format.

Neutrino source
The required source key allows the user to provide a description of the incident neutrino energy spectrum using a JSON object. This object must always contain at least two keys: neutrino and type.
The neutrino key represents the neutrino species to be used as the projectile. This may be specified using either an integer PDG code or a string, as shown in table 3.
The type key may be used to request one of several built-in neutrino spectra or to indicate to MARLEY that a user-defined spectrum should be used. In the latter case, users should typically not apply any cross section weighting themselves to the neutrino spectrum, as MARLEY will weight the incident spectrum using the appropriate reaction cross section at runtime (see section 4.1.3). In unusual cases where this behavior is not desirable, the user may disable cross section weighting by following the procedure described in section 6.7.3.
The following paragraphs describe each of the allowed options for the type key, together with any additional parameters needed for each source type.

Quasiparticle Random Phase Approximation
Neutrino PDG code string ν e 12 "ve" ν e −12 "vebar" ν µ 14 "vu" ν µ −14 "vubar" ν τ 16 "vt" ν τ −16 "vtbar" Table 3: Allowed values for the neutrino key in the neutrino source specification 6.5.1. Monoenergetic A neutrino source whose type is mono or monoenergetic always produces a neutrino with a fixed energy E ν . This energy is specified in MeV using the energy key. Lines 19-23 of listing 4 provide an example configuration for a monoenergetic ν e source with E ν = 15 MeV.

Fermi-Dirac distribution
The fermi-dirac or fd source type emits neutrinos whose energies are drawn from a Fermi-Dirac distribution with temperature T (MeV) and pinching parameter η. The incident neutrino spectrum is given by where C is a normalization constant chosen so that The form for the pinched Fermi-Dirac distribution used by MARLEY is based on that of the "Livermore model" of the supernova neutrino energy spectrum [158].

"Beta fit" spectrum
The beta-fit or bf source type produces neutrinos with energies drawn from the spectrum The values of the parameters E min ν , E max ν , E ν , and β are specified using the keys Emin, Emax, Emean, and beta, respectively. Omitting the beta key will yield the default value β = 4.5. The other three parameters must be given explicitly. Negative β values are unphysical and will cause MARLEY to halt and issue an error message. Listing 6 shows an example "beta fit" source specification in which the parameter values E min ν = 0 MeV, E max ν = 60 MeV, E ν = 15 MeV, and β = 3 have been chosen.
The "beta fit" distribution has been employed in a recent theoretical study [159] of supernova neutrino detection and is equivalent to the widely-used "Garching model" parameterization of supernova neutrino spectra [160]. The sole difference between the two parameterizations is in the convention for the pinching parameter. The Garching model employs a parameter α and makes the substitution β → α + 1 in eq. (87). For α = β − 1, the distributions used by the two approaches are identical.

Muon decay-at-rest
The decay-at-rest or dar source type emits neutrinos according to a Michel spectrum for a muon decaying at rest in the lab frame. For a ν e source (corresponding to the decay µ + → e + + ν e +ν µ ), the energy spectrum is given by [161] φ where m µ is the muon mass. Aν µ source has the spectrum Production ofν e and ν µ (from µ − → e − +ν e +ν µ ) may also be simulated using a dar neutrino source. Because the decay kinematics fully determine the shape of the Michel spectrum, no additional input parameters are needed for this source type. Listing 7 shows an example configuration for a muon decay-at-rest ν e source.
Neutrinos from muon decay-at-rest produced at facilities like the Oak Ridge Spallation Neutron Source (SNS) [162] and the Fermilab Neutrinos at the Main Injector (NuMI) beamline [36] offer a valuable opportunity to study neutrino-nucleus scattering at tens-of-MeV energies with a terrestrial source.

Histogram
The histogram or hist source type allows the user to define the incident neutrino spectrum as a histogram with one or more energy bins. The n bins are specified by their lower edges E ν (k) and weights W (k) for k ∈ {1, 2, . . . , n}. Within bin k, incident neutrino energies are distributed uniformly on the half-closed interval [E ν (k), E ν (k + 1)). Weights that do not sum to unity will be automatically renormalized by the code. The keys E bin lefts and weights are used to provide JSON arrays of E ν (k) and W (k) values, respectively. The sizes of these arrays must be equal and are used by MARLEY to determine the bin count n. The energy bin edges E ν (k) must be given in increasing order. The upper edge of the final energy bin E ν (n + 1) must also be specified using the key Emax.
Listing 8 shows an example histogram source specification in which three equally-spaced bins are defined between 7 and 10 MeV.

Grid
The grid source type emits neutrinos with energies drawn from a tabulated probability density function φ(E ν ). To define this function, a grid of at least two monotonically increasing E ν values is specified as a JSON array via the energies key. The corresponding φ(E ν ) values (which need not be normalized by the user to integrate to unity) must also be given in an array of the same size using the prob densities key.

Rule
Interpretation "const" histogram-like "lin" linear in both E ν and φ(E ν ) "linlog" linear in E ν , logarithmic in φ(E ν ) "loglin" logarithmic in E ν , linear in φ(E ν ) "log" logarithmic in both E ν and φ(E ν ) Table 4: Allowed values of the rule key for the grid neutrino source type.
To evaluate φ(E ν ) at energies between the grid points, MARLEY employs an interpolation rule specified by the rule key. The default "lin" rule (linear in both E ν and φ) is typically most useful, but MARLEY is also capable of interpolating logarithmically along one or both axes. A "const" interpolation rule is also available which treats φ(E ν ) as a step function at each energy grid point. Table 4 lists all allowed values of the rule key together with their interpretations. Listing 9 shows an example grid source specification which defines a symmetric triangular distribution between 10 and 20 MeV.

TH1 and TGraph
If MARLEY's optional interface to the ROOT [3] data analysis framework is enabled (see section 5), then two additional neutrino source types, th1 and tgraph, are allowed. The th1 type takes a one-dimensional ROOT histogram object (of C++ type TH1) as input and converts it into a histogram neutrino source using MARLEY's native representation of the distribution. 25 The tgraph source type is similar, except that it converts a ROOT TGraph object into a grid neutrino source.
To use either of the th1 or tgraph source types, the user must provide the name of an input ROOT file (including any needed path specification) using the tfile key. The string value given for the namecycle key, which must also be provided, will be used by MARLEY in a call to TDirectoryFile::Get(const char* namecycle) in order to retrieve the object of interest from the input ROOT file. Units of MeV should always be used for neutrino energies when preparing a TH1 or TGraph object for use with MARLEY.
Listing 10 shows an example th1 source specification which uses a neutrino energy histogram named "MyHist" stored in a ROOT file in the working directory called my root file.root.

Executable settings
The optional executable settings key is used to provide a JSON object that controls the marley command-line executable. However, the entirety of the executable settings will be ignored if the configuration file is used to initialize MARLEY in a different context (e.g., from within a Geant4 application that links to the MARLEY shared libraries, see section 8).
Within the executable settings object, the events key is used to specify the number of events that the marley executable should generate before terminating. Because the MARLEY JSON parser expects the associated value to be an integer literal, using scientific notation to specify the number of events (e.g., 1e4 instead of 10000) is not currently allowed. If the events key is omitted, a default value of 1000 will be assumed.
The output key in the executable settings object is used to describe zero or more output files that will receive the generated events. The corresponding value should be a JSON array containing one JSON object per output file. Each object in the array should define the following keys: file The name of the output file (with any needed path specification).
format The format to use when storing events in the output file. Valid values for this key are "ascii", "hepevt", "json", and (if MARLEY's ROOT interface is active) "root". Descriptions of each of the allowed output file formats are given in section 7.2.
mode The approach that the marley executable should use if the output file is not initially empty. For the ASCII and HEPEVT formats, valid values are "overwrite" (erase any previously existing file contents) and "append" (continue output immediately after any existing file contents). For the JSON and ROOT formats, valid values are "overwrite" and "resume". If the "resume" mode is chosen, the generator will restore its previous state from an incomplete run (e.g., a run that was interrupted by the user by pressing ctrl+C) that was saved to the output file and continue from where it left off. If the "resume" mode is selected for more than one output file, then MARLEY will halt and print an error message.
Two additional keys may be used in specific contexts: force Boolean value used only for the "overwrite" mode. If it is true, then the marley executable will not prompt the user before overwriting existing data. If this key is omitted, a value of false is assumed.
indent Nonnegative integer value used only for the JSON output format. It represents the number of spaces that should be used as a tab stop when pretty-printing the JSON output. If this key is omitted, all unnecessary whitespace will be suppressed. This default behavior results in the most compact JSON-format output files.
If the output key in the executable settings object is omitted entirely, then the default configuration output : [ { file : " events . ascii " , format : " ascii " , mode : " overwrite " } ] will be used. Lines 26-32 of listing 4 define an example executable settings object that includes the default settings described in this section.

Less commonly-used parameters
Although the job configuration file shown in listing 4 provides usage examples for all of the key-value pairs typically needed to define a MARLEY simulation, a number of additional parameters are recognized by the configuration parser. 26 These parameters, which are intended for advanced users or for addressing unusual situations, are documented in the remainder of this section.

Logging
The singleton marley::Logger class provides rudimentary support for writing diagnostic messages to zero or more output streams including standard output, standard error, and plaintext log files. This is done in a prioritized way by configuring a logging level for each output stream. Messages are likewise categorized by a logging level. In order of severity, the logging levels recognized by marley::Logger are debug, info, warning, error, and disabled, with the last available for output streams but not for messages. When a message with logging level MSG LEVEL is passed by MARLEY to the Logger, the message will be ignored by all streams whose logging levels are more severe than MSG LEVEL. The message will be written to all other configured streams.
Nearly all output written to the terminal by the marley executable is handled by the Logger class. The main exception is the status display discussed in section 5.4 (see listing 1).
The default behavior of the Logger class may be adjusted via the optional log key in the job configuration file. This key is used to specify an array of JSON objects with one element per output stream. Each object in the array may define the following keys: file A string giving the name of a text file that will receive the output logging messages. If the string value is "stdout" ("stderr"), then the messages will be written to standard output (standard error) rather than a text file.
level A string specifying the logging level that should be used as a threshold for writing messages to the stream. Valid values are "debug", "info", "warning", "error", and "disabled".
overwrite A boolean value indicating whether new messages should be appended to the end of the output file (false) or whether any previously existing file contents should be erased (true).
This key is ignored when writing to standard output and standard error. A default value of false is assumed when this key is not present.
If the log key is omitted from the job configuration file, then MARLEY will use default settings equivalent to log : [ { file : " stdout " , level : " info " } ] 6.7.2. Disabling simulation of nuclear de-excitations As described in section 2, MARLEY models neutrino-nucleus scattering as proceeding via a prompt 2 → 2 reaction and a subsequent sequence of zero or more binary decays of the remnant nucleus. For certain calculations, 27 it may be convenient to simulate events in which only the primary 2 → 2 scattering process is considered. This behavior may be enabled in the job configuration file by using the optional do deexcitations key together with a boolean value. If the accompanying value is true, then the MARLEY event loop will proceed normally, and nuclear de-excitations will be simulated. If it is false, then the de-excitation loop will be skipped.

Disabling cross-section weighting of the neutrino spectrum
By default, MARLEY samples the reacting neutrino energy E ν for each event from a probability density function P (E ν ) ∝ σ(E ν ) φ(E ν ), where σ(E ν ) is the abundance-weighted total cross section and φ(E ν ) is the incident neutrino spectrum (see section 4.1.3). For some use cases, it may be desirable to sample E ν directly from φ(E ν ) without any cross-section weighting. This may be achieved by adding the optional weight flux key, which takes a boolean value, to the contents of the source JSON object (section 6.5). If the weight flux value is set to true, then E ν will be sampled from the usual probability density function P (E ν ). If the weight flux value is set to false, then E ν will be drawn directly from φ(E ν ).

Orbital angular momentum and multipolarity cutoffs
When computing differential decay widths for nuclear fragment (γ-ray) emission from a compound nucleus, MARLEY truncates the infinite sum over orbital angular momenta (multipolarities) that appears in eq. (29) (eq. (49)) by imposing a cutoff value max (λ max ). The user may modify the default value ( max = λ max = 5) by including two keys in the job configuration file. The fragment lmax key is used to specify a nonnegative integer which will be used as the maximum orbital angular momentum max to consider in nuclear fragment emission. Similarly, the gamma lmax key is used to specify a positive integer which will serve as the maximum multipolarity λ max considered in γ-ray emission.

Coulomb correction factor
As described in section 2.1.2, MARLEY accounts for final-state Coulomb effects in charged-current neutrino-nucleus interactions by the inclusion of a Coulomb correction factor F C in the expression for the differential cross section. By default, this factor is computed according to eq. (27), which uses the smaller of the two corrections given by the Fermi function (eq. (18)) and the modified effective momentum approximation (MEMA, eq. (26)). The MARLEY prescription for the Coulomb correction factor may be controlled in the job configuration file by use of the coulomb mode key. The following string values are allowed for this key: "none" No Coulomb corrections will be applied (F C = 1 in all cases).
"Fermi" The Fermi function defined in eq. (18) will always be used to compute the Coulomb correction factor (F C = F Fermi ).
"Fermi-EMA" The Fermi function will be compared to the EMA, and the approach that yields the smaller correction (i.e., the F C value closer to unity) will be used. The Coulomb correction factor F C is calculated as in eq. (27) with the substitution F MEMA → F EMA .
"Fermi-MEMA" Similar to "Fermi-EMA", except that the MEMA is used instead of the EMA. The Coulomb correction factor F C is calculated as in eq. (27). This is MARLEY's default approach.

Status update interval
As discussed in section 5.4, the marley executable prints a text-based status display at the bottom of the screen. By default, this display is updated after each set of n update = 100 events has been generated. The user may modify the value of n update by adding the status update interval key to the executable settings object in the job configuration file (see section 6.6). Only positive integer values of n update may be specified for this key.

Manually specifying a maximum value of the neutrino energy PDF
To select a reacting neutrino energy for each event, MARLEY relies on rejection sampling from the probability density function given in eq. (73). In situations where automatic estimation of the global maximum P max of this PDF is found to be inadequate, MARLEY will print a set of messages to alert the user (see section 5.4.1). These messages will include an improved estimate of P max (see line 6 of listing 2) which may be manually supplied by the user in the job configuration file as a floating-point value associated with the energy pdf max key. Doing so will disable numerical estimation of P max in favor of using the specified value, which may help to resolve some rejection sampling problems. xsec dump steps Used to set the number of equally-spaced points at which the total cross section should be evaluated as a function of projectile kinetic energy. Only positive integer values are allowed. If this key is absent, then mardumpxs will use a default value of 10,000.

Interpreting the output
The neutrino scattering events generated by the marley command-line executable may be saved to disk in four distinct output formats. Section 6.6 explains how one or more of these output formats may be selected by the user in the job configuration file. Following a brief discussion of the particle numbering scheme used by MARLEY in section 7.1, documentation for each of the output file formats is provided in section 7.2. Section 7.3 presents a C++-based API which may be used to parse all of the standard output formats via a unified interface. Section 7.4 provides guidance on how histograms of quantities of interest from the generated events may be converted into easily interpretable physics distributions, e.g., differential cross sections and event rates.

Particle numbering scheme
Like many other modern physics event generators, MARLEY identifies particle species using a standard set of integer codes defined by the Particle Data Group (PDG) [163, sec. 44, pp. 661-664]. Each kind of particle is assigned a unique positive integer, and the corresponding antiparticle is assigned a negative integer with the same absolute value. These PDG codes are used by MARLEY both internally and in output files. Table 5 provides the interpretation of the most common PDG codes that appear in MARLEY events. In general, a nucleus with proton number Z and mass number A is labeled with the PDG code P ID = 10,000 Z + 10 A + 1,000,000,000 .
The only two exceptions to the rule given in eq. (90) are for a free neutron and 1 H, which are represented by the codes 2112 and 2212, respectively.

Output file formats
In the current version of MARLEY, generated events may be saved to disk using four possible file formats: ASCII The native text-based format for MARLEY events.
HEPEVT A legacy format for interfacing event generators with each other and with other software. Despite its relative age and emergence in the specific context of QCD event generation for Large Electron-Positron Collider experiments [164], compatibility with this format is maintained in many modern high energy physics software tools. Examples include the HepMC3 [165] event record library and (via its TextFileGen module) the LArSoft toolkit discussed in section 8.3.
JSON In addition to serving as the input language for MARLEY job configuration files (with slight extensions to the standard grammar, see section 6 for details), JSON is also available as an output file format.
ROOT If MARLEY's interface to ROOT is active, a binary representation of the full marley::Event objects may be saved to disk in an output TTree. A simplified "flat" form of the ROOT output that can be read without recourse to the MARLEY shared libraries may also be produced using the marsum command-line tool (see section 7.3.5).

ASCII file format
An ASCII-format output file begins with the line where Ni (Nf) is the number of particles in the initial (final) state. The three remaining fields in the event header report properties of the recoiling nucleus immediately following the prompt 2 → 2 scattering reaction. The quantity Ex is the excitation energy (MeV), twoJ is an integer equal to two times the total spin, 28 and P is a single character representing either positive (+) or negative (-) parity. Following the event header, each of the Ni initial-state particles is described by a single line of the form

PDG Etot Px Py Pz M Q
where PDG is the PDG code (see section 7.1) identifying the particle species and Etot, Px, Py, and Pz are the Cartesian components of the particle 4-momentum (in MeV). The particle mass M (MeV) and (net) electric charge Q (in units of the elementary charge) are also included in each line. To complete the event record, Nf lines describing the final-state particles appear in the same format used for the initial-state particles.
To preserve full numerical precision when converting back and forth between marley::Event objects held in memory and ASCII-format event files, all floating point numbers are output in scientific notation with the required number of base-10 digits needed to uniquely represent all distinct values of the C++ type double. 29 However, this level of precision is not enforced by the code when reading events as input from an ASCII file. Listing 11 shows an example MARLEY output file in ASCII format. Line 1 gives the value of 5.984 × 10 −19 MeV −2 = 2.330 × 10 −40 cm 2 for the flux-averaged total cross section. Line 2 begins the record for the first event, which involves a transition to the 40 K nuclear level with an excitation energy of 3.797 MeV above the ground state and spin-parity 1 + . Lines 3-4 describe the initial state: a 10 MeV ν e traveling in the +z direction toward a 40 Ar atom at rest. Lines 5-8 describe the final-state particles: a 5.2 MeV electron, the recoiling 40 K ion, and two de-excitation γ-rays with energies of 1.54 and 2.26 MeV. The second and final event in the file, which appears on lines 9-14, involves a ν e -40 Ar collision which produces an electron, a 39 K ion, and a neutron.

HEPEVT file format
For the sake of brevity, only those aspects of the HEPEVT format needed to interpret MARLEY output are discussed in this section. A full description of the HEPEVT standard is available in ref. [164, pp. 327-330].
A HEPEVT-format output file consists of one or more text-based event records. Each of these records begins with the header

NEVHEP NHEP
where NEVHEP is the event number (untracked by MARLEY and thus always set to zero) and NHEP is the number of particles in the event. The header is followed by NHEP lines, each representing a single particle. These have the format where ISTHEP is an integer code identifying the particle status and IDHEP is the PDG particle ID code. In agreement with the HEPEVT standard, MARLEY uses status code 1 for the finalstate particles and 3 for the initial-state particles, the latter of which are not considered part of the event history [164]. The JMOHEP1, JMOHEP2, JDAHEP1, and JDAHEP2 entries record the indices j (1 ≤ j ≤ NHEP) of particles in the event record that correspond to the first mother, second mother, first daughter, and last daughter of the current particle, respectively. These indices are set to zero in cases where they do not apply (e.g., a particle which has not decayed will have JDAHEP1 = JDAHEP2 = 0). Entries PHEP1 through PHEP3 record the x, y, and z components of the particle 3-momentum, while PHEP4 gives the total energy and PHEP5 gives the particle mass (all in GeV). Entries VHEP1 through VHEP3 store the x, y, and z positions of the particle production vertex (mm), and VHEP4 gives the production time (mm/c).
Because MARLEY currently treats all nuclear de-excitations as instantaneous and does not perform any particle tracking, VHEP1 through VHEP4 are always identically zero in HEPEVT output files. Intermediate de-excitation steps are also not currently stored in the event record, so JMOHEP1, JMOHEP2, JDAHEP1, and JDAHEP2 are also identically zero in most cases.
In addition to the initial-and final-state particles, MARLEY adds a dummy particle with ISTHEP = 11 to each HEPEVT event record. All data fields are zero for this particle except for (1) JMOHEP1, which contains the nuclear spin multiplied by two, (2) JMOHEP2, which reports the parity of the nucleus as an integer, (3) PHEP4, which gives the excitation energy of the nucleus (MeV), and (4) PHEP5, which records the flux-averaged total cross section (see section 7.4) in units of MeV −2 / atom. As is the case for the ASCII format, the excitation energy, spin, and parity values refer to the nuclear state that existed immediately after the initial 2 → 2 scattering reaction.
Each output file includes two top-level keys: events, which is associated with an array of event objects, and gen state, which stores a JSON object representation of the generator state at the moment that the file was created.
Each element of the events array is a JSON object containing five key-value pairs. The first three of these, Ex, twoJ, and parity, provide the excitation energy (MeV), two times the total spin, and the parity of the final nucleus after the initial 2 → 2 scattering reaction but before any de-excitations have taken place. The other two keys, initial particles and final particles, are used store arrays of particles represented as JSON objects. Each particle object defines the following keys: (1) charge: the (net) electric charge in units of the elementary charge, (2) pdg: the PDG code identifying the particle type, (3) E: the total energy, (4) px: the x momentum component, (5) py: the y momentum component, (6) pz: the z momentum component, and (7) mass: the particle mass. The particle 4-momentum components and mass are all given in natural units (MeV).
The gen state JSON object includes several key-value pairs. The config key refers to a JSON object which reproduces the full contents (except for comments) of the job configuration file used to generate the events. The event count, flux avg xsec, and seed keys label the total number of events generated at the time the file was written, the flux-averaged total cross section (MeV −2 / atom, see section 7.4), and the integer random number seed used to initialize the event generator. A final key, generator state string, records a string representation of the internal state of the std::mt19937 64 object used by MARLEY to obtain pseudorandom numbers.
An example MARLEY output file in JSON format (example.json) is included in the supplemental materials.

ROOT file format
If MARLEY has been built against the appropriate shared libraries from the ROOT data analysis toolkit (see section 5), then the marley command-line executable may also produce output files in the standard ROOT compressed binary format.
Within a ROOT-format output file, access to the generated events is managed by an instance of the ROOT TTree class, which provides a hierarchical data structure for storing many objects belonging to the same C++ type. In general, a TTree may own one or more branches (represented by the TBranch class), each of which owns one or more leaves (represented by TLeaf). Branches may be read from a file independently of one another, allowing for efficient access to elements of a large dataset stored in a suitably-organized TTree [166].
The ROOT-format output files generated by MARLEY contain a single TTree called MARLEY event tree. This TTree contains a single branch called event, which stores one marley::Event object per tree entry. Although direct access to the events is possible by manipulating the MARLEY event tree, use of the simplified C++ API described in section 7.3 is recommended instead.
In addition to the MARLEY events themselves, four pieces of metadata are stored in a ROOTformat output file: In listing 13, an example C++ function is shown which retrieves all of these metadata objects from a file called events.root and prints their contents to standard output.
1 // Standard library includes 2 # include < iostream > 3 # include < string > 4 5 // ROOT includes 6 # include " TFile . h " 7 # include " TParameter . h " TFile event_file ( " events . root " , " read " ) ; The ROOT output files described above may be converted into an alternative "flat" format which is readable by ROOT without a need for the MARLEY shared libraries. This is done by running the marsum command-line tool described in section 7.3.5.

C++ analysis API
Although the file format descriptions given in section 7.2 should be sufficient to enable processing of MARLEY events using any programming language, a C++ API allowing manipulation of events stored in any of the standard output formats has been developed for the convenience of users. The API is usable within compiled code as well as via the interactive C++ interpreters included with ROOT. 30

The EventFileReader class
Programmatic access to MARLEY event records stored in a file is provided by the marley::Event FileReader class. The constructor of this class takes as its sole argument a std::string containing the name (including any needed path specification) of the file to be parsed. Events may be read one-by-one from the file using the stream extraction operator >>, which returns a boolean value indicating whether a new event was successfully loaded.
Listing 14 shows the source code for examples/executables/minimal/efr.cc, an example C++ program that illustrates the recommended usage of the EventFileReader class. Line 12 creates a new EventFileReader object that will read MARLEY events from a file whose name is given as the first command-line argument when the program is run. The while loop on lines 19-21 streams all events in the file to standard output, where they will be printed in ASCII format. Line 15 contains a call to the EventFileReader member function flux averaged xsec, which returns the flux-averaged total cross section (see section 7.4) used to generate the events in units of 10 −42 cm 2 / atom. Passing a boolean value of true to this function will cause it to return the cross section in natural units (MeV −2 / atom) instead.
1 # include < iostream > 2 3 # include " marley / Event . hh " 4 # include " marley / EventFileReader . hh " To build the efr.cc example program, it is necessary to link the executable to the MARLEY shared library. This is most easily done using the marley-config script with the syntax described in section 5.5.1.
The EventFileReader class can parse MARLEY output files written in the ASCII, HEPEVT, and JSON formats. If MARLEY has been built against ROOT (see section 5), then support for the ROOT output format in addition to the others is available via the RootEventFileReader class. This class is derived from EventFileReader and implements an identical user interface. Making the replacement EventFileReader → RootEventFileReader on lines 4 and 12 of listing 14 will yield a program that is capable of reading events from any MARLEY output file. Compiling this modified program requires linking to the shared library containing the MARLEY interface to ROOT. The build command given in section 5.5.1 will automatically include the necessary compiler flags for linking to ROOT if MARLEY was successfully built with ROOT support.

Accessing information from an event record
As discussed in section 2, MARLEY conceives of each scattering event as consisting of a primary 2 → 2 reaction possibly followed by nuclear de-excitations. The particle content of the primary reaction may be written in the general form where particles a, b, c, and d, are respectively termed the projectile, target, ejectile, and residue. Where a mass difference exists, MARLEY chooses the projectile (ejectile) to be the lighter of the two particles in the initial (final) state. Otherwise, the choice is arbitrary. All four-vector components stored in a MARLEY event record are expressed in the laboratory frame, i.e., the rest frame of the target. A full MARLEY event record is represented in the code itself by the marley::Event class. Instances of this class own two vectors of pointers to marley::Particle objects, with one (initial particles ) describing the initial state and the other (final particles ) describing the final state. Within the initial particles vector, the first and second elements always correspond to the projectile and target, respectively. A similar ordering (ejectile followed by residue) is enforced for the first two elements of the final particles vector. If present, elements of final particles beyond the second correspond to nuclear de-excitation products, which are listed in the order that they were emitted. If nuclear de-excitations were enabled in the simulation (as they are by default, see section 4.1.7), then the Particle object for the residue stores its properties after it has reached the ground state.
The Event class provides the following member functions for accessing Particle objects stored in the event record: projectile(), target(), ejectile(), residue() Returns a reference to a specific particle involved in the primary 2 → 2 reaction get initial particles() Returns a reference to the owned vector of initial particles get final particles() Returns a reference to the owned vector of final particles initial particle count() Returns the number of initial-state particles in the event final particle count() Returns the number of final-state particles in the event initial particle(size t idx) Returns a constant reference to the initial-state particle at position idx final particle(size t idx) Returns a constant reference to the final-state particle at position idx The Particle class stores a 4-momentum (represented as a C-style array of type double[4]) together with variables representing particle properties. It implements the following member functions for data retrieval: All of these functions return a value of type double except for pdg code(), which returns an int.
In addition to Particle objects, the Event class also stores information about the state of the residue immediately after its creation by the prompt 2 → 2 reaction. Access to this information is provided by these member functions: Ex() Initial residue excitation energy (double, MeV), measured with respect to its own ground state twoJ() Two times the initial residue spin (int) parity() Initial residue parity (marley::Parity) The Parity object returned by Event::parity() provides a type-safe representation of a parity value ±1. It may be converted to a boolean (true ↔ +1, false ↔ −1) or integer representation via an explicit cast. A Parity object may also be converted to a char value (+ ↔ +1, -↔ −1) via the member function to char(). Full technical documentation for the Event, Particle, and Parity classes is available on the official MARLEY website [149]. Similar HTML documentation may be automatically generated for offline use from the MARLEY source code using Doxygen [169]. After Doxygen has been installed, one may build the MARLEY API documentation files by executing make doxygen from within the build/ folder. The generated documentation may then be viewed using any web browser by opening the file docs/ build/html/doxygen/index.html

Executable example: marprint
The marprint program (examples/executables/marprint.cc) included with MARLEY provides a detailed example of how the C++ analysis API may be used in compiled code. After instantiating either an EventFileReader or a RootEventFileReader object (depending on whether MARLEY was built against ROOT), the program prints all information stored in each event to standard output using a human-readable format. Usage examples for many of the functions listed in section 7.3.2 are provided within the marprint.cc source file. After sourcing the setup marley.sh script, one may build the marprint example program via the commands The marprint executable takes the names of one or more files containing MARLEY events as commandline arguments. For example, one may print all of the events stored in the files /home/events1.ascii and /home/events2.json by executing marprint / home / events1 . ascii / home / events2 . json Listing 15 shows the printout generated by marprint while parsing the first event given in the example ASCII-format output file shown above (listing 11). *** Event 0 has 2 initial particles and 4 final particles . *** The residual nucleus initially had excitation energy 3.79748 MeV and spin -parity 1+ Initial particles particle with PDG code = 12 has total energy 10 MeV , 3 -momentum = (0 MeV , 0 MeV , 10 MeV ) , mass = 0 MeV , and charge = 0 times the proton charge . particle with PDG code = 1000180400 has total energy 37224.

Example ROOT macros
As a second example of the MARLEY analysis API, several ROOT macros, i.e., programs intended to be run using the interactive C++ interpreter distributed with ROOT, are provided in the examples/macros/ folder. Although RootEventFileReader and the other MARLEY classes are fully compatible with version 6 of ROOT, version 5 lacks support for C++11 features, which are used extensively by MARLEY. To work around this problem, MARLEY provides a class called MacroEvent FileReader which is usable within macros written for both versions 5 and 6 of ROOT. This class provides the same interface to MARLEY events as RootEventFileReader, but the use of C++ syntax that is incompatible with ROOT 5 has been hidden from the interpreter.
In order to correctly interact with MARLEY classes, ROOT requires precompiled dictionaries to be loaded at runtime. These are stored in the MARLEY ROOT shared library. To avoid the need for users to manually load the MARLEY class dictionaries at the start of each ROOT interpreter session, a startup script called mroot is placed in the build/ folder whenever MARLEY is successfully built against ROOT. After sourcing the setup marley.sh setup script (see section 5.3), issuing the command mroot will start the interactive ROOT interpreter and automatically load the MARLEY class dictionaries. For version 5 of ROOT, only the Event, Particle, Parity, and MacroEventFileReader classes (all defined within the marley namespace) may be used within an mroot session. For ROOT 6, all MARLEY classes will be available if the appropriate header files are loaded via #include statements.
The README.md file in the examples/macros/ folder gives brief descriptions of each of the example ROOT macros together with usage instructions. Users are encouraged to adopt these macros as a starting point for implementing their own calculations.

"Flat" ROOT files: the marsum utility
Although they are expected to be suitable for many applications, the standard ROOT-format output files generated by MARLEY (see section 7.2.4) have two important limitations: (1) They are only fully readable in environments in which both ROOT and MARLEY are installed, and (2) Each Event object to be analyzed must be loaded from disk in its entirety.
To address these limitations, MARLEY provides a command-line utility called marsum, which takes as input one or more files containing MARLEY events stored in any of the standard output formats. The marsum program creates a new file in which many quantities of interest from the input events have been saved as individual branches of a ROOT TTree. Following a successful build of MARLEY with ROOT support (see section 5), the marsum executable will be present in the build/ folder. Assuming that the setup marley.sh script has already been sourced, one may execute the command marsum myout . root / home / events1 . root / home / events2 . ascii to create a new file called myout.root in the working directory. This file will contain a single ROOT TTree named mst (for "MARLEY summary tree" 31 ) with one entry for each event present in the two input files (/home/events1.root and /home/events2.ascii). The mst TTree will contain the following branches: The "flat" file created in this way will be readable by ROOT without the need to load the MARLEY class dictionaries. Individual branches may also be loaded and plotted (e.g., via the TTree::Draw member function) without the need to manipulate the event as a whole.

Converting event distributions to physics quantities
The theoretical predictions made by an event generator like MARLEY may most usefully be compared to competing calculations and experimental data in the form of cross sections and event rates. To see how these quantities may be obtained from a set of simulated events, first note that the expression given in eq. (73) for the probability density P (E ν ) of the energy E ν of a reacting neutrino may be rewritten in the form is the total neutrino flux and is the (abundance-weighted) flux-averaged total cross section. Let x denote an arbitrary, continuouslydistributed observable that is computable from a MARLEY event record. Then, for a single event, the probability P j that x falls within the jth bin x ∈ [x j , x j+1 ) is given by where is the conditional probability density of x at fixed neutrino energy, dσ(E ν )/dx is the total differential cross section with respect to x, and is the flux-averaged total differential cross section. A Monte Carlo estimator for the average value of this quantity in the jth bin may be obtained via where ∆x j ≡ x j+1 − x j is the bin width and f j = n j /N is the ratio of the n j events that fall within the jth bin to the total number N of simulated events. Recognizing that n j follows a binomial distribution allows for an estimate of the associated Monte Carlo statistical uncertainty (standard deviation) via The results from eqs. (98) and (99) may readily be extended to multiple dimensions. For a Monte Carlo calculation of an n-dimensional differential cross section, simply let ∆x j denote the product of the n widths of the jth bin in the n-dimensional phase space. For a discrete observable, similar expressions may be used with the bin width ∆x j omitted. The flux-averaged partial cross section σ j to produce events fulfilling a criterion j (e.g., involving emission of a single neutron) is estimated from the simulation results via σ j ≈ σ f j , where f j = n j /N is the fraction of simulated events satisfying the criterion. The statistical uncertainty of this estimator is approximated by the expression on the right-hand side of eq. (99) with the substitution ∆x j → 1. Figure 3 shows the results of two example calculations of physics observables performed using MARLEY events and the procedure outlined in section 7.4. The histograms shown in the left and right panels of fig. 3 were computed using independent simulations of N = 2 × 10 6 events each. To produce the left-hand plot, coherent elastic neutrino-nucleus scattering (CEvNS) on 40 Ar was simulated forν µ produced by µ + decays at rest (see section 6.5.4). A histogram describing the distribution of the kinetic energy T f of the recoiling final-state nucleus was prepared from the simulated events using a uniform bin width ∆T f = 1.5 keV. The event counts n j from each bin were renormalized according to eq. (98) to yield a Monte Carlo estimator for the mean value of the flux-averaged differential cross section in the jth bin:

Examples
Here σ = 26.69 × 10 −40 cm 2 / 40 Ar is the MARLEY prediction for the CEvNS flux-averaged total cross section. Equation (100) gives the expression used to obtain the content of each histogram bin shown in the left-hand plot of fig. 3. Based on eq. (99), an estimate of the statistical uncertainty associated with each bin was also calculated but is small on the scale of the plot. A similar procedure was used to obtain the distribution shown in the right-hand plot of fig. 3, but the quantity of interest is T e , the kinetic energy of the electron produced in the charged-current reaction 40 Ar(ν e , e − ) 40 K * . The simulation of this process was carried out using the reaction input file ve40ArCC Bhattacharya1998.react (see section 6.4.1). The incident neutrino spectrum was defined using the Fermi-Dirac source shown in listing 5, which was treated as a toy model of the time-integrated ν e flux at Earth produced by a core-collapse supernova. The flux-averaged differential cross section was converted into a differential event rate via where a total time-integrated flux 32 of Φ = 1.0 × 10 11 ν e / cm 2 was assumed. There are n targets = 1.5 × 10 32 atoms in 10 kt of pure 40 Ar. 32 The quoted value is a rough estimate for a core-collapse supernova at 10 kpc from Earth.

Energy-dependent total cross sections: mardumpxs
One of the observables predicted by MARLEY that allows the most straightforward comparison to other theoretical models is the (abundance-weighted) total cross section as a function of neutrino energy σ(E v ). While a Monte Carlo estimate of this quantity may be computed from generated events 33 using the approach described in section 7.4, σ(E ν ) is exactly calculable given only the neutrino target composition (see section 4.1.6) and the information stored in the reaction input file(s) (see section 4.1.4). For the convenience of users, an example C++ program called mardumpxs (examples/executables/mardumpxs.cc) is provided with MARLEY that produces tables of σ(E ν ). Assuming that the setup marley.sh script has already been sourced (see section 5.3), one may build mardumpxs by executing the commands cd $ { MARLEY }/ build make mardumpxs The mardumpxs executable takes the name of an output file followed by the name of a MARLEY job configuration file (see section 6) as command-line arguments. For example, the command mardumpxs xsec_table . txt config . js will write a table of σ(E ν ) values to the output file xsec table.txt after configuring MARLEY according to the settings given in config.js. Each line of the output file will have the format

KE XSec
where KE is the projectile kinetic energy 34 in MeV and XSec is the abundance-weighted total reaction cross section σ(E ν ) (see section 4.1.6) in 10 −42 cm 2 / atom. By default, mardumpxs assumes a ν e projectile and produces a table of 10,000 cross section values using a regularly-spaced energy grid between 0 and 100 MeV. This behavior may be altered using mardumpxs-specific keys in the job configuration file, as described in section 6.7.8.

Interfacing MARLEY with external tools
While it is hoped that MARLEY's capabilities as a standalone software package will be beneficial to the low-energy neutrino physics community, interfacing the generator with external codes has the potential to greatly extend its usefulness. This is particularly true for neutrino detection experiments, which typically rely on end-to-end simulations of neutrino production, neutrino interactions in and around the detector, final-state particle propagation, and the detector electronics response in order to interpret their measurements. For some applications, simply passing MARLEY output files as input to the next stage of a multi-step simulation may be satisfactory. In other contexts, direct calls to MARLEY functions by a client code may be more appropriate.
Section 8.1 presents the recommended approach to generating MARLEY events within an external C++ application, which involves use of the marley::JSONConfig and marley::Generator classes. An example application of this kind, marg4, is discussed in section 8.2. The marg4 program produces MARLEY events and tracks the final-state particles through a simple geometry using the popular Geant4 [100, 101] particle transport code. Section 8.3 then briefly discusses the MARLEY interface included in the LArSoft [171][172][173][174] toolkit.

C++ event generation API
The core functionality of MARLEY is encapsulated for use by external applications in the form of the marley::Generator class. While a Generator object may be instantiated and configured without reference to a file, doing so is not recommended in most situations. Instead, users are encouraged to construct Generator objects indirectly by means of the create generator member function of the marley::JSONConfig class. The constructor of JSONConfig takes a single string argument containing the name of a MARLEY job configuration file. The contents of this file are parsed and used to initialize a Generator object during a subsequent call to JSONConfig::create generator. With the exception of the parameters described in section 6.6 (which are unique to the marley command-line executable), all configuration options listed in section 6 will be recognized and respected by the JSONConfig class. Once a fully-initialized Generator object has been created, a single neutrino scattering event may be simulated and returned as a marley::Event object by calling the create event member function. A new event will be generated with each call to this function.
Listing 16 shows the source code for examples/executables/minimal/evgen.cc, an example C++ program that uses the recommended event generation API. On line 13, a JSONConfig object called cfg is constructed using settings from the job configuration file /home/config.js. The settings parsed by cfg are then used to create a Generator object called gen on line 14. A simple event loop is defined on lines 16-19 and iterates until ten events have been generated (line 17) and printed in ASCII format to standard output (line 18). The evgen.cc example program may most easily be compiled by using the marley-config script as described in section 5.5.1.
1 // Standard library includes 2 # include < iostream > 3 4 // MARLEY includes 5 # include " marley / Event . hh " 6 # include " marley / Generator . hh " 7 # include " marley / JSONConfig . hh " 8 9 constexpr int NUM_EVENTS = 10; In cases where handling of ROOT-dependent configuration file parameters (see, e.g., section 6.5.7) is desired, the RootJSONConfig class should be used instead of JSONConfig . This amounts to making the replacement JSONConfig → RootJSONConfig on lines 7 and 13 of listing 16. Assuming that MARLEY has been built with ROOT support, the compilation command given in section 5.5.1 will automatically link to the required ROOT libraries.

Interfacing MARLEY with Geant4
The C++-based Geant4 software framework provides a large suite of tools for simulating particle propagation through matter. In the main function of a typical Geant4 application, a G4RunManager object is constructed and used to drive the simulation. Before the simulation can begin, the G4RunManager must be initialized with pointers to three objects, each of which is derived from a distinct abstract base class defined by Geant4. The three required classes used to initialize the run manager are G4VUserDetectorConstruction Defines the geometry and material composition of the world volume (and zero or more subvolumes) through which the simulated particles will be tracked G4VUserPhysicsList Defines the particle species and physics processes to be included in the simulation G4VUserPrimaryGeneratorAction Defines a member function, GeneratePrimaries, which will be used to populate each new event with the starting locations, momenta, etc. of all primary particles to be tracked In addition to these required classes, pointers to objects that instantiate optional user action classes may also be registered with the run manager. These allow user-defined functions to be executed at various stages of the simulation. Further details about general Geant4 application development are available in ref. [175].

The marg4 Geant4 application
The folder examples/marg4/ contains the source code for marg4, an example Geant4 application that uses MARLEY events as a source of primary particles. The world volume is defined by the DetectorConstruction class and consists of a single uniform sphere of liquid argon with a radius of 10 m and centered on the origin. A built-in Geant4 physics list suitable for MeV-scale particle transport, QGSP BIC HP, is constructed using a factory method in examples/marg4/src/marg4.cc. As a trivial example of a user action, an EventAction object is used to print the current event count to standard output at the beginning of every hundredth event.
An example of a direct interface between MARLEY and Geant4 is provided by the MarleyPrimary GeneratorAction class, which is derived from G4VUserPrimaryGeneratorAction. The constructor of this class takes a single std::string argument containing the name of a MARLEY job configuration file. This file name is used to create either a JSONConfig or a RootJSONConfig object, with the latter being chosen if the USE ROOT preprocessor macro is defined (see section 5.5). The member variable marley generator , which is a Generator object, is then initialized using the parameters from the configuration file.
Listing 17 shows the GeneratePrimaries member function defined by the MarleyPrimary GeneratorAction class. This function is called once during initialization of each Geant4 event.
On line 4, a new G4PrimaryVertex object is created at the spacetime origin. All primary particles which are associated with it will begin their Geant4 trajectories at the same 4-position. After a single marley::Event object is created on listing 17, each of its final-state particles is converted into a new G4PrimaryParticle by the loop defined on lines 14-28. Line 27 associates each fully-initialized G4PrimaryParticle with the primary vertex defined previously, and line 33 adds the completed primary vertex to the current G4Event object. Propagation of the primary particles obtained from the MARLEY event is simulated by Geant4 after the GeneratePrimaries function returns.

Building and running marg4
Building the example marg4 program requires the geant4-config utility (included as part of a standard Geant4 installation) to be present on the system's executable search path. Additionally, use of the QGSP BIC HP physics list mentioned above requires installation of the data files belonging to the Geant4 Neutron Data Library (G4NDL) [176,177]. These files, which are available for download from the Geant4 website (https://geant4.web.cern.ch/support/download), are required for high-precision (HP) tracking of low-energy neutrons by Geant4. The use of Geant4's HP treatment of neutron transport is strongly recommended for propagation of neutrino-induced neutrons from MARLEY events.
Assuming that the setup marley.sh script has already been sourced (see section 5.3), the marg4 executable may be built against Geant4 via the commands cd $ { MARLEY }/ build make marg4 The marg4 executable takes the number of desired events followed by the name of a MARLEY job configuration file as its command-line arguments. For example, invoking the program via marg4 500 / home / myconfig . js will simulate 500 Geant4 events. Each of these events will contain a single primary vertex populated with the final-state particles from one MARLEY event. The MARLEY events will be generated using the job configuration given in the file /home/myconfig.js.
The marg4 program is provided as a simple usage example for the API described in section 8.1. As such, it does not produce any output other than logging messages from both MARLEY and Geant4. Interested users are encouraged to copy and modify the marg4 source code to meet their needs. Information from the marley::Event objects themselves may be accessed using the member functions described in section 7.3.2. The functions needed to extract quantities of interest from the Geant4 simulation are documented in ref. [175].

LArSoft interface
The liquid argon software toolkit (LArSoft) provides a flexible simulation, reconstruction, and analysis framework designed for liquid argon time projection chamber (LArTPC) experiments.
The LArSoft source code is hosted on GitHub and split into multiple respositories, each of which provides a particular type of functionality. The larsim respository (https://github.com/ LArSoft/larsim) includes, among other things, interfaces to external physics event generators for neutrino interactions (e.g., GENIE) and for other processes (e.g., the cosmic-ray generators CORSIKA [185] and CRY [186]). Version 6.04.00 of LArSoft was the first 35 to include a direct larsim interface to MARLEY contributed by the present author. 36 Ever since the initial version of the interface was added, MARLEY has been included as part of the standard LArSoft distribution. A detailed description of the LAr-Soft interface to MARLEY, which is more sophisticated than the examples given in section 8.1 and section 8.2, is beyond the present scope. However, brief descriptions of the relevant C++ classes defined in the larsim source code (see the larsim/EventGenerator/MARLEY/ subfolder of the larsim repository) are given below. All four of these classes are defined within the evgen namespace used by LArSoft for event generation.
ActiveVolumeVertexSampler Used to sample neutrino vertex positions uniformly over the active volume(s) of a detector. The approach used by this class is only suitable for simulations in which (1) the detector is uniformly illuminated by the incident neutrino flux, and (2) neutrino interactions occurring outside the detector active volume(s) are not of interest.
MarleyGen Creates neutrino scattering events in the native LArSoft format using the Active VolumeVertexSampler and MARLEYHelper classes MARLEYHelper Implements the low-level interface between MARLEY and LArSoft. At the beginning of a simulation job, this class initializes a marley::Generator object by converting a LArSoft configuration given in the Fermilab Hierarchical Configuration Language (FHiCL) [188] to the JSON-like format (see section 6) used by MARLEY. With the exception of the "executable settings" described in section 6.6, all MARLEY job configuration file parameters are available for use via FHiCL.
This class also provides a member function (create MCTruth) which generates a MARLEY event and stores a representation of it in an instance of the simb::MCTruth class. The simb ::MCTruth class is used by LArSoft as a generator-agnostic event record.
MarleyTimeGen Similar to MarleyGen, but provides experimental support for generating events using a time-dependent neutrino spectrum Further technical details about these and other LArSoft classes are available in the LArSoft Doxygen documentation (https://nusoft.fnal.gov/larsoft/doxsvn/html/index.html). 35 Ref. [187] incorrectly identifies the initial LArSoft version as 6.03.00. This number refers instead to the corresponding version of the larsim subpackage. 36 The first MARLEY version to be included in LArSoft was 0.9.5, an August 2016 public beta release.

Prospects for future development
This work describes a new Monte Carlo event generator, MARLEY, suitable for simulating tensof-MeV neutrino scattering on complex nuclei. The release of MARLEY 1.2.0, the first version of the code to be documented with the present level of detail, represents a significant milestone. However, active, open-source development of MARLEY is ongoing, and contributions from the community are encouraged. To coordinate development efforts, those interested in contributing improvements to MARLEY are asked to contact the author at the earliest opportunity.
Potential avenues for future development of MARLEY include the following: • Preparation of additional reaction input files. In addition to new scattering modes (inelastic NC andν e CC) for 40 Ar, the addition of input files for several other nuclides will likely be of immediate interest. These include but are not limited to 1. 12 C and 16 O, both of which make small but non-negligible contributions to the total event rate in hydrocarbon (e.g., NOvA [189]) and water Cherenkov (e.g., Super-Kamiokande [190]) supernova neutrino detectors 2. The stable lead isotopes (especially 208 Pb), which serve as the target material for the HALO [47,48] detector 3. The stable iron and copper isotopes (especially 56 Fe and 63 Cu), which, together with lead, are under study by the COHERENT [35] experiment in an effort to understand low-energy neutrino-induced neutron production • Inclusion of forbidden transitions in the MARLEY neutrino-nucleus cross section model. These are currently neglected via the allowed approximation. Recent theoretical calculations using a Continuum Random Phase Approximation (CRPA) approach [191,192] indicate that the forbidden contributions can have a significant impact on both the total cross section and the kinematic distributions of the outgoing lepton.
• Addition of new keys to the job configuration file format to allow variations of model parameters, e.g., those used by the nuclear optical model. With some related enhancements to other parts of the code, e.g., storage of the full de-excitation history in the MARLEY event record, these variations could be used to assess theoretical uncertainties on MARLEY predictions via event reweighting. An approach of this kind is commonly used by accelerator neutrino experiments. Documentation of the GENIE event reweighting framework, for example, is given in ref. [170, ch. 9, pp. 129-154].
• Improvements to the MARLEY treatment of nuclear de-excitations. These could include -Realistic angular distributions for the de-excitation products (instead of the isotropic emission that is currently assumed) -Storage of particle emission times in the event record. Related updates could also be made to the nuclear structure data format to allow measured discrete level lifetimes to be listed.
-Competition of internal conversion with γ-ray emission -Pre-equilibrium emission of nuclear fragments -Modeling of neutrino-induced fission • Creation of new subclasses of marley::Reaction to simulate events for projectiles other than neutrinos. Processes likely to be of interest include electron-nucleus scattering and beyond the Standard Model reactions like nuclear absorption of MeV-scale dark matter [193].
• Implementation of a tool to facilitate comparisons of MARLEY model predictions with cross section measurements from low-energy neutrino experiments. This could potentially be accomplished by adding low-energy datasets and a MARLEY interface to the NUISANCE [194] framework used by the accelerator neutrino community.
• Simulation of events within a detector geometry in which both the target material and the incident neutrino flux are spatially non-uniform. This may most easily be achieved by interfacing MARLEY with flux and geometry navigation drivers present in an existing generator (see ref. [170, sec. 6] for a description of GENIE's). A "community based" version of these tools, independent of any particular generator, has also been proposed [195].

Acknowledgements
I am grateful to Myung-Ki Cheoun for providing the QRPA Gamow-Teller matrix elements which are used to compute 40 Ar(ν e , e − ) 40 K * cross sections at high excitation energies (see section 6.4.1 and ref. [108]). I also thank the authors of the TALYS code for sharing their nuclear level data files under the GNU General Public License.
Vishvas Pandey provided helpful feedback on a draft of this paper, and I thank Robert Svoboda, Ramona Vogt, and Michael Mulhearn for their comments on the PhD thesis [187] which served as the first formal description of MARLEY. I am also indebted to Sam Hedges and Erin Conley for their tests of the MARLEY 1. A. Nuclear structure data file format The MARLEY nuclear structure data files mentioned in section 4.1.5 contain tables of nuclear energy levels and γ-ray branching ratios for one or more nuclides. Each table follows a whitespacedelimited format that begins with the header Z A num_levels in which Z is the proton number, A is the mass number, and num_levels is the number of tabulated nuclear levels. The header is followed by data blocks for each level in order of increasing excitation energy. Each block begins with the line

Ex twoJ Pi num_gammas
where Ex is the level excitation energy (MeV), twoJ is two times the level spin (the factor of two allows half-integer spins to be represented by an integer), and Pi is the parity (denoted by the character + or -). The data block for a level terminates with a set of lines describing each available γ-ray transition (for a total of num_gammas transitions). Each of these lines has the format E_gamma RI Lf_index where E_gamma is the energy of the emitted γ-ray (MeV) and RI is the relative intensity of the transition. Although the relative intensities in the official MARLEY structure data files are normalized to sum to unity, this is not required. The parameter Lf_index gives the position of the final nuclear level accessed by the transition, with Lf_index = 0 corresponding to the first listed level (presumably the ground state). Listing 18 shows the table for 43 Cl which appears in the structure data file data/structure/ Cl.dat included with MARLEY. Five nuclear levels appear in the table, each (apart from the ground state) having one listed γ-ray transition. Although the lines giving the nuclide header, level descriptions, and γ-ray descriptions are shown here with different indentations to improve readability, this is not required by the file format.

B. Reaction input file format
The reaction input files mentioned in section 6.4 provide information needed for MARLEY to configure cross section calculations according to eq. (5) (for neutrino-nucleus scattering) or eq. (55) (for neutrino-electron elastic scattering). Any line of a reaction input file that begins with a # character will be treated as a comment and ignored by the parser. The first line of the file that is not a comment contains a header of the form

ProcessType NucPDG
where ProcessType is an integer code representing the reaction mode and NucPDG is a nuclear PDG code (see section 7.1) identifying the target nuclide involved in the initial state. Table B.6 lists the allowed values of ProcessType in the second column together with their corresponding reaction modes and isospin operators. In the first column, the symbol A is used as a stand-in for any target nucleus. The third column lists the elements of the enumerated type used in the source code to represent each reaction mode.
For reaction input files describing neutrino-electron elastic scattering (ProcessType = 3), only the header shown above is required to be present. If the user wishes to enable simulation of this
process for multiple atomic targets, additional values of NucPDG may be included on subsequent lines. Listing 19 shows an example reaction input file which instructs MARLEY to simulate neutrinoelectron elastic scattering for the three stable isotopes of argon: 36 Ar, 38 Ar, and 40 Ar. For neutrino-nucleus reaction modes, each line of the input file following the header is used to specify the value of a reduced matrix element describing a transition to a particular final nuclear level. These lines are whitespace-delimited and have the format where Ex is the excitation energy of the final level (MeV, measured with respect to the ground state of the residual nucleus), B is the corresponding value of either B(F) or B(GT), and type is an integer code that represents whether B should be interpreted as a value of B(F) (type = 0) or of B(GT) (type = 1). The matrix elements must be listed in order of increasing Ex in the reaction input file. If this is not the case, then MARLEY will halt with the error message [ ERROR ]: Invalid reaction dataset . Level energies must be unique and must be given in ascending order .
when the file is used. In most cases, the matrix element values given in the reaction input file should be computed using the full expressions shown in eq. (6) and eq. (7). This includes, e.g., any assumed value of the axial-vector coupling constant g A . The sole exception is the Fermi matrix element for neutralcurrent scattering, which should be listed with a factor of Q 2 W /4 removed. For an NC transition to the ground state of a spin-zero nucleus, B(F ) = g 2 V Q 2 W /4, but B = g 2 V = 1. Listing 20 shows an excerpt from the reaction input file data/react/ve40ArCC Bhattacharya 1998.react, which tabulates nuclear matrix elements for CC neutrino-argon scattering. Two B(GT) values (0.90 and 1.50) are given, both for final states between 2 and 3 MeV above the 40 K ground state.