A critical point in the distribution of lepton energies from the decay of a spin-1 resonance

TL;DR

This paper analyzes the energy distribution of leptons from a spin-1 resonance decay, revealing a critical point at half the resonance mass that can aid in precise mass measurements, especially for W bosons.

Contribution

It identifies a non-analytic point in the lepton energy distribution at half the resonance mass and shows how finite width effects smooth this feature, providing a new method for resonance mass estimation.

Findings

Half the resonance mass is a special, non-analytic point in the lepton energy distribution.

Finite resonance width smooths singularities, enabling more accurate mass measurements.

The approach reduces dependence on production and decay dynamics for resonance mass determination.

Abstract

We consider a spin-$1$ resonance produced with an arbitrary spectrum of velocities and decaying into a pair of massless leptons, and we study the probability density function of the energy of the leptons in the laboratory frame. A special case is represented by the production of $W$ bosons in proton-proton collisions, for which the energy of the charged lepton from the decaying $W$ can be measured with sufficient accuracy for a high-precision measurement of $M_W$. We find that half of the resonance mass is a special value of the lepton energy, since the probability density function at this point is in general not analytic for a narrow-width resonance. In particular, the higher-order derivatives of the density function are likely to develop singularities, such as cusps or poles. A finite width of the resonance restores the regularity, for example by smearing cusps and poles into local stationary points. The quest for such points offers a handle to estimate the resonance mass with much reduced dependence on the underlying production and decay dynamics of the resonance.

Tables

Table 1. Table 1: The statistical uncertainty (in MeV) and the scale and PDF uncertainty on x ^ 1 subscript ^ 𝑥 1 \hat{x}_{1} and x ^ 3 subscript ^ 𝑥 3 \hat{x}_{3} , separately for W + superscript 𝑊 W^{+} , W − superscript 𝑊 W^{-} , and their combination, without (top rows) or with (bottom rows) acceptance requirements, obtained from a Monte Carlo simulated sample of p ​ p → W ± ​ X → 𝑝 𝑝 superscript 𝑊 plus-or-minus 𝑋 pp\to W^{\pm}X , W ± → μ ± ​ ν μ → superscript 𝑊 plus-or-minus superscript 𝜇 plus-or-minus subscript 𝜈 𝜇 W^{\pm}\to\mu^{\pm}\nu_{\mu} events, corresponding to an integrated luminosity of 1.9 1.9 1.9 fb -1 . The relative uncertainty on σ scale subscript 𝜎 scale \sigma_{\rm scale} and σ PDF subscript 𝜎 PDF \sigma_{\rm PDF} is about 30% and 10%, respectively. The relative statistical uncertainty on σ stat subscript 𝜎 stat \sigma_{\rm stat} is estimated using psuedo-experiments to be about 5% and 25% for x ^ 1 subscript ^ 𝑥 1 \hat{x}_{1} and x ^ 2 subscript ^ 𝑥 2 \hat{x}_{2} , respectively.

	$W^{+}$			$W^{-}$			$W^{\pm }$
	$\frac{{\sigma }_{\mathrm{stat}}}{\mathrm{MeV}}$	$\frac{{\sigma }_{\mathrm{scale}}}{{\sigma }_{\mathrm{stat}}}$	$\frac{{\sigma }_{\mathrm{PDF}}}{{\sigma }_{\mathrm{stat}}}$	$\frac{{\sigma }_{\mathrm{stat}}}{\mathrm{MeV}}$	$\frac{{\sigma }_{\mathrm{scale}}}{{\sigma }_{\mathrm{stat}}}$	$\frac{{\sigma }_{\mathrm{PDF}}}{{\sigma }_{\mathrm{stat}}}$	$\frac{{\sigma }_{\mathrm{stat}}}{\mathrm{MeV}}$	$\frac{{\sigma }_{\mathrm{scale}}}{{\sigma }_{\mathrm{stat}}}$	$\frac{{\sigma }_{\mathrm{PDF}}}{{\sigma }_{\mathrm{stat}}}$
${\hat{x}}_1$	120	22%	39%	110	38%	40%	80	10%	49%
${\hat{x}}_3$	230	5%	6%	370	18%	4%	180	17%	6%
${\hat{x}}_1$	80	14%	44%	90	41%	39%	60	31%	52%
${\hat{x}}_3$	210	6%	10%	350	11%	5%	180	6%	10%

Equations

$$E$$

$$h(y)=\frac{1}{\pi }\frac{\Delta }{(y-1{)}^2+{\Delta }^2},$$

$$z\in \left[\gamma -\sqrt{{\gamma }^2-1},\gamma +\sqrt{{\gamma }^2-1}\right]$$

$$c^{\ast }=\frac{1}{\sqrt{{\gamma }^2-1}}\left(\frac{E}{E^{\ast }}-\gamma \right)\to dE=\sqrt{{\gamma }^2-1}{E}^{\ast }d{c}^{\ast },$$

$$\frac{1}{\sigma }\frac{d\sigma }{d{c}^{\ast }}=\frac{3}{8}\left[\left(1+\frac{A_0(\gamma )}{2}\right)+A_4(\gamma )c^{\ast }+\left(1-\frac{3}{2}{A}_0(\gamma )\right){c^{\ast }}^2\right],$$

$$f_0=\frac{A_0}{2},f_L=\frac{1}{4}\left(2-A_0\pm A_4\right),f_R=\frac{1}{4}\left(2-A_0\mp A_4\right),$$

$$\frac{1}{\sigma }\frac{d\sigma }{dE}(E|\gamma ,E^{\ast })=\frac{1}{\sigma }\frac{d\sigma }{d{c}^{\ast }}\left|\frac{d{c}^{\ast }}{dE}\right|=$$

$$\frac{3}{8E^{\ast }}\frac{1}{\sqrt{{\gamma }^2-1}}\left[\left(1+\frac{A_0}{2}\right)+A_4\left(\frac{E/E^{\ast }-\gamma }{\sqrt{{\gamma }^2-1}}\right)+\left(1-\frac{3}{2}{A}_0\right){\left(\frac{E/E^{\ast }-\gamma }{\sqrt{{\gamma }^2-1}}\right)}^2\right],$$

$$\frac{d\sigma }{dx}(x|\gamma ,y)=$$

$$\times I\left(\gamma -\sqrt{{\gamma }^2-1}\le \frac{x}{y}\le \gamma +\sqrt{{\gamma }^2-1}\right),$$

$$f(x)=$$

$$\times \frac{3}{8}\left[\frac{1+\frac{1}{2}{A}_0(\gamma )}{{\left({\gamma }^2-1\right)}^{\frac{1}{2}}}+\frac{A_4(\gamma )}{\left({\gamma }^2-1\right)}\left(\frac{x}{y}-\gamma \right)+\frac{1-\frac{3}{2}{A}_0(\gamma )}{{\left({\gamma }^2-1\right)}^{\frac{3}{2}}}{\left(\frac{x}{y}-\gamma \right)}^2\right]$$

$$g(1+\kappa )=g^{(0)}+g^{(1)}\kappa +\dots ,$$

$$A_{0,4}(1+\kappa )=A_{0,4}^{(0)}+A_{0,4}^{(1)}\kappa +\dots .$$

$$\frac{1}{2}\left(x+\frac{1}{x}\right)=1+\frac{{ϵ}^2}{2}+\mathcal{O}({ϵ}^4),{\left({\gamma }^2-1\right)}^{-\frac{k}{2}}\approx 2^{-\frac{k}{2}}{\kappa }^{-\frac{k}{2}}$$

$$f(1+ϵ)={\int }_{{ϵ}^2/2}^{\delta }d\kappa \text{mbox}{Pol}_2\left(ϵ;\kappa ,g^{(k)},A_{0,4}^{(k)}\right)+{\int }_{\delta }^{+\infty }d\kappa \text{mbox}{Pol}_2\left(ϵ;\kappa ,g(\kappa ),A_{0,4}(\kappa )\right),$$

$${\int }_{{ϵ}^2/2}^{\delta }d\kappa \left(g^{(0)}+g^{(1)}\kappa \right)\left[1+\frac{1}{2}\left(A_0^{(0)}+A_0^{(1)}\kappa \right)\right]2^{-\frac{1}{2}}{\kappa }^{-\frac{1}{2}}=$$

$$=-g^{(0)}\left(1+\frac{A_0^{(0)}}{2}\right)|ϵ|-\frac{1}{6}\left[\frac{g^{(0)}}{2}{A}_0^{(1)}+g^{(1)}\left(1+\frac{A_0^{(0)}}{2}\right)\right]|ϵ{|}^3+K_{\delta }+\mathcal{O}({ϵ}^5)$$

$${\int }_{{ϵ}^2/2}^{\delta }d\kappa \left(g^{(0)}+g^{(1)}\kappa \right)\left(A_4^{(0)}+A_4^{(1)}\kappa \right)2^{-1}\left(-1+ϵ{\kappa }^{-1}\right)=$$

$$=\frac{1}{2}\left[g^{(0)}{A}_4^{(0)}\ln \delta +\delta \left(g^{(0)}{A}_4^{(1)}+g^{(1)}{A}_4^{(0)}\right)+\frac{1}{4}{g}^{(1)}{A}_4^{(1)}{\delta }^2\right]ϵ+$$

$$-g^{(0)}{A}_4^{(0)}ϵ\ln |ϵ|+\frac{g^{(0)}}{4}{A}_4^{(0)}{ϵ}^2-\frac{1}{4}\left(g^{(0)}{A}_4^{(1)}+g^{(1)}{A}_4^{(0)}\right){ϵ}^3+K_{\delta }^{\prime }+\mathcal{O}({ϵ}^3)$$

$${\int }_{{ϵ}^2/2}^{\delta }d\kappa \left(g^{(0)}+g^{(1)}\kappa \right)\left[1-\frac{3}{2}\left(A_0^{(0)}+A_0^{(1)}\kappa +\dots \right)\right]2^{-\frac{3}{2}}{\kappa }^{-\frac{3}{2}}\left({\kappa }^2-2ϵ\kappa +{ϵ}^2\right)=$$

$$=-\sqrt{2}{g}^{(0)}\left(1-\frac{3}{2}{A}_0^{(0)}\right){\delta }^{\frac{1}{2}}ϵ+g^{(0)}\left(1-\frac{3}{2}{A}_0^{(0)}\right)|ϵ|+g^{(0)}\left(1-\frac{3}{2}{A}_0^{(0)}\right)|ϵ|ϵ$$

$$\left[-\frac{{\delta }^{-\frac{1}{2}}}{\sqrt{2}}{g}^{(0)}\left(1-\frac{3}{2}{A}_0^{(0)}\right)+\frac{1}{\sqrt{2}}\left(-\frac{3}{2}{g}^{(0)}{A}_0^{(1)}+g^{(1)}\left(1-\frac{3}{2}{A}_0^{(0)}\right)\right){\delta }^{\frac{1}{2}}\right]{ϵ}^2+$$

$$-\frac{1}{2}\left[-\frac{3}{2}{g}^{(0)}{A}_0^{(1)}+g^{(1)}\left(1-\frac{3}{2}{A}_0^{(0)}\right)+\frac{1}{6}{g}^{(0)}\left(1-\frac{3}{2}{A}_0^{(0)}\right)\right]|ϵ{|}^3+K_{\delta }^{\prime \prime }+\mathcal{O}({ϵ}^3)$$

$$f(1+ϵ)$$

$$C$$

$$E$$

$$G$$

$$H$$

$$E_0=\sqrt{E_{-}{E}_{+}}.$$

$$ϵ\to {ϵ}^{\prime }=z-1=\frac{x}{y}-1\equiv \frac{1+ϵ}{y}-1,$$

$$f(1+ϵ)\approx \int \frac{dy}{y}h(y)f\left(\frac{1+ϵ}{y}\right).$$

$$\frac{d}{dϵ}{\left[\int \frac{dy}{y}h(y)\left|\frac{1+ϵ}{y}-1\right|\right]}_{ϵ=0}=\int \frac{dy}{y^2}h(y)\text{mbox}sign(1-y)=\mathcal{O}({\Delta }^2),$$

$$\underset{\gamma \to 1^{+}}{\lim }g(\gamma )$$

$$=\underset{\gamma \to 1^{+}}{\lim }{\int }_0^{|\mathbf{q}{|}^2}d{q}_{\mathrm{T}}^2\frac{d^2\sigma }{d{q}_{\mathrm{T}}^{2}d{q}_z^2}=\underset{\gamma \to 1^{+}}{\lim }{\int }_0^{|\mathbf{q}{|}^2}d{q}_{\mathrm{T}}^2\frac{1}{2E\sqrt{|\mathbf{q}{|}^2-q_{\mathrm{T}}^2}}\frac{d^2\sigma }{d{q}_{\mathrm{T}}^{2}dy}$$

$$\approx {\left[\frac{d^2\sigma }{d{q}_{\mathrm{T}}^{2}dy}\right]}_0\frac{|\mathbf{q}|}{M}{\int }_0^{1}d\zeta \frac{1}{\sqrt{1-\zeta }}\propto \sqrt{{\gamma }^2-1},$$

$$\frac{d\sigma }{d{q}_{\mathrm{T}}^2}(q_{\mathrm{T}}|y)={\sigma }_0+{\sigma }_1{q}_{\mathrm{T}}^2+\dots$$

Full text

1 Introduction

The problem of estimating the mass $M$ of a resonance that partially decays into undetectable particles often arises in collider experiments. For example, it occurs when some of the decay products of the resonance interact too weakly with the detector to produce a signal, or when they are measured with insufficient precision. If the kinematics of the collision event can be closed by other means, for example by using energy-momentum conservation, the problem has an obvious solution, otherwise it is under-constrained.

The loss of information due to the unobserved particles, which prevents $M$ from being unambiguously determined on an event-by-event basis, can be statistically recovered if the dymanics of both production and decay of the resonance are known. When such a prior knowledge is available, the probability density function (p.d.f.) of the visible particle momenta $\{\bm{p}_{\ell}\}$, denoted by $\sigma^{-1}d\sigma/d\{\bm{p}_{\ell}\}$, can be computed by marginalizing the unobserved degrees of freedom. This marginal p.d.f. depends on the unknown resonance mass through the kinematics of the visible decay products. In general, the multi-dimensionality of the observable space makes the analytic calculation of this function of paramount complexity. The problem is then best tackled by the use of Monte Carlo (MC) simulations of the process of interest, resulting in a discrete set of MC templates $\sigma^{-1}_{\rm MC}\,d\sigma^{\rm MC}(M_{i})/d\{\bm{p}_{\ell}\}$ generated at different trial values of $M$. With these templates at hand, a numerical evaluation of the likelihood function of the data is possible, and the standard theory of likelihood-based estimators can be then used for estimating the unknown mass James . By construction, such an approach is model-dependent, as it relies on theoretical assumptions (in fact, the complete $S$-matrix for the process of interest) for relating $\sigma^{-1}d\sigma/d\{\bm{p}_{\ell}\}$ to $M$. There are indeed examples where model uncertainties represent the limiting factor to the experimental accuracy. The determination of the $W$ boson mass at hadron colliders represents perhaps the most remarkable of such cases WMassCDF ; WMassD0 ; WMassATLAS .

An alternative approach, which allows the aforementioned limitation to be partly overcome, consists in exploiting singularities in the phase-space of the visible observables Kim , i.e. special points where the tangent plane to the phase-space manifold is aligned with one of the invisible particle directions. The position of such pointed features in the spectra of kinematic variables can be related to the unknown mass, or, more generally, to combination of masses when there is more than one resonance in the decay chain Han . Besides being ideally independent from the details of the underlying dynamics, the main advantage of the phase-space singularity method is that a multi-dimensional problem is recast into a search for striking features, like sharp edges or cusps, on univariate distributions. A study of the phase-space singularity method in the context of the $W$ boson mass measurement at hadron colliders has been documented in Ref. Rujula . Not surprisingly, the optimal of such singularity variables is highly correlated with the usual transverse mass which, being a function of the transverse hadronic recoil, is affected by other well-known sources of experimental uncertainty WMassATLAS .

Motivated by the need of reducing the model-dependence in the measurement of the $W$ boson mass without having to rely on the hadronic recoil, we will concentrate hereafter on the special case of a spin-1 resonance that decays into a pair of massless leptons, of which one is assumed to be measured with high precision, whereas the other is undetected. It has been pointed out in Ref. Franceschini that a two-body decay kinematics of this type features an obvious, yet subtle, invariance under boosts. Indeed, the mass of the mother particle plays a special role in the distribution of energy $E$ of the visible daughter particle. In particular, it can be proved that $M/2$ is a local maximum of the energy distribution $\sigma^{-1}d\sigma/dE$, if the mother particle is produced unpolarized. In this case, one can just measure $M$ by locating the point in the observed energy spectrum featuring the largest density. An application of this technique in the context of the top-quark mass measurement has been documented in Ref. CMS-top .

The argument leading to the identity $\mbox{argmax}[d\sigma/dE]=M/2$ relies on the assumption that the resonance is unpolarized. Instead, we would like to be as agnostic as possible with respect to the mechanism of production and decay of the resonance. In this spirit, we will study the mathematical properties of the p.d.f. of the lepton energy in full generality by deriving exact results in the approximation of a narrow-width resonance. Strictly speaking, any unstable resonance has a finite width $\Gamma>0$. In practice, the latter has to be compared with the experimental resolution $\sigma_{E}$ on the visible particle energy, which sets the minimum granularity at which differential properties of the p.d.f. $\sigma^{-1}d\sigma/dE$ can be defined. The case $\Gamma/\sigma_{E}\ll 1$, is then mathematically equivalent to treating the resonance in the narrow-width approximation. We will then validate the results against selected toy examples of production and decay dynamics. The results of this study motivate the usage of stationary points in the higher-order derivatives of the energy p.d.f., in particular of the second derivative, as an estimator of the resonance mass. Finally, we will study this method in the context of the $W$ boson mass measurement at the LHC using a MC simulation of the reaction $pp\to W^{\pm}X$, $W^{\pm}\to\mu^{\pm}\nu_{\mu}$ in proton-proton ($pp$) collisions at a center-of-mass energy $\sqrt{s}=13$ TeV. Quantitative estimates of the statistical and of the dominant theoretical uncertainty affecting the newly proposed method of measurement are also provided.

2 Kinematics in the laboratory frame

Let $E$ ($E^{*}$) be the lepton energy in the laboratory (center-of-mass) frame, and $c^{*}\equiv\cos\theta^{*}$ the cosine of the polar angle in the center-of-mass frame with respect to the mother particle velocity $\bm{\beta}$ in the laboratory. We also define $E_{0}=M/2$ and introduce the dimensionless parameters $x=E/E_{0}$, $y=E^{*}/E_{0}$, and $z=E/E^{*}=x/y$. The lepton energy in the laboratory is related to $E^{*}$ and $c^{*}$ via a Lorentz transformation that depends only on the boost factor $\gamma=(1-\beta^{2})^{-\frac{1}{2}}$, with $\beta=|\bm{\beta}|$, namely:

[TABLE]

The distribution of energies in the center-of-mass frame is assumed to be described by a Breit-Wigner function:

[TABLE]

where $\Delta=\Gamma/2M$. Since we are mostly concerned with narrow-width resonances, i.e. $\Delta\ll 1,$ we can safely neglect the fact that the function $h$ in Eq. (2) should be truncated at $y=0$ to prevent the center-of-mass energy from becoming negative. In fact, Eq. (2) coincides with the more correct relativistic Breit-Wigner distribution PDG only when $y\approx 1$ (although it is somehow simpler for the calculations to use the non-relativistic version of Eq. (2), the results presented here do not depend on this assumption). Finally, we remark that this p.d.f. converges to the Dirac delta function in the limit $\Delta\to 0$.

From Eq. (1), the domain of $z$ is found to be:

[TABLE]

where the relation $\gamma^{2}\beta^{2}=\gamma^{2}-1$ has been used. If $\gamma\neq 1$, Eq. (1) can be inverted yielding:

[TABLE]

which implies a linear relation between $c^{*}$ and $E$ at a fixed value of $\gamma$ and $E^{*}$.

In the center-of-mass frame of a spin-1 resonance decaying to a pair of spin-$1/2$ particles, the cosine of the polar angle of the decaying lepton with respect to a given quantization axis is described by a p.d.f. of the form Mirkes :

[TABLE]

where the angular coefficients $A_{0,4}$ have been introduced as arbitrary dimensionless functions of the boost factor $\gamma$. The $A_{0}$ coefficient controls the fraction of longitudinal polarization ($f_{0}$) and satisfies the requirement $0\leq A_{0}\leq 2$, whereas $A_{4}$ regulates the fractions of left ($f_{L}$) and right ($f_{R}$) transverse polarization. For a pure $V-A$ interaction, the angular coefficients are related to the polarization fractions $f_{0,L,R}$, relative to the direction of flight of the resonance, by the linear relations:

[TABLE]

where the choice of sign depends on the lepton charge. Special cases of Eq. (6) are the values $(A_{0},A_{4})=(0,\pm 2)$, which corresponds to a purely left/right polarized resonance, and $(2/3,0)$, which corresponds to an unpolarized resonance. By combining Eqs. (4)-(5), we obtain the conditional p.d.f. of $E$:

[TABLE]

where the explicit dependence of the angular coefficients on $\gamma$ has been omitted for simplicity. Multiplying both sides of Eq. (7) by the constant $E_{0}$, we obtain:

[TABLE]

where $I(\cdot)=1$ if the argument is true, [math] otherwise. The p.d.f. of $x$ can be now obtained by marginalizing both $\gamma$ and $y$. We assume $\gamma\sim g(\gamma)$ independently of $y$, which is usually appropriate for a narrow-width resonance. Under this assumption, we can write:

[TABLE]

The p.d.f. $g$ is positive-definite and normalized to unity: $\int_{1}^{\infty}d\gamma\,g(\gamma)=1$. We first consider the case that $g$ is an analytic function everywhere, in particular at $\gamma=1$ (the alternative case will be discussed later). Under this assumption, it can be replaced by its Taylor series centered at $\gamma=1$:

[TABLE]

where $\kappa\equiv\gamma-1\geq 0$. Likewise, we assume that $A_{0,4}(\gamma)$ are analytic at $\gamma=1$ such that:

[TABLE]

We now move to study the behavior of $f$ when $x\approx 1$. To this purpose, we expand the right-hand side of Eq. (9) in terms of a small parameter $\epsilon$, such that $x=1+\epsilon$. In this limit, we have:

[TABLE]

where $k$ is an integer.

2.1 The narrow-width approximation

We first consider the case of a narrow-width resonance, i.e. we set ${h}(y)=\delta(1-y)$. After integrating-out $y$, the right-hand side of Eq. (9) can be rewritten symbolically as:

[TABLE]

where $\mbox{Pol}_{2}(\epsilon\,;\,\cdot)$ stands for a second-order polynomial in $\epsilon$. In Eq. (13), the integration region has been split into two disjoint intervals: $[\epsilon^{2}/2,\delta]$, where the cut-off $\delta$ is sufficiently small that the approximations in Eq. (10)-(11) are valid to first order, and the complementary interval $[\delta,+\infty]$. The first integral provides the contribution inside a neighborhood of $x=1$ from the phase-space region $\gamma\approx 1$, i.e. when the decaying particle is almost at rest; the second integral accounts for the contribution stemming from larger boost values. By virtue of the spin-1 assumption, the integrand function within both integrals is a quadratic polynomial in $\epsilon$, hence it has vanishing derivatives beyond the second order. We can now compute explicitly the first integral at the right-hand side of Eq. (13). After a straightforward integration, we get:

[TABLE]

where $K_{\delta}$ are constants that depend only on the cut-off $\delta$. By rearranging the various terms in Eq. (14), we finally obtain:

[TABLE]

where the constants $A,\ldots,H$ are independent of $\epsilon$. There are terms in this expansion which are not analytic at $\epsilon=0$. They are are proportional to the constants:

[TABLE]

As a consequence of Eq. (15), the higher-order derivates of $f$ can develop various types of singularity at $x=1$: kinks or cusps (from $|\epsilon|$ terms), discontinuities (from $\mbox{sign}(\epsilon)$ terms), delta functions (from the derivatives of the latter), and poles (from the derivative of the $\epsilon\ln|\epsilon|$ term). This non-regular behavior should not come as a surprise: even if $g$, $A_{0}$, and $A_{4}$ are smooth functions, the transformation in Eq. (4) becomes singular in the limit $\gamma\to 1^{+}$. When convoluted with a continuous spectrum of boosts, this primordial singularity is weighted by an infinitesimal cross section $g(1)dx$, but still percolates to the final p.d.f., in a way that depends on how the phase-space $(\gamma,\,c^{*})$ gets populated. We anticipate that the appearance of a singularity in a strict mathematical sense is a consequence of treating the resonance in the narrow-width approximation. Within this approximation, however, its existence is a robust result, as discussed later.

The nature of such singularity is, to some extent, akin to the phase-space singularity discussed in Ref. Kim . Indeed, for a fixed value $\gamma>1$, the variable $x$ has two wall singularities associated with the decay of the visible particle collinear or anti-collinear with the velocity of the resonance. These configurations correspond to edges of the phase-space. When $\gamma=1$, a singularity of higher degree appears because the dimensionality of the phase-space shrinks from a line to a point. The singularity studied here has, however, a richer phenomenology compared to the algebraic singularity of Ref. Kim , because it depends not just on the geometry of the phase-space manifold, but also on how the dynamics of production and decay distributes events across the phase-space.

It is interesting to consider some limiting cases of Eq. (16). As expected, an unpolarized resonance gives rise to a p.d.f. of the form $f(x)=A+C|x-1|+\mathcal{O}((x-1)^{2})$, implying that $x=1$ is a local maximum of the density (in particular, it is a cusp if $g^{(0)}>0$ and a stationary point otherwise). This is in agreement with the result obtained in Ref. Franceschini . Whenever the boost p.d.f. has a vanishing amplitude in the neighborhood of $\gamma=1$, i.e. $g^{(k)}=0$ for the first $k$ derivatives, the coefficients in Eq. (16) are also vanishing. A special case is when there is a minimum momentum threshold on the production of the resonance, such that $g(\gamma)\equiv 0$ for $\gamma\leq\gamma_{\rm thr.}$ (in this case, all the derivatives at $x=1$ are formally zero). The expansion of Eq. (15) thus contains only terms of order $\epsilon^{k}$, with $k=0,1,2$: in the neighborhood of $x=1$, the energy p.d.f. is proportional to a parabola. Equation (3) can be then used to express $E_{0}$ in terms of the lower ($E_{-}$) and upper ($E_{+}$) edges of the interval in which $f(x)\propto\mbox{Pol}_{2}(x)$, as:

[TABLE]

We now briefly consider the possibility that either of the functions in the integrand is not analytic at $\gamma=1$, such that the Taylor expansion of Eq. (10)-(11) are not defined. To fix the ideas, we consider the case $g(\gamma)\sim(\gamma-1)^{\alpha}$, with $\alpha>0$, for which $\gamma=1$ is a cusp point. As we will see later, this choice of p.d.f. finds at least one remarkable physical application. In this case, Eq. (14) gets modified by the appearance of terms like $|\epsilon|^{2\alpha+m}$, where $m=1,2,3,\ldots$ is an integer, which, for arbitrary values of $\alpha$, gives rise to the same phenomenology of non-regularity on $f$.

2.2 Finite-width effects

The regularity of laboratory energy p.d.f. is restored by integrating over a continuous spectrum of center-of-mass energies. Indeed, in the limit $\gamma\to 1^{+}$, the laboratory frame coincides with the center-of-mass frame, i.e. $f(x)=h(x)$, which is a smooth function. In the case of a finite width, Eq. (13) applies with the replacement:

[TABLE]

where $\epsilon$ is again defined as $\epsilon=x-1$. Equation (13) thus becomes:

[TABLE]

Consider for example a term like $|\epsilon^{\prime}|$ in the expansion of $f$. Upon integration over $y$, its first derivative calculated at $x=1$ gives:

[TABLE]

hence the new minimum/maximum of $f$ gets displaced from $x=1$ by an amount of $\mathcal{O}(\Delta^{2})$. Notice that if the $k$-th order derivative has a kink such that $|\lim_{x\to 1^{+}}f^{(k)}|\neq|\lim_{x\to 1^{-}}f^{(k)}|$, the integration over $y$ smears this singularity into a stationary point whose position depends not just on $\Delta$, but also on $g$ and $A_{0,4}$, which determine the left and right slopes of $f^{(k)}$. In the latter case, nothing can be said a priori about $E_{0}$, unless that it must be close to the stationary point. For example, in the case of a symmetric kink in the second-order derivative, such a displacement is of $\mathcal{O}(\Gamma^{2}/M)$ and could be in principle subtracted from the measured stationary point, if $\Gamma$ were also known. In general, by knowing the width of the resonance and by relying on some theory prior on the coefficients of Eq. (15), the measured stationary point can be calibrated to recover an unbiased estimator of $E_{0}$.

2.3 An explicit example

We now discuss the case of $W$ boson production at the LHC which allows us to specialize some of the generic formulas derived before.

By using the fact that $d\gamma^{2}=d(|\bm{q}|^{2}/M^{2})$, we can write

[TABLE]

with $\zeta=q_{\rm T}/|\bm{q}|$. Hence, $g(1)=0$. Furthermore, since $g\approx\sqrt{\gamma^{2}-1}$, the boost spectrum is not analytic at $\gamma=1$. The finiteness of the double-differential cross section $\left[d^{2}\sigma/dq^{2}_{\rm T}dy\right]_{0}$, where $y$ is the rapidity of the $W$ boson, is a general result that arises from the small transverse momentum behavior of cross sections in hard processes Petronzio :

[TABLE]

Similarly, the differential cross section $d\sigma/dy$ at $y=0$ is finite because it is proportional to the product of the partonic densities evaluated at $x_{1,2}=M/\sqrt{s}$, where $\sqrt{s}$ is the proton-proton center-of-mass energy. By using a complete Monte Carlo simulation of this reaction, discussed in Sec. 4, we also find that the $\sim(\gamma-1)^{\frac{1}{2}}$ rise is limited to a tiny region of phase-space, typically $\gamma-1\lesssim 5\times 10^{-4}$, corresponding to the region $|\bm{q}|\lesssim 4$ GeV, i.e. where the differential cross section in $q_{\rm T}$ is rapidly growing. For $\gamma$ values in excess of about $10^{-3}$, the boost p.d.f. is well approximated by a power law of the form $g(\gamma)\sim(\gamma-p_{0})^{-p_{1}}$, with $p_{0}\approx 0.9$ and $p_{1}\approx 0.8$.

The angular coefficients are dimensionless functions of $|\bm{q}|/M$ and $y$ encoding the average polarization of $W$ bosons produced in hadron collisions as a function of the $W$ boson kinematics Mirkes . They are simultaneously determined by the partonic density functions (PDFs) of the proton and by emission of additional QCD radiation. Furthermore, they depend on the reference frame, i.e. they are not rotation-invariant. The coefficients $A_{0,4}$ of Eq. (5) are calculated in a particular rest frame of the $W$ boson defined by a boost along the velocity $\bm{\beta}$, which we will refer to as the helicity frame. For fixed and small values of $\gamma$, their values are determined by the average of $A_{0,4}$ over all momenta $\bm{q}$ defining the surface of a sphere of radius $|\bm{q}|$, such that $\gamma\approx 1+\frac{|\bm{q}|^{2}}{2M^{2}}$:

[TABLE]

Similarly to Eq. (2.3), one can easily verify that all directions $\bm{q}/|\bm{q}|$ are equally likely in the limit $|\bm{q}|\to 0$:

[TABLE]

where the right-hand side does not depend on the direction $\Omega$. When $|\bm{q}|$ is small, the direction of the quark and antiquark in any rest frame of the $W$ boson remain almost anti-parallel. They both carry spin parallel or anti-parallel to their respective momentum. While averaging over the full solid angle, their directions, which are fixed in the laboratory frame, move isotropically in the helicity frame around the quantization axis defined by $\bm{q}$ itself. The net result must be a uniform distribution in $c^{*}$, i.e. $A^{(0)}_{0}=2/3$ and $A^{(0)}_{4}=0$, or, using Eq. (6), $f_{L}=f_{R}=f_{0}$. This is only true when $|\bm{q}|\to 0$. For small finite momenta, a tiny polarization is produced. The amount by which a polarization is built by the misalignment of the quark directions can be estimated looking at a particular configuration where $y=0$ and $|\bm{q}|=q_{\rm T}$. In this case, the center-of-mass frame is related to the Collins-Soper (CS) frame CS by a rotation of $\pi/2$. From a MC simulation, we find that the longitudinal polarization in the CS frame is built at a pace of about $A_{0}(q_{\rm T})\approx 2\times(q_{\rm T}/M_{W})^{2}$, or equivalently $A_{0}(\gamma)\approx 4\gamma$. The values of $(f_{L},f_{R},f_{0})$ in the helicity frame is thus perturbed by an amount of similar size from the values $(\frac{1}{4},\frac{1}{4},\frac{1}{2})$ in the limit $|\bm{q}|\to 0$. When $|\bm{q}|$ grows, large values of $q_{\rm T}$ in Eq. (23) become increasingly unlikely and the cross section favors a longitudinal motion with $|q_{z}|\gg q_{\rm T}$. In this latter case, a net transverse polarization is built as a consequence of the PDF ratio $q/\bar{q}$ growing fast at large rapidities. Again using a MC simulation, we find an empirical slope $dA_{4}/dy\approx\pm 0.3$ for $q_{\rm T}=0$, where the sign depends on the charge of the $W$ boson. Given that $dy=d\ln\gamma$, we also have $A_{4}^{(1)}\approx 0.3$.

We notice that the same argument applies to the case of a proton-antiproton collider, for which $W$ bosons produced almost at rest in the laboratory frame are preferentially polarized in the direction of the antiproton. When integrating over the full solid angle, however, the average polarization in the helicity frame vanishes.

2.4 Discussion

We now summarize the results obtained so far. When considering the two-body decay of a spin-1 resonance, the probability density function $f$ describing the laboratory energy of any of the two daughter particles, assumed to be massless spin-$\frac{1}{2}$ particles, is in general not analytic at $x=1$, when the narrow-width approximation for the resonance is made. In particular, the derivatives of $f$ are likely to be non-derivable, discontinuous, or divergent at that point, depending on the boost factor p.d.f. and on the polarization of the resonance. The appearance of a local maximum of the density at $x=1$ is in general a fortuitous occurrence. A pole at $x=1$ in the first derivative is associated with the presence of a non-zero transverse polarization at rest. Conversely, cusps in the second derivative appear quite naturally as a result of terms of the form $|x-1|^{3}$ in the expansion of $f$ around $x=1$, which do not compete with $(x-1)^{3}$ terms from higher boost factors. This is a general result that only depends on the spin-1 assumption for the resonance. The condition for which a cusp is generated amounts to $|G|>|F|$, as defined in Eq. (15). If $g^{(0)}=0$, this is satisfied if $A_{4}^{(0)}=0$ (otherwise the coefficient $E$ in Eq. (16) would be non-zero, giving rise to a pole in $f^{(1)}$). If instead $g^{(0)}>0$ and $A_{4}^{(0)}=0$, the condition for developing a cusp in $f^{(2)}$ is:

[TABLE]

which is satisfied if $A_{4}$ is a slowly varying function of $\gamma$. We remark that even if Eq. (25) were not satisfied, kinks at $x=1$ for at least one among $f^{(0)}$, $f^{(1)}$, and $f^{(2)}$ will be present, so that a divergence in the higher-order derivatives will eventually show up.

The following search algorithm is then proposed. For simplicity, we assume the laboratory energy $E$ to be normalized to a constant $E_{0}^{\prime}$, playing the role of a trial mass. We define $x=E/E_{0}^{\prime}$ and set $f^{(0)}\equiv f$ for consistency of notation. Then:

if the resonance is known to be unpolarized, then define $\hat{x}_{1}=\mathrm{argmax}[f^{(0)}]$ and stop, else compute $f^{(1)}$; 2. 2.

if a pole or cusp in $f^{(1)}$ is found, then define such point $\hat{x}_{2}$ and stop, else compute $f^{(2)}$; 3. 3.

if $f^{(2)}=\mathrm{const.}$ over a range $[x_{-},x_{+}]$, then define $\hat{x}_{3}=\sqrt{x_{-}x_{+}}$ and stop. Else: if there is a cusp, define such point $\hat{x}_{3}$ and stop; 4. 4.

if no such points exist, then compute $f^{(k)}$, with $k\geq 3$, and continue searching for a singularity $\hat{x}_{k+1}$ .

The mass estimator is then defined as $\hat{M}=2\hat{x}_{k}E_{0}^{\prime}$. When a broad distribution of energies in the center-of-mass frame is accounted for, the analyticity of $f$ is restored. In particular, poles and cusps are regularized into local stationary points. These points are in general displaced from $x=1$ by an amount that vanishes in the limit $\Delta\to 0$. Furthermore, since there may be a multiplicity of such stationary points, a prior on $M$ will be in general needed to disambiguate among them and for an ultimate calibration of the estimator. The determination of the unknown resonance mass is then recast as a univariate optimization problem, in a way that decouples from the details of the underlying production and decay dynamics to the extent that the resonance width can be neglected.

3 Numerical examples

The predictions of Eq. (15) have been verified numerically for selected choices of the functions $g$, $A_{0}$, and $A_{4}$. The three following functional forms for the boost factor p.d.f. have been studied:

• •

$g_{\rm exp}(\gamma)\propto(\gamma-1)\exp^{-(\gamma-1)}$. This function is analytic in $\gamma=1$, and is chosen as the prototype of a p.d.f. with $g^{(0)}=0$.

• •

$g_{\rm pow}(\gamma)\propto(\gamma-0.9)^{-0.8}$. This function is analytic in $\gamma=1$, but this time $g^{(0)}>0$. The numerical values of the coefficients are somehow tuned on the empirical boost distribution for $W$ bosons production at the LHC when $\gamma-1$ is in excess of about $10^{-3}$, see Sec. 2.3.

• •

$g_{\rm sqrt}(\gamma)\propto(\gamma-1)^{\frac{1}{2}}$. This function is chosen as the prototype of a p.d.f which is not analytic in $\gamma=1$. In particular, it is finite for $\gamma\to 1^{+}$, but its first derivative is infinite.

For sake of numerical precision, the integration over the boost factors is restricted to the range $\gamma\in[1,3]$. Both $g_{\rm exp}$ and $g_{\rm pow}$ are not integrable, so strictly speaking they cannot be interpreted as probability density functions. However, they can still provide a good approximation of physical densities for small values of $\gamma$. For all three functions, the following values for $(A_{0},A_{4})$ have been studied:

$(0,0)$, corresponding to an equal left and right polarization; 2. 2.

$(\frac{2}{3},0)$, corresponding to an unpolarized resonance; 3. 3.

$(0,1)$, corresponding to partial transverse polarization; 4. 4.

$(\mbox{tanh}\left[4(\gamma-1)\right],\mbox{tanh}\left[0.3(\gamma-1)]\right)$, corresponding to a resonance which has equal left and right polarization at rest, and then it acquires both a longitudinal and a transverse polarization as the boost factor increases. The choice of numerical constants is somehow inspired by the case of $W$ production as discussed in Sec. 2.3.

The resulting probability density functions $f^{(0)}$ are shown in Fig. 1-3 together with their first $f^{(1)}$ and second $f^{(2)}$ derivatives. The latter are estimated from finite differences of $f^{(0)}$ over an equally-spaced mesh of $x_{i}$ values:

[TABLE]

where $d$ is the mesh size.

Figure 1 shows the results for $g_{\rm exp}$ for each choice of the angular coefficients. Since $g^{(0)}=0$, we have $C=E=H=0$, as for Eq. (16). Apart from case $2)$, where $x=1$ is also a local maximum of $f^{(0)}$, the first derivative does not vanish in general at $x=1$. However, the presence of a term like $|\epsilon^{3}|$ in Eq. (15) induces the presence of a cusp in the second derivative.

For $g_{\rm pow}$, which has $g^{(0)}>0$, additional sources of non-analyticity are present in Eq. (15), clearly visible in Figure 2. In case $1)$, the Taylor expansion of $f^{(0)}$ contains a term of the form $E\epsilon|\epsilon|$ with $E\neq 0$, hence the second derivative receives a contribution from a step-function centered at $x=1$. In case $2)$, $C\neq 0$ so that $x=1$ is a cusp: the first and second derivatives are thus locally proportional to a step-function and a delta function, respectively. In case $3)$, $H<0$, so that the first order diverges to $+\infty$ like $\ln|\epsilon|$ when $\epsilon\to 0$, whereas the second order derivative goes like $1/\epsilon$. Case $4)$ is qualitatively similar to the first.

Whenever $g$ or any of the two angular coefficients are not analytic at $\gamma=1$, like for $g_{\rm sqrt}$, Eq. (15) does not apply. The general appearance of step functions and poles in $x=1$ is however unchanged, as shown by Figure 3. In particular, a term of the form $\epsilon^{2}\ln|\epsilon|$ stems from the last but one line of Eq. (13). The choice of $g_{\rm sqrt}$ is, however, special since, for an unpolarized resonance, it provides an analytic p.d.f.:

[TABLE]

In this last case, the mass estimator would be provided by $\mbox{argmax}[f]$.

4 $W$ bosons at the LHC

A special case of the problem studied in Sec. 2 is represented by $W$ bosons produced in hadron-hadron collisions and decaying into a lepton-neutrino pair. For the purpose of studying this particular process, a sample of proton-proton collision events at $\sqrt{s}=13$ TeV simulating the $pp\to W^{\pm}X$, $W^{\pm}\to\mu^{\pm}\nu_{\mu}$ reaction has been generated with NLO QCD accuracy using the MG5_aMC@NLO amc event generator interfaced with Pythia8 for the parton shower pythia . The NNPDF3.0 nnpdf set is used to simulate the proton PDFs. A total of about $84$ millions of events are generated, with a fraction of negative weights such that the effective number of events is reduced by roughly a factor of two compared to the case of unweighted events. Given that the cross section for $W\to\mu\nu_{\mu}$ production is about $20.5$ nb at $\sqrt{s}=13$ TeV FEWZ , the simulated sample used for this study has the same statistical power of a sample of collision events corresponding to $1.9$ fb*-1* of integrated luminosity.

The natural width of the $W$ boson is $\Gamma_{W}\approx 2.08$ GeV PDG , corresponding to a value $\Delta\approx 10^{-2}$ in Eq. (2). This is not negligible on the scale of a high-precision measurement of $M_{W}$, which targets a relative accuracy on the mass as low as $10^{-4}$ WMassCDF ; WMassD0 ; WMassATLAS . Hence, an ultimate calibration of the estimator is required to meet this level of accuracy. In Sec. 2.3, it was found that the boost factor p.d.f. for $W$ bosons produced in proton-proton collisions can be roughly approximated by a power law $g\sim(\gamma-1)^{\alpha}$: for small values of $\gamma-1$, i.e. $\lesssim 5\times 10^{-4}$, we have $\alpha\approx 0.5$ and the $W$ boson is almost unpolarized in the helicity frame; for higher boost values, $\alpha\approx-0.8$, and a net polarization is built, ultimately dominated by a particular transverse mode. From the numerical simulations of Fig. 2-3, we could thus expect $x=1$ to be a local minimum of $f^{(2)}$. Indeed, the rising edge of $g$ populates only the region $|x-1|\lesssim 5\times 10^{-4}$, where it provides a smooth function $f^{(2)}$, similarly to the rightmost panel in the second row of Fig. 3. For larger boosts, $f^{(0)}$ should resemble more closely the plots in Fig. 2, featuring a deep minimum of $f^{(2)}$ at $x=1$. The whole picture is then smeared by the finite width of the $W$ boson.

Figure 4 shows the binned density $f^{(0)}$ obtained from the simulated sample of events. The first and second derivatives are estimated bin-wise in the same fashion as Eq. (26). A deep minimum in the histogram of $f^{(2)}$ at $x$ values close to unity is clearly visible. Interestingly, $x=1$ is also close to be a global maximum of $f^{(0)}$, a result that qualitatively recalls the last toy example in Fig. 2, where the boost p.d.f. and the angular coefficients were indeed tuned on the values extracted from the Monte Carlo simulation.

4.1 Detector acceptance and final-state radiation

For the case of $W$ boson production and decay at the LHC, two effects break the mathematical hypotheses assumed to derive Eq. (8): the presence of acceptance selection requirements, which are unavoidable in experiments at hadron colliders, and the emission of final-state photon radiation (FSR) from the charged lepton. Both affect the center-of-mass dynamics, albeit in different ways as discussed below.

When the detector coverage is incomplete, the harmonic polynomials that depend on $\phi^{*}$, which are themselves proportional to $\sin\theta^{*}$ and $\sin 2\theta^{*}$ Mirkes , don’t average exactly to zero in some regions of the phase-space, thus adding spurious terms to the $c^{*}$ expansion of the decay angle distribution. Furthermore, the detector acceptance requirements, being based only on the kinematics of the visible decay products, affect the lepton reconstruction efficiency differently depending on the kinematics of the $W$ boson. The overall result is to modify the angular coefficients by boost-dependent efficiency factors $\rho_{0,4}(x\,|\,\gamma)$ in Eq. (9), which are in general non-trivial functions of $x$. Here, we will study the effect of selection requirements realistic for general-purpose experiments like ATLAS ATLAS and CMS CMS , namely $|\eta|\leq 2.5$ and $p_{\rm T}\geq 25$ GeV, where $\eta$ and $p_{\rm T}$ are the muon pseudorapidity and transverse momentum, respectively. In the MC simulation we find these cuts to have an efficiency of about $77\%$ for $W^{+}$ and $84\%$ for $W^{-}$ for events with lepton energy $E\approx M_{W}/2$.

The emission of FSR by the charged lepton perturbs the center-of-mass dynamics. The overall effect of such perturbation can be thought of as the convolution of the original harmonic polynomials with a smearing kernel, which introduces infinite harmonics in $c^{*}$. Furthermore, the emission of extra particles (photons and lepton-pairs) reduces the center-of-mass energy available for the muon and thus primarily affects the visible energy spectrum by an overall downward shift. This process is well-known Calame though, so that it could be in principle unfolded at the detector level to recover a pure QCD description of the final-state kinematics. For the purpose of studying this process, we will consider both an unrealistic scenario, where the charged leptons do not undergo photon radiation (pre-FSR leptons), and a realistic scenario where a QED-shower of the muons is simulated by the Pythia8 MC (bare leptons).

4.2 The search for a stationary point

We now consider the problem of finding the stationary points of the higher-order derivatives of $f$. The rightmost histogram in Fig. 4 clearly shows that a local minimum of $f^{(2)}$ is present at $x\approx 1$. The estimator of such point is, however, not uniquely defined. We won’t address here the problem of finding the statistically optimal of such estimators. Instead, we decide to define the estimator implicitly as the root of a conveniently chosen function of the data. To this purpose, we first approximate the density $f$ with a polynomial function of degree $D$ centered at $x=1$:

[TABLE]

The coefficients $c_{n}$ in Eq. (28) are determined from a fit to the simulated data by means of an analytic $\chi^{2}$ method. We then define the stationary points of the $i$-th order derivative implicitly as the roots of the $(i+1)$-th order derivative. The latter are determined numerically by using Halley’s root-search method Halley , a variant of the classical Newton method. Statistical uncertainties on the coefficients of the polynomial fit are propagated to the roots $\hat{x}_{i}$ by means of pseudo-data, resulting in 68% confidence level (CL) intervals. This approach has a twofold advantage: it regularizes the statistical bin-by-bin fluctuations by the use of smooth functions and it allows for an analytic evaluation of the derivatives at any point $x$.

The energy spectrum is provided as a histogram with 100 MeV large bins. The central value of each bin is normalized to the constant $E_{W}=M_{W}/2$ to yield the dimensionless variable $x$. The fit is performed in the interval $E\in[36.2,44.3]$ GeV, corresponding to invariant masses of the $W$ bosons in a window of about $\pm 4\Gamma_{W}$ around $M_{W}$. Such range is large enough to provide acceptable fits with $D=4$, which is the minimum degree to define a unique root of the third derivative $\hat{x}_{3}$. We notice that this way of estimating the roots $\hat{x}_{i}$ is quite sensitive to border effects related to the choice of the fit range: since Eq. (28) is only a local approximation of the density, discrepancies between the true spectrum and $f^{(0)}$ at the edges of the fit range tend to pull more strongly the coefficients associated with the large powers of $n$, which in turn affect more strongly the roots of the higher order derivatives. The bias associated with the choice of fit range will be eventually reabsorbed as part of the calibration procedure discussed in the next section.

For the sake of comparison, the root of the first derivative $\hat{x}_{1}$, which corresponds to a local maximum of $f^{(0)}$, is also studied. Positive and negative muon events are first considered separately. Since the two samples of events provide consistent results, they are ultimately combined to maximize the statistical accuracy of the analysis. The result of the fit to the simulated data is shown in Fig. 5, together with the first, second, and third derivative of the fitted polynomial function. As expected, the second derivative features a local minimum around $x=1$ identified by the root $\hat{x}_{3}$ of the third derivative.

4.3 Calibration curve

The calibration of the $\hat{x}_{i}$ estimator is determined by reweighting the same MC sample to different values of $M_{W}$. The fit is then repeated for each mass-reweighted sample and new roots $\hat{x}_{i}$ are computed, resulting in a calibration curve $\hat{M}_{W}=\hat{M}_{W}(\hat{x}_{i})$. Figure 6 shows such curves separately for muons in the full phase-space and for muons within the detector acceptance as defined in Sec. 4.1. The points $\hat{x}_{i}$ are then interpolated through a linear regression.

The response of $\hat{x}_{3}$ to a change of $M_{W}$ is found to be linear to better than 1%. This fact is reassuring and confirms that $\hat{x}_{3}$ is indeed a good estimator of $M_{W}$. For comparison, the root of the first derivative $\hat{x}_{1}$, and the mean value $x_{\mu}$ in the same range of $x$ values considered in the fit, are also reported as a function of $M_{W}$. The former is found to have a good linear response but a larger offset compared to $\hat{x}_{3}$ (800 MeV against 100 MeV). Instead, the mean value $x_{\mu}$ is very mildly related to $M_{W}$, which makes it a rather poor estimator of the mass. This is however mainly an artifact of considering a narrow range of $x$ values: as illustrated by the first panel of Fig. 5, the function $f^{(0)}$ in the neighborhood of $x=1$ is a concave function: a tiny shift $\delta x$ of the peak position does not change the mean of the distribution to first-order in $\delta x$. It is also interesting to study the response of the three estimators to a restriction of the lepton phase-space. This is shown in the right panel of Fig. 6. Both ${x}_{\mu}$ and $\hat{x}_{1}$ are found to be significantly affected by acceptance selection requirements, i.e. their values change compared to the full-acceptance case by more than their statistical uncertainty. On the contrary, $\hat{x}_{3}$ is more stable, changing by less than one standard deviation.

Finally, an identical analysis is repeated considering bare leptons instead of pre-FSR leptons. The results are shown in Fig. 7. Besides an overall shift of about $200$ MeV, which can be ascribed to the loss of energy drained away by FSR Calame , the linear response of $\hat{x}_{1}$ and $\hat{x}_{3}$ to $M_{W}$ is found to be preserved.

4.4 Residual model dependence

The residual model-dependence of $\hat{x}_{i}$ on the production and decay dynamics will be eventually incorporated as a systematic uncertainty on the calibration curve. For the purpose of assessing the level of such relic model-dependence, the same MC sample has been reweighted to nine different sets of values for the renormalization and factorization scales ($\mu_{R}$, $\mu_{F}$) amc and to one-hundred MC replicas of the same PDF fit nnpdf . In the first case, the maximal deviation of the roots compared to the value obtained for the central choice of scales is symmetrized and used as a proxy of the perturbative scale uncertainty on $\hat{x}_{i}$. In the second case, the RMS of the distribution of roots is considered as systematic uncertainty related to the imperfect knowledge of the PDFs.

The results are shown in Table 1. The variance of the weights used to modify the simulated sample contributes to these uncertainties. The latter is assessed by means of pseudo-experiments where the data in the nominal sample are randomized to account for the extra uncertainty added by the reweighting. We find the additional statistical fluctuation introduced by the reweighting to contribute to $\hat{x}_{i}$ by an amount corresponding to about 30% (10%) of the Poisson uncertainty on the fitted value for $i=1$ ($3$). The reduction in statistical uncertainty after applying the acceptance selection requirements, most noticeable for $\hat{x}_{1}$, is due to the change of the functional form of the energy spectrum. Indeed, if the selection efficiency were independent of $E$, the effect of the acceptance requirements would be to enlarge $\sigma_{\rm stat}$ by some 10% due to the reduced event yield. Instead, the acceptance selection efficiency increases with $E$: the net effect is to sculpt the energy spectrum in such a way that the curvature of $f^{(0)}$ around $x=1$ increases, thus reducing the statistical uncertainty on the position of the maximum.

We can now summarize the results of this study as follows:

The root $\hat{x}_{3}$ is less affected than $\hat{x}_{1}$ by changes in the modeling of $W$ production and decay dynamics induced by different choices of perturbative scales and PDFs. 2. 2.

The root $\hat{x}_{1}$ shows a significant systematic uncertainty, i.e. larger than the Poisson fluctuation introduced by the reweighting. This corroborates the observation that $\hat{x}_{1}$ depends more than $\hat{x}_{3}$ on the underlying dynamics, as also deduced by its larger sensitivity to the acceptance requirements. 3. 3.

The residual scale and PDF uncertainty on $\hat{x}_{3}$ is consistent with the Poisson fluctuation introduced by the reweighting; within the statistical accuracy of this study, there is no indications of residual systematic bias. 4. 4.

A statistical-only uncertainty on $\hat{x}_{3}$ corresponding to $15$ MeV uncertainty on $M_{W}$ could be reached with about 300 fb*-1* of LHC data, which is within the reach of the Run 3 of the LHC.

4.5 Background

The analysis discussed in the previous paragraphs neglects the presence of background sources. Since the measurement of a stationary point $\hat{x}_{i}$ relies on a local description of the energy p.d.f. in the neighborhood of $x=1$, any background with a non-flat p.d.f. has the potential to bias the mass estimator. For the case of $W$ boson production at hadron colliders WMassCDF ; WMassD0 ; WMassATLAS , three major background processes should be considered: $i)$ multi-jet production, where the muon comes from hadron decays within a jet, $ii)$ Drell-Yan production of a muon pair, where one of the muons escapes detection, and $iii)$ top quark ($t$) production. The latter two cases are the least harmful. Indeed, neutral Drell-Yan events pass the event selection criteria to the extent that one of the muon is emitted with either soft $p_{\rm T}$ or large values of $|\eta|$. In turn, this condition preferentially selects events where the intermediate $Z/\gamma^{*}$ boson is produced with a finite boost. As for Eq. (15) with $g^{(k)}=0$, this implies that the $E$ spectrum of the selected muons in the neighborhood of $M_{Z}/2$ must have a flat second-order derivative. Similarly, muons in top quark events come from the decay of boosted $W$ bosons, since, in the rest frame the decaying $t$ quark, the $W$ boson recoils againts a $b$ quark with $|\bm{p}_{W}^{*}|\approx 0.4\,m_{t}$. For multi-jet production, these arguments do not hold and a detailed data-driven estimation of the functional form should be performed. However, we remark that the analysis discussed here is robust against changes in muon acceptance, as observed in Sec. 4.1. Since the multi-jet background is reducible by either tighter identification criteria on the muon or by stricter requirements on the missing energy or transverse mass in the event, we do expect room for optimization in case the functional form of this background were found to be measurable with only limited accuracy.

4.6 Outlook

A more refined analysis of the residual theoretical uncertainties would require the simulation of a much larger data sample and a careful treatment of other model effects, like non-perturbative physics, mixed QCD-QED corrections, etc. (see e.g. Ref. Alioli for a recent review). Similarly, experimental uncertainties from the backgrounds, the bias due to the choice of a fixed-order polynomial to fit the data, the impact of the lepton energy scale uncertainty, etc., should be thoroughly assessed. This is beyond the scope of this work. The study presented here confirms that a stationary point in the second derivative of the lepton energy density is a good estimator of $M_{W}$ and that it is robust against changes of the underlying $W$ boson production and decay dynamics, detector acceptance requirements, and the emission of photon radiation.

5 Conclusions

We have considered the two-fermion decay of a spin-1 resonance of mass $M$, and analyzed the lepton energy distribution in the laboratory frame in full generality, i.e. regardless of the underlying production and decay dynamics of the resonance. In particular, we have studied the analyticity of the probability density in the neighborhood of $M/2$. We find that the density at this point is not analytic for a narrow-width resonance. In particular, we have studied the conditions for which a singular point appears in the higher-order derivatives of the density, and found that the second derivative is likely to develop a cusp or a pole. Exact formulas have been derived under the assumption that the distribution of boost factors $\gamma$ and the polarization of the resonance are described by regular functions of $\gamma$. The formulas have been qualitatively validated with toy examples of production and decay of a narrow-width resonance. When a finite width of the resonance is accounted for, the regularity is restored such that cusps or poles are smoothed into local stationary points potentially displaced from $M/2$. The size of such displacement depends on the width of the resonance, but partially also on the production dynamics, thus requiring an ultimate calibration. The quest for stationary points in the higher-order derivatives of the energy density function is thus advocated as a way to estimate $M$ with possibly limited knowledge of the underlying production and decay dynamics of the resonance. A special case is represented by the production of $W$ bosons in proton-proton collisions, which has been studied on a Monte Carlo simulation of this process, assuming LHC-like conditions on the proton beams. As expected, a stationary point in the second derivative is found close to $M_{W}/2$. The robustness of this point as an estimator of $M_{W}$ has been studied by considering the effect of detector acceptance requirements, the emission of final-state radiation, changes of the perturbative calculation of scattering amplitudes for $W$ production, the proton PDFs, and the input $W$ boson mass. Interestingly, such a mass estimator features a good linearity, a small bias, and is rather resilient to changes in the lepton acceptance and in the modeling of the $W$ boson production dynamics. An ultimate assessment of the residual model-dependence is left for future work.

Acknowledgements.

We would like to thank Roberto Franceschini for a fruitful discussion during the CMS Italian National Meeting held in Piacenza in 2017, which triggered this work. We are also indebted to Nico Kleijne for reading the manuscript and double-checking the formulas.