Singularity variables for missing energy event kinematics

We discuss singularity variables which are properly suited for analyzing the kinematics of events with missing transverse energy at the LHC. We consider six of the simplest event topologies encountered in studies of leptonic W -bosons and top quarks, as well as in SUSY-like searches for new physics with dark matter particles. In each case, we illustrate the general prescription for finding the relevant singularity variable, which in turn helps delineate the visible parameter subspace on which the singularities are located. Our results can be used in two different ways — first, as a guide for targeting the signal-rich regions of parameter space during the stage of discovery, and second, as a sensitive focus point method for measuring the particle mass spectrum after the initial discovery.


Introduction
Events with missing transverse energy (MET) at the Large Hadron Collider (LHC) are of great interest to both theory and experiment. On the experimental side, the MET is a very challenging object, and a great amount of effort has gone into the proper calibration of the detector and in the evaluation of its missing energy performance in both ATLAS [1] and CMS [2]. On the theoretical side, events with MET are likely to hold the key to understanding some of the great unsolved puzzles of the Standard Model (SM). For example, the dark matter problem greatly motivates searches for new physics beyond the Standard Model (BSM) with dark matter candidates [3], which would escape the detector leaving a MET signature [4]. Similarly, the flavor problem provides impetus for focusing on the third generation in the SM, where the top quark, the bottom quark and the tau lepton all have decay channels with invisible particles (neutrinos) in the final state. Finally, the W -bosons, JHEP04(2020)027 whose leptonic decays exhibit a classic MET signature, have long been considered promising probes of the electroweak symmetry breaking sector [5,6], and more recently have rounded out the suite of Higgs discovery channels [7,8]. Therefore, the sound theoretical understanding of the event kinematics in MET events should be a high priority.
The fundamental problem with MET events is the incomplete information about the final state, since the energies ε i and momenta q i of the invisible particles (neutrinos or dark matter particles) are not measured. At the same time, there is partial information available in the form of the energies e j and momenta p j of the visible final state particles, which typically are the (approximately massless) leptons, photons and/or QCD jets, but may also be massive 1 reconstructed visible particles like a W -boson, a Z-boson or a Higgs boson. In the spirit of the simplified models approach [9], in what follows we shall remain agnostic about the underlying physics which might give rise to a particular event topology, and instead shall focus on the salient features of its kinematics as represented by phase space singularities. Correspondingly, we shall also ignore any secondary dynamical effects such as spin correlations, helicity suppressions, etc., since they do not affect the singular phase space features which we are interested in.
With this backdrop, we are ready to introduce the three relevant (and related) questions regarding MET events which will be addressed in this paper.
• Are there any singular features in the visible phase space of a given event topology?
• What are the relevant kinematic variables which best describe such singular features?
• What do measurements of such features tell us about the underlying mass spectrum?
Let us briefly motivate and comment on each question before proceeding to the main analysis in the following sections.

Visible phase space singularities
In this paper, we define a singularity in the visible parameter space { p j } as a point where the event number density formally becomes infinite. 2 The origin of such singularities is very well understood [10][11][12][13]: they arise in the process of projecting the allowed region in the full phase space { q i , p j } (which does not exhibit any singularities) onto the visible subspace { p j }. Similar to the phenomenon of caustics in optics, astrophysics [14] or accelerator physics [15], singularities are formed at points where the visible projection onto { p j } of the allowed phase space in { q i , p j } gets folded. Mathematically this is expressed as the reduction in the rank of the Jacobian matrix of the coordinate transformation from the relevant set of kinematic constraints to { q i } (alternatively, from the generator-level event parameters to the visible space { p j }), which is why such singularities are sometimes known as Jacobian peaks. 1 For this reason, in our analysis below we shall try to retain the visible particle masses mj ≡ e 2 j − p 2 j as arbitrary whenever possible. 2 In the literature, kinematic endpoints, cusps and kinks are sometimes also referred to as singularities [10], even though they are defined in terms of suitable derivatives of the event number density. In what follows, we shall adopt the more narrow definition of a phase space singularity as stated in the text above. JHEP04(2020)027

Singularity variables
The visible parameter space { p j } may be parametrized simply in terms of the visible momentum components, as indicated, but it also allows infinitely many alternative reparametrizations, involving, e.g., angles, rapidities, invariant masses, etc. An interesting question then is which of those re-parametrizations is "the best". In his pioneering paper [10], Kim proposed to construct an optimized one-dimensional kinematic variable, called a singularity coordinate, which is defined in terms of the measured visible momenta { p j }, plus possibly some mass parameters {M k }. The latter set includes the masses ε 2 i − q 2 i of the invisible particles in the final state, as well as the masses of any intermediate resonances in the event topology. As the name suggests, the singularity coordinate is designed to capture the singular behavior and so must satisfy the following criteria put forth in [10]: (i) it must vanish at the singularity locations; (ii) its direction must be perpendicular to the singularity hypersurface in the observable phase space { p j }; (iii) events which are equally far away from the singularity should produce the same value.
Unfortunately, Kim's paper went largely unnoticed -in the past ten years, there have been very few explicitly worked out examples of practical significance, with the exception of two follow-up investigations by De Rujula and Galindo, who introduced and studied several different versions of a singularity coordinate for the case of single W production [11] and h → W W [12]. Our goal in this paper is to expand the set of worked out examples, on occasion taking the opportunity to point out connections to other results in the literature which have been obtained by different means. For completeness, we shall also review and further expand on the two case studies in [11,12].
In parallel with the one-dimensional approach of a singularity variable, we shall also try to analyze and visualize the singularity hypersurface from a multi-dimensional perspective. 3 For example, consider the classic supersymmetry (SUSY) example of a single decay chain proceeding through three successive two body decays. The relevant visible parameter space is three-dimensional and can be parametrized by the pairwise invariant masses of the final state visible particles. The singularity is then found at the two-dimensional surface boundary of the allowed region [16][17][18][19]. In order to maximize the sensitivity, the experimental analysis should target this two-dimensional surface; this can be done either directly [20][21][22] (e.g., using Voronoi tessellations of the data [23,24]) or by means of a suitable one-dimensional singularity coordinate [25]. parameters, and finally, a set {M k } of mass parameters, which, depending on the circumstances, may or may not be known a priori. As already mentioned, {M k } would typically include the masses m i = ε 2 i − q 2 i of the invisible particles in the final state. Now, if one is willing to assume that the invisible particles are all neutrinos, as is usually 4 the case with studies of SM signatures, those masses can be safely set to zero. On the other hand, in studies of BSM signatures the invisible particles are some new dark matter particles whose masses are a priori unknown and should be explicitly retained in (1.1). Also included in the set {M k } are the masses of intermediate resonances, which are known if the resonance is a SM particle and unknown otherwise.

Using singular features for particle mass measurements
The constraint (1.1) can be viewed in different ways. First, the function g can be regarded as a bona fide singularity coordinate in the sense of ref. [10]. Second, given a choice of masses {M k }, the constraint (1.1) defines the locus of points where the signal event density becomes singular (by construction), while the behavior of the background event density is typically unremarkable. Therefore, this locus of points is precisely the region of phase space which should be targeted (with suitable selection cuts) in an analysis aimed at a discovery [25]. Finally, we can turn the last argument around and instead of studying the phase space { p j } for a given choice of mass parameters {M k }, we can study the mass parameter space {M k } for given points in phase space as sampled by the events in the data. Correspondingly, for each event, eq. (1.1) can be viewed as a constraint on the allowed values of the mass parameters for which the current event would be located on a singularity hypersurface. In a companion paper [13], this idea was formulated as a new mass measurement method, 5 called the "focus point method", which was inspired by previous related work in [29][30][31][32][33]. The analysis in [13] was illustrated with one specific event topology, dilepton tt events, and we shall now show that it is also applicable for the remaining 5 event topologies considered here.
The paper is organized as follows. In section 2.1 we first introduce the six event topologies to be studied in this paper, along with our notation and conventions, as well as some details on our simulations. Then in section 2.2 we flesh out the general method for deciding whether a singularity feature exists, and if so, for deriving the corresponding constraint (1.1). In the next five sections 3-7, we illustrate the general method for each individual event topology. In the process, we shall sometimes rederive some existing results in the literature, albeit in a different, universal and perhaps more intuitive way. We hope that the expert reader will enjoy seeing the underlying commonality between those familiar results, as well as appreciate the novelty of the others. The novice reader is perhaps best advised to first read refs. [10][11][12][13] and be ready to consult a standard reference like [34,35] when necessary. Finally, in section 8 we present our conclusions and outlook for future studies. 1 Figure 1. The event topologies studied in this paper. Diagrams (a-c) in the top row represent single decay chains, each terminating in a single invisible particle with mass M 0 and 4-momentum q = (ε, q), while diagrams (d-f) in the bottom row represent pair-production processes leading to two decay chains and two invisible particles with 4-momenta q = (ε, q) and Q = (E, Q), respectively. p i = (e i , p i ) and P i = (E i , P i ) label the 4-momenta of visible final state particles of mass m i , while M k with k = 1, 2, 3 label the masses of on-shell intermediate resonances.

Notations and setup
The six event topologies to be considered in this paper are shown in figure 1, where, in accordance with the notation of [36], the letters p and q are reserved for denoting momenta of visible and invisible particles, respectively. Diagrams (a-c) in the top row represent single decay chains, each terminating in a single invisible particle with mass M 0 and 4-momentum q = (ε, q), which is taken to be on-shell: 6 ε 2 − q 2 = M 2 0 . The decay chains in figure 1(a-c) differ by the number of two-body decays -here we have limited ourselves to at most three (in figure 1(c)), but the generalization to four or more is straightforward [22]. Similarly, p i = (e i , p i ), with i = 1, 2, 3, label the 4-momenta of visible final state particles of mass m i , with e 2 i − p 2 i = m 2 i . The masses M k with k = 1, 2, 3 label the masses of on-shell intermediate resonances. The diagrams (d-f) in the bottom row of figure 1 represent pairproduction processes leading to two decay chains and thus to two invisible particles with 4-momenta q = (ε, q) and Q = (E, Q), respectively. The 4-momenta of the visible particles JHEP04(2020)027 Value (in GeV) 700 800 1000 1300 in the second decay chain are similarly capitalized: The six diagrams in figure 1 cover a host of interesting physics processes at the LHC, both within and outside the SM. For any given such process, there will be a corresponding choice for the masses m i , M 0 and M k . Consider, for example, single top decay as an illustration for figure 1(b). In that case, m 1 is the mass of a lepton (electron or muon), m 2 is the mass of the reconstructed b-jet, M 0 is the neutrino mass, M 1 is the Wboson mass and M 2 is the top quark mass. Any other process will dictate a similar choice of mass parameters; and of course, BSM processes may involve dark matter particles and intermediate resonances with a priori unknown masses M 0 and M k . In order to keep everything on the same footing, in our numerical examples we shall use the same mass spectrum for all six processes as shown in table 1.
Events will be generated at parton level with MadGraph [38]. Intermediate resonances will be decayed by phase space, i.e., ignoring any spin correlations. This assumption will not impact our results, since they are based on purely kinematics arguments and are thus independent of the model dynamics. For clarity of the presentation, and to better see the singularity structures, we will not employ any detector simulation and will keep the intermediate resonances strictly on-shell. Of course, in a real experiment there will be some smearing due to both the detector resolution and the finite particle widths; the extent of those effects depends on the particular realization of the topology in question and is beyond our scope here. The visible final state particles, which are typically leptons and QCD jets, will be taken to be massless. As we shall see below, in the case of figure 1(a), the presence of initial state radiation (ISR) is important and actually quite beneficial for the analysis. This is why ISR will be simulated for that case, while in the remaining event topologies the ISR will not be modelled. Finally, we shall ignore the combinatorial problem arising in figure 1(f), since it has already been addressed in the literature [39][40][41], and for simplicity we will assume that the visible particles have been properly assigned to the two decay branches.

The basic idea
We are now ready to outline the general method for deriving a singularity variable [10][11][12][13]. The first step is to collect all kinematic constraints on the invisible momenta q (and Q, if JHEP04(2020)027 applicable). In general, the kinematic constraints fall into several categories: • On-shell relations. These exist for the final state invisible particles as well as for the intermediate resonances: (2.3b) Notice that the equations (2.2b) and (2.3b) are only applicable in the case of the bottom row diagrams in figure 1.
• Missing transverse momentum constraint. Momentum conservation in the transverse plane constrains the sum of the transverse momenta of the invisible particles: (2.4) • Additional constraints due to the specific experimental setup. Additional constraints may arise at specific experimental facilities, e.g., at a lepton collider, where the center-of-mass (CM) energy and momentum of the initial state are known.
Combining all of these constraints, we obtain a set of N C equations 7

JHEP04(2020)027
In what follows, we shall focus on situations where the number of constraints N C is just enough so that one can solve (2.5) and (2.6) for the invisible momenta in terms of the mass parameters {M k }. In other words we shall always have 8 Note that this does not imply that the final state kinematics is fully solved -we are just trading one set of unknowns, the components of the invisible momenta, for another: the masses {M k } of the intermediate resonances and of the invisible final state particles. In other words, we are still dealing with an undetermined problem, in which we are not able to compute the exact momenta of the invisible particles in the event.
In order to illustrate the basic idea, it is sufficient to consider the simpler version (2.6) of the kinematic constraints -the same argument goes through for any number of invisible particles in the final state, as long as (2.8) is in effect. Consider some particular solutioñ q µ of (2.6): A singularity atq is obtained when at least one of the directions of the local tangent plane to the full phase space is aligned with an invisible momentum direction [10]. This means that we can make an infinitesimal change δq in the unmeasured invisible 4-momentum components while continuing to satisfy the original system, i.e. The condition (2.8) ensures that the Jacobian matrix is a square matrix. The existence of non-trivial solutions δq µ for the system (2.11) is guaranteed if the determinant of D αµ vanishes: Note that the left-hand side of this equation is a function of only visible momenta {p} and mass parameters {M k }, since the invisible momenta can be eliminated via eqs. (2.9). Comparing to (1.1), we see that the left-hand side of (2.13) can be taken to be the desired singularity coordinate. 9 Repeating the same argument for the case with two invisible particles, i.e., starting from eq. (2.5), leads us to a similar condition ∂f α ∂(q, Q) = 0. (2.14)

JHEP04(2020)027
From a mathematical point of view, eqs. (2.13) and (2.14) are simply the reduced rank conditions leading to a critical point. In our case their importance lies in the fact that the distribution of the corresponding kinematic variable is singular at the critical point, which is why refs. [11,12] also referred to these reduced rank conditions as "singularity conditions".
In the subsequent sections we shall explore the implication of (2.13) or (2.14) for the event topologies of figure 1. As outlined in the introduction, for each event topology we shall focus on the following three issues: • Derivation of the relevant singularity coordinate.
• Delineation of the signal-rich regions of the visible phase space, i.e., where the signal density becomes singular. In doing so, we shall be careful to use the symmetries of the problem in order to maximally reduce the dimensionality of the observable phase space without washing out any singular kinematic features.
• Demonstration of the focus point method for mass measurements proposed in [13].

Single decay chain, one two-body decay
In this section we shall consider the single two-body decay diagram from figure 1(a). We shall revisit and expand the discussion in ref. [11], which showed that the singularity coordinate in this case is nothing but the usual transverse mass variable m T [43,44]. By now, the transverse mass is one of the standard kinematic variables, which has been widely used in precision measurements of the W -boson mass [45][46][47] as well as in new physics searches for W resonances [48][49][50][51]. The new element in our discussion will be the role of the focus point method of [13] and its connection to the kink method for mass measurements [52][53][54][55][56].

Derivation of a singularity coordinate
Let us begin by listing the four kinematic constraints for figure 1(a): or after a simple rearrangement,
The latter is nothing but the equal rapidity condition which is known [36] to provide the link between transverse invariant mass variables (like the transverse mass m T , the Cambridge m T 2 and others) and their respective 3+1 dimensional analogues [57][58][59][60][61][62][63][64][65][66][67][68]. In order to cast the singularity condition in the desired form (1.1), we must eliminate the invisible 4-momentum q by using four out of the five equations appearing in eqs. (3.2) and (3.5). For example, using (3.2a), (3.2c), (3.2d) and (3.5), one can find the components of the 4-vector q = (ε, q x , q y , q z ) as Substituting the result (3.6) into the remaining fifth equation (3.2b), we obtain the final singularity condition explicitly in the form where in the left-hand side we recognize the transverse mass variable m T written as If we introduce the transverse energies (3.11) Note that while we have succeeded in eliminating the unknown components of the invisible 4-momentum q, there still remains one a priori unknown parameter, namely, the mass M 0 of the invisible particle, which enters through the transverse energy ε T . The singularity condition (3.7) can then be rewritten in the very compact form This confirms that, up to the additive constant M 1 , the relevant singularity variable for the event topology of figure 1(a) is indeed the transverse mass m T . Furthermore, eq. (3.12) shows that the singularity occurs at the mass M 1 of the parent particle, i.e.
However, there is an important subtlety -the latter statement is true only if we have made the correct choice for the invisible mass parameter M 0 entering the definition (3.11).
In general, the true value of M 0 is a priori unknown, and we have to adopt a certain ansatz for it (denoted in the following byM 0 ) in order to compute the transverse mass m T from (3.11). 10 Once the ansatzM 0 differs from the true value M 0 , the existence of a singularity in the m T (M 0 ) distribution is generally not guaranteed, and the singular feature in the dN/dm T distribution predicted by eqs. (3.12) and (3.13) will be washed out. This is the key observation behind the focus point method for mass measurements proposed in [13]. There the method was illustrated for the case of the dilepton tt event topology of figure 1(f). In section 3.3 below we shall show that it is also applicable to the simple event topology of figure 1(a) as well. But first let us understand the phase space geometry behind the singularity condition (3.12).

The phase space geometry of the singularity condition
In general, the parent particle is produced inclusively, and the event depicted in figure 1(a) would also contain additional visible particles due to initial state radiation (ISR), or decays upstream to the parent particle. Let us denote the total transverse momentum of these additional visible particles as P ISR T , which can be measured in the detector. Then, by 10 In some sense, the situation here is analogous to the well-known behavior of the kinematic endpoint which is obtained by considering the whole sample of events and finding the maximum value of mT . Since mT (M0) ≤ M1 by construction, the measured endpoint value m max T (M0) can be interpreted as the corresponding massM1 of the parent particle for this choice ofM0:

JHEP04(2020)027
definition, the missing transverse momentum / P T is due to the invisible recoil against both p 1T and P ISR T : (3.14) 3.2.1 The case of no upstream visible momentum: P ISR In order to gain some intuition, it is instructive to first consider the simpler case of no upstream visible momentum, P ISR T = 0. In that case, (3.14) reduces to / P T = − p 1T and the transverse mass (3.11) can be written simply as Note that in this case m T is a function of only one degree of freedom, namely, the magnitude of p 1T . This can be easily understood in terms of the symmetries of the problemp 1z does not enter due to the invariance under longitudinal boosts, while the direction of p 1T is irrelevant due to the azimuthal symmetry. The singularity condition (3.12) can then be written as The set of points satisfying this relation belong to a circle in the p 1T plane which is centered at the origin and has radius equal to where λ(x, y, z) ≡ x 2 +y 2 +z 2 −2xy −2xz −2yz. Such points were categorized as "extreme" in [13], since they delineate the boundary of the allowed phase space, where the solutions for the invisible momenta become degenerate. This is pictorially illustrated 11 in figure 2, where the top left panel corresponds to the current case of P ISR T = 0. We plot the event number density (as indicated by the color bar) in the p 1T plane, which we choose to parametrize as (p T , p T ⊥ ), where p T (p T ⊥ ) is the component in the direction along (orthogonal to) P ISR T [69,70]. We see that the allowed phase space in the p 1T plane is indeed a circle, and furthermore, that the maximal event density is found along the circumference of the circle, in agreement with (3.16). With the mass spectrum from table 1, eq. (3.17) predicts the radius of the circle to be p max 1T = 93.75 GeV, which is confirmed in the top left panel of figure 2. Since the extreme events along the circumference of the circle have the same value of p 1T = p max 1T , they will also share the same value of m T , regardless of the choice of test massM 0 , see eq. (3.15). This means that for any value ofM 0 , the m T (M 0 ) distribution will continue to exhibit a singularity at For simplicity, in figure 2 the visible particle is assumed massless (m1 = 0) in which case (3.17) reduces to  We see that in all three cases, the m T distribution has a very sharp singularity at its upper kinematic endpoint (3.19). However, this situation is rather atypical -we shall see below that, in general, when we use the wrong value for the invisible mass parameterM 0 , the JHEP04(2020)027 singularity will be washed out. The reason why it persists here is that the singularity coordinate (3.15) is parametrized by a single degree of freedom, p 1T .

3.2.2
The case with non-zero upstream visible momentum: P ISR We are now in position to discuss the more general case of non-vanishing upstream visible momentum, P ISR T = 0. Using (3.14), the transverse mass formula (3.11) becomes [36] Since P ISR T breaks the azimuthal symmetry, the transverse mass is now a function of two visible momentum degrees of freedom, namely both the magnitude and the direction of p 1T in the transverse plane. We can parametrize the latter degree of freedom by the angle ϕ measured with respect to the direction defined by P ISR T , in which case the doubly projected transverse components of p 1T used in figure 2 are given by As before, let us find the locus of points which satisfy the singularity condition (3.12), now in the presence of non-zero P ISR T . For simplicity, let us only focus on massless visible particles, m 1 = 0, which is an excellent approximation for leptons and jets. In that case, the singularity condition (3.12) reads which can be solved to give the location of the singularity surface in parametric form (see the left panel in figure 4) As a consistency check, we see that eq. (3.24) reduces to (3.18) while a and b are the semi-major and semi-minor axes, respectively:  is always oriented along the positive p T axis, the recoil of the mother particle is in the negative p T direction, which explains the increasing preference for negative p T values as P ISR T gets larger (see also the right panel in figure 4). At the same time, the (doubly) transverse components p T ⊥ are unaffected by the boost of the mother particle, and the maximal p T ⊥ value in each panel stays the same. This is also reflected in the fact that the semi-minor axis (3.27) of the ellipse remains constant, independent of P ISR T . In deriving the equation of the singularity surface (3.24), we have achieved our main goal for this subsection. From here on, how the result (3.24) will be used in practice, depends on the specific purpose of the experimental analysis. If the aim is a discovery of signal events with the event topology of figure 1(a), one should study the distribution of events in the three-dimensional space (p T , p T ⊥ , P ISR T ) depicted in the left panel of figure 4, where the signal-rich regions will be found at the locations singled out by eq. (3.24). If, on the other hand, the goal is to measure the mass spectrum, the singularity surface contains all the kinematic information to do that as well, and the masses M 0 and M 1 can be extracted from a parameter fit to eq. (3.24).
The fitting procedure will have to account for the different available statistics at different values of P ISR T . Operationally this can be accomplished as follows (following on an idea from ref. [56]). One can select a subset of events with (approximately) the same value of P ISR T , and plot them in the p 1T = (p T , p T ⊥ ) plane as in figure 2. The signal events will exhibit an overdensity along the singularity ellipse given by (3.24), as also illustrated in the right panel in figure 4. Then, by fitting to (3.24), one can find the distances from JHEP04(2020)027 the focus to the two vertices of the ellipse, i.e., the maximum value of p 1T in the direction along P ISR T and in the direction opposite to P ISR T : These two relations can be easily inverted and solved for M 1 and M 0 : (3.31) This demonstrates that the two measurements (3.28)-(3.29) are sufficient to determine the masses M 1 and M 0 [56]. The procedure can be repeated for different ranges of P ISR T , as long as there is sufficient statistics to reconstruct the singularity ellipse and from there extract the values of p max 1T (ϕ = 0) and p max 1T (ϕ = π).

The focus point method
The discussion in the previous section 3.2 revealed that for the event topology of figure 1(a), the location of the singularity is nicely exhibited as a two-dimensional surface in a threedimensional space of observables (p T , p T ⊥ , P ISR T ), as shown in the left panel in figure 4. But is it possible to simplify matters further, e.g., by projecting onto an observable space of even lower dimensionality, while retaining all singular features? As an extreme example, is it possible to define a single kinematic variable whose distribution will capture all of the singular behavior, over the whole surface parametrized by eq. (3.24)? In section 3.2.1 we saw that for the special case of P ISR T = 0 this is possible, and the relevant kinematic variable was the transverse mass m T (M 0 ) (regardless of the choice of test massM 0 , see figure 3), or alternatively, the magnitude p 1T of the transverse visible momentum. However, the subsequent discussion in section 3.2.2 makes it clear that if we wish to continue using m T (M 0 ) in the more general case of P ISR T = 0, we run into a problem -the parametrization of the singularity surface (3.24) involves the mass spectrum, and in particular the mass M 0 of the invisible particle. Therefore, as implied in the singularity condition (3.12), the transverse mass m T (M 0 ) will continue to be the relevant singularity coordinate, but only for the correct choice of the test massM 0 = M 0 , since only in that case the singularity surface (3.24) is a surface of constant m T . If we make a wrong choice forM 0 , which is different from the true mass M 0 , the singularity surface is not a surface of constant m T , and therefore, the singular behavior observed in the three-dimensional picture of the left panel in figure 4 will tend to be washed out in the one-dimensional m T (M 0 ) distribution. These observations are precisely the motivation behind the focus point method for mass measurement [13], which can also be applied to the event topology of figure 1(a) considered here, with the intent of measuring the two masses M 1 and M 0 . The method is illustrated in figures 5 and 6, where, in order to avoid overcrowding the plot, we use just a handful of events (in this case ten, chosen at random). The main idea is for each event to delineate the allowed region in the hypothesized mass parameter space (M 0 ,M 1 ), which would lead to viable solutions for the invisible momentum q, given the kinematic constraints (3.2). It is well known that for a given test massM 0 , the transverse mass m T (M 0 ) provides the lowest kinematically allowed value for the parent mass M 1 [31,36,53], therefore the boundary of the allowed region in our case will be given simply by the functioñ After superimposing these kinematic boundaries from many different events as in figures 5 and 6, the true values of M 0 and M 1 are revealed by the location of the focus point of the kinematic boundaries [13]. This is indeed what is seen in figures 5 and 6 -even with the low statistics of just 10 events, a focus point emerges near the red "+" symbol marking the true mass point (M 0 , M 1 ) = (700 GeV, 800 GeV). In an actual experimental analysis, the available statistics is expected to be much larger, perhaps as much as several orders of magnitude, thus there should be no problem observing this focus point. Note that when JHEP04(2020)027 plotted in the (M 0 ,M 1 ) plane as in figure 5, the lines tend to be parallel to each other, and their crossing is difficult to trace with the naked eye without zooming in on the relevant region near the red cross. This is why in figure 6 we have replotted the same data, only now replacingM 1 on the y-axis with the more relevant combination of masses (M 2 1 −M 2 0 )/M 1 which enters the analytical formulas (see, e.g., eqs. (3.18) and (3.25)-(3.27)). As seen in figure 6, this has the benefit of removing the dependence on the overall scaleM 0 , which allows us to better concentrate on the relative differences exhibited by lines corresponding to different events.
As observed in figures 5 and 6 (and supported by the previous discussion in section 3. be even more beneficial to showcase our method, we limit ourselves here to these more typical values expected to come from initial state radiation (which has a steeply falling P T spectrum) or from decays upstream, where the typical P T is governed by the mass splitting of the new particles, which is likely to be around the M 1 mass scale.
Let us first discuss the upper left panels corresponding to the case of P ISR T = 0 which was the subject of section 3.2.1. These plots demonstrate that in the absence of any P ISR T , different values ofM 0 are indistinguishable. In the plots, this is indicated by the fact that the lines stay roughly parallel to each other, and there are no line crossings at all. Notice that quite a few events in the figure are close to being extreme, i.e., their lines pass very close to the true mass point marked with the red cross. This is simply a consequence of the singularity condition (3.12) -since the ten events entering the plot were picked at random, it is more likely that they belong to the region where the event density exhibits a (singular) peak. Correspondingly, their lines appear to be "bunched up", in the sense that their m T values forM 0 = M 0 tend to be very similar, being so close to the true mass M 1 of the parent particle. Then, as we vary the test massM 0 away from the true value M 0 , the lines in the upper left panels continue to stay bunched up, confirming the presence of a singularity in the m T (M 0 ) distribution for all other values ofM 0 as well (recall figure 3); the exact location of the singularity as a function ofM 0 is given by eq. (3.19).
Fortunately, the situation changes completely in the presence of non-zero P ISR T , as illustrated in the remaining five panels in each of figures 5 and 6. We use the same 10 events as before, only now they have been boosted accordingly in order to generate the desired P ISR T . For illustration, let us only focus on the lower right panels with the largest P ISR T = 800 GeV, where the effects are easiest to see. We observe that the lines for the nearextreme events are still bunched up in the vicinity of the red + symbol, but significantly diverge away from it in the region of either very low or very high values ofM 0 (this is especially easy to see with the parametrization used in figure 6). In other words, the crossing pattern of the lines has formed a focus point in the vicinity of the true mass point, and therefore, finding this focus point is a way to find the true masses [13]. Note that Both of those mass measurement methods are illustrated in figures 7 and 8, which are the analogues of figures 5 and 6, only this time we are using the full data sample and, following [13], we represent the density of curves as a heatmap where the color corresponds to the fraction of events 12 whose lines pass through a square bin of width 2 GeV.  test massM 0 , but the true value of M 0 remains unknown. This is why, such events with low P ISR T should not be used in the analysis for measuring M 0 with the focus point method. Before concluding this section, let us discuss the connection between the focus point method presented here and the kink method for mass measurements [52][53][54][55][56]. The two methods are closely related -in fact, figures 5-8 also nicely illustrate the kink method itself, where one instead tries to measure the maximal possible value of m T as a function of the test massM 0 . In other words, the kink method is essentially targeting the upper boundaries of the colored regions in figure 7, which for P ISR T = 0 exhibit a kink 13 at M 0 = M 0 . In contrast, the focus point method is targeting the highest line density bin on the plot, and thus is designed to take full advantage of the available statistics -note that the boundaries of the colored regions in figures 7 and 8 are defined by just a handful of events and the extraction of the boundary (from kinematic endpoints or otherwise) is statistics-limited. On the other hand, as illustrated in figures 5 and 6, events which do not contribute to the boundary may still pass close to the focus point and thus usefully contribute to the focus point method. Of course, the kink occurs precisely at the location of the highest line density bin, 14 thus the two methods in principle give the same results, as confirmed in figures 5-8.

Single decay chain, two successive two-body decays
In this section we shall discuss the event topology of figure 1(b) which is a cascade decay involving three new particles with masses M 0 , M 1 and M 2 . Correspondingly, there are three on-shell conditions which can be rewritten as

2c)
At this point, we have N C = 3 constraints and N q = 4 unknowns which are the components of q µ . Thus in order to obtain the necessary match (2.8), we can do one of two thingseither add an additional constraint, e.g., in the form of a measurement of one of the / P T components (but not the other), or reduce the number of unknowns by going to 2 + 1 dimensions, where N q = 3. Both of these options are of mostly academic interest, so for concreteness we choose the latter, which has the added benefit of somewhat simpler math 13 The kink is easier to see with the alternative parametrization of figure 8. 14 This connection was present, but overlooked in the existing literature, e.g., note the resemblance of figures 6(c) and 6(d) in ref. [54] to our figures 7(a) and 7(b), respectively.

JHEP04(2020)027
(involving 3 × 3 instead of 4 × 4 matrices). Consequently, for the remainder of this section, we shall be working in 2 + 1 dimensions, where the momenta have only transverse and no longitudinal components, i.e., p i = (e i , p ix , p iy ) and q i = (ε, q x , q y ). The next step is to construct the Jacobian matrix (2.12) for the set of constraints (4.2): The singularity condition (2.13) becomes 15 ε q x q y e 1 p 1x p 1y e 2 p 2x p 2y = 0. (4.5) Now, in order to rewrite the singularity variable appearing on the left-hand side in terms of observable momenta only, we just need to eliminate the invisible momentum components, for example, using eqs. (4.2). However, a more straightforward approach is to note that and combine the last two equations as By rewriting the singularity condition in this form, we have not only eliminated any reference to the invisible momentum components, but we also verified explicitly that the event-wise kinematic information only enters through the dot product p 1 · p 2 , or equivalently, through the invariant mass m 12 of the two visible particles since the latter is related to p 1 · p 2 as m 2 Eq. (4.7) suggests that instead of taking the whole determinant as the singularity coordinate for this event topology, a much simpler choice would be either p 1 · p 2 or m 12 , both of 15 Note that it can be equivalently written in a more compact form as

JHEP04(2020)027
which have the important advantage of avoiding the necessity of an ansatz for the a priori unknown masses M 0 , M 1 and M 2 . For definiteness here we shall choose to work with p 1 ·p 2 , but m 12 will be equally good. 16 The locations of the singularities in the p 1 · p 2 distribution can be found from eq. (4.7), which leads to a quadratic equation and correspondingly two solutions. In order to simplify the formulas, let us again focus only on the case of massless visible particles, i.e., m 1 = m 2 = 0, when (4.7) becomes leading to the quadratic equation also happen to be the two kinematic endpoints of the p 1 · p 2 distribution. The actual shape of the distribution according to pure phase space is given by (4.12) The distribution of the singularity variable p 1 ·p 2 is illustrated in figure 9. In each panel, we plot the theoretical prediction (4.12) superimposed on the result from our numerical simulations. The spectrum was chosen as in table 1: M 0 = 700 GeV, M 1 = 800 GeV and M 2 = 1000 GeV, which according to (4.11b) gives (p 1 · p 2 ) max = 42, 187.5 GeV 2 . We see that, as expected, the distribution develops a sharp singularity at each end. These singularities are quite striking when viewed on a linear scale (as in the left panel). In order to better see the shape of the distribution in the intermediate p 1 · p 2 range, in the right panel of figure 9 we replotted the same data using a log scale for the y-axis.
Having derived the singularity coordinate for this case as either p 1 · p 2 or m 12 and identified the locations (4.11) of the two singularities, we have accomplished two of the three stated goals at the end of section 2. The last goal, showcasing the focus point method of ref. [13] for mass measurements, is not applicable in this case, since the singularity variable 16 The reader should keep in mind that we are working in 2 + 1 dimensions here, and it is only in that case that m12 is a singularity coordinate and its distribution has singularities. In 3 + 1 dimensions, the distribution of m12 has no singularities and is simply proportional to m12 (in other words, the distribution of m 2 12 is flat).  Probability density (log scale) Expectation Simulation Figure 9. Distribution of the singularity variable p 1 · p 2 on a linear scale (left panel) and on a log scale (right panel). With the mass spectrum from table 1, eq. (4.11) predicts the locations of the two singularities to be at (p 1 · p 2 ) min = 0 and (p 1 · p 2 ) max = 42, 187.5 GeV 2 . The black solid line shows the theoretical prediction (4.12) for the shape of the distribution while the histogram shows the result from the numerical simulation.
is constructed out of visible momenta only, with no reference to any hypothesized mass parameters. The situation is analogous to the one already encountered in section 3.2.1 -there we saw that when P ISR T = 0, the singularity variable for the event topology in figure 1(a) can be taken to be simply the p T of the visible particle, and a singularity occurs for any choice of test mass as shown in the upper left panels of figures 5-8. The best one can do in our case here, therefore, is to obtain one constraint on the three masses M 0 , M 1 and M 2 from the measurement of the upper kinematic endpoint (4.11b).

Single decay chain, three successive two-body decays
In this section we shall work out the classic SUSY decay chain of three successive two-body decays depicted in figure 1(c). This cascade involves four new particles with masses M 0 , M 1 , M 2 and M 3 , which in SUSY are typically identified with the lightest neutralino, a charged slepton, the second-to-lightest neutralino, and a squark, respectively. The same decay chain also pops up in other phenomenological models of new physics such as Universal Extra Dimensions with KK-parity [71][72][73], Little Higgs with T -parity [74][75][76][77], etc. In any given such new physics scenario, there is a definite assignment for the spins of the new particles, however in order to stay as general as possible, we shall ignore any spin correlations (which are typically rather small anyway [78][79][80]) and let the particles decay according to pure phase space. Unlike the toy example of the previous section, here and for the rest of the paper we shall work in 3 + 1 dimensions, thus a single invisible particle with 4-momentum q contributes N q = 4 unknown degrees of freedom.
We begin by listing the four on-shell conditions for the event topology of figure 1(c):

JHEP04(2020)027
which can be rewritten in analogy to (4.2) as Since we already have N C = 4 constraints for N q = 4 unknowns, the condition (2.8) is already met and we can proceed with the derivation of the Jacobian matrix (2.12) and the singularity condition (2.13) 17 In order to obtain a singularity variable in terms of the visible momenta only, we need to eliminate the invisible momentum components with the help of (5.2). However, the same task is more easily accomplished with the determinant trick used in the previous section: after noting that we can combine the last two equations into Det D = 0 ⇐⇒ Det DηD T = 256 17 Note the alternative compact notation in analogy to (4.4):

JHEP04(2020)027
With the use of (5.2), the dot products of momenta appearing in (5.6) can now easily be traded for the relevant masses. Once again, the result simplifies significantly in the case of massless visible particles, i.e., p 2 0 = 0, (5.8) where m ij are the pair-wise invariant masses of the three (massless) visible particles: Up to the numerical prefactor of 256, the left-hand side of (5.8) is precisely the ∆ 4 variable introduced in [18,34]. We have thus rederived from first principles 18 the well-known fact that ∆ 4 is the relevant singularity variable for the event topology of figure 1(c), and that the locations of singularities are those where ∆ 4 = 0. Eq. (5.8) also confirms that the relevant observable phase space is only three-dimensional, 19 and can be conveniently parametrized with the pair-wise invariant masses as Given the ubiquity of the decay chain of figure 1(c) in SUSY and elsewhere, it is not surprising that the properties of the allowed region within this invariant mass phase space have been extensively studied in the literature. Therefore, rather than reproducing previously published work, here we shall only state the results most relevant to the current discussion, and for further details we refer the reader to the corresponding literature.
• The shape of the allowed region in phase space. The allowed region in the visible invariant mass space (5.10) is compact, and is bounded by the (closed) two-dimensional surface defined by eq. (5.8). Three-dimensional plots of the allowed region can be found in figure 1 of [16], figure 9 of [19] and on pages 568-572 of the TASI lectures [17].
• Density enhancement on the boundary. The phase space singularities occur on the two-dimensional boundary of the allowed region, i.e., the condition ∆ 4 = 0 (which is equivalent to (5.8)) defines both the boundary of the allowed region as well as the singularity locations [18]. Since the allowed region is three-dimensional, it is difficult to visualize this enhancement unless one looks at two-dimensional slices through the allowed region -such plots can be found in figure 8 of [20] and figure 12 of [21]. 18 Alternative derivations leading to ∆4 = 0 as the defining condition for the boundary of the allowed phase space can be found in [16,17]. The fact that the event number density is in addition singular at that boundary was later emphasized in [18]. 19 The visible momenta p1, p2 and p3 parametrize a 9-dimensional phase space, but the constraints (5.2) are invariant under the 6-parameter Lorentz group, leaving only 3 relevant visible degrees of freedom.

JHEP04(2020)027
• The shape of the one-dimensional distribution of the singularity variable. The differential distribution of ∆ 4 is known analytically. In terms of the unit-normalized variable q ≡ ∆ 4 /∆ max 4 it is given by [20] dN dq = arcsin( The sharp peak at ∆ 4 = 0 can be used for discovering such new physics signal over the smooth SM background, as discussed in [25]. • Mass measurements and the focus point method. Since computing the singularity variable ∆ 4 requires an ansatz for the mass spectrum, the focus point method for mass measurements [13] is in principle applicable, and one would be looking for a peak in the 4-dimensional parameter space of {M 0 , M 1 , M 2 , M 3 }. However, the relevant observable parameter space (5.10), being only three-dimensional, is already simple enough so that in practice it may be easier to just perform a four-parameter fit to the boundary of the allowed region, as demonstrated in [21].
6 Two decay chains, each with one two-body decay In this section we shall simultaneously address the two event topologies shown in figures 1(d) and 1(e). The latter is known as the "antler" topology [81] and has been previously discussed in the context of both hadron colliders [12,[81][82][83][84][85][86][87][88][89][90] and lepton colliders [91,92]. At the same time, the diagram of figure 1(d) is extremely common, and may represent many processes, including but not limited to W -pair production in the SM, squark or slepton production in SUSY, etc. This diagram also has been extensively discussed in the literature -at lepton colliders [93][94][95], and especially at hadron colliders, where it offers a formidable challenge, despite its apparent simplicity (for reviews, see [27,28]). In fact, it was precisely the diagram of figure 1(d) which initially motivated a large number of now popular kinematic variables and techniques, including the Cambridge M T 2 variable [96,97] and its variants [70,[98][99][100], the contransverse mass M CT [101], the M CT 2 variable [102,103], the MAOS method [104,105], and many others. Given the similarities in the diagrams of figures 1(d) and 1(e), in this section we shall discuss them in one go by assuming the 4-momentum vector of the initial state P = (P 0 , P) to be fixed. This is certainly true at lepton colliders, where the kinematics of the initial state is completely known: P = (E CM , 0, 0, 0), where E CM is the beam CM energy, which in the case of the antler diagram of figure 1(e) can be tuned to be equal to the mass M 2 of the intermediate resonance, so that P becomes P = (M 2 , 0, 0, 0). At hadron colliders, we will primarily focus on the antler diagram of figure 1(e), for which the 4-momentum of the initial state can be written as P = ( M 2 2 + P 2 z , 0, 0, P z ), since in general the resonance of mass M 2 will be produced with some non-zero longitudinal momentum P z , whose size will be governed by the value of M 2 and the parton distribution functions (pdfs) of the initial state partons. For completeness, we shall initially retain P z in our formulas, but in the end, following [12], we shall take the "gluon collider" approximation P z ≈ 0, which can be JHEP04(2020)027 justified in cases where the pdfs of the initial state partons are the same (or similar) and are fast-falling functions, e.g., as in the SM process H → W + W − [12].

Derivation of a singularity coordinate
The event topologies of figures 1(d) and 1(e) have the following constraints in common: where in the last two equations we have assumed that there is no P ISR T accompanying our event topology, hence / P T = − p 1T − P 1T . Since we are treating the initial state kinematics as fixed, we can add two more relations representing energy conservation and longitudinal momentum conservation, respectively with some fixed P 0 and P z as discussed above. Eqs. (6.1e)-(6.2b) can now be used to eliminate the 4-momentum Q as Substituting this into eqs. (6.1b) and (6.1d), and again limiting ourselves for simplicity to the case of massless visible particles, p 2 1 = P 2 1 = 0, we can rewrite the four remaining constraints (6.1a)-(6.1d) as From here we can compute the Jacobian matrix (2.12) Using (6.4) to eliminate q, after a couple of row and column manipulations, this simplifies to where the left-hand side is precisely the desired singularity variable ∆ antler for the case of figures 1(d) and 1(e). Eq. (6.8) indicates that the relevant phase space of kinematic observables is three-dimensional 20 and can be parametrized, e.g., as {P · p 1 , P · P 1 , p 1 · P 1 } . (6.9) Just like the singularity variable ∆ 4 from section 5, the singularity variable ∆ antler depends not only on the phase space (6.9), but also requires an ansatzM 1 andM 0 for the two mass parameters, i.e., ∆ antler (M 1 ,M 0 ). If the ansatz is correct (M 1 = M 1 andM 0 = M 0 ), then the singularity condition (6.8) guarantees that the distribution of ∆ antler exhibits a singularity peak at ∆ antler (M 1 , M 0 ) = 0. (6.10) This is demonstrated in the left panel of figure 10, which shows the one-dimensional unit-normalized distribution of ∆ antler (M 1 , M 0 ) at a lepton collider with E CM (= M 2 ) = 3000 GeV, and with the mass spectrum from table 1, M 1 = 800 GeV and M 0 = 700 GeV. As expected, there is a sharp peak at ∆ antler = 0. We also note that with our conventions the values of ∆ antler are negative -this can be traced back to eq. (6.6) and the fact that with the correct choice of mass parameters D is guaranteed to be real and thus Det DD T > 0, while Det η = −1.
However, the mass spectrum may not always be known a priori, thus one would like to reduce the number of mass ansatze as much as possible. To this end, we can follow the idea of the transverse mass m T from section 3, which takes onlyM 0 as an input and then uses the corresponding singularity condition (3.7) to define the singularity variable. In our case, we can use (6.10) to define implicitly an alternative singularity variable M antler (M 0 ) as the solution to the equation to two possible solutions, both of which are entered in the plot. The larger of the two solutions, which we label M max antler , comprises the orange histogram, while the smaller one, labelled M min antler , makes up the blue histogram, and the two histograms are individually normalized to 1. As expected, the right panel of figure 10 exhibits the presence of a sharp peak at the correct mass of the parent particle, M 1 = 800 GeV. It is interesting to note that both solutions M min antler and M max antler are contributing to the singularity -from below and from above, respectively. Another noteworthy feature of the plot is that the orange and blue histograms do not overlap at all -in fact, they meet at the true value of the mass M 1 = 800 GeV (which is also the location of the singularity peak). This can be seen even more clearly in figure 11, which shows the two-dimensional distribution of events as a heatmap in the plane of (M max antler , M min antler ).

JHEP04(2020)027
The heatmap in figure 11 reveals an overdensity of events at both M min antler = M 1 (the horizontal blue-shaded band) and M max antler = M 1 (the vertical blue-shaded band). This confirms that both solutions for M antler play a role in forming the singularity peak observed in the right panel of figure 10. Note also the absence of any events with M min antler > M 1 and M max antler < M 1 , which implies that the following hierarchy is always true: In other words, the true parent mass M 1 is always located between the two found solutions for M antler , and furthermore, M 1 is the upper kinematic endpoint of M min antler (the blue histogram in the right panel of figure 10) and at the same time it is also the lower kinematic endpoint of M max antler (the orange histogram in the right panel of figure 10). This can be understood as follows. The two solutions for M antler obtained from the singularity condition (6.11) represent the two values for the trial massM 1 where the number of solutions for the invisible momentum q changes [13], in this case between 0 and 2. Since the true mass M 1 will always give valid solutions for q, it belongs to the allowed interval forM 1 with 2 solutions, which is sandwiched between the two disallowed regions with 0 solutions.
In summary, the discussion in this subsection (and figure 10 in particular) shows that both ∆ antler and M antler are valid singularity variables, albeit the latter has the added advantages of having a clear physical meaning and being of mass dimension 1 only.

The phase space geometry of the singularity condition
Having derived the singularity variables for the event topologies of figures 1(d) and 1(e), we can now discuss the geometry of the singularity surface in the relevant observable phase space. The latter can be parametrized as in (6.9), but the equation of the singularity surface (6.8) can be written more compactly if we use an alternative set of observables X = 2(P · p 1 + P · P 1 ), (6.13a) Y = 2(P · p 1 − P · P 1 ), (6.13b) Z = 4(P · p 1 )(P · P 1 ) P 2 − 2(p 1 · P 1 ) (6.13c) which reduces (6.8) to the constraint This equation describes a closed surface whose cross-sections at fixed Z are ellipses in the (X , Y) plane. As discussed earlier, for illustration purposes, we shall now fix the momentum of the initial state as P = (M 2 , 0, 0, 0), which can be viewed as a lepton collider running at a CM energy E neglects the longitudinal momentum of the heavy s-channel resonance [12]. In that limit, the variables (6.13) become In the last line we recognize the quantity M C introduced in [101], which is invariant under contra-linear (back-to-back) boosts. We can use this boost invariance to bring the two intermediate particles of mass M 1 to their corresponding rest frames (along with their decay products), which then allows us to express the quantity M C as (1 + cos θ * ), (6.16) where θ * is the angle between p 1 and P 1 after the respective boosts. With the help of (6.15) and (6.16), the equation of the singularity surface (6.14) becomes The singularity surface defined by this equation is pictorially illustrated in figures 12 and 13. Figure 12 is analogous to the plot in the left panel of figure 4 from section 3.2.2, where the space of relevant observables was three-dimensional as well. In figure 12, we show two different views of the singularity surface when, as suggested by eq. (6.15), the observable phase space is parametrized as (e 1 , E 1 , M 2 C ). 21 Figure 13, on the other hand, shows a series of plots in analogy to figure 2. Each individual panel depicts the allowed range of values for the energies e 1 and E 1 of the two visible particles for a given fixed value of M C , or equivalently, for a fixed value of cos θ * , 21 The plots in figure 12 can be contrasted to the plots in figure 3 of ref. [12], which were done with the parametrization p1T , P1T , p 1T · P 1T since the two are related by eq. (6.16). For this figure we prefer to work with cos θ * (and equally spaced fixed values for it) since the distribution in cos θ * is uniform, thus each panel in figure 13 has the same total number of events. The event number density in the (e 1 , E 1 ) plane is indicated by the color bar. We see that, as expected, the events tend to cluster on the singularity boundary, whose shape is an ellipse of varying eccentricity depending on the value of cos θ * . The distortion in the shape of the elliptical boundary can be easily tracked and understood with the help of eq. (6.17). Consider, for example, the case of cos θ * = +1 shown in the upper left panel of figure 13. This implies that θ * = 0 and the first term in the left hand side of (6.17) dominates, which in turn implies the linear relation = 351.6 GeV seen in the plot. As the value of cos θ * decreases, the elliptical boundary becomes less eccentric, and for cos θ * = 0, i.e., θ * = π 4 , it eventually becomes a circle, as seen in the middle plot of the middle row. As the value of cos θ * decreases further, the ellipse begins to stretch along the orthogonal e 1 = E 1 direction, and for cos θ * = −1 it simply becomes the line e 1 = E 1 . We use 20 randomly chosen events in the left panel, in which the red "+" symbol marks the location of the true masses, and 10,000 events to produce the heatmap in the right panel, where the color indicates the fraction of events whose solvability boundaries pass through a given 2 × 2 GeV bin.

The focus point method
Before concluding this section, we shall demonstrate that the focus point method of ref. [13] applies to the event topologies of figures 1(d) and 1(e) as well. The point is that the singularity variable ∆ antler derived in section 6.1 requires a mass ansatzM 1 andM 0 as input. The need for such mass ansatze was once considered undesirable, but, as more recent studies have shown, it is precisely the dependence on the test masses that opens the door to new methods for extracting useful information about the mass spectrum, case in point being the kink method for measuring M 0 [52][53][54][55][56].
In our case, the ansatz forM 1 andM 0 is needed to provide the needed number of kinematic constraints, which would allow us to solve for the invisible momenta. However, not all choices ofM 1 andM 0 will lead to real solutions. Each event will thus delineate a viable region in the (M 0 ,M 1 ) mass parameter space. The idea of the focus point method is to superimpose the boundaries of the allowed regions selected by different events, as illustrated in the left panel of figure 14. The plot shows the solvability boundaries defined by ∆ antler (M 1 ,M 0 ) = 0 (6.18) for 20 randomly chosen events, where for better visualization, we have rescaled the y-axis as was previously done in figures 6 and 8. We see that even with just a handful of events, the solvability boundary curves tend to focus near the true mass point, marked with the red "+" symbol. With a lot more statistics, we obtain the heatmap shown in the right panel of figure 14 (contrast to the analogous heatmaps seen in figure 8). The heatmap clearly identifies the singularity peak which is situated at the true values of the masses, thus establishing the viability of the method.
7 Two decay chains, each with two successive two-body decays For completeness, in this section we shall review the final event topology from figure 1, namely, the dilepton tt event topology in figure 1(f), which was also the one used in ref. [13] JHEP04(2020)027 to introduce and illustrate the idea of the focus point method for mass measurements. Correspondingly, we shall not repeat the analysis of ref. [13] here, and simply refer the readers interested in mass measurements aspects to that paper. Here we shall focus more narrowly on the derivation of a singularity coordinate for that case, following the general method outlined in section 2.2 and illustrated with the examples from the previous sections. Let us begin by listing the kinematic constraints for the event topology from figure 1(f): 1b) (q + p 1 ) 2 = M 2 1 , (7.1c) (Q + P 1 ) 2 = M 2 1 , (7.1d) (q + p 1 + p 2 ) 2 = M 2 2 , (7.1e) (Q + P 1 + P 2 ) 2 = M 2 2 , (7.1f) q x + Q x = / P x , (7.1g) q y + Q y = / P y , (7.1h) which can be rewritten as The Jacobian matrix (2.12) for the system (7.2) is where the dependence on the mass ansatz enters through the solutions for the invisible momentaq andQ from (7.3). We can now take the singularity variable for the event topology of figure 1(f) to be ∆ tt (M 2 ,M 1 ,M 0 ) ≡ Det D tt (M 2 ,M 1 ,M 0 ). (7.5) As indicated, its computation again requires an ansatz for the masses, but this is precisely the property which makes it relevant for mass measurements, since the singularity condition ∆ tt (M 2 , M 1 , M 0 ) = 0 (7.6) holds only if we use the true mass spectrum [13]. The distribution of ∆ tt (M 2 , M 1 , M 0 ) is shown in figure 15. One subtlety of the ∆ tt computation is that there can be multiple (up to four) solutions (7.3) for the invisible momenta. In making figure 15, we made sure that each event contributes equally to the plot, by entering the result for each solution with a weight 1/N s , where N s is the total number of solutions found in that event. 22

JHEP04(2020)027
Having derived the singularity variable for the event topology of figure 1(f), our next task would have been to illustrate the singularity surface in the relevant observable phase space, similarly to figure 4 (left panel) and figure 12. Unfortunately, the observable phase space in this case is nine-dimensional, 23 and we shall not attempt to visualize it here. The third and final task, the demonstration of the focus point method, was already accomplished in ref. [13].

Conclusions and outlook
In this paper we outlined the general prescription for deriving a singularity variable for a given event topology with missing energy, i.e., where some of the final state particles are invisible in the detector. We then illustrated the procedure with several common event topologies shown in figure 1. We started with the case of a single two-body decay in section 3 and re-derived the well known result that the distribution of the transverse mass m T has a Jacobian peak. In the subsequent sections, we demonstrated that similar Jacobian peak features are present in the distributions of the relevant kinematic variables for the remaining five event topologies. We also identified, parametrized and studied the shapes of the singularity surfaces in the appropriate visible phase spaces. In some special circumstances (see sections 3.2.1 and 4) the singularity variable can be computed directly in terms of the available kinematic information, without the need for any additional hypothesized inputs. However, more often than not, the singularity variable depends on the masses of the intermediate resonances and/or the masses of the invisible final state particles, and thus its computation requires a mass ansatz. If the event topology is applied to a SM process, the input masses will be known, but when applied to a BSM process, the mass ansatz is a priori unknown. However, this can be used to our benefit -it is precisely this dependence on the mass ansatz that makes the focus point method for mass measurements possible, as explicitly demonstrated in section 3.3 and 6.3. 24 The main advantage of the focus point method is that it maximally utilizes the singularity structures in phase space. As a consequence, the true masses are identified as a kinematic peak instead of a kinematic endpoint -endpoints are more difficult to observe experimentally, once we include the finite widths and detector resolution effects [109]. It is also worthwhile to contrast the focus point method to the polynomial method [110][111][112][113]. While both methods use the same type of kinematic constraints, the latter requires a larger set of constraints in order to avoid the need for a mass ansatz.
The ideas presented in this paper may find immediate application in a large number of LHC analyses targeting the event topologies of figure 1.
• Standard Model measurements. Precision studies of SM processes typically require the identification of a specific event topology, e.g., figure 1(a) for W production, 23 The momenta of the four visible particles are parametrized by 4 × 3 = 12 degrees of freedom, three of which can be removed by an azimuthal rotation and separate z-boosts for each of the two decay chains in figure 1(f). 24 For an application of the focus point method to the event topologies of figures 1(c) and 1(f), see refs. [21] and [13], respectively. JHEP04(2020)027 figure 1(b) for top quark decay, figure 1(d) for W pair-production, and figure 1(f) for dilepton tt production. The corresponding singularity variables are ideal for tagging such events, and may also be used as input features for event selectors based on machine learning.
• New physics searches. The Jacobian peaks in the distributions of the singularity variables can be used to discover new physics processes over the SM background [25]. Likewise, the peak structures in the heatmaps constructed in the focus point method can also be used for discovery of new physics, with the added advantage that background processes will not develop fake peaks 25 in the signal regions.
In this paper we mostly focused on the case (2.8) when the number of unknowns N q matches the number of kinematic constraints N C . However, the under-constrained case N C < N q is also worth investigating, e.g., following the analysis of ref. [12]. All of these topics are being pursued in a future publication [42].