Das ist der HAMMER: Consistent new physics interpretations of semileptonic decays

Precise measurements of $b\to c\tau\bar\nu$ decays require large resource-intensive Monte Carlo (MC) samples, which incorporate detailed simulations of detector responses and physics backgrounds. Extracted parameters may be highly sensitive to the underlying theoretical models used in the MC generation. Because new physics (NP) can alter decay distributions and acceptances, the standard practice of fitting NP Wilson coefficients to SM-based measurements of the $R(D^{(*)})$ ratios can be biased. The newly developed HAMMER software tool enables efficient reweighting of MC samples to arbitrary NP scenarios or to any hadronic matrix elements. We demonstrate how HAMMER allows avoidance of biases through self-consistent fits directly to the NP Wilson coefficients. We also present example analyses that demonstrate the sizeable biases that can otherwise occur from naive NP interpretations of SM-based measurements. The HAMMER library is presently interfaced with several existing experimental analysis frameworks and we provide an overview of its structure.


Introduction
Precision analyses of semileptonic b-hadron decays typically rely on detailed numerical Monte Carlo (MC) simulations of detector responses and acceptances. Combined with the underlying theoretical models, these simulations provide MC templates that may be used in fits, to translate experimental yields into theoretically well-defined parameters. This translation though can become sensitive to the template and its underlying theoretical model, introducing biases whenever there is a mismatch between the theoretical assumptions used to measure a parameter and subsequent theoretical interpretations of the data.
Such biases are known to arise in the analyses of semileptonic decays of b hadrons, in particular, for the measurements of the CKM element |V cb |, and the ratio of semitauonic vs. semileptonic decays to light leptons, where H b,c denote b-and c-flavor hadrons. To avoid this, the size of these biases need to be either carefully controlled when experiments quote their results arXiv:2002.00020v1 [hep-ph] 31 Jan 2020 by reversing detector effects, or they can be avoided by using dedicated MC samples for each theoretical model the measurement is confronted with. In this paper we present the newly developed tool, Hammer (Helicity Amplitude Module for Matrix Element Reweighting), designed expressly for the latter purpose. Semitauonic b hadron decays have long been known to be sensitive to new physics [1][2][3][4][5][6][7], and were first constrained at LEP [8]. At present, the measurements of the R(D ( * ) ) ratios show about a 3σ tension with SM predictions, when the D and D * modes are combined [9]. In the future, much more precise measurements of semitauonic decays are expected, not only for the B → D ( * ) τν channels, but also for the not yet studied decay modes, Λ b → Λ c τν, B s → D ( * ) s τν, as well as involving excited charm hadrons in the final state.
All existing measurements of R(D ( * ) ) rely heavily on large MC simulations to optimize selections, provide fit templates in discriminating kinematic observables, and to model resolution effects and acceptances. Both the τ and the charm hadrons have short lifetimes and decay near the interaction point and measurements rely on reconstruction of the ensuing decay cascades. To reconstruct the decay products, often complex phase space cuts and detector efficiency dependencies come into play, and the measurement of the full decay kinematics is impossible due to the presence of multiple neutrinos. In addition, depending on the final state, a significant downfeed with similar experimental signatures from misreconstructed excited charm hadron states can be present. Isolation of semitauonic decays from other background processes and the light-lepton final states, then requires precise predictions for the kinematics of the signal semitauonic decay. 1 Often the limited size of the available simulated samples, required to account for all these effects, constitutes a dominant uncertainty of the measurements, see e.g. [10][11][12].
In the literature on the R(D ( * ) ) anomaly, it has become standard practice to reinterpret the experimental values of R(D ( * ) ) in terms of NP Wilson coefficients, even though all current ratio measurements were determined assuming the SM nature of semitauonic decays. However, NP couplings generically alter decay distributions and acceptances. Therefore, they modify the signal and possibly background MC templates used in the extraction, and thus affect the measured values of R(D ( * ) ). This may introduce biases in NP interpretations: preferred regions and best-fit points for the Wil-son coefficients can be incorrect; an instructive example of this is provided in Sec 2.3.
Consistent interpretations of the data with NP incorporated requires dedicated MC samples, ideally for each NP coupling value considered, which would permit directly fitting for the NP Wilson coefficients. This approach is sometimes referred to as 'forward-folding', and is naively a computationally prohibitively expensive endeavour. Such a program is further complicated because none of the MC generators current used by the experiments incorporate generic NP effects, nor do they include state-of-the-art treatments of hadronic matrix elements.
In this paper we present a new software tool, Hammer, that provides a solution to these problems: A fast and efficient means to reweight large MC samples to any desired NP, or to any description of the hadronic matrix elements. Hammer makes use of efficient amplitude-level and tensorial calculation strategies, and is designed to interface with existing experimental analysis frameworks, providing detailed control over which NP or hadronic descriptions should be considered. The desired reweighting can be implemented either in the event weights or in histograms of experimentally reconstructed quantities (both further discussed in Sec. 3). The only required input are the event-level truth-fourmomenta of existing MC samples. Either the event weights and/or histogram predictions may be used, e.g., to generate likelihood functions for experimental fits. Some of the main ideas of Hammer were previously outlined in Refs. [13,14].
In Sec. 2 we demonstrate the capabilities of Hammer by performing binned likelihood fits on mock measured and simulated data sets, that are created using the Hammer library, and corrected using an approximate detector response. In Sec. 3 a brief overview of the Hammer library and its capabilities are given. Section 4 provides a summary of our findings. Finally, Appendix A provides a detailed overview of the Hammer application programming interface.

New physics analyses
We consider two different analysis scenarios: 1. In order to explore what biases may arise in phenomenological studies if NP is present in Nature, we perform an illustrative R(D ( * ) ) toy measurement. This involves carrying out SM fits to mockups of measured data sets, that are generated for several different NP models. The recovered R(D ( * ) ) values are then compared to their actual NP values.
2. To demonstrate using a forward-folded analysis to assess NP effects without biases, we carry out fits to (combinations of) NP Wilson coefficients themselves, with either the SM or other NP present in the mock measured data sets.
The setting of these analyses is a B-factory-type environment. We focus on leptonic τ decays, but the procedures and results in this work are equally adaptable to the LHCb environment, and other τ decay modes or observables. We emphasize that the derived sensitivities shown below are not intended to illustrate projections for actual experimental sensitivities per se. Such studies are better carried out by the experiments themselves.

MC sample
The input Monte Carlo sample used for our demonstration comprises four distinct sets of 10 5 events: one for each of the two signal cascades B → D(τ → eνν)ν, B → (D * → Dπ)(τ → eνν)ν and for the two background processes, B → Deν and B → (D * → Dπ)eν. These are generated with EvtGen R01-07-00 [15], using the Belle II beam energies of 7 GeV and 4 GeV. The second B meson decay, often used for identifying or 'tagging' the BB event and constraining its kinematic properties, are not included in the current analysis for simplicity, but can be incorporated in a Hammer analysis straightforwardly.
In each cascade, the b → clν vertices are generated equidistributed in phase space ("pure phase space") instead of using SM distributions. This reduces the statistical uncertainties that can otherwise arise from reweighting regions of phase space that are undersampled in the SM to NP scenarios in which they are not. 2

Reweighting and fitting analysis
Hammer is used to reweight the MC samples into twodimensional 'NP generalized' histograms (see Sec. 3), with respect to the reconstructed observables |p * | and m 2 miss , the light lepton momentum in the B rest frame and the total missing invariant mass of all neutrinos, respectively. Both variables are well-suited for separating signal from background decays involving only light leptons. In the cascade process of the leptonic τ decay in 2 For an actual experimental analysis one would instead use Hammer to reweight SM MC samples. The correct statistical uncertainty of the reweighting can be incorporated, using weight squared uncertainties computed by the library. This information could be used, e.g., to adaptively generate additional pure phase space MC in undersampled regions. B → D ( * ) τ ν, the signal lepton carries less momentum than the lepton from prompt B → D ( * ) ν decays. Similarly, the missing invariant mass of B → D ( * ) ν decays peaks strongly near m 2 ν 0, in contrast to B → D ( * ) τ ν in which the multiple neutrinos in the final state permit large values of m 2 miss . The B → D ( * ) processes are reweighted to the BLPR form factor parametrization [16], which includes predictions for NP hadronic matrix elements using HQET [17][18][19][20] Charged particles are required to fall in the Belle II angular acceptance of 20 • and 150 • , and leptons are required to have a minimum kinetic energy of 300 MeV in the laboratory frame. An additional event weight is included to account for the slow pion reconstruction efficiencies from the D * → Dπ decay, based on an approximate fit to the pion reconstruction efficiency curve from BaBar data [10,21]. The analysis assumes that the second tagging B meson decay was reconstructed in hadronic modes, such that its fourmomentum, p Btag , is accessible. In conjunction with the known beam four-momentum p e + e − , the missing invariant mass can then be reconstructed as m 2 miss ≡ (p e + e − −p Btag −p D ( * ) −p ) 2 , and the four-momentum of the reconstructed lepton can be boosted into the signal B rest frame. A Gaussian smearing is added to the truth level m 2 miss with a width of 0.5 GeV 2 to account for detector resolution and tagging-B reconstruction. No additional correction is applied to |p * |. Higher dimensional histograms including the reconstructed q 2 and the D * → Dπ helicity angle may also be incorporated, but are omitted here for simplicity.
Hammer can be used to efficiently compute histograms for any given NP choice. The basis of NP operators is defined in Table 1, with respect to the Lagrangian where Γ X(Y ) is any Dirac matrix and c XY is a Wilson coefficient. We shall generally write explicit Wilson coefficients as c XY = S qXlY , V qXlY , T qXlY , where the S, V , T denotes the Lorentz structure, and X, Y = L, R denotes the chirality. In this simplified analysis, we assume that NP only affects the b → cτ ν decays, and not the light-lepton modes.
In order to carry out Wilson coefficient fits, we wrap the Hammer application programming interface with a gammaCombo [22] compatible class. This allows one to use Hammer's efficient reweighting of histogram bins to generate the relevant quantities required to calculate a likelihood function. We then carry out a fully twodimensional binned likelihood fit in |p * | and m 2 miss , assuming Gaussian uncertainties. The fit uses 12×12 bins  Table 2 The Asimov data set components. The fractions were motivated by Refs. [10,21].
with equidistant bin widths for |p * | ∈ (0.2, 2.2) GeV and m 2 miss ∈ (−2, 10) GeV 2 . The fits determine either R(D ( * ) ), or the real and imaginary parts of Wilson coefficients. The preferred SM coupling is determined simultaneously, in order to remove explicit dependence on |V cb |.
We construct an Asimov data set [23] assuming the fractions and total number of events in Table 2, following from the number of events in Ref. [10,21]. In the scans, the total number of events corresponds to an approximate integrated luminosity of 5 ab −1 of Belle II collisions. We assume events are reconstructed in two categories targeting B → D τν and B → D * τν. A fit for the real and imaginary parts of a single Wilson coefficient plus the (real) SM coupling thus has 2 × 12 × 12 − 3 = 285 degrees of freedom.
A sizable downfeed background from D * mesons misreconstructed as a D is expected in the B → D τν channel via both the B → D * τν and B → D * ν de-cays. This is taken into account by partitioning the simulated B → D * τ ν and B → D * ν events into two samples: One with the correct m 2 miss = (p B −p D * −p ) 2 and the other with the misreconstructed m 2 miss = (p B −p D − p ) 2 , which omits the slow pion. This downfeed reduces the sensitivity for the case that NP couplings induce opposite effects on the B → Dτν versus B → D * τν total rates or shapes. In addition to semileptonic processes, we assume the presence of an irreducible background from secondaries (i.e., leptons from semileptonic D meson decays), fake leptons (i.e., hadrons that were misidentified as leptons) and semileptonic decays from higher charm resonances (i.e., D * * states). The irreducible background is modeled in a simplified manner by assuming 10 background events in each of the 12×12 bins, totaling overall 1440 events per category. Figure 1 shows the impact on the fit variables of three benchmark models that we use to investigate the effects of new physics: i) The R 2 leptoquark model, which sets S qLlL 8 T qLlL (see, e.g., Refs. [24,25]); ii) A pure tensor model, via T qLlL ; iii) A right-handed vector model, via V qRlL .
For the ratio plots in Fig. 1, we fix the NP Wilson coefficients to specific values to illustrate the shape changes they induce in |p * | and m 2 miss . The R 2 leptoquark model and tensor model exhibit sizable shape changes. The right-handed vector model shows only an overall normalization change for B → D τν, with no change in shape compared to the SM, because the axial-vector B → D hadronic matrix element vanishes by parity and angular momentum conservation. For B → D * , both vector and axial vector matrix elements are nonzero, so that introducing a right-handed vector current leads to shape and normalization changes. Figure 2 shows the projections of the constructed Asimov data set, as well as the distributions expected for the three NP models. The latter have the same couplings as those shown in Fig. 1.

R(D ( * ) ) biases from new physics truth
Many NP analyses and global fits to the R(D ( * ) ) measurements -together with other potentially templatesensitive observables, including q 2 spectra -have been carried out by a range of phenomenological studies (see, e.g., Refs. [24][25][26][27][28][29][30][31][32]). As mentioned above, the standard practice has been to fit NP predictions to the worldaverage values of R(D ( * ) ) (and other data) to determine confidence levels for allowed and excluded NP couplings. However, because the R(D ( * ) ) measurements use SM-based templates, and because the presence of     NP operators can strongly alter acceptances and kinematic distributions, such analyses can lead to incorrect best-fit values or exclusions of NP Wilson coefficients.
To illustrate such a bias, we fit SM MC templates to NP Asimov data sets, that are generated with Hammer for two different NP 'truth' benchmark points: the 2HDM Type II with S qRlL = −2, corresponding to tan β/m H + 0.5 GeV −1 ; and the R 2 leptoquark model with S qLlL = 8 T qLlL = 0.25+0.25 i. (These models and couplings are for illustration and are arbitrary choices; our goal here is only to demonstrate the generic biases that can occur.) We replicate the fit of all existing measurements, allowing the normalizations of the D and D * modes (and the light leptonic final states) to float independently, without imposing e.g. their predicted SM relationship. This fit leads to a best-fit ellipse in the R(D) -R(D * ) plane.
In Fig. 3 we show the recovered values, R(D ( * ) ) rec , obtained from this procedure, and compare them to the actual predictions of the given NP truth benchmark point, R(D ( * ) ) th . For ease of comparison, we normalize the R(D ( * ) ) values against the SM predictions for R(D ( * ) ) SM . For both NP models, the recovered ratios from fitting the Asimov data set exclude the truth R(D ( * ) ) th values at 4σ. This illustrates the sizeable bias in the measured R(D ( * ) ) values that ensues from carrying out fits with an SM template, if instead NP actually contributes to the measurements.
We also show in Fig. 3 the equivalent bias arising from a naïve fit of the R(D ( * ) ) NP prediction that attempts to recover the complex Wilson coefficient. This is done by parametrizing R(D ( * ) ) th = R(D ( * ) )[c XY ], and fitting this expression to the recovered R(D ( * ) ) rec values. Explicitly, one calculates CLs in the Wilson coefficient space via the 2 degree of freedom chi-square The resulting best fit Wilson coefficient regions similarly exclude the truth values.
Thus, the allowed or excluded regions of NP couplings determined from fits to the R(D ( * ) ) measurements must be treated with caution, as these fits do not include effects of the NP distributions in the MC templates. Similarly, results of global fits should be interpreted carefully when assessing the level of compatibility with specific NP scenarios.

New physics Wilson coefficient fits
Instead of considering observables like R(D ( * ) ), for phenomenological studies to be able to properly make interpretations and test NP models, experiments should provide direct constraints on NP Wilson coefficients themselves. For example, this could be done with simplified likelihood ratios that profile out all irrelevant nuisance parameters from, e.g., systematic uncertainties or information from sidebands or control channels, or by other means.  As an example, we now use Hammer to perform such a fit for the real and imaginary parts of the NP Wilson coefficients, using the set of three NP models in Sec. 2.2 as templates. These are fit to the same two truth benchmark scenarios as in Fig. 4: a truth SM Asimov data set; and a truth Asimov data set reweighted to the 2HDM Type II with S qRlL = −2. Figure 4 shows in shades of red the 68%, 95% and 99% confidence levels (CLs) of the three NP model scans of SM Asimov data sets. For the SM truth benchmark, the corresponding best fit points are always at zero NP couplings. The derived CLs then correspond to the expected median exclusion of the fitted NP coupling under the assumption the SM is true.
We further show in shades of yellow the same fit CLs for the 2HDM truth benchmark Asimov data set. These latter fits illustrate a scenario in which NP is present, but is analyzed with an incomplete or incorrect set of NP Wilson coefficients. Depending on the set of coefficients, we see from the ∆χ 2 of the best fit points that the new physics might be obfuscated or wrongly identified. This underlines the importance for LHCb and Belle II to eventually carry out an analysis with the full set of operators listed in Table 1.

The Hammer library
In this section we present core interface features and calculational strategies of the Hammer library. Details Fig. 4 The 68%, 95% and 99% CL allowed regions of three models under consideration, from fitting the SM (red) and 2HDM type II (yellow and with S qRlL = −2) Asimov data sets. of the code structure, implementation, and use, can be found in the Hammer manual [33]; here we provide only an overview.

Reweighting
We consider an MC event sample, comprising a set of events indexed by I, with weights w I and truth-level kinematics {q} I . Reweighting this sample from an 'old' to a 'new' theory requires the truth-level computation of the ratio of the differential rates Historically, the primary focus of the library is reweighting of b → c ν semileptonic processes, often in multistep cascades such as B → D ( * , * * ) (→ DY ) τ (→ Xν)ν. However, the library's computational structure is designed to be generalized beyond these processes, and we therefore frame the following discussion in general terms, before returning to the specific case of semileptonic decays.

New Physics generalizations
The Hammer library is designed for the reweighting of processes via theories of the form where O α are a basis of operators, and c α , are SM or NP Wilson coefficients (defined at a fixed physical scale; mixing of the Wilson coefficients under RG evolution, if relevant, must be accounted for externally to the library). We specify in Table 1 the conventions used for various b → c ν four-Fermi operators and other processes included in the library. The corresponding process amplitudes may be expressed as linear combinations c α A α . They may also be further expressed as a linear sum with respect to a basis of form factors, F i , that encode the physics of hadronic transitions (if any). 3 In general, then, an amplitude may 3 In all b → c processes currently handled by Hammer (see Table 3 for a list) the form factors are functions of or equivalently of the dimensionless kinematic variable, be written in the form in which {s} are a set of external quantum numbers and {q} the set of four-momenta. 4 The object A αi is an NP-and FF-generalized amplitude tensor. In the case of cascades, relevant for B → D ( * , * * ) (→ DY ) τ (→ Xν)ν decays, the amplitude tensor may itself be the product of several subamplitudes, summed over several sets of internal quantum numbers. The corresponding polarized differential rate in which the phase space differential form dPS includes on-shell δ-functions and geometric or combinatoric factors, as appropriate.
The outer product of the amplitude tensor, defined as W ≡ AA † , is a weight tensor. The object (7) is independent of the Wilson coefficients: Once this object is computed for a specific {q} -an event -it can be contracted with any choice of NP to generate an event weight. Similarly, on a patch of phase space Ω -e.g., the acceptance of a detector or a bin of a histogram -the marginal rate can now be written as The Wilson coefficients factor out of the phase space integral, so that the integral itself generates a NPgeneralized tensor. After it is computed once, it can be contracted with any choice of NP Wilson coefficients, c α , thereafter.
The core of Hammer's computational philosophy is based on the observation that this contraction is computationally much more efficient than the initial computation (and integration). Hence efficient reweighting is achieved by -Computing NP (and/or FF, see below) generalized objects, and storing them; with four velocities v = p H b /m H b and v = p Hc /m Hc . For decays with multi-hadron final states, such as the τ → nπ, n ≥ 3, the form factors are also dependent on multiple invariant masses of the final state hadrons. Thus, b → cτ ν decays followed by with subsequent hadronic τ decays involve at least two separate sets of hadronic functions at the amplitude level. 4 The momenta of an event passed to the library must all be specified in the same frame. The choice of frame is arbitrary.
-Contracting them thereafter for any given NP (and/or FF) choice to quickly generate a desired NP (and/or FF) weight.

Form factor generalizations
Similarly to the NP Wilson coefficients, it is often desirable to be able to vary the FF parameterizations themselves. For instance, one might contemplate variations along the error eigenbasis of a fit to the FF parameters, or FF parametrizations that are linearized with respect to a basis of parameters, such as the BGL parametrization [36][37][38] in B → D ( * ) ν. To this end, an FF parametrization with a parameter set {µ} can be linearized around a (best-fit) point, {µ 0 } so that where 'a' is one or more variational indices and e a is the variation. In the language of the error eigenbasis case, F i,a is the perturbation of F i in the ath principal component e a of the parametric fit correlation matrix.
Defining ξ a ≡ (1, e a ) and Φ i,a+1 ≡ (F i , F i,a ), so that Eq. (9) becomes then the differential rate with U an NP-and FF-generalized weight tensor. The ξ a are independent of {q} and factor out of any phase space integral just as the Wilson coefficients do. That is, an integral on any phase space patch, One may thus tensorialize the amplitude with respect to Wilson coefficients and/or FF linearized variations, to be contracted later with with NP or FF variation choices (the latter within the regime of validity of the FF linearization). Hereafter, the ξ a are referred to as 'FF uncertainties' or 'FF eigenvectors' following the nominal fit correlation matrix example.

Rates
In certain cases, it is also useful to compute and fold in an overall ratio of rates Γ old /Γ new , or the rates themselves, Γ new,old , may be required. For example, if the ? MC sample has been initially generated with a fixed overall branching ratio, B new , one might wish to enforce this constraint via an additional multiplicative factor B old /B new . The different components computed by Hammer are then: (i) The NP-and/or FF-generalized tensor for (dΓ new I /dPS)/(dΓ old I /dPS), via Eq. (11), noting that the denominator carries no free NP or FF variational index. (The ratio r I is then itself generally at least a rank-2 tensor.); (ii) The NP-and/or FF-generalized rate tensors Γ old, new , which need be computed only once for an entire sample. (These rates require integration over the phase space, which is achieved by a dedicated multidimensional Gaussian quadrature integrator.)

Primary code functionalities
The calculational core of Hammer computes the NP or FF generalized tensors event-by-event for any process known to the library (see Table 3 for a list), and as specified by initialization choices (more detail is provided in Sec. 3.6) and specified form factor parametrizations. This core is supplemented by a wide array of functionalities to permit manipulation the resulting NP-and FFgeneralized weight tensors as needed. This may include binning -equivalent to integrating on a phase space patch -the weight tensors into a histogram of any desired reconstructed observables, and/or it may include folding of detector simulation smearings, etc. Such histograms have NP-and FF-generalized tensors as bin entries, and we therefore call them generalized or tensor histograms. Once such NP-and FF-generalized tensor objects are computed and stored, contraction with NP or FF eigenvector choices permits the library to efficiently generate actual event weights or histogram bin weights for any theory of interest. The architecture of Hammer is designed around several primary functionalities: 1. Provide an interface to determine which processes are to be reweighed, and which (possibly multiple) schemes for form factor parametrizations are to be used. This includes handling for (sub)processes that were generated as pure phase space. 2. Parse events into cascades of amplitudes known to the library, and compute their corresponding NPand/or FF-generalized amplitude or weight tensor, as well as the respective rate tensors, as needed. 3. Provide an interface to generate histograms (of arbitrary dimension), and bin the event weight tensors -i.e., r I w I , as in Eq. (3) -into these histograms, as instructed. This includes functionality for weightsquared statistical errors, functionality for generation of ROOT histograms, as well as extensive internal architecture for efficient memory usage. 4. Efficiently contract generalized weight tensors or bin entries against specific FF variational or NP choices, to generate an event or bin weight. This includes extensive internal architecture to balance speed versus memory requirements. 5. Provide interface to save and reload amplitude or weight tensors or generalized histograms, to permit quick reprocessing into weights from precomputed or 'initialized' tensor objects.
Examples of the implementation of these functionalities are shown in many examples provided with the source code.

Code flow
A Hammer program may have two different types of structure: An initialization program, so called as it runs on MC as input, and may generate Hammer format files; or a analysis program, which may reprocess histograms or event weights that have already been saved in an initialization run. Pertinent details of the elements of the application programming interface mentioned below are provided in Appendix A, with more details in the Hammer manual. An initialization program has the generic flow: By contrast, an analysis program (from a previously initialized sample, stored in a buffer) has the generic flow: 1. Create a Hammer object and specify the input file. 2. Load or merge the run header -include or forbid specifications, FF schemes, or histograms -with loadRunHeader (after initRun). One may further declare additional histograms to be compiled (from saved event weight data) via addHistogram. (e) Fill histograms with event weights via processEvent.

Conclusions
Precision measurements of b → cτν decays require large Monte Carlo samples, which incorporate detailed simulations of detector responses and physics backgrounds. The limited statistics due to the computational cost of these simulations are often a leading systematic uncertainty in the measurements, and it is prohibitively expensive to generate fully simulated MC samples for arbitrary NP models or descriptions of hadronic matrix elements. In this paper we described the Hammer library, and illustrated its utility. Hammer allows the fast and efficient reweighting of existing SM (or phase-space based) MC samples to arbitrary NP models. In addition, Hammer can be used to change form factor parametrizations and/or incorporate uncertainties from form factors into experimental measurements. Hammer provides a computationally fast way for binned fits to generate predictions, and we implement a demonstrative forwardfolding fit to constrain NP Wilson coefficients using this feature. Such a fit should be carried out by experimental collaborations in future measurements to provide reliable constraints on NP contributions in semileptonic b → cτν decays. The results will allow people outside the collaborations to make correct interpretations of the data, which has not been possible to date without potentially sizeable biases. To demonstrate this latter point, we carried out toy NP analyses using SM fits to NP Asimov data sets, and showed that sizeable biases can indeed occur. Hammer is open source software and we are looking forward to the experimental results and interpretations it will enable. A: Core elements of the Application Programming Interface The user interface of the Hammer library provides four main classes: the Hammer class itself; the Process and Particle classes, used to create events; and the IOBuffer class used for saving and loading precomputed objects. A schematic of the architecture of Hammer is shown in Fig. 5.
In the following we describe various core parts of the Hammer Application Programming Interface (API), with many more details available in the code manual. The library itself is implemented in C++, along with a Python3 wrapper of the API; we will consider here the C++ interface only. This discussion is ordered by scope, rather than the typical code flow. Further details can be found in the Hammer manual [33].

A.1: Building processes and events
A typical decay cascade is contained in the library by the Process class; an event may contain multiple Process instances as e.g., is the case for a signal plus tag B-B pair. Each cascade may be simply represented in graphical terms as a 'process tree', as shown in Fig. 6: Each particle in the cascade is assigned an index, and each decay vertex is represented as a map from a parent index, to the indices of all its daughters. Hammer assembles the process tree through two methods Process::addParticle and Process::addVertex. The former adds a Particle class object -a momentum and a PDG code -to a container of particles; the latter fills the map of each parent index to its daughters for each decay vertex.

A.2: Specifications
The Hammer library contains an interpreter that maps a string representation of a vertex -a vertex string -to all possible charge conserving processes allowed by the specified particle names. The interpreter uses the syntax that particle names are parsed by a capital letter: the full list of names is provided in the manual. (For example the vertex string "D*DPi" is interpreted as all twelve possible D * → Dπ vertices, while "D*+DPi" is interpreted as only the D * + → D + π 0 , D * + → D 0 π + , and (the heavily CKM suppressed) D * + →D 0 π + decay.) The decay processes to be reweighed by Hammer are specified via Hammer::includeDecay, which takes a single vertex string or vector of vertex strings {V 1 , V 2 , . . . , V n } as an argument, and may be invoked multiple times. Each includeDecay specification is inclusive and permits any process tree whose full set of vertices contains all of {V 1 , V 2 , . . . , V n }. The boolean logic applied by includeDecay is AND between each vertex string element, and OR between separate invocations of includeDecay.
The Hammer library allows the user to specify multiple form factor 'schemes' to be used in reweighting. A form factor scheme is a set of FF parameterization choices for each hadronic transition involving form factors, and is labelled by a 'scheme name'. These schemes are set by the method Hammer::addFFScheme, which takes a scheme name plus a map from hadronic string representation to FF parametrization. The hadronic string follows the same syntax and uses the same particle symbols as for vertex strings. For example, ham.addFFScheme("Scheme1", {{"BD", "BLPR"}, {"BD*", "BLPR"}}; ham.addFFScheme("Scheme2", {{"BD", "BGL"}, {"BD*", "CLN"}}); declares two different FF schemes, choosing BLPR for both B → D and B → D * form factors in "Scheme1", and a mixture of schemes for "Scheme2". Separate histograms and event weights are generated for each scheme name. The list of available FF parametrizations are provided in Table 3. The hadronic strings are charge sensitive, hence, e.g., {"B+D", "BLPR"} versus {"B0D", "CLN"} assigns two different FF parametrizations to charged and neutral B → D decays. Specification of the form factor schemes used to generate the MC sample, i.e., the denominator or input form factors, must be specified in order for Hammer to be able to generate the reweighting tensors. These schemes are specified by the method Hammer::setFFInputScheme.
Units of the input MC sample may/should be specified via Hammer::setUnits, for instance ham.setUnits("MeV"). The default is GeV.
The Hammer library permits the user to declare particular vertices, in either the denominator or numerator amplitude, to be evaluated as pure phase space. This is achieved by the method Hammer::addPurePSVertices, which takes a set of string vertices as an argument, and an optional enum WTerm taking values COMMON (default), changes the two BGL outer function parameters from their default settings. (Note that the YAML key for the relevant FF class has a "to" inserted in the hadronic transition, e.g., "BtoDBGL", rather than "BDBGL", to make it clear we are identifying settings for a particular class -the B → D BGL class -and not a process.) By default the library computes the total rate (or looks up a partial width) the first time each unique vertex is encountered in a run. This behavior may be disabled, e.g., if the required integration is multidimensional and time consuming, via ham.setOptions("ProcessCalc: {Rates: false}").
Various other additional specification features are provided by the library, including e.g. specialization of Wilson coefficients to a particular global choice in the event tensor weight calculations. Specifications may also be declared through a card interface, as shown in demo...card.cc example programs provided with the source code.

A.3: Histogramming
Histograms of arbitrary dimensionality may be created by the Hammer library. In general, histogram bins contain event weight tensors (or direct products of them if there multiple processes in the event). 5  The first declaration creates a histogram set each with 20 × 15 bins, no under/overflow, binned uniformly over the respective ranges 3-12 and 0-2.5 (in appropriate units). With reference to the above addFFScheme example, this histogram set contains one histogram for each combination of either "Scheme1" or "Scheme2" with each unique B → D decay cascade. The second declaration similarly creates a set of 3 × 2 histograms with nonuniform bins and additional under/overflow bins. Filling of histograms for a specific event is performed by Hammer::fillEventHistogram, which takes the histogram name and the values of the observables corresponding to each histogram dimension. For example, ham.fillEventHistogram("q2VsEmu", {4., 0.5}) fills the appropriate bin element for the "q2VsEmu" histograms belonging to the event being processed, and fills the relevant histograms for each declared FF scheme name. (Invocations of fillEventHistogram must occur before Hammer::processEvent, discussed in Sec. A.4 below.) If fillEventHistogram is not invoked for a particular histogram for a particular event, the events tensor weight is not added to the histogram. When the under/overflow bool is set to false, events outside the bin ranges are ignored by fillEventHistogram.
Computation of the weight-squared uncertainties is off by default. This may be enabled globally via the options setting ham.setOptions("Histos: {KeepErrors: true}"). However, for computational speed and/or memory efficiency, it may be instead enabled or disabled for individual histograms via Hammer::keepErrorsInHistogram, which takes the name of the histogram as an argument, and a bool. For instance ham.keepErrorsInHistogram("q2VsEmu", true); enables weight-squared computation for this particular histogram. This method should be invoked before initRun.
Various additional histogramming methods are provided by the library, that enable histogram compression, projection, and Wilson coefficient or FF special-ization. These permit reduction of memory requirements or speed enhancements, and are detailed in the manual.

A.4: Processing
An event may contain multiple instances of Process, in order to account for the fact that a single event may feature, e.g., two B decay processes. The Event class is initialized by Hammer::initEvent(), which may take an optional initial event weight (this can also be set by Hammer::setEventBaseWeight). Process instances are added by Hammer::addProcess(proc) which also returns a hash ID of the process. If the process is not allowed according to the includeDecay or forbidDecay specifications, the returned hash ID is zero, and the process is not added to the relevant Event containers.
Once a process is added, it is automatically initialized, which chiefly involves: calculating the signatures of each vertex in the decay cascade; identifying the various subamplitudes making up the cascade, as well as relevant form factor parametrizations and vertex decay rates, for both the numerator/output and denominator/input; and calculating the total rate for the vertex (this is done only the first time each unique vertex is encountered, i.e., only once per run per unique vertex and per FF scheme). The amplitude tensors and form factors are not computed, however, until the invocation of Hammer::processEvent.
Once all processes are added and relevant histograms (if any) have been denoted to be filled, the weights are actually computed and added to the histogram (if any) bins by invocation of processEvent. A pseudo-example on a single event with a set of processes might look like Once processEvent is completed, the event weight may be retrieved by Hammer::getWeight, that takes the FF scheme name. For instance ham.getWeight("Scheme1") computes the currently loaded event weight for the currently specified WCs and FFs. The latter may be set via Hammer::setWilsonCoefficients and setFFEigenvectors.
The method setWilsonCoefficients takes a string that identifies which operator WCs are being set, and either a vector of the WC values or a map. The default WC settings are the SM. A typical example of the usage of this method is ham.setWilsonCoefficients("BtoCTauNu", {{"S_qLlL", 1.}, {"T_qLlL",0.25}}); where the first argument specifies b → cτ ν four-Fermi WCs are being set, and the second argument is a std::map<std::string, std::complex<double>> of each WC to its desired value. The full list of WCs and their definitions is supplied in the manual. An optional third argument is the WTerm enum, that declares whether the evaluation should be applied to the numerator and/or denominator (numerator by default). As an alternative, one may instead pass as second argument a std::vector<std::complex<double>>, corresponding to the ordered basis {"SM", "S_qLlL", "S_qRlL", "V_qLlL", "V_qRlL", "T_qLlL", "S_qLlR", "S_qRlR", "V_qLlR", "V_qRlR", "T_qRlR"}, with the conventions for these WCs shown in Table 1.
It is important to note that the setWilsonCoefficients method, when taking a std::map, produces incremental settings changes. A subsequent invocation ham.setWilsonCoefficients("BtoCTauNu", {{"S_qLlL", 0.5}}) will result in S_qLlL = 0.5 and T_qLlL = 0. 25 would contract the bin weights with the specified NP Wilson coefficients (and FF eigenvectors, if any) for each histogram in the "q2VsEmu" histogram set with FF scheme "Scheme2", and then combines them together into a single final histogram. This contracted histogram output histo is a (row-major) flattened vector of BinContents structs. This struct has members sumWi, sumWi2 and n for sum of weights, sum of squared weights and number of events in the bin, respectively. (By contrast, the method getHistograms extracts all histograms of a specific name and scheme, producing a map of event hash IDs to histogram for all available "q2VsEmu" histograms with the FF scheme "Scheme2".) Integrated rates or partial widths for a specific vertex may be retrieved via Hammer::getRate. The vertex is specified via either a vertex string, or the parent and daughter PDG codes, plus an FF scheme. (Partial widths are returned in the units specified by both return the partial width for the B 0 → D * − µ + ν vertex, using the form factor parameterization specified in "Scheme2", and whatever WCs or FF uncertainties have been specified. (The getRate method is charge conjugate sensitive, so the vertex string must specify sufficient charges to make the vertex charge unique. For example, writing just "B0D*MuNu" would correspond to not only B 0 → D * − µ + ν, but also the very heavily suppressed B 0 → D + µ −ν .) The method getDenominatorRate similarly returns the partial width according to the specified denominator/input FF parametrization chosen in setFFInputScheme, and the denominator/input WCs or FF eigenvectors.

A.7: Multithreading
The library has the ability to perform lock-free parallelization of the getHistogram(s) and getWeight evaluations. This requires use of the thread local methods setWilsonCoefficientsLocal and setFFEigenvectorsLocal to set the desired WC or FF uncertainties. These ...Local methods take the same syntax as setWilsonCoefficients and setFFEigenvectors, but with different behaviour: They do not set the values incrementally from the current settings, but always increment from the SM and zero FF uncertainties, respectively. Global values of the WCs or FF variations are unaffected by the ...Local methods.

A.8: Saving
Hammer provides the ability to store header settings, generated event weights, histograms, and/or rates in binary buffers for later retrieval and reprocessing. These buffers are built on the cross-platform serialization library flatbuffers 6 : The buffer structs Hammer::IOBuffer and Hammer::RootIOBuffer permit writing/reading of Hammer internal objects using C++ binary files and ROOT trees, respectively.

The methods
Hammer::saveRunHeader, saveEventWeights, saveRates, saveHistogram may be used to save: specification settings (like includeDecay etc; the process weight(s) of the event currently loaded in memory (this should be invoked inside an event loop, after processEvent); the computed rates; and histograms. Each of these methods returns a IOBuffer, which can be stored as sequential records in the buffer via an ostream operator. For example, outFile << ham.saveRunHeader(); outFile << ham.saveHistogram("q2VsEmu"); writes the declared run header, with all its settings, into an IOBuffer and passes it as a record into the buffer, and then does the same for the histogram A histogram is always saved sequentially as a definition record then the histogram data record. The saveHistogram method may optionally take additional arguments -such as an FF scheme namein order to save only part of an entire histogram set; see the manual for further details. Saving a buffer in ROOT format is achieved by passing the IOBuffer output of the save... methods into a RootIOBuffer, that may then be stored in a ROOT TTree. Explicit implementations of this functionality are provided in various demo...root.cc example programs.
It is the responsibility of the user to curate the logic and order under which a buffer is saved and then read. Once an object is loaded, it behaves just as the originally computed instance. So one may invoke getHistogram for a reloaded histogram as described in Sec. A.5.
Loading a buffer in ROOT format is achieved by reading the RootIOBuffer stored in a TTree into an IOBuffer that can be passed to the load... methods. Explicit implementations of this functionality are provided in various demo...root.cc example programs.
In order to permit parallelization of initialization runs, the load... methods accept an additional bool, to specify whether to merge the buffer contents with existing objects in memory (true), or overwrite them (false, default). Merging of histograms occurs if two histograms are loaded with a matching name. This merging is additive for histograms in each histogram set with the same FF scheme and hash IDs, and otherwise results in the new unique histograms being appended to the existing histogram set. (If one wishes instead to overwrite a histogram one may instead first invoke removeHistogram, and then reload the desired components of the histogram set.) The methods loadEventWeights and loadRates behave similarly. For weights (rates) with matching hash IDs, merging permits appending of process weights (rates) computed with new form factor schemes to the process weights (rates). Finally, loadRunHeader permits merging of two sets of header specifications into their union. More details are provided in the manual.