Model-independent approach for incorporating interference effects in collider searches for new resonances

The presence of large-mass resonances in the data collected at the Large Hadron Collider would provide direct evidence of physics beyond the Standard Model. A key challenge in current resonance searches at the LHC is the modelling of signal--background interference effects, which can severely distort the shape of the reconstructed invariant mass distribution relative to the case where there is no interference. Such effects are strongly dependent on the beyond the Standard Model theory that must be considered as unknown if one aims to minimise any theoretical bias on the search results. In this paper, we describe a procedure which employs a physically-motivated, model-independent template functional form that can be used to model interference effects, both for the characterisation of positive discoveries, and in the presentation of null results. We illustrate the approach with the example of a scalar resonance decaying into a pair of photons.


Introduction
With the discovery of the Higgs boson at the Large Hadron Collider (LHC) in 2012 [1,2], the experimental evidence for the particle spectrum of the Standard Model (SM) is seemingly complete. In addition to this, there is a fast-growing body of successful comparisons between accurate SM predictions and the corresponding LHC data, for a very large number of observables measured in a significantly diverse set of production processes. Nevertheless, and regardless of this satisfactory phenomenological framework, theoretically the SM is almost universally understood to be the low-energy manifestation of a theory whose validity extends up to the Planck scale. It is hoped that glimpses of such a theory could be found at the TeV scale, which is what motivates the large variety of searches for physics beyond the Standard Model (BSM) at the LHC and at future colliders.
A common feature of BSM physics is the existence of new resonances, whose discovery and characterisation can, for example, be achieved by studying the invariant mass distributions of their decay products. A statistically meaningful quantification of a discovery, or indeed a null result, can be inferred from experimental data after the selection of a benchmark model. In the absence of strong interference effects between the resonant signal and the relevant SM backgrounds, it is straightforward to take a lineshape for the signal (which depends on both the mass and width of the resonance), convolve it with the known detector resolution, and perform a fit to the observed data with both the signal and background invariant mass distributions. In the case of positive results in the data, this allows for the extraction of the discovery significance; for null results, one can extract exclusion limits on the product of the cross section for the resonance production and the branching ratio for decay into the final state of interest. Needless to say, both of these results will depend crucially on the underlying theoretical assumptions; in other words, the same data might lead to different results had different benchmark models been adopted.
Unfortunately, neglecting interference effects is a pragmatic compromise rather than a well-motivated assumption. The details of signal-background interference are highly dependent on the new physics that gives rise to the resonance, which is unknown a priori. A scalar being produced via gluon-fusion through a fermion loop, for example, interferes with existing gluon-induced SM processes that produce the same final state as the scalar, with the precise effects depending on both the mass of the scalar and the masses of any fermions that can run in the loop. Different models will generate different patterns of interference, the ultimate effect of which is to change the production rate of the resonance, whilst distorting the invariant mass distribution of the decay products in such a way as to change the apparent mass of the scalar resonance [3,4,5,6]. The potential presence of these effects complicates both the interpretation of a new discovery in a resonance search, and the presentation of null results in the form of cross-section times branching ratio limits, which are not well-defined in the case of interference.
In this work, we adopt an alternative point of view and present a practical approach for incorporating interference effects in resonance searches in a model-independent way. The key idea is that although the actual lineshape (of the resonance invariant mass) depends on the unknown parameters of an unknown physics model, the space of its possible functional forms is largely dictated by general Quantum Field Theory arguments. Thus, we employ a physically-motivated functional form that is capable of describing the distortions of the lineshape encountered in the presence of signal-background interference, and illustrate how LHC-experiment fitting procedures can be modified to use this functional form in the presentation of both positive and null results. We demonstrate the technique using an assumed model of a scalar resonance produced via gluon fusion and decaying to pairs of photons, but the approach easily generalises to other models and final states. A modelindependent approach to resonance searches has also recently been presented in ref. [7], but we note that the method we propose is vastly different; in particular, ref. [7] relies heavily on a Fourier representation of non-periodic functions, which need not be introduced in this paper.
This paper is structured as follows. Preliminary considerations are first presented in section 2. We then introduce a general, model-independent functional form in section 3. A benchmark signal model is presented in section 4, which will serve as our assumed choice of a scenario that exists in Nature. In section 5, we demonstrate that the general functional form of section 3 is able to describe the physics of the benchmark signal model of section 4, and assume a generalisability of its description to other signal models due to the wide range of behaviour covered. In section 6, we make our tests more realistic by using fully-generated Monte Carlo (MC) samples of the signal, interference and background diphoton invariant mass distributions for the benchmark model, along with simulated detector effects. Finally, we conclude in section 7.

Preliminary considerations
The most straightforward way to present both positive and null search results is that of working in the context of a given BSM theory; an approach of this type is, by construction, a top-down one. While statistically clean, top-down procedures have two main drawbacks. Firstly, they often have to be repeated, even if the datasets are unchanged, whenever a different theoretical model is chosen. Secondly, the details of how BSM theories are treated in such procedures are under the control of the experimental collaborations, which, among other things, renders it difficult for theorists to assess how tweaking different aspects of the models might improve, worsen, or otherwise affect the search results.
For these reasons, it is interesting to consider the opposite viewpoint, namely that of a bottom-up approach in which data are manipulated, and the search results presented using the fewest possible number of theoretical assumptions. This is the goal of the present paper. More specifically, a model-independent functional form for describing the lineshape of a resonance and its interference with the background is employed, and the search results presented as allowed or forbidden regions in the space of parameters relevant to such a form. The idea is that if experimental results are given in this way, any theoretical model can be quickly checked to be compatible or incompatible with the data by means of a simple computation whose results are expressed in terms of the same parameters.
In order to simplify the approach we are proposing, a number of assumptions need to be made. In particular: 1. We consider one resonance at a time; if several resonances are present, they must be sufficiently well separated for the procedure to work independently for each of them.
2. We work with the invariant mass of the resonance, which can be reconstructed by means of the four-momenta of the decay products.
3. A single partonic process is responsible for the signal-background interference pattern.
There is at least one implication of item #3 that requires an immediate explanation. The overarching understanding is that we presently have a solid confidence in the SM, as well as in the correctness of the theoretical tools that are used to simulate both SM and BSM physics processes with a good control on the systematics. Thus, in the entirely realistic possibility that the background to the search proceeds through more than one partonic channel (H → γγ being a chief example of this situation), the channels that are not interfering must be subtracted from the data prior to the fitting procedure that we shall describe below. This operation will contribute to the overall systematics of the procedure we are proposing.

Model-independent template functional form
Given the assumptions listed in section 2, let us denote by m and Γ the mass and width, respectively, of the resonance whose characteristics we seek to determine. We write the amplitude for the partonic process that features the signal-background interference as follows:Ā where we have denoted by q 2 the resonance virtuality 1 , and by h (with 1 ≤ h ≤ N ) the label of the helicity configurations. Loosely speaking, one can identify the first and the second term on the r.h.s. of eq. (3.1) with the "signal" and "background" contributions, respectively. Indeed, the amplitudeB is by construction the one relevant to the production process of interest when all BSM effects are neglected; we remark that we find it convenient to work with B h , rather than directly withB h , owing to the fact that its canonical dimensions are equal to those of S h . We write the complex numbers S h and B h by making their dependences on complex phases explicit, as follows: Thus, the square of the amplitude in eq. (3.1) is: In order to make the forthcoming discussion as transparent as possible, we assume that only one helicity configuration exists, i.e. N = 1; later, we shall consider the case N > 1.
We simplify our notation accordingly, by dropping the index h wherever it appears. With this assumption, the amplitude squared of eq. (3.5), when multiplied by the flux and phase-space factors, is the differential cross section for the signal plus background plus signal-background interference; henceforth, we shall refer to this quantity as to the "full" cross section. By computing its ratio over its analogue stemming from eq. (3.2) (which is thus the background-only cross section), flux and phase-space factors mutually cancel, and we obtain what follows: where the "even" and "odd" dimensionless functions E(q 2 ) and O(q 2 ), respectively, are: with: In the vicinity of the resonance mass, q 2 m 2 , by neglecting all dynamical effects, i.e. by replacing E(q 2 ) with E(m 2 ) and O(q 2 ) with O(m 2 ), eq. (3.7) exhibits the wellknown interference pattern of pure kinematical origin. Namely, the functional form in q 2 is a linear combination of a Breit-Wigner (BW henceforth), which is even under the (q 2 − m 2 ) → (m 2 − q 2 ) transformation, and of a BW times a (q 2 − m 2 ) factor, which is odd under the said transformation. This does not imply that the functions E(q 2 ) and O(q 2 ) are even and odd, respectively. However, we do expect that in a neighbourhood of m 2 the kinematical effects be dominant over the dynamical ones. We can formalise this statement by re-writing eq. (3.7) by Taylor-expanding E(q 2 ) and O(q 2 ) around q 2 = m 2 : having denoted: Equation (3.11) gives us the first opportunity to discuss the bottom-up approach introduced in section 2. One first truncates the Taylor expansion to some order K, i.e. one writes: The terms of O((q 2 − m 2 ) K+1 ) and higher are then discarded, and the parameters that appear in eq. (3.14), namely the following ones: m , Γ , a 0 , . . . a K , (3.15) have to be regarded as parameters to be determined by a fit to the data. The results of such a fit will be compared with the theoretical predictions for the same set of parameters (bar for m and Γ, which must be considered as inputs to theoretical simulations). We point out that, for any given choice of K, the values of the parameters of eq. (3.15) emerging from fitting eq. (3.14) to the data will differ from those one would obtain if one had retained all orders in the Taylor expansion as is done in eq. (3.11), even in the ideal case of infinite statistics. This is because the fit based on eq. (3.14) will tend to compensate for the lack of the missing higher-order terms by suitably adjusting the fit parameters, which can happen rather effectively (i.e. without changing significantly the quality of the fit) if the fitting range in q 2 is chosen in an appropriate manner. It is obvious that, by progressively enlarging such a range, the fit quality will degrade, and eventually lead to an unstable procedure. We shall comment at length on this point in the following, and show that the flexibility in choosing the fitting range is an effective self-diagnostic tool.
The set in eq. (3.15) constitutes a convenient choice since, for any given K, it allows one to include all of the information resulting from the fit in a minimal number of parameters. On the other hand, a possible drawback associated with it is the fact that the parameters a k do not have one-to-one relationships with quantities that emerge directly from matrixelement computations, such as the Taylor coefficients of R(q 2 ) and φ(q 2 ). However, one can express the former parameters in terms of the latter ones. By exploiting eqs. (3.8), (3.9), and (3.12) we obtain, for K = 2 (which will be our default choice henceforth): where, analogously to eq. (3.12), we have defined: (3.19) Thus, after having determined the values of the a k parameters, one solves eqs. (3.16)-(3.18) for the Taylor coefficients of the R(q 2 ) and φ(q 2 ) functions. There are two issues with the procedure. Firstly, the system of eqs. (3.16)-(3.18) is underconstrained: there are more unknowns than equations, the more so the larger K. This implies that the solutions can not be given as central values plus uncertainties for each parameter, but rather as allowed hyperplanes in the space of parameters. For example, the set of possible solutions for R (0) and R (1) will sketch out a band in the R (0) , R (1) plane, with finite width due to uncertainties and the effect of projecting out the remaining parameters. Secondly, the system of eqs. (3.16)-(3.18) is in any case not easy to solve, particularly owing to the presence of trigonometric functions whose argument is φ(q 2 ). This problem can be alleviated by solving directly for the sine and cosine of φ(q 2 ), which is what the notation of eqs. (3.16)-(3.18) already implicitly suggests. While this implies that the system of equations is even more underconstrained, it is mostly an academic issue: in fact, we shall see that it is inevitable in the realistic case of multiple helicity configurations.
An alternative, and much more practical, procedure is that of regarding the parameters on the r.h.s. of eqs. (3.16)-(3.18) directly as fit parameters. This implies employing the r.h.s. of those equations in the fitting template of eq. (3.14), thereby replacing eq. (3.15) with: (3.20) The comparison of eq. (3.15) with eq. (3.20) renders it manifest the first issue discussed above: there are more parameters in the latter set than in the former one. Conversely, since the parameters of eq. (3.20) are more directly related to physical quantities (or rather, to quantities that naturally emerge in theoretical computations), it is possible to reduce their number by means of physics considerations. For example, we expect the complex phases to be more slowly-varying than the absolute values of the amplitudes, and thus we may neglect the q 2 dependence of the former. This implies trimming eq. (3.20) down to: Furthermore, rather than regarding eq. (3.14) as emerging from the Taylor expansion of eq. (3.7), one can start from Taylor-expanding the functions that appear on the r.h.s. of eqs. (3.8) and (3.9), and then replacing the resulting expressions into eq. (3.7), discarding consistently the terms of orders higher than those stemming from the original expansions.
As an explicit example, we consider again our default case K = 2, where eq. (3.21) reads as follows: Conversely, by using a first-order Taylor expansion for R (and by still neglecting the q 2 dependence of the complex phases) the fit parameters are: The corresponding template functional form is the same as in eq. (3.14), with the a 0 , a 1 , and a 2 coefficients given in eqs. We can repeat here the comment made after eq. (3.15). Namely, the simplifying assumptions that lead from eq. (3.20) to eq. (3.22) and thence to eq. (3.23) imply that the values of the parameters that are common to these three sets will in general be different in the three cases. However, at variance with the case of the a k parameters, such differences will not necessarily be small even in fits of comparable good quality, since the system is underconstrained: thus, the individual parameter has more latitude to accommodate for the neglected terms than any of the a k ones. The trigonometric parameters c (0) φ and s (0) φ will give the clearest example of this behaviour.
Furthermore, at variance with the case of eq. (3.15) which is unambiguously determined once K is chosen, the sets in eqs. (3.20)-(3.23) differ from each other owing to considerations stemming from their underlying physics meaning, which is more direct than for eq. (3.15). While this is an appealing characteristic, it must be kept in mind that the parameters in eqs. (3.20)-(3.23) are still not measurable quantities. Thus, the considerations mentioned above must be subject to a level of scrutiny that is deeper than that relevant to the parameters of eq. (3.15); we shall further this point in section 5.
Before closing this section, we return to considering the case of multiple helicity amplitudes, i.e. we work with N > 1 and start from eq. (3.5). The amplitudes squared relevant to the full and background-only cross sections are: It is then a matter of simple algebra to show that eqs. (3.11)-(3.23) are unchanged, provided that the function R(q 2 ) is defined as follows: and that the functions cos φ(q 2 ) and sin φ(q 2 ) are replaced by c φ (q 2 ) and s φ (q 2 ), respectively, where: (3.28) While eqs. (3.27) and (3.28) imply that: and therefore that both c φ (q 2 ) and s φ (q 2 ) can indeed be seen as the cosine and the sine of an angle, in general this is not the same angle. This fact, which has been anticipated before, is what forces one to treat the Taylor coefficients of c φ (q 2 ) and s φ (q 2 ) as independent fit parameters, as we have done in eqs. (3.20)-(3.23). We conclude this section by summarising our fit setup. In what follows we will show results based on two template fits: T R and T a corresponding to the fit parameters listed in eq. (3.22) and eq. (3.15) respectively, the latter with K = 2. For completeness, the two sets of parameters are: which, along with eq. (3.14) (and substitutions similar to those of eqs. (3.16), (3.17), and (3.18) in the case of the T R set), define the two functional forms that we will use. We point out again that the results obtained by employing either of the T R or T a sets constitute alternative descriptions of the same underlying physics; their different characteristics can be exploited depending on the emphasis of the specific new-physics search or modelling.

Benchmark physics model
The discussion of the previous section is general, model-independent, and can be used to interpret the results of any discovery (or null result), even without making assumptions of an underlying resonant physics model. In order to test that this model-independent functional form is indeed appropriate for a realistic physics analysis, we now consider a benchmark model that contains a new Higgs-like CP -even scalar (spin-0) resonance produced via gluon fusion and decaying to two photons. The leading-order (LO) Feynman diagram of the process that we consider X f g γ g γ q g γ g γ Figure 1.
Left: the gg → X → γγ signal process. The resonance X is a CP -even spin-0 particle. The f denotes heavy virtual fermions. Right: the leading order gg → γγ interfering SM background, with circulating quarks q.
is displayed on the left panel of fig. 1. The LO diagram relevant to the SM process that interferes with this signal is shown on the right panel. The two template forms of the previous section will be tested against the invariant mass distribution of the γγ pairs produced in these interactions, which are assumed to represent the model chosen by Nature.
For testing purposes, it is useful to have semi-analytic descriptions of the signal process, the dominant SM background, and the expected resonant-background interference to facilitate the generation of very high-resolution distributions of the diphoton invariant mass. To distinguish these descriptions from the contents of the previous section, we will refer to them collectively as the physics model (PM) functional form. We now discuss the signal, the background, and their interference in turn.

Signal model
In the assumed benchmark model, the interaction of the scalar resonance, X, with gluons is mediated by heavy fermion loops, and can therefore be described by the effective interaction: where G µν is the gluon field strength tensor. Its decay into photons is described by the dimension-5 operator: with A µν the electromagnetic field strength tensor. These effective interactions can be used to compute the amplitude of the production and decay of the scalar resonance. The differential cross section with respect to the diphoton invariant mass is given by: where L gg (q) is the gluon-gluon luminosity function, and A S (q 2 ) is the signal amplitude, which can be written as:  Figure 2. A fit of the signal differential cross section functional form to a histogram of the MC events generated at the m X = 400 GeV, Γ X /m X = 5% point.
with A ggX and A Xγγ the amplitudes of the production loop and decay vertex respectively, and f BW the BW function: where m X and Γ X are the mass and the width of the resonance, X, respectively. Both the effective production and decay vertices contribute a factor to the amplitude with a simple q-dependence, A Xγγ/XGG ∝ q 2 . Thus, using eq. (4.3), we posit that the diphoton invariant mass distribution of our chosen signal can be described by [8]: where f s is a proportionality factor involving all other q-independent factors. An approximation of the gluon luminosity lineshape was extracted using APFEL [9] and the NNPDF2.3 set of parton luminosity functions (PDFs) [10] to leading order: where q is expressed in GeV, and E CM = 13 TeV is the centre-of-mass energy corresponding to LHC Run 2 specifications. A sample of MC signal-only events was generated using MadGraph5_aMC@NLO, interfaced with Pythia8 for a parton shower simulation [11,12,13]. The Higgs characterisation (HC) framework [14] was used for an implementation of our assumed resonant model. Events were generated for a resonant mass m X = 400 GeV, with a width Γ X /m X = 5%. Using the ROOT data analysis framework [15], the analytic form of eqs. (4.6) and (4.7) was tested against a binned histogram of the diphoton invariant mass distribution, for events at the generator level (i.e. those obtained without performing any detector simulation). The result is presented in fig. 2. The values extracted for m X and Γ X agree well with the inputs selected, and a good description of the data is found, with χ 2 /ndf ≈ 1. 5 Figure 3.
A fit of the background template functional form to the histogram of generated background gg → γγ events.
clear that the analytic form of eqs. (4.6) and (4.7) indeed provide a description equivalent to that of the HC model, but which will be quicker to run, and is unaffected by statistical fluctuations.

Background parametrisation
The interfering SM gg → γγ diagram at the LO is shown in fig. 1 (right). While it is possible to proceed in a similar manner as in the signal case in order to obtain a description of the background contribution to the diphoton invariant mass distribution, the physics of the background is not the main interest of our study. We adopt the methodology of experimental collaborations and employ a template functional form for an ad hoc description of the background differential cross section [16]: where E CM = 13 TeV. The parameters of this template consist of a normalisation constant, f b , and exponents, A and B.
The description of eq. (4.8) was tested against a sample of background gg → γγ events generated using MadGraph5_aMC@NLO. The result of its fit to the events is presented in fig. 3; a good description of the simulated background diphoton invariant mass distribution is found, with χ 2 /ndf ≈ 0.9.

Interference between signal and background
Using eqs. (4.6) and (4.8), we can write the interference contribution to the full differential cross section as follows: where c φ X and s φ X are analogous to the quantities of eqs. (3.27) and (3.28), but specifically defined under the assumed benchmark model.  For a heavy resonance that decays via an effective contact interaction, and in the limit of infinite fermion masses for the loop-induced resonant production and background interaction, the phase difference φ h (q 2 ) (eq. (3.6)) vanishes for all interfering helicity amplitudes. Thus, to generate interference-only event samples corresponding to non-trivial phase differences, we have modified the signal amplitude appearing within the HC framework by means of the replacement S → S × e iθ , such that the value chosen for the artificial phase θ then corresponds to the phase difference between the signal and background helicity amplitudes, θ ≡ φ h .
To verify that we can extract the expected phase from interference-only event samples generated in this way, we note that since the same phase difference is defined for each of the interfering helicity configurations, eq. (4.9) can be written in the following form: where f i is defined as: with S h and B h the helicity amplitudes of the signal and background respectively, as defined in eq. (3.1). The four free parameters of eq. (4.10) are m X , Γ X , θ, and f i , where the invariant mass dependence of f i is neglected. Eq. (4.10) was tested against an interference MC sample generated with m X = 400 GeV, Γ X = 20 GeV, and a complex phase θ = −3π/4 ≈ −2.36. The result of a fit to the event sample is shown in fig. 4. The values extracted for m X , Γ X , and θ are in good agreement with those used to generate the events, and a good description of the interference differential cross section lineshape is found, with χ 2 /ndf ≈ 1.2.

Template fits to physics model toys
In this section, we shall test the model-independent functional form of eq. (3.14) against the chosen benchmark physics model. These tests will be performed using toy samples obtained from the analytic description of the benchmark model, rather than from a Monte Carlo generator, in order to achieve very fine resolutions of the invariant mass distributions.
We consider the template functional form as expressed in terms of the T R or T a parameter sets, of eqs. (3.30) and (3.31) respectively, and compare and contrast the results that we obtain for the two different template fits. By design, the parameters of the functional form should be extracted through a fit of the ratio of full to (interfering) background-only differential cross sections, so as to reduce the bias of a result on the particular partonic luminosity assumed for the data. However, in the current study, we shall assume that the background is known exactly (i.e. without any uncertainties associated with its description). In this case, it is equivalent to perform the fits in terms of absolute differential cross sections, using a description that follows from eq. (4.3): where F B (q) is the background-only differential cross section (eq. (4.8)), and the ratio of amplitudes squared is given by eq. (3.14). Binned likelihood fits are performed using the MultiNest implementation of the nested sampling algorithm [17]. We employ a log-likelihood function defined assuming independent Gaussian random variables for each bin: where y i and σ i respectively represent the content and uncertainty of the i th bin, and Y (q i ; Θ) denotes the fit function evaluated at the central bin value q i (in the case of template fits, Y is given by eq. (5.1), with Θ corresponding to either the T a or T R parameter sets).

Construction of Asimov toys
The analytic description of the PM invariant mass distribution, given by the sum of eqs. (4.6), (4.8) and (4.9), is used to construct the Asimov toy histograms that will be used in our tests. Asimov datasets contain no statistical fluctuations [18], and are constructed by setting the content of the k th bin equal to the value of the input distribution, evaluated at the central invariant mass value, q k , of that bin. Bin uncertainties are given by: where A(q) denotes the generating distribution of the Asimov (in this case, the analytic PM description), and N total is the total number of events assumed for the dataset. The square-rooted contribution arises from regular counting statistics, while the rightmost factor accounts for the fact that bin contents do not correspond to a number of events. Thus, the toys perfectly represent the PM invariant mass distributions used to generate them, and provide a clean test-bed for the general template. Ten input PM mass points are considered, with values ranging from m X = 400 GeV to m X = 1300 GeV in increments of 100 GeV. Input widths are set equal to Γ X /m X = 5%,  Figure 5. Visualisation of Asimov datasets, generated using the PM functional form at various mass points with Γ X /m X = 5%, for the two different sets of input c φ X , s φ X , and f s . and two different sets of values for the remaining PM parameters are chosen as follows: with a notable difference being the "height" of the signal, f s , which is two orders of magnitude smaller in set 2. The collection of Asimov histograms constructed for these inputs is shown in fig. 5, for 2 GeV bin widths and N total = 10 million events across the 100-1600 GeV invariant mass range. We note that this number of diphoton events is approximately of the same order as that collected in Run 2 of the LHC. Finally, we point out that, as can be inferred from fig. 5, the contribution of the signal to the full cross section in the high-mass tail is larger than that of the background in the adopted PM. While this might not be the case in actual physics scenarios, we stress that it does not have any implications on the procedure we present in this paper, owing to the limited fit window that is central to the latter; we shall comment explicitly on this matter in section 5.2.

Test of fit windows
The general template we consider retains only terms up to O((q 2 −m 2 ) 2 ) in its form. Thus, it is expected that its description of an invariant mass distribution will deteriorate if one considers masses too far away from the resonance peak. In this section, we fit a range of invariant mass windows centred on the true mass for each toy. Window widths from 2w = 40 GeV to 2w = 400 GeV are chosen, in increments of 40 GeV, such that fit windows are equal to m X ± w. From such a collection of fit results, one can (approximately) determine an appropriate range of fit windows by noting that −2 log L(Θ) follows a χ 2 ν distribution, where for a given fit ν is the number of fitted bins minus one for each free parameter. A 1σ cut-off in the fit quality, for example, can then be determined by evaluating the quantile function of the χ 2 ν distribution, Q χ 2 ν (p), at p ≈ 0.68. If a fit returns a best-fit χ 2 smaller than this value, then we can conclude that the general functional form is able to provide a good description of the data, to within 1σ, over the corresponding mass range.
The result of fitting the T a parameters to PM toys for all of the masses and windows chosen is presented in fig. 6. The left and right panels are obtained using the set 1 and set 2 input PM parameter values of eq. (5.4), respectively. The half-width of the fit window, w, and the input PM mass, m X , are reported on the x and y axes. For each combination of these quantities, the z-axis is represented as colour-coded values according to the scale depicted on the right of the two panels, corresponding to the fit quality in terms of the ratio of best-fit to 1σ cut-off chi-squares, χ 2 bf /χ 2 1σ , with the latter obtained from a χ 2 distribution of appropriate dimensions. Thus, larger z-axis values indicate poorer fits, with a value of one marking the 1σ boundary. Note that the same colour code is used in the two results, but corresponds to different z-axis scales. We do not show the analogous results in terms of the T R parameters, as both parametrisations yield identical lineshapes.
Considering each input mass point separately, we find the trend of decreasing fit quality with increasing fit window size, as is expected. The results of set 1 show a much faster deterioration in contrast to those of set 2, with the latter indicating a good fit for every window tested. Nevertheless, provided that a sufficiently small fit window is chosen, these results show that the general parametrisation is capable of correctly characterising a wide range of physical lineshapes.

Profile likelihood contours in template parameter space
Once a suitable choice for the fit window has been made, aided by results akin to those of fig. 6 or otherwise, the next step of an analysis is to extract a result in terms of the general parameters. As an example, let us refer to the particular instance of fig. 6 corresponding to the set 1 input PM parameters, m X = 700 GeV, Γ X /m X = 5%, and the m X ± 40 GeV fit window. A ratio of chi-squares that is close to zero is found for this configuration, indicating a very good fit of the toy data. Fig. 7 visualises the T a parameter space of this result as a collection of two-dimensional profile likelihood ratios. Each point in a parameter plane corresponds to the ratio of local to global maximum likelihoods: the former is found by profiling over the remaining T a parameters (i.e. by allowing them to adopt values that maximise the likelihood), while the latter corresponds to the overall best-fit likelihood. The white contours represent the 1σ and 2σ confidence boundaries, and a red circle marks the PM input point (m X , Γ X ) in the m, Γ plane. The expected input mass and width values are recovered in the fit to be well within the 1σ confidence level, and we find tightly-constrained contours with very  mild correlations between the T a parameters. Fig. 8 shows the analogous result in terms of the T R parameter set, for a selection of possible two-parameter planes. The m, Γ contour is essentially identical to that in fig. 7, which confirms the consistency of the two m and Γ determinations obtained by means of the T R and T a sets. The R (i) , c φ parameter planes are of more interest: we find large, relatively flat regions of high likelihood, with highly non-trivial degeneracies between the parameters. These are a sharp contrast to the neat solutions of T a space. As we previously mentioned in section 3, such a behaviour is the expected consequence of the system being underconstrained by the T R parameters, and affirms that the results obtained using this parametrisation cannot meaningfully be presented as a set of parameter values with associated uncertainties. Furthermore, flat regions in the fit-parameter space may induce larger uncertainties in cases less ideal than those constituted by Asimov datasets.
From the perspective of a fit, c  The toy corresponds to the m X = 700 GeV, Γ X /m X = 5%, and set 1 input parameters. Twodimensional profile likelihood plots are presented, with 1σ and 2σ boundaries outlined in white. Where applicable, input PM parameter values are indicated with a red circle. Grey contours represent the 1σ and 2σ regions that correspond to physical solutions, c  In these results, the physical contours are found to be fully contained within the enlarged contours. However, this is not necessarily the case in general, since the underconstrained nature of the parametrisation implies that c (0) φ and s (0) φ values can vary to compensate for the lack of higher order terms in the functional form. This can lead to solutions that strongly prefer unphysical regions of the parameter space.
Given that the only reason one would favour the T R over the T a parameters lies in the more immediate physical interpretation of the former, even at the cost of more severe degeneracies, results presented in terms of T R parameters should preferably be consistent with physical constraints. Nonphysicality of the trigonometric parameters can be used as a diagnostic tool: results that do not admit any physical solutions must be discarded. On the other hand, this implies that as long as physical regions are not rejected by the fit, it is safe to present the result in which the physical considerations are imposed on c φ solutions outside of their physical ranges occur more readily with the set of eq. (3.23) (which emerges from a strict first-order expansion of the function R(q 2 )), and for this reason, we suggest the baseline T R fitting procedure to involve the second-order set of eq. (3.30).
We conclude the section by remarking that both of the T R and T a parameter sets can accommodate background-only solutions, and as such, can also be employed for the description of null signals. The T a parameters provide the simplest demonstration of a null result: the solution corresponding to a k = 0 for all k is the background-only one. The analogous solution in T R parameter space is less simple, owing to the complicated correlations between its parameters; however, its advantage lies in that they admit the space of possible signal and interference contributions that sum to an apparent backgroundonly distribution. We present an example of a null result in fig. 9, corresponding to a fit of the T a parameters to a background-only Asimov histogram over the 100-1600 GeV mass range. Note that the result in the m, Γ plane is included for completeness, but is largely meaningless (and can thus be neglected) in the case of a null result.

Physics model closure test
For a theorist to be able to test their specific physics model given a set of results presented in terms of either the T R or T a parameters, there are two possible strategies. Either one computes the T R or T a parameters directly from one's chosen theoretical model, and compares the predictions to those measured by the experiments; or one transforms the T R or T a experimental results into the parameters of one's physics theory. While the former approach is more straightforward in that it follows the same procedure we have outlined so far, we demonstrate in this section that the latter one is viable, too. We adopt again the case corresponding to m X = 700 GeV and Γ X /m X = 5% as an example; by using the results of figs. 7 and 8, we show that the expected (i.e. input) PM parameter values can be recovered.
We begin by using the MultiNest (see eq. (5.2)) sampling outputs corresponding to the T a and T R fits to construct "results" histograms using their respective functional forms. Such histograms are defined over an invariant mass range equal to the fitted invariant mass window (in this case, 660-740 GeV), with the same binning as that of the fitted dataset  (2 GeV). The points in parameter space satisfying were extracted from the MultiNest files, where ν = 5 for the T a parameters (Q χ 2 5 (0.68) ≈ 5.86) and ν = 7 for the T R ones (Q χ 2 7 (0.68) ≈ 8.14). These are the points that lie within the 1σ confidence region of the respective ν-dimensional parameter spaces. For a given histogram bin, the range of functional form values prescribed by these points can thus be interpreted as the 1σ confidence interval of the bin. In this way, we construct histograms, corresponding to a diphoton invariant mass distribution, that represent the band of T a and T R lineshapes agreeing with the original PM toy to within 1σ. Fig. 10 shows the profile likelihood ratios resulting from fits of the analytic PM description to the results histograms constructed as detailed above. The main set of results correspond to the fit of the T a results histogram, with dotted white contours showing the result in the T R case. Since the T a and T R results both describe the same underlying physics, the two sets of contours are very similar, as expected; their slight difference arises due to the greater dimensionality of the T R parameter space, which implies larger bin uncertainties being defined for its results histogram, in accordance with eq. (5.5). For all of the parameter planes visualised, we find the corresponding PM inputs (eq. (5.4), set 1)  Profile likelihood plots of the PM parameters from a fit to T a and T R results histograms, corresponding to input m X = 700 GeV, Γ X /m X = 5% and set 1 parameter values. The fit is performed over the 700 ± 40 GeV mass window. 1σ and 2σ contours are outlined in solid white (T a ) and dotted white (T R ). Input PM points are marked with a red circle.
to lie within the high-likelihood regions obtained. Thus, the procedure described above can indeed be used to construct histograms representative of physical data given a set of results in terms of the general parameters, and thence constrain the parameter space of the physics model of interest.
While we have utilised an analytic description of the physics model in our demonstration, in the case that such a form is not readily available one can generate a Monte Carlo histogram for each point in the parameter space of the model one wishes to test, and compare that to the histograms derived from the T a or T R fit results. This comparison will define a likelihood for the fit of the new model parameters.

Test of discovery significance
The results obtained in the previous sections used knowledge of the true (input) resonance mass in selecting fit windows. However, in a blind search, there is the lower-level question of whether or not a signal is even present. Thus, some preliminary analysis is required to identify candidate signal regions.
For this purpose, a background-only hypothesis is commonly used to quantify the significance of a discovery. The test statistic of interest is given by with L 0 the likelihood of a background-only description against the data, and L(Θ) the likelihood of a model that includes both signal and background, with model parameters Θ (in this case, the T a or T R parameter sets). The q 0 statistic is thus a measure that compares the likelihoods of the null (background-only) and alternative (i.e. including a signal) hypotheses; given a particular dataset and the observed value, q 0,obs , one can calculate the p-value, that subsequent data will exhibit an incompatibility with the background-only hypothesis to an equal or greater degree. Here, f (q 0 |0) represents the probability distribution function for q 0 conditional on the null hypothesis being true, and as such is always non-negative; a larger value of q 0,obs thus yields a smaller p 0 , and indicates a greater disagreement with the background-only hypothesis. An alternative metric that quantifies the statistical incompatibility of the backgroundonly assumption against the observed dataset is the local significance: where Q G (p) is the quantile function for the standard Gaussian distribution evaluated at probability p. For the test statistic of eq. (5.6), the local significance is given by [18]: The common threshold required for a claim of discovery in particle physics is Z 0 = 5.
For the purposes of a background-only hypothesis test, the entire invariant mass region of the search is typically fitted; thus, the T R and T a parameter sets (eqs. (3.30) and (3.31)), which are accurate only within a restricted mass range about the signal peak, are technically inadequate. However, we note that one can still achieve the goal of identifying prospective signals using the second-order functional forms, since the test only seeks to compare the background-only description with one that includes a signal, and not to extract accurate estimations of the parameters of the latter. To prevent an underestimation of the local significance, one should define sufficiently large prior volumes for the parameters during a fit, so as to give them the freedom needed to compensate for the lack of higher-order terms.
We demonstrate the procedure using an Asimov dataset generated from the physics model, with the input parameters as we have adopted in the previous sections (m X = 700 GeV, Γ X /m X = 5%, and set 1 parameter values). The functional form of eq. (5.1) was used as the model alternative to the null hypothesis, and was fitted over the full (100-1600 GeV) invariant mass range of the toy. We note that since parameter estimation is not of import, one can choose to do this using either the T a or T R parameters; in the interest of computational efficiency, we recommend the T a parameter set be employed, as a fit of its parameters typically converges more quickly in comparison to the T R one. The prior volumes for m and Γ were uniformly partitioned into 100 separate regions over the 100-1600 GeV and 0-160 GeV ranges, respectively. A scan over all of the partitions was performed: for each iteration, the maximum likelihood was obtained by profiling over the remaining parameters, and eqs. (5.6) and (5.10) were then used to calculate the local significance, in units of standard deviations (σ).
The results are visualised in fig. 11. A distinct region of high local significance can be identified at the mass value expected from the input value chosen. The width is constrained less precisely by this result, but in a blind search one would only require the former in order to presuppose the existence and approximate mass of a resonance. The remaining general parameters (including the width) can then be more precisely obtained by following our procedure of performing a dedicated fit of the model-independent functional form within a suitable fit window (see in particular sections 5.2 and 5.3).

Procedure for a general resonance search
Our tests of the T R and T a functional forms have demonstrated their flexibility in fitting a wide range of physical distributions, provided one makes a suitable choice for the invariant mass window over which the fit is performed.
A summary of our results, and an outline of the procedure that we propose for a general resonance search, is as follows: 1. Perform a background-only hypothesis test by scanning over the m, Γ plane of the T a parameters, and produce a plot of the local discovery significance. This can also be performed with the T R parameters, but we recommend the T a for simplicity and computational reasons. Apart from the mass and width, all of the other T a (or T R ) parameters are to be profiled over to yield the maximum significance for each iteration of the scan.
2. (a) If notably high significance regions exist: perform fits of the template functional form to the data in these regions over mass windows of increasing size, using either of the T R or T a parameter sets. A suitable mass window that yields a good fit quality should be identified from the results of these fits. Then, extract the general parameters (T R or T a ) from a fit over the chosen mass window. If this is performed using the T R parameters, one should confirm that unphysical solutions are not strongly preferred for c φ . Fit results are to be presented in the form of 2-dimensional profile likelihood contours in the general parameter space.
(b) If notably high significance regions do not exist: perform a fit of the template functional form to the data over the entire mass region. We recommend the T a parameter set for such a fit, due to its simpler interpretation in the case of null signals: one should simply find a k ≈ 0 for each k. Fit results are to be presented in the form of 2-dimensional profile likelihood contours in the general parameter space.
We will also reiterate our preliminary considerations: a single resonance should contribute to the invariant mass range considered, and if multiple regions of high significance are found in step #1, they must contribute to sufficiently separated regions of the data, such that independent analyses can be conducted for each. Furthermore, it is assumed that the data contains only backgrounds produced in the same partonic channel as the signal (and hence generating an interference). Any additional backgrounds yielding the same final state but induced by a different partonic process from the signal should be subtracted from the data prior to the analysis.

Incorporation of detector effects in template forms
The results presented thus far have been obtained under the assumption of statistically perfect datasets, whilst also neglecting event reconstruction effects stemming from the limitations of particle detectors, which will affect the events observed in any physical experiment. To account for these effects, one can choose one of two approaches. The first option is to obtain an approximation of truth-level data using unfolding algorithms, to which the procedure, as summarised in section 5.6, can then be applied. The alternative approach is to incorporate a parametrisation of the detector reconstruction effects within the template descriptions, so that they can be fitted directly to reconstructed data.
In this section, we study the latter approach. We work with MC event samples generated using MadGraph5_aMC@NLO and the HC model, as described in section 4, coupled with a fast simulation of the CMS detector using the Delphes 3.4.1 framework [19]. In accordance with current CMS trigger requirements [20], we designate photons passing a p T > 60 GeV selection cut as candidates for the diphoton pair produced in hard scattering events.
We begin by describing a method of parametrising detector reconstruction effects, before demonstrating that it holds for the reconstructed MC event samples. In applying the procedure to the template functional form, several assumptions are made to simplify the method. It is important to note that these assumptions, while expected to be reasonable, may not be generalisable to all experimental scenarios. As such, the content of this section should not be seen as a comprehensive guideline for an analysis, but rather as a procedure that could be adopted, or adapted, if the assumptions are deemed reasonable. Since our procedure for doing so is based on that used to mimic the full simulation of the ATLAS and CMS experiments, we note that it should therefore be applicable within the LHC experimental collaborations.

Description of detector effects through convolution
The standard procedure of an analysis is to model the detector-smeared invariant mass lineshape, R(q), as the convolution of the truth-level description, T (q), with a detector resolution function, DR(q) [16]: This convolution is assumed to hold for the separate modelling of smeared signal, interference, and background distributions, although for the latter one would typically opt for the alternative of obtaining a direct parametrisation instead. If R(q) is normalised to unity, one can then write the reconstructed differential cross section as: where ε(q) is the (generally q-dependent) reconstruction efficiency, and σ the (truth-level) total integrated cross section. We choose to parametrise DR(q) as a Double-Sided Crystal Ball (DSCB) function, which comprises of a Gaussian core with power law tails: for t > α high , The parameters µ DSCB and σ DSCB are the mean and standard deviation of the Gaussian core. α low and α high determine the mass point at which the Gaussian transitions into power law distributions, with exponents n low and n high respectively. N DSCB is a factor that normalises the DSCB to unity. The parametrisation of the detector resolution function as a DSCB function (whose parameters potentially depend on the invariant mass) is largely inspired by the procedure adopted in ATLAS analyses [21]. The exact description of the DSCB parameters is obtained using narrow-width signal samples generated at a range of mass points; such samples assume sharply-peaked truth-level distributions, T (q) ∼ δ(q − m X ), and it thus follows from eq. (6.1) that the reconstructed narrow-width signal distributions should be well described by the chosen detector resolution function. Given the parametrised DSCB resolution function obtained in this way, a reconstructed lineshape can then be modelled according to eqs. (6.1) and (6.2), for any appropriate truth-level lineshape T (q), and for DR(q) corresponding to the DSCB function with parameters evaluated at the resonance mass, q = m X .
That being said, we point out that a more sophisticated, and fully correct, approach would be that of modelling detector effects by means of a convolution analogous to that of eq. (6.1), but where DR(q) is replaced by a detector response function that depends on several variables (e.g. the photon transverse momenta in addition to the pair invariant mass), all of which would be integrated over. Such a convolution would account for both the resolution and the efficiency contributions to the q-dependence manifest in the observed reconstructed lineshape. Nevertheless, we will adopt the simpler parametrisation of eqs. (6.1) and (6.2) for our current study.

Verification of the convolution description
In this section, we show that eqs. (6.1) and (6.2) can be used to describe reconstructed large-width signal, background-only, and interference-only events. We demonstrate that the description works separately for each of these components, using events that are reconstructed from the corresponding truth-level samples of figs. 2, 3 and 4. In each case, the truth-level description that enters the convolution is taken to be that of eqs. (4.6) (in conjunction with eq. (4.7)), (4.8), and (4.9), respectively, with all parameter values fixed to those reported from their corresponding fits of the truth-level data. A parametrisation of the DSCB function was extracted by using narrow-width signal samples, with its parameter values evaluated at the known resonance mass (m X = 400 GeV), and a numerical value for N DSCB computed and implemented accordingly to ensure normalisation of the DSCB to unity.
Requiring a minimum photon p T during the reconstruction of diphoton events introduces a distortion (dependent on the physics model) in the invariant mass lineshape that is most prominent at the lower end of its spectrum. Thus, to model the lineshape of a reconstructed invariant mass distribution in its entirety, one requires accurate knowledge of the reconstruction efficiency and how it varies with the invariant mass. However, since this effect is expected to manifest mainly in the low mass region, one can also posit that the approximation of a q-independent efficiency should suffice to describe reconstructed lineshapes above some invariant mass threshold. In such a case, only a single free parameter remains in the parametrised description, namely that of an overall normalisation factor corresponding to the reconstruction efficiency of the diphoton events, ε. Note that in general this will differ for signal, interference and background-only events, since the rejection of events below a hard p T threshold removes a different fraction of events from each component, owing to their different photon p T distributions.
We test this assumption by performing one-parameter fits over the 300-1000 GeV invariant mass range. Results are presented in fig. 12. The ratio of the reconstructed (i.e. fitted) function and the corresponding truth-level distribution is visualised beneath each plot; in particular, we note that in the case of the background, this ratio varies slowly with the invariant mass. We find that the simple convolution provides a good description of the reconstructed lineshape for all of the components, and the approximation of q-independent efficiencies has not prevented a good fit being found in each result. There is approximately a 15% discrepancy, at most, between the efficiencies reported for each component, which stems from the p T selection. In a realistic analysis, and particularly as one begins to consider TeV-scale resonances, this difference will diminish and most likely become subdominant in comparison to other sources of uncertainty. For this reason, we expect that it will be valid firstly to adopt a q-independent reconstruction efficiency approximation, and secondly, to neglect the difference between component efficiencies in realistic search Fits to reconstructed signal, interference and background distributions, using a convolution of their respective truth-level analytic functions with a DSCB detector resolution parametrisation. The fits are performed over the 300-1000 GeV invariant mass window. The free parameter in each of the fits is ε, the event reconstruction efficiency. The ratio of the best-fit function and its corresponding truth-level distribution is plotted beneath each result. scenarios.

Convolving the template functional form
We now proceed to apply the convolution method of incorporating detector effects to the general functional form. Note that the template description of eq. (3.14) is formulated as a ratio of full to background-only differential cross sections. However, we will again treat the (reconstructed) background as a quantity that is known exactly, i.e. with zero associated uncertainty. As was remarked in section 5, this implies that it is equivalent to perform fits of the full differential cross sections, and in this case, it allows us to work with the quantity that is consistent with eqs. (6.1) and (6.2).
To further simplify the procedure, we make the following additional assumptions, based on the discussion of the previous section: 1. the reconstruction efficiencies of the various components are approximately equal, 2. the dependence of the efficiency on the invariant mass is negligible; 3. the convolution of the background with the detector resolution function produces negligible change in its shape: where F truth B (q) and F reco B (q) denote the truth and reconstructed backgrounds, respectively.
Along with the assumption of a perfect reconstructed background description, 3 these items constitute the assumptions underpinning the procedure we describe in this section. It is important to point out that these assumptions need only hold within the invariant mass window over which a fit is to be performed.
Under the first assumption, one can exploit the linearity of the convolution operator to convolve the full reconstructed differential cross section, and obtain what follows: where the rightmost equality is made using eq. (5.1), with an according equality between F truth B (q) and the F B (q) of eq. (5.1). The second and third assumptions can then be used to write eq. (6.6) in terms of only the reconstructed background: In this way, the efficiency factor is absorbed into the description of the background, and one obtains a parametrisation that can be expressed without pertaining to any additional free parameters: a fit to reconstructed events can be performed using the same set of parameters as those in the truth scenario, namely those of the T R or T a sets. Note that returning to the ratio regime (cf. eq. (3.14)) is a simple matter of dividing this equation through by the reconstructed background differential cross section. A few comments are in order here. Firstly, eq. (6.7) stems from eq. (6.6) if one can parametrise the background lineshape with the same functional form before and after a convolution with detector effects. If the effect of the convolution is not entirely negligible, the parameters relevant to the two scenarios will generally be different. One might wonder, then, if detector effects are not double counted on the r.h.s. of eq. (6.7). We posit that this is not the case, since these must largely cancel in the ratio F reco B (q)/ε (if that were not the case, the approximate equality between the leftmost and rightmost sides of eq. (6.5) simply could not hold). Secondly, if the convolution with detector effects is so significant for background lineshapes that eq. (6.5) cannot be correct, one needs to instead use eq. (6.6); this is acceptable, but it entails an increased dependence on theoretical predictions, which would likely increase the overall systematics. Finally, in the case of significant departure from all of the assumptions made, it might be preferable to pursue the alternative procedure altogether, namely that of fitting eq. (3.14) directly to (the ratio of) unfolded datasets.
We perform a straightforward test of eq. (6.7) by fitting it to a sample of reconstructed events, and comparing the extracted result to that of the corresponding truth-level fit. An agreement between these two sets of results implies that eq. (6.7) has correctly characterised the underlying physics via a fit of the reconstructed distribution. Following the discussion of the previous section, we conduct this test using MC samples generated at a larger resonance mass of m X = 800 GeV, with width Γ X /m X = 5%, to diminish the impact of approximating equal selection efficiencies between the signal, interference, and background events. The comparison is demonstrated using the T R parameter set.
To determine a suitable mass window for estimating the parameters, preliminary fits were performed using ROOT [15] across a range of invariant mass windows centred on the true resonance mass. The left panel of fig. 13 shows the ratio of best-fit (χ 2 bf ) to 1σ cutoff (χ 2 1σ ) chi-squared values against the fit windows tested. Unlike the results of fig. 6, a monotonically increasing ratio with the fit window is not seen. This is due to the stochastic   Figure 13. Left: the variation of χ 2 bf /χ 2 1σ against the width of the fit window, for fits of the convolved second-order template functional form to reconstructed MadGraph5_aMC@NLO events. Right: A fit of the convolved second-order template functional form to reconstructed MadGraph5_aMC@NLO events over the 300-1200 GeV mass range.
nature of MC samples, in contrast to Asimovs that perfectly capture the underlying theory. 4 To prevent random fluctuations from yielding volatile results, a sufficiently large fit window will thus be required; for the current case, the ratio of chi-squares reaches a minimum at w = 250 GeV, though the distribution effectively plateaus past w ≈ 150 GeV and does not worsen significantly for larger windows. This suggests a resonance small enough that a second-order approximation sufficiently characterises the entire range of its invariant mass lineshape. For samples with larger signals, one would expect to instead see a relationship of the fit quality with a clear minimum in w. The right panel of fig. 13 shows that the T R parameter set is indeed able to fit the reconstructed distribution well over a large 300-1200 GeV mass range. However, for the purpose of parameter estimation, it is still advisable to choose a restricted fit window to avoid the possibility of parameter values drifting to compensate for the missing higher order terms of the functional form of eq. (3.14).
Choosing w = 150 GeV, we perform a fit of the T R parameters over the 650-950 GeV mass window using MultiNest. A fit of the truth-level template was also performed on the truth events. A selection of resulting profile likelihood distributions in T R parameter space are presented and compared in fig. 14. The main set of results, represented by the colour gradient and white contours, correspond to the fit to reconstructed events, while truth-level contours are drawn in grey. A red circle marks the expected mass and width in the corresponding parameter plane.
The reconstructed fit is able to recover the input mass and width values within 1σ confidence. We find the truth and reconstructed-level results to agree well, with the white and grey contours largely overlapping in all of the parameter planes. The largest discrepancies are seen in the R (0) and R (2) parameters; the reconstructed fit yields slightly larger values for these, with an approximately 10% and 30% increase for R (0) and R (2) , respectively, when comparing the truth-level contours to the reconstructed results. This behaviour arises from approximating equal reconstruction efficiencies for the signal, interference and background components. For the current benchmark point, where ε s ε i ε b , the ap-proximation introduces a small enhancement to R(q 2 ), which manifests most noticeably in the zeroth and second order expansion coefficients in the current set of results. However, as we have noted in section 6.2, such effects are expected to become negligible in realistic searches. Our results thus show that the description of eq. (6.7) can be used to extract the general parameters by means of a fit directly to reconstructed events.

Conclusion
The top-down approach typically employed in the presentation of LHC resonance search results, while constituting a straightforward procedure in the scenario where there is strong motivation to believe a priori in a particular model, becomes less ideal for reporting the findings of a general search for BSM resonances. In the latter case, it is desirable to present results in a form that is sufficiently generalisable to any BSM models, such that a given theorist is able to relate them a posteriori to their particular model of interest.
In this paper, we have proposed a procedure that is apt to the presentation of search results in the form advocated above. This procedure assumes that data can be organised as the differential distribution in the invariant mass of the would-be resonance, e.g. as is reconstructed by means of the four-momenta of its decay products, and is based on employing a functional form for the lineshape of such a distribution that stems from general quantum field theory considerations and is fully model-independent. A definite advantage of the procedure is that it allows the characterisation of the resonance and of its interference pattern with the Standard Model background, given in terms of the free parameters that enter the lineshape; the method is shown to work in the case of null search results as well.
The core of the procedure is described in section 3, and its results in terms of the lineshape parameters are presented in section 5.3 using as a test case the simple physics model that is summarised in section 4. While this is the essence of the method we have proposed, we have also discussed various ancillary techniques that supplement it, and which can either be adopted as they are, or replaced by alternative mechanisms if the latter will be deemed more convenient in certain working conditions. More specifically, the lineshapeparameter determination works well if one knows with a good approximation the value of the pole mass of the resonance. In section 5.6 we have shown how a preliminary determination of such a value can be achieved; however, nothing prevents experiments from pursuing a different determination strategy, whose result can then be used in the context of our characterisation procedure. Likewise, when experimental results are presented in terms of the lineshape parameters, any theorist can check the compatibility of their models with the data by computing the same parameters by following exactly the same procedure as is done by experiments. However, in section 5.4 we have also discussed an alternative approach, which employs the lineshape parameters extracted by the experiments to determine allowed regions in the parameter space of the theoretical model one wishes to test. Finally, while the procedure we advocate would typically be applied to unfolded data, we have shown in section 6 how, under certain conditions, it can be employed directly on raw data, by means of a suitable and relatively simple description of detector effects.
We have explicitly shown how the general lineshape functional form can be parametrised by means of two different parameter sets that are characterised by different features-one being more tightly constrained but with a less direct physical interpretation, the other being more closely connected with an underlying theoretical description but liable to have flat directions in the parameter space. Ultimately, the choice of which set to adopt depends  Two dimensional profile likelihood plots of the T R parameters, for a fit of reconstructed diphoton events over the 650-950 GeV invariant mass window. The grey boundaries correspond to 1σ and 2σ contours at the truth level. The physical constraint −1 ≤ c on the emphasis that the analysers will want to give to their searches and/or tests. In view of that, we conclude by pointing out that the two parameters sets we have discussed constitute minimal options that can be systematically extended if necessary, by following the methodology presented in section 3.