Detecting kinematic boundary surfaces in phase space: particle mass measurements in SUSY-like events

We critically examine the classic endpoint method for particle mass determination, focusing on difficult corners of parameter space, where some of the measurements are not independent, while others are adversely affected by the experimental resolution. In such scenarios, mass differences can be measured relatively well, but the overall mass scale remains poorly constrained. Using the example of the standard SUSY decay chain $\tilde q\to \tilde\chi^0_2\to \tilde \ell \to \tilde \chi^0_1$, we demonstrate that sensitivity to the remaining mass scale parameter can be recovered by measuring the two-dimensional kinematical boundary in the relevant three-dimensional phase space of invariant masses squared. We develop an algorithm for detecting this boundary, which uses the geometric properties of the Voronoi tessellation of the data, and in particular, the relative standard deviation (RSD) of the volumes of the neighbors for each Voronoi cell in the tessellation. We propose a new observable, $\bar\Sigma$, which is the average RSD per unit area, calculated over the hypothesized boundary. We show that the location of the $\bar\Sigma$ maximum correlates very well with the true values of the new particle masses. Our approach represents the natural extension of the one-dimensional kinematic endpoint method to the relevant three dimensions of invariant mass phase space.


Introduction
The dark matter problem is currently our best experimental evidence for the existence of new particles and interactions beyond the Standard Model (BSM). A great number of ongoing experiments are trying to discover dark matter particles through direct [1] or indirect detection [2]. In principle, dark matter particles could also be produced in high energy collisions at the Large Hadron Collider (LHC) at CERN, providing a complementary discovery probe in a controlled experimental environment [3].
Since the dark matter particles must be stable on cosmological timescales, in many popular BSM models they carry some conserved quantum number. The simplest choice is a Z 2 parity, which is known as R-parity in models with low energy supersymmetry (SUSY) [4], Kaluza-Klein (KK) parity in models with universal extra dimensions (UED) [5], Tparity in Little Higgs models [6], etc. As a result, the dark matter particles are necessarily produced in pairs: either directly, or in the cascade decays of other, heavier BSM particles [7]. The prototypical such cascade decay is shown in Fig. 1, in which a new particle D undergoes a series of two-body decays, terminating in the dark matter candidate A, which Figure 1: The typical cascade decay chain under consideration in this paper. Here D, C, B and A are new BSM particles, while the corresponding SM decay products are: a QCD jet j, a "near" lepton ± n and a "far" lepton ∓ f . This chain is quite common in SUSY, with the identification D =q, C =χ 0 2 , B =˜ and A =χ 0 1 , whereq is a squark,˜ is a slepton, andχ 0 1 (χ 0 2 ) is the first (second) lightest neutralino. In what follows we shall quote our results in terms of the D mass m D and the three dimensionless squared mass ratios R CD , R BC and R AB defined in eq. (1.6).

Introduction
SUSY is a primary target of the LHC searches for new physics beyond the Standard Model (BSM). In SUSY models with conserved R-parity the superpartners are produced in pairs and each one decays through a cascade decay chain down to the lightest superpartner (LSP). If the LSP is the lightest neutralinoχ 0 1 , it escapes detection, making it rather difficult to reconstruct directly the preceding superpartners and thus measure their masses and spins. In recognition of this fact, in recent years there has been an increased interest in developing new techniques for mass  and spin  measurements in such SUSY-like missing energy events.
Roughly speaking, there are three basic types of mass determination methods in SUSY 1 . In this paper we concentrate on the classic method of kinematical endpoints [1]. Following the previous SUSY studies, for illustration of our results we shall use the generic decay chain D → jC → j ± n B → j ± n ∓ f A shown in Fig. 1. Here D, C, B and A are new BSM particles with masses m D , m C , m B and m A . Their corresponding SM decay products are: a QCD jet j, a "near" lepton ± n and a "far" lepton ∓ f . This decay chain is quite common in SUSY, with the identification D =q, C =χ 0 2 , B =˜ and A =χ 0 1 , whereq is a squark,˜ is a slepton, andχ 0 1 (χ 0 2 ) is the first (second) lightest neutralino. However, our analysis is not limited to SUSY only, since the chain in Fig. 1 also appears in other BSM scenarios, e.g. Universal Extra Dimensions [77]. For concreteness, we shall assume that all three decays exhibited in Fig. 1 are two-body, i.e. we shall consider the mass hierarchy (1.1) 1 For a recent study representative of each method, see Refs. [43,47,49].
-2 - is a difficult problem, which has been attracting a lot of attention over the last 20 years (for a review, see [8]). The main challenge stems from the fact that the momentum of particle A is not measured, so that the standard technique of directly reconstructing the new particles as invariant mass resonances does not apply. Instead, one has to somehow infer the new masses (1.1) from the measured kinematic distributions of the visible SM decay products.
In the decay chain of Fig. 1, the SM decay products are taken to be a quark jet j and two leptons, labelled "near" n and "far" f . This choice is motivated by the following arguments: • At a hadron collider like the LHC, strong production dominates, thus particle D is very likely to be colored. At the same time, the dark matter candidate A is neutral, therefore the color must be shed somewhere along the decay chain in the form of a QCD jet. Here we assume that this "color-shedding" occurs in the D → C transition 1 , since one expects the strong decays of particle D to be the dominant ones.
• The presence of leptons among the SM decay products in Fig. 1 is theoretically not guaranteed, but is nevertheless experimentally motivated. First, leptonic signatures have significantly lower SM backgrounds and thus represent clean discovery channels. 1 We note that in principle one can test this assumption experimentally, e.g. by constructing suitably defined on-shell constrained M2 variables corresponding to the competing event topologies [9], or by studying the shapes and the correlations for the invariant mass variables considered below [10]. Such an exercise is useful, but beyond the scope of this paper.
Second, the momentum of a lepton is measured much better than that of a jet, therefore the masses (1.1) will be measured with a better precision in a leptonic channel (as opposed to a purely jetty channel). Finally, if the SM decay products in Fig. 1 were all jets, in light of the arising combinatorial problem [11][12][13][14][15][16], we would have to resort to sorted invariant mass variables [17,18], whose kinematic endpoints are less pronounced and thus more difficult to measure over the SM backgrounds.
• From a historical perspective, the best motivation for considering the decay chain of Fig. 1 is that it is ubiquitous in SUSY, where D represents a squarkq, C is a heavier neutralinoχ 0 2 , B is a charged slepton˜ and A is the lightest neutralinoχ 0 1 , which escapes the detector and leads to missing transverse energy / E T . In the two most popular frameworks of SUSY breaking, gravity-mediated and gauge-mediated, the combination of (a) specific high scale boundary conditions, and (b) renormalization group evolution of the soft SUSY parameters down to the weak scale, leads to just the right mass hierarchy for the decay chain of Fig. 1 to occur. In the late 1990's and early 2000's, this prompted a flurry of activity on the topic of mass determination in such "SUSY-like" missing energy events. Soon afterwards, it was also realized that the decay chain of Fig. 1 is not exclusive to supersymmetry, but the same final state signature also appears in other models, e.g. minimal UED [19] and littlest Higgs [20].
To date, a large variety of mass measurement techniques for SUSY-like events have been developed. Roughly speaking, they all can be divided into two categories.
• Exclusive methods. In this case, one takes advantage of the presence of two decay chains in the event (they are often assumed identical) and the available / E T measurement. Several approaches are then possible. For example, in the so-called "polynomial methods" one attempts to solve explicitly for the momenta of the invisible particles in a given event, possibly using additional information from prior measurements of kinematic endpoints [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35]. 2 Alternatively, utilizing information from both branches, one could introduce suitable transverse 3 variables whose distributions exhibit an upper kinematic endpoint indicative of the parent particle mass [48][49][50][51][52][53][54][55][56][57][58][59][60][61][62][63]. In the latter case, one still retains a residual dependence on the unknown dark matter particle mass m A , which must be fixed by some other means, e.g. via the kink method [64][65][66][67][68][69][70][71] or by performing a sufficient number of independent measurements [53,57]. While they could be potentially quite sensitive, these exclusive methods are also less robust, since they rely on the correct identification of all objects in the event, and are thus prone to combinatorial ambiguities, the effects from / E T resolution, initial and final state radiation, underlying event and pileup, etc.
2 For long enough decay chains, the polynomial methods are able to solve for the invisible momenta, even without additional experimental input and without a second decay chain in the event. If the decay chain of Fig. 1 contained an additional two-body decay to a visible particle, just 5 events are sufficient for solving the event kinematics [22,27]. 3 Transversality is not strictly necessary, in fact it may even be beneficial to work with 3 + 1-dimensional variants of those variables [36][37][38][39][40][41][42][43][44][45][46][47].
• Inclusive methods. In this case, one focuses on the decay chain from Fig. 1 itself, disregarding what else is going on in the event. Using only the measured momenta of the visible SM decay products, i.e., the jet and the two leptons, one could form all possible invariant mass combinations 4 , namely m , m j n , m j f , and m j , measure their respective upper kinematic endpoints and use them to solve for the four input parameters (1.1). As just described, this approach is too naive, as it overlooks the remaining combinatorial problem involving the two leptons n and f . Since "near" and "far" cannot be distinguished on an event by event basis, the variables m j n and m j f are ill defined. This is why it has become customary to redefine the two jet-lepton invariant mass combinations as 5 to invert and solve for the input mass parameters (1.1). This procedure constitutes the classic kinematic endpoint method for mass measurements, which has been successfully tested for several SUSY benchmark points [77][78][79][80][81][82][83][84].
However, despite its robustness and simplicity, the kinematic endpoint method still has a couple of weaknesses. As we show below, taken together, they essentially lead to an almost flat direction in the solution space, thus jeopardizing the uniqueness of the mass determination. The first of these two problems is purely theoretical -it is well known that in certain regions of the parameter space (1.1) the four measurements (1.5) are not independent, but obey the relation [81] m max (1.6) In practice, this means that the measurements (1.5) fix only three out of the four mass parameters (1.1), leaving one degree of freedom undetermined. In what follows, we shall choose to parametrize this "flat direction" with the mass m A of the lightest among the four new particles D, C, B and A. One can then use, e.g., the first three of the measurements 4 In general, one is not limited to Lorentz-invariant variables only, e.g., recently it was suggested to study the peak of the energy distribution as a measure of the mass scale [72][73][74][75]. 5 A more recent alternative approach is to introduce new invariant mass variables which are symmetric functions of m j n and m j f , thus avoiding the need to distinguish n from f on an event per event basis [76].
in (1.5) and solve uniquely for the three heavier masses m D , m C , and m B , leaving m A as a free parameter. We list the relevant inversion formulas in Appendix A. The obtained one-parameter family of mass spectra will satisfy the three measured kinematic endpoints m max ll , m max jll , and m max jl(lo) by construction. What is more, in parameter space regions where eq. (1.6) holds, the family (1.7) will also obey the fourth measurement of m max jl(hi) , so that the four measurements (1.5) will be insufficient to lift the m A degeneracy in (1.7).
These considerations beg the following two questions, which will be addressed in this paper.
1. In the remaining part of the parameter space, where (1.6) does not hold and m max jl (hi) provides an independent fourth measurement, how well is the m A degeneracy lifted after all? With the explicit examples of Sections 3 and 4 below, we shall show that although in theory the additional measurement of m max jl(hi) determines the value of m A , in practice this may be difficult to achieve, since the effect is very small and will be swamped by the experimental resolution.
2. In the region of parameter space in which (1.6) holds, what additional measurement should be used, and how well does it lift the degeneracy? In the existing literature, the standard approach is to consider the constrained 6 distribution m jll(θ> π 2 ) , which exhibits a useful lower kinematic endpoint m min jll(θ> π 2 ) [85,86]. In what follows, we shall therefore always supplement the original set of 4 measurements (1.5) with the additional measurement of m min jll(θ> π 2 ) to obtain the extended set so that in principle there is sufficient information to determine the four unknown masses. Even then, we shall show that the sensitivity of the additional experimental input m min jll(θ> π 2 ) to the previously found flat direction (1.7) is very low. First of all, it is already well appreciated that the measurement of m min jll(θ> π 2 ) is very challenging, since 6 The distribution m jll(θ> π 2 ) is nothing but the usual m jll distribution taken over a subset of the original events, namely those which satisfy the additional dilepton mass constraint In the rest frame of particle B, this cut implies the following restriction on the opening angle θ between the two leptons [86] θ > π 2 , thus justifying the notation for m jll(θ> π 2 ) .
in the vicinity of its lower endpoint, the shape of the signal distribution is concave downward, which makes it difficult to extract the endpoint with simple linear fitting, and one has to use the whole shape of the m jll(θ> π 2 ) distribution [87]. Secondly, as we shall show in the examples below, the variation of the value of m min jll(θ> π 2 ) along the flat direction (1.7) can be numerically quite small, and therefore the sensitivity of the added fifth measurement along the flat direction (1.7) is not that great.
Either way, we see that the known methods for lifting the degeneracy of the flat direction (1.7) will face severe limitations once we take into account the experimental resolution, finite statistics, backgrounds, etc. [88,89] Thus the first goal of this paper will be to illustrate the severity of the problem, i.e. to quantify the "flatness" of the family of solutions (1.7). For this purpose, we shall reuse the study points from Ref. [89], which at the time were meant to illustrate discrete ambiguities, i.e. cases where two distinct points in mass parameter space (1.1) accidentally happen to give mathematically identical values for all five measurements (1.8). Here we shall extend those study points to a family of mass spectra (1.7) which give mathematically identical values for the first three 7 of the measurements (1.8), and numerically very similar values for the remaining measurements.
Having identified the problem, the second goal of the paper is to propose a novel solution to it and investigate its viability. Our starting point is the observation that the signal events from the decay chain in Fig. 1 populate the interior of a compact region in the (m j n , m , m j f ) space, whose boundary is given by the surface S defined by the constraint [17,90,91] S : which, for convenience, is written in terms of the unit-normalized variableŝ We note that both the shape and the size of the surface S depend on the input mass spectrum (1.1), i.e., S(m A , m B , m C , m D ), and this dependence is precisely what we will be targeting with our method to be described below.
In its traditional implementation, the kinematic endpoint method is essentially 8 using the kinematic endpoints (1.2) of the one-dimensional projections of the signal population onto each of the three axes m j n , m and m j f , as well as onto the "radial" direction m j = m 2 j n + m 2 + m 2 j f . This approach is suboptimal because it ignores correlations and misses endpoint features along the other possible projections. The only way to guarantee that we are using the full available information in the data is to fit to the three-dimensional boundary (1.9) itself [92,93], which will be the approach advocated here. As previously 7 And sometimes four, if we are in parts of parameter space where (1.6) holds. 8 The fact that one has to use m jl(lo) and m jl(hi) in place of m j n and m j f does not change the gist of the argument. observed in [92] (and extended to a broader class of event topologies in [94]), most of the signal events are populated near the phase space boundary (1.9), on which the signal number density ρ s formally becomes singular. This fact is rather fortuitous, since it implies a relatively sharp change in the local number density as we move across the phase space boundary, even in the presence of SM backgrounds (with some number density ρ b , which is expected to be a relatively smooth function). Thus, we need to develop a suitable method for identifying regions in phase space where the gradient of the total number density ρ ≡ ρ b + ρ s is relatively large, and then fit to them the analytical parametrization (1.9) in order to obtain the best fit values for the four new particle masses (1.1).
The first step of this program was already accomplished in our earlier paper [93], building on the idea originally proposed in [95] for finding "edges" in two-dimensional stochastic distributions of point data. Ref. [95] suggested that interesting features in the data, e.g., edge discontinuities, kinks, and so on, can be identified by analyzing the geometric properties of the Voronoi tessellation [96] of the data. 9 The volume v i of a given Voronoi cell generated by a data point at some location r i provides an estimate of the functional value of the number density ρ at that location, Therefore, in order to obtain an estimate of | ∇ρ( r)|, we can construct variables which compare the properties of the Voronoi cell and its direct neighbors. Among the different options investigated in Refs. [95,98], the relative standard deviation (RSD),σ i , of the volumes of neighboring cells, was identified as the most promising tagger of edge cells. The RSD was defined as follows. Let N i be the set of neighbors of the i-th Voronoi cell C i , with volumes, {v j }, for j ∈ N i . The RSD,σ i , is now defined bȳ where we have normalized by the average volume of the set of neighbors, Subsequently, in [93] we showed that this procedure for tagging edge cells can be readily extended to three-dimensional point data, as is the case here. The end result of the method was a set of Voronoi cells which have been tagged as "edge cell candidates" since their values ofσ i were above the chosen threshold [93]. With the thus obtained set of edge cells in hand, it appears that we are in a perfect position to perform a mass measurement, simply by finding the set of values for {m A , m B , m C , m D } which maximize the overlap between our tagged edge cells and the hypothesized surface S. We have checked that this approach 9 We note the existence of efficient codes for finding Voronoi tessellations in the form of the qHULL algorithms [97]. Wrappers that allow the use of these algorithms in many frameworks also exist, and in this work we use a private Python code to compute the geometric attributes of the Voronoi cells.
indeed works and gives a reasonable estimate of the true mass spectrum. However, here we prefer to suggest a slightly modified alternative, which accomplishes the same goal, but with somewhat better precision. The problem with fitting to a subset of the original data set (namely the set of Voronoi cells which happened to pass theσ i cut) is that we are still throwing away useful information, e.g., the Voronoi cells which barely failed the cut. In spite of formally failing, those cells are nevertheless still quite likely to be edge cells. Thus, in order to retain the full amount of information in our data, we prefer to abandon this "cut and fit" approach, and instead design a global variable which is calculated over the full data set. The only requirement is that the variable is maximized (or minimized, as the case may be) for the true values of the masses {m A , m B , m C , m D }.
In order to motivate such a variable, consider for a moment the case when the function ρ( r) is known analytically, then let us investigate the (normalized) surface integral for some arbitrary trial 10 values (m A ,m B ,m C ,m D ) of the unknown masses (1.1). The meaning of the quantity (1.14) is very simple: it is the average gradient of ρ( r) over the chosen surfaceS. We expect the dominant contributions to the integral to come from regions where the gradient is large, and we know that the gradient is largest on the true phase space boundary S(m A , m B , m C , m D ), defined in terms of the true values of the particle masses. However, if our choice for (m A ,m B ,m C ,m D ) is wrong, the integration surfaceS will be far from the true phase space boundary S, and those large contributions will be missed. The only way to capture all of the large contributions to the integral is to haveS coincide with the true S, and this is only possible if in turn the trial masses are exactly equal to the true particle masses. This suggests a method of mass measurement whereby the true mass spectrum is obtained as the result of an optimization problem involving the quantity (1.14).
Of course, in our case the analytical form of the integrand | ∇ρ( r)| is unknown, but we can obtain a closely related quantity using the Voronoi tessellation of the data. Following [93,95], we shall utilize the RSDσ i defined in (1.12), which has been shown to be a good indicator of edge cells, and replace the integrand in (1.14) as In other words, the gradient estimator 11 function g( r) is defined so that it is equal to the RSDσ i of the Voronoi cell C i in which the point r happens to be. Eqs. (1.14) and (1.15) 10 From here on trial values for the masses will carry a tilde to distinguish from the true values of the masses which will have no tilde. Correspondingly,S stands for a hypothesized "trial" boundary surface (1.9) obtained with trial values of the mass parameters. 11 Note that g( r) is not supposed to be an approximation for | ∇ρ( r)|, the crucial property for us is that the two functions peak in the same location.
suggest that the variable which we should be maximizing is It obviously depends on our choice of trial masses (m A ,m B ,m C ,m D ), and as argued above, we expect the maximum ofΣ to occur for the correct choice This hypothesis will be tested and validated with explicit examples below in Sections 3 and 4. The paper is organized as follows. In the next Section 2 we shall first review the well known formulas for the one-dimensional kinematic endpoints (1.8) and introduce the corresponding relevant partitioning of the mass parameter space into domain regions. In the next two sections we shall concentrate on the two most troublesome regions, (3,2) and (3,1), where the problematic relationship (1.6) holds. We shall pick one study point in each region, then study how well our conjecture (1.17) is able to determine the true mass spectrum. In principle, (1.17) involves optimization over 4 continuous variables, which is very time consuming (additionally, we have to perform the integration in the numerator of (1.16) by Monte Carlo). This is why for simplicity we choose to illustrate the power of our method with a one-dimensional toy study along the problematic flat direction (1.7). In particular, for each of our two study points we shall assume that the first three kinematic endpoints m max ll , m max jll , and m max jl(lo) are already measured, leaving us only the task of determining the remaining degree of freedom m A along the flat direction defined in (1.7). Correspondingly, we shall consider the whole family of mass spectra (1.7) which passes through a given study point. This family will eventually take us into the neighboring parameter space regions, including the third potentially problematic region, namely (2, 3), in which (1.6) is satisfied. For each family, we shall perform the following investigations • As a warm up, we shall first illustrate that for each of the three distributions, m ll , m jll , and m jl(lo) , the endpoint along the flat direction is the same (as expected by construction).
• We shall then investigate the variation of the kinematic endpoint of the m jl(hi) distribution along the flat direction (1.7). The endpoint value m max jl(hi) is expected to be constant in regions (3,2), (3,1) and (2,3), so the main question will be, how much does it vary in the remaining parameter space regions.
• We shall similarly investigate the variation of the lower kinematic endpoint m min along the flat direction (1.7). Together with the previous item, this will serve as an illustration of the main weakness of the classic kinematic endpoint method for mass measurements.
• Then we shall illustrate the distortion of the kinematic boundary surface (1.9) along the flat direction (1.7). The size of the distortion will be indicative of the precision with which one can hope to perform the mass measurement (1.17) using the kinematic boundary surface in phase space.
• Finally, we shall perform the fitting (1.17) along the flat direction parameterized bym A . We shall show results in two cases: (a) when the background events are distributed uniformly in m 2 phase space, and (b) when the background is coming from dilepton tt events.
We shall summarize and conclude in Section 5. Appendix A contains the inversion formulas needed to define the flat direction (1.7).
2 Endpoint formulas and partitioning of parameter space

Notation and conventions
Following [89], we introduce for convenience some shorthand notation for the mass squared ratios Note that in (2.1) there are only three independent quantities, which can be taken to be the set {R AB , R BC , R CD }. To save writing, we will also introduce convenient shorthand notation for the five kinematic endpoints as follows Note that these represent the kinematic endpoints of the mass squared distributions 12 .
In the next two sections we shall use the three endpoint measurements m max ll , m max jll , and m max jl(lo) to fix m D , m C and m B , leaving m A as a free parameter. Another way to think about this procedure is to note that the parameter space (1.1) can be equivalently parametrized as Then, the endpoint measurements of m max ll , m max jll , and m max jl(lo) can be used to fix the ratios R CD , R BC and R AB (see Appendix A), leaving the overall mass scale undetermined and parametrized by m A .

Endpoint formulas
The kinematical endpoints are given by the following formulas: Finally, the endpoint m min jll(θ> π 2 ) introduced earlier in the Introduction, is given by

Partitioning of the mass parameter space
One can see that the formulas (2.5-2.7) are piecewise-defined: they are given in terms of different expressions, depending on the parameter range for R CD , R BC and R AB . This divides the {R CD , R BC , R AB } parameter subspace from (2.3) into several distinct regions, illustrated in Fig. 2. Following [81], we label those by a pair of integers (N jll , N jl ). As already indicated in eqs. (2.5-2.7), the first integer N jll identifies the relevant case for m max jll , while the second integer N jl identifies the corresponding case for (m max jl(lo) , m max jl(hi) ). One can show that only 9 out of the 12 pairings (N jll , N jl ) are physical, and they are all exhibited For the purposes of this paper, only six of those regions will be in play, and we have color-coded them as follows: region (3, 1) in red, region (4, 1) in blue, region (3,2) in cyan, region (4, 2) in yellow, region (4,3) in magenta, and region (2, 3) in green.
within the unit square of Fig. 2. In what follows, an individual study point within a given region (N jll , N jl ) will be marked with corresponding subscripts as P N jll N jl .
Using (2.4), (2.5) and (2.7), it is easy to check that the "bad" relation (1.6), which can be equivalently rewritten in the new notation as b = a + d, (2.12) is identically satisfied in regions (3,1), (3,2) and (2,3) of Fig. 2. Therefore, as already discussed, in these regions one would necessarily have to rely on the additional information provided by the measurement of the e endpoint (2.11).
Before concluding this rather short preliminary section, we direct the reader's attention to the color-coding in Fig. 2, where we have shaded in color six of the parameter space regions: region (3, 1) in red, region (4, 1) in blue, region (3,2) in cyan, region (4, 2) in yellow, region (4,3) in magenta, and region (2, 3) in green. It will turn out that the two families of mass spectra considered in the next two sections will visit the six color-shaded regions. For the benefit of the reader, in the remainder of the paper we shall strictly adhere to this color scheme -for example, results obtained for a study point from a particular region will always be plotted with the color of the respective region: study points in region (3,1) are red, study points in region (4, 1)

Kinematical properties along the flat direction
In this section we shall study the flat direction (1.7) in mass parameter space which is generated by a study point P 31 from region (3, 1) (the same study point was used in [89] for a slightly different purpose). Table 1 lists some relevant information for the study point P 31 : the input mass spectrum (1.1), the corresponding mass squared ratios (2.1), and the predicted kinematic endpoints (1.8), also reminding the reader of the alternative shorthand notation (2.2). As discussed in the Introduction, starting from the point P 31 , we can follow a one-dimensional trajectory (1.7) through the parameter space (2.3) so that everywhere along the trajectory the prediction for the three endpoints a, b and c is unchanged (see Fig. 6 below). This trajectory is illustrated in Fig. 3, where we show its projections onto the three planes (R BC , R AB ) (left panel), (R AB , R CD ) (middle panel) and (R BC , R CD ) (right panel). The lines in Fig. 3 are parametrized by the continuous test mass parameterm A . For any given fixed value ofm A , the trajectory in Fig. 3 predicts the test values for the The lines are colored in accordance with the coloring convention for the regions depicted in Fig. 2. The red square marks the original study point P 31 from Table 1, while the circles denote the other three study points from Table 1: P 41 in region (4, 1) (blue circle), P 43 in region (4, 3) (magenta circle), and P 23 in region (2, 3) (green circle). other three mass parameters, namelym B ,m C andm D . This is shown more explicitly in  Fig. 2. Initially, as we move away from point P 31 (marked with the red square in Fig. 3), we are still within the red region (3,1), and the trajectory is therefore colored in red and parametrically given by eqs. (A.8-A.10). As the value ofm A is reduced from its nominal value (236.6 GeV) at the point P 31 , the mass spectrum gets lighter and eventually we reachm A = 0, where (the red portion of) the trajectory terminates at R AB = 0, R BC 0.67 and R CD 0.58. If, on the other hand, we start increasingm A from its nominal P 31 value, the spectrum gets heavier, and we start approaching the neighboring region (4,1). Eventually, at aroundm A ∼ 3600 GeV, the trajectory crosses into region (4, 1) and thus changes its color to blue. This transition is illustrated in the left panel of Fig. 5, where we plot the mass squared ratios R AB , R BC and R CD (solid lines), together with some other relevant quantities (dotted lines). In particular, the boundary between regions (3, 1) and (4, 1) is given by the relation R AB = R BD , see (A.2) and (A.29). We can see that crossover more clearly in the insert in the left panel of Fig. 5, where the line color changes from red to blue as soon as the R BD (dotted) line crosses the R AB (solid) line.
Once we are in region (4, 1), we follow the blue portion of the trajectory in Fig. 3, which is parametrically defined by eqs. (A.36-A.38). We choose a representative study point for region (4, 1) as well -it is denoted by P 41 and listed in the third (blue shaded) column of Table 1. The corresponding mass spectrum is clearly very heavy, but is nevertheless perfectly consistent with the three measured endpoints a, b and c, as shown in Fig. 6. As seen in Fig. 3, the blue portion of the mass trajectory appears headed for the point (R AB , R BC , R CD ) = (1, 1, 1), which is indeed reached in the limit ofm A → ∞, without ever entering into the neighboring region (1, 1) 13 . Fig. 3 reveals that the mass family (1.7) through our study point P 31 includes a segment which starts at (R AB , R BC , R CD ) = (0, 0.67, 0.58) and ends at (R AB , R BC , R CD ) = (1, 1, 1), visiting regions (3, 1) and (4, 1). Since the actual study point P 31 belongs to this segment, in what follows we shall refer to it as "the true branch". However, Fig. 3 also shows that there is an additional disconnected segment of the mass trajectory through the green region (2, 3) and the magenta region (4,3). In the following, we shall refer to this additional  Table 1: red for P 31 , blue for P 41 , magenta for P 43 and green for P 23 . segment as "the auxiliary branch". Note that this terminology is introduced only for clarity and should not be taken too literally -as far as the measured endpoints a, b and c are concerned, all points on the true and auxiliary branches are on the same footing, since the experimenter would have no way of knowing a priori which is the true branch and which is the auxiliary branch. This is why we have to seriously consider points on the auxiliary branch as well. We choose two representative study points, which are listed in the last two columns of Table 1: point P 43 belongs to the magenta region (4, 3), while point P 23 is in the green region (2, 3). As shown in Fig. 3 To summarize our discussion so far, we have imposed the three endpoint measurements a, b and c on the four-dimensional parameter space (1.1), reducing it to the one-dimensional parameter curve depicted in Figs. 3 and 4. The curve consists of two branches which visit four of the colored regions in Fig. 2, and we have chosen one study point in each region. The four study points are listed in Table 1, and their predicted invariant mass distributions from the ROOT phase space generator [99] are shown in Fig. 6: m in the left panel, m j (lo) in the middle panel and m j in the right panel. By construction, for any points along the mass trajectory (1.7), and in particular for the four study points from Table 1, these distributions share common kinematic endpoints. Furthermore, as Fig. 6 reveals, the shapes of most distributions are also very similar, which makes it difficult to pinpoint our exact location along the mass trajectory (1.7). This is why in the remainder of this section, we shall focus on the question, what additional measurements may allows us to discriminate experimentally points along the two branches in Figs. 3 and 4, and in  The same as Fig. 7, but for the endpoint √ e and the corresponding distribution m jll(θ> π 2 ) . particular distinguish between the four study points in Table 1.
One obvious possibility is to investigate the remaining kinematic endpoints d and e, which are analyzed in Figs. 7 and 8, respectively. The left panels show the theoretical predictions for the kinematic endpoints √ d = m max j (hi) and √ e = m min jll(θ> π 2 ) along the flat direction (1.7) as a function ofm A , while the right panels exhibit the corresponding invariant mass distributions for each of our four study points from Table 1.
Let us first focus on Fig. 7 which illustrates them A dependence of the m j (hi) distribution and its kinematic endpoint √ d. As we have already discussed, in regions (3, 1) and (2, 3) the additional measurement of √ d is not useful, since it is not independent -the value of d is predicted by the relation (2.12), as confirmed by the left panel in Fig. 7, where the red and green dotted lines representing those two regions are perfectly flat and insensitive tom A . However, this still leaves open the possibility that in the remaining two regions, namely (4, 1) and (4,3), the measurement of the d endpoint will be able to lift the degeneracy and determine the value of m A , since, at least in theory, d is a non-trivial function ofm A , see (A.39b) and (A.63b). Unfortunately, Fig. 7 demonstrates that this is not the case in practice -them A dependence is extremely weak, and the endpoint value for √ d only changes by a few tens of MeV asm A is varied over a range of several TeV! This lack of sensitivity is the reason why we have been referring to the family of mass spectra (1.7) as a "flat direction" in mass parameter space. Clearly, due to the finite experimental resolution, an endpoint measurement with a precision of tens of MeV is not feasible, the anticipated experimental errors at the LHC are significantly higher, on the order of a few GeV [100].
It is instructive to understand this lack of sensitivity analytically, by studying, e.g. the mathematical expression (A.39b) for d which is relevant for region (4, 1). Figs. 4 and 7 already showed that region (4, 1) occurs at large values ofm A , where the spectrum is relatively heavy -on the order of several TeV. At the same time, the measured parameter inputs into (A.39b), namely the endpoints a, b and c, are all on the order of several hundred GeV. This suggests an expansion in terms of 1/m A as Using (A.39b), we get the expansion coefficients to be Interestingly, the numerical value of K 0 is extremely close to b − a: Since K 0 is the leading order prediction for d, (3.4) implies that even in region (4, 1), the relation (2.12) will still hold to a very good approximation -any deviations from it will be 1/m A suppressed. We can formalize this observation by introducing the value m

(b)
A which the parameterm A takes when the mass trajectory (1.7) crosses the boundary between regions (3, 1) and (4, 1). Using the continuity of the function d(a, b, c,m A ), we can write where the left-hand side is the value of d in region (3, 1) which is given by (A.11), while the right-hand side is the value of d as predicted by the Taylor expansion (3.1) in region (4, 1). Eliminating K 0 from (3.5), we can rewrite the expansion (3.1) in the form which manifestly shows that the deviations from the relation (2.12) are 1/m A suppressed. One can check that the sign of the K 1 coefficient (3.3) is positive, then (3.6b) explains why d is a decreasing function ofm A in region (4, 1), as observed in the left panel of Fig. 7.
Starting from (A.63b), one can repeat the same analysis for the magenta portion of the auxiliary branch which is located in region (4, 3). As the left panel of Fig. 7 shows, the conclusions will be the same -the d endpoint is still given approximately by the "bad" relation (2.12), and the corrections to it are tiny and 1/m A suppressed. The right panel in Fig. 7 explicitly demonstrates that the variation of the d endpoint along the flat direction is unnoticeable by eye even with perfect resolution, large statistics and no background. The shapes of the m j (hi) distributions are also very similar. As a result, we anticipate that the additional measurement of the d kinematic endpoint and the analysis of the associated m j (hi) distribution will not help much in lifting the degeneracy of the flat direction (1.7).
We now turn to the discussion of the fifth and final kinematic endpoint, e, illustrated in Fig. 8. The left panel now shows a more promising result -the variation along the flat direction is much larger than what we saw previously in Fig. 7. This is especially noticeable for the auxiliary branch, where the prediction for √ e can vary by as much as 17 GeV, suggesting that one might be able to at least rule out some portions of it. At the same time, the variation of √ e along the true branch is only 4 GeV, once again making it rather difficult to pinpoint an exact location along the true branch. Unfortunately, these theoretical considerations are dwarfed by the experimental challenges in measuring the e endpoint, as suggested by the right panel of Fig. 8. Unlike the other four kinematic endpoints, e is a lower endpoint (a.k.a. "threshold"), which places it in a region where one expects more background. More importantly, the signal distribution is very poorly populated near its lower endpoint -the vast majority of signal events appear sufficiently far away from the threshold, and the measurement will suffer from a large statistical uncertainty. This casts significant doubts on the feasibility of this measurement -in previous studies, the √ e endpoint was either the measurement with the largest experimental error from the fit (on the order of 10 GeV [79]), or one could not obtain a measurement for it at all [81].
One could hope to improve on the precision by utilizing shape information [101], but this introduces additional systematic uncertainty, since the background shape and the shape distortion due to cuts has to be modeled with Monte Carlo. Being mindful of the challenges involved with the measurement of the e endpoint, in this paper we shall look for an alternative method for lifting the degeneracy along the flat direction. Our proposal is to study the shape of the kinematic boundary (1.9), which is a two-dimensional surface in the three-dimensional space of observables As a proof of principle, we first illustrate the change in the shape of the surface (1.9) as we move along the flat direction. Our results are shown in Fig. 9 (for the true branch) and in Fig. 10 (for the auxiliary branch). Following [93], we visualize the surface (1.9) by showing a series of two-dimensional slices in the (m 2 j (lo) , m 2 j (hi) − m 2 j (lo) ) plane, where the slight modification of the "y-axis" was done in order to avoid wasted space on the plots due to the unphysical areas with m j (lo) > m j (hi) . Each slice is taken at a fixed value of m 2 , starting from a very low value (10 GeV 2 ) and going up all the way until the kinematic endpoint (m max ) 2 = 20, 976 GeV 2 . The red solid lines in Fig. 9 correspond to Figure 10. The same as Fig. 9, but for the auxiliary branch going through regions (2, 3) and (4, 3). The dashed lines represent points withm A = 100 GeV (black),m A = 500 GeV (green), m A = 2000 GeV (blue) andm A = 6000 GeV (yellow). For reference, we also show the case of the true mass spectrum for point P 31 (red solid lines), although P 31 does not belong to the auxiliary branch.
the nominal case of the study point P 31 . In each panel, the signal events will be populating the areas delineated by these red solid lines. As pointed out in [92], the density of signal events is enhanced near the phase space boundary, i.e. signal events will cluster close to the solid red lines; this property can be incorporated into the algorithm for detecting the surface boundary [93]. It is worth noting that in general, each panel contains two signal populations, which arise from the reordering (1.3-1.4) [89]. As we vary the value of m 2 , the shape of the red solid lines changes in accordance with eq. (1.9), which follows from simple phase space considerations. However, the main purpose of Figs. 9 and 10 is to check how much the shape is modified relative to the nominal case of P 31 when we vary the value ofm A along the flat direction (1.7). The dashed lines in Fig. 9 show results for several representative values ofm A along the true branch:m A = 0 (black),m A = 100 GeV (gray),m A = 500 GeV (green),m A = 1000 GeV (blue),m A = 2000 GeV (yellow) andm A = 5000 GeV (magenta). We observe noticeable shape variations, especially at low to intermediate values of m 2 , which bodes well for our intended purpose of measuring the value of m A . Fig. 9 aids in visualizing why sensitivity is lost when performing onedimensional projections. Consider, for example the variable m j (lo) . The top two rows of Fig. 9 show that asm A is varied along the flat direction, the boundary contours are being stretched vertically, which does not have any effect on the m j (lo) endpoint. Later on, when the events are projected vertically on the m j (lo) axis to obtain the m j (lo) distribution seen in the middle panel of Fig. 6, the effects from this vertical stretching tend to be washed out and the resulting m j (lo) distributions have very similar shapes. Fig. 10 shows the analogous results for the auxiliary branch. Once again, the red solid lines represent the study point P 31 , while the dashed lines correspond to four values ofm A :m A = 100 GeV (black),m A = 500 GeV (green),m A = 2000 GeV (blue) and m A = 6000 GeV (yellow). This time the shape variation along the flat direction is much more significant compared to what we saw in Fig. 9. This observation agrees with our expectation based on Fig. 8 that points on the auxiliary branch behave quite differently from our nominal study point P 31 , especially at lowm A .

A toy study with uniformly distributed background
In the remainder of this section we shall illustrate our proposed method for mass measurement with two exercises. In each case, we shall assume that the standard set of onedimensional kinematic endpoints (1.2) has already been well measured and used to reduce the relevant mass parameter space (1.1) to the flat direction (1.7) parametrized by the test massm A for the lightest new particle A. This is done only for simplicity -in principle, our method would also work without any prior information from endpoint measurements, but by using those, we are reducing the 4-dimensional optimization problem in (1.17) to the much simpler one-dimensional optimization problem (3.8) wherem B (m A ),m C (m A ), andm D (m A ) are the masses of particles B, C and D along the flat direction. Our main emphasis here is on demonstrating the advantages of our method relative to the method of kinematic endpoints. In Section 3.1 we already showed that while the method of kinematic endpoints does a good job in reducing the unknown mass parameter space (1.1) to the flat direction (1.7), it does a poor job of lifting the degeneracy along the flat direction. Thus, if we can show that our method can perform the remaining mass measurement along the flat direction, we will have accomplished our goal. In order to make contact with our previous studies in [93], we begin with a simple toy exercise where in addition to the signal events from the cascade decay in Fig. 1, we also consider a certain number of background events, which we take to be uniformly distributed in the mass squared space of observables (3.7). While the assumption of uniform background density is unrealistic, such an exercise is nevertheless worth studying for several reasons. First, our method is completely general and applies in any situation where we have a decay of the type shown in Fig. 1, while to correctly identify the relevant backgrounds, we must be a lot more specific -we need to fix the signature, the type of production mechanism (which determines what else is in the event), the cuts, etc. In order to retain generality, we choose to avoid specifying those details and instead we generate background events by pure Monte Carlo according to a flat hypothesis. Second, as shown in [93], a uniform background distribution is actually a pretty good approximation to more realistic backgrounds resulting, e.g., from dilepton tt events (compare to the results in Section 3.3 below). Finally, our method is attempting to detect a discontinuity in the measured event density caused by a signal kinematic boundary, so the exact shape of the background distribution is not that important, as long as it is smooth and without any sharp kinematic features.
In order to detect the exact location of the kinematic boundary, we shall be computing the quantityΣ defined in (1.16) along the flat direction (1.7), i.e. (3.9) We shall perform several versions of the exercise, with varying levels of signal-to-background. For this purpose, we vary the ratio of signal to background events inside the true "samosa" surface S(m A , m B , m C , m D ): where V S is the volume inside the samosa S(m A , m B , m C , m D ), while ρ s and ρ b are the signal and background event densities from Section 1, respectively. In this exercise, we shall fix the overall normalization by choosing N B = 1000 background events inside S. Our main result is shown in Fig. 11 There are several important lessons from Fig. 11: • Viability of the method. We see that in each panel, the maximum ofΣ is obtained for a value ofm A which is close to the true value m A = 236.6 GeV. This validates our conjecture 15 , eq. (1.17), and proves the viability of our method.
• Precision of the method. Of course, we did not recover exactly the input value for m A , but in each case, came relatively close. Each panel of Fig. 11 contains an insert which zooms in on the region near the peak, which is sampled more finely.  perfect, it may be instructive to compare the theoretical boundary for the input study point P 31 to the boundary surface found by the fit. This is illustrated in Fig. 12, where in analogy to Figs. 9 and 10 we show two-dimensional slices at fixed m 2 of the Voronoi tessellation of the data for the case of S/B = 3. The Voronoi cells are color coded by their value ofσ i defined in (1.12). As in Figs. 9 and 10, the red solid line in each panel is the expected signal boundary for the nominal case of point P 31 . We notice that the cells with the highest values ofσ i are indeed clustered near the nominal boundary, in agreement with the results from Refs. [93,95]. On the other hand, the boundary delineated by the black dashed lines in Fig. 12 corresponds to the best fit value ofm A = 280 GeV, which was found in the upper left panel in Fig. 11. The difference between the solid red and black dashed contours in Fig. 12 is essentially a measure of the resolution of our method.
• Elimination of the auxiliary branch and the largem A tail of the true branch. One very positive piece of news from Fig. 11 is that the whole auxiliary branch has very low values forΣ which makes it easy to rule it out -one can see that no point on the auxiliary branch was ever in contention for the top spot. Similar comments, albeit to a lesser extent, also apply to the long tail along the true branch at largem A . In particular, region (4, 1) seems to be ruled out, as well as the largem A portion of region (3,1). In effect, the range of possible values form A along the flat direction has been significantly narrowed down to a small interval within a few tens of GeV of the true value m A .
• The adverse effect of the background. Comparing the different panels in Fig. 11, we see that as we make S/B smaller, the difference between the true and auxiliary branch is reduced, but the auxiliary branch is still disfavored. As for the true branch, the peak near m A still persists, even in the case when the data is dominated by background events. This is not surprising, since the background distribution is relatively smooth, so that in the background-dominated regions of phase space there aren't too many Voronoi cells with large values ofσ i , which could adversely affect the fit.

A study with tt dilepton background events
We are now in position to repeat the exercise from Section 3.2, with signal events from D+A associated production and background taken from dilepton tt events, which represent the main background to the signature from Fig. 1 Fig. 1, and considered the associate production of a squarkq with the lightest neutralinoχ 0 1 , namely pp →qχ 0 1 [103,104]. Since each tt background event contains two jets, there is a two-fold ambiguity in the jet selection. We will use both possible pairings, so that each background event will contribute two entries to our data. Of course, we do not know a priori how many of those entries will end up inside the nominal boundary surface S(m A , m B , m C , m D ), which is why we have to use a slightly different normalization from Sec 3.2. We shall fix the number of signal events to N S = 3000, and then we shall consider several values 16 for the number of dilepton tt events: N B = {3000, 4000, 5000, 6000}. From Monte Carlo we then find that these choices correspond to S/B = {1.52, 1.14, 0.91, 0.76} inside the S boundary, see (3.10).
Our main result is shown in Fig. 13, which is the analogue of Fig. 11 for this case. Once again, we find that the functionΣ(m A ) is maximized in the vicinity ofm A = m A = 236.6 GeV. The fitting procedure is illustrated in Fig. 14, which is the analogue of Fig. 12. The red solid lines show the boundary contours for the nominal value of m A = 236.6 GeV, while the black dashed lines are for the best fit value ofm A = 280.0 GeV, which was found in the top left panel of Fig. 13. Fig. 13 shows that once again, our procedure has disfavored the whole auxiliary branch and narrowed down the range of viable values ofm A to a few tens of GeV around the nominal value m A .

A case study in region (3, 2)
In this section, we shall repeat the analysis from Section 3, only this time our nominal study point, from now on labelled as P 32 , will be chosen within the cyan region (3, 2) of Fig. 2. Recall that the problematic relation (2.12) was satisfied in three of the colored regions in Fig. 2, namely (2, 3), (3,1) and (3,2). The former two regions, (2, 3) and (3, 1), 16 The anticipated signal-to-background ratio is model-dependent. In this sense, SUSY may not be the best case for discovery, since other scenarios, e.g., UED [105][106][107], have higher signal cross-sections. were already visited by the mass trajectory studied in Section 3, thus here for completeness we will also illustrate the case of region (3, 2).

Kinematical properties along the flat direction
The mass spectrum for the study point P 32 and the corresponding mass squared ratios and kinematic endpoints are shown in the cyan-shaded column of Table 2. Point P 32 was used previously in Ref. [89] as an example of a discrete two-fold ambiguity, while here it serves to define a flat direction (1.7) in mass parameter space. This flat direction is illustrated in Figs. 15 and 16, which are the analogues of Figs. 3 and 4, respectively. According to Figs. 15 and 16, the mass trajectory now goes through five of the six colored regions in Fig. 2: there is a true branch through the red region (3, 1), the cyan region (3,2) and the the yellow region (4, 2), as well as an auxiliary branch through the yellow region (4, 2), the magenta region (4, 3) and the green region (2, 3). As in Section 3, we choose one representative study point in each of these regions. The four additional study points, P 31 , P 42 , P 43 and P 23 , are also listed in Table 2, and their columns are shaded with the color of their respective regions in Fig. 2. The flat direction depicted in Fig. 15 is again parametrized by the trial valuem A for the mass of the lightest new particle A. However, as seen in Fig. 16, this time the allowed range form A does not extend all the way tom A = 0, and instead the true and auxiliary branch meet inside the yellow region (4, 2) around at the lowest valuem A ∼ 89 GeV.
The transitions between two neighboring regions along the flat direction can be under-  Table 2. The same as Table 1, except now the starting point is a point (P 32 ) from region (3,2). Figure 15. The same as Fig. 3, but for the flat direction generated by point P 32 from Table 2. stood from Fig. 17, which plots the mass squared ratios R AB , R BC and R CD (solid lines) and several other quantities (dotted lines) which are relevant for defining the regions from   By construction, all five study points from Table 2 predict identical values for the three kinematic endpoints a, b and c. This is demonstrated in Fig. 18, which is the analogue of Fig. 6, but for the five study points from Table 2. As before, the distributions in Fig. 18 are color-coded according to our color conventions from Fig. 2. The nominal input study point P 32 is represented by the solid line, while the dotted lines mark the other four study points. Given that the five study points look very similar on Fig. 18, we now focus on the remaining two distributions, m j (hi) and m jll(θ> π 2 ) , which are investigated in Figs. 19 and Figure 18.
The analogue of Fig. 6, but for the five study points exhibited in Table 2. The distributions are color-coded according to our color conventions for the regions in Fig. 2. Figure 19. The analogue of Fig. 19, but for the flat direction defined in Fig. 15  By now, we should not be surprised by the extreme flatness of the curves exhibited in the left panel of Fig. 19. The mass trajectory from Fig. 15 passes through all three of the regions where the endpoint d is not an independent quantity, but is fixed by the relation (2.12) and is therefore strictly independent ofm A . In the remaining two regions, (4, 2) and (4, 3), Fig. 19 shows a maximal deviation of only 2 GeV from the prediction (2.12). Taken together, the left panels of Figs. 7 and 19 justify our terminology of the mass trajectory (1.7) as a "flat direction" in mass parameter space.
On the other hand, the left panel of Fig. 20 shows a much more significant variation of the kinematic threshold variable √ e along the flat direction. The total variation is on the order of 15 GeV, which is of the same order as our previous result in Fig. 8. However, as we already discussed in Section 3, the measurement of √ e presents significant experimental distributions shown in the right panel of Fig. 20. This motivates searching for alternative methods for lifting the degeneracy along the flat direction. As already discussed in Section 3, one such method is to track the deformation of the shape of the kinematic boundary (1.9) along the flat direction. The effect is illustrated in Figs Fig. 9, we observe a much larger variation in the shape of the kinematic boundary in the present case, which promises good prospects for the mass measurement exercise to follow. The results shown in Fig. 22 for the auxiliary branch are also quite good. This should not come as a surprise, since the exercise in Section 3 already indicated that the auxiliary branch has a different kinematic behavior, as reflected in the shape of the phase space boundary.

A toy study with uniformly distributed background
In the remainder of Section 4 we shall repeat the two exercises from Sections 4.2 and 4.3, only this time we shall use P 32 as our input study point, and perform the measurement along the corresponding flat direction described in Figs. 15 and 16.
First we consider the case of uniformly distributed (in mass squared) background events, and proceed to evaluate the quantityΣ along the flat direction. As before, we fix N B = 1000 and then vary the signal-to-background ratio inside the boundary surface S.  In all four cases, the auxiliary branch is disfavored, as it always gives low values forΣ, while the true branch is restricted to a very narrow region near the true mass spectrum.

A study with tt dilepton background events
Our final task will be to repeat the P 32 exercise with dilepton tt events as was done in Section 3.3. As before, we fix the number of signal events N S = 3000 and then consider several values for the number of background events: N B = {3000, 4000, 5000, 6000}. In each case, we compute the functionΣ(m A ) along the flat direction of Fig. 15. The results are shown in Fig. 25, which has the same qualitative behavior as Fig. 23. TheΣ values for the Figure 23. The same as Fig. 11, but now taking point P 32 as input and measuring along the flat direction depicted in Fig. 15. auxiliary branch tend to be low, and the branch is disfavored. The global peak ofΣ(m A ) is again found in the vicinity of the right answer (for N B = {3000, 4000, 5000, 6000}, the peak is atm A = {116, 125, 125, 125} GeV), and the largem A tail of the true branch is also disfavored. One final consistency check is provided by Fig. 26, which shows a comparison of the kinematic boundaries for the nominal study point P 32 with m A = 126.5 GeV (red solid lines), and the boundaries for the best fit valuem A = 125 GeV (black dashed lines).

A detector level study
In this paper, we introduced the new Voronoi-based method for mass measurement as a proof of principle, and showed that at the parton level it does reasonably well in the two examples considered so far in Sections 3 and 4. Before concluding, we would like to also test the method in the presence of detector effects (this subsection) and combinatorics (see Sec. 4.5 below). For this purpose, we first repeat the exercise from Section 4.3, only this time we account for the finite detector resolution by smearing the jet energies with the typical hadronic calorimeter resolution and electromagnetic calorimeter resolution in CMS [108], with the energy measured in GeV. Smearing of the muon momenta is done according to the "Full System" values in Fig. 1.5 of [108]. The result of the fitting exercise is shown in Fig. 27. We see that the peak structure in the vicinity of the correct mass value (126.5 GeV) is preserved, but somewhat degraded due to the detector resolution effects. Figure 25. The same as Fig. 13, but using study point P 32 as input.

D-pair production and combinatorics effects
Our proposed mass measurement method uses a single decay chain like the one depicted in Fig. 1. In that sense, the method is inclusive and model-independent, since it does not depend on what else is going on in the event. In particular, the method is equally applicable when particle D is produced singly, in pairs, or in association with another object. Nevertheless, a well-motivated and widely studied class of models are the SUSYlike dark matter scenarios in which all particles A, B, C and D carry negative parity under the additional Z 2 symmetry. In that case, D has to be produced in association with another negative parity object. If A is a neutral dark matter candidate, then D must carry color, therefore D pair-production is strong and may dominate the inclusive cross-section for D production.
The presence of a second D decay chain in the event can either be a blessing or a curse. If both D particles decay the same way, as in Fig. 1, we can attempt to double our statistics by considering the second decay chain as well. However, this comes at a cost, since now we have to face the combinatorial problem of associating the different reconstructed final state objects to one of the two decay chains. First, there is a two-fold ambiguity of associating each jet to the correct side, and furthermore, there is an additional two-fold ambiguity in the case when all leptons are the same flavor. The simplest approach would be to consider all possible combinations and use the resulting data set for building the Voronoi tessellation, then proceeding with the fitting of the boundary surface as before. The result from this exercise is shown in Fig. 28 for the case of 100 (left panel) and 500 (right panel) signal events. We see that despite the pollution from wrong combinatorics, the peak inΣ is still very well visible, and found in the right location.
Our final study is reserved for the most challenging example so far -a case with a severe combinatorics problem, in the presence of SM (tt) background events. For signal, let us again consider D-pair production, only this time let all three of the decay products of the second D particle be QCD jets. Thus, each one of our signal events has 2 leptons and 4 jets, and picking the correct jet becomes a difficult task. Now, instead of using all possible  combinations, we design a preselection cut in order to improve our chances of capturing the correct jet pairing. For this purpose, we consider the four possible jet-lepton-lepton combinations and compute the corresponding three-body invariant masses. Next, we rankorder these four values [17,18] and eliminate from further consideration the two jets which correspond to the two largest jet-lepton-lepton invariant masses, since those jets are very likely to come from the decay chain opposite the two leptons. The remaining two jets still cause a two-fold ambiguity, which we handle as in Fig. 28: by simply plotting both Figure 29. The same as the lower two plots in Fig. 25, but for signal events where particles D are pair-produced and one of them decays as in Fig. 1, while the other decays to 3 jets and particle A.
combinations. The end result from the analysis is shown in Fig. 29. As in the example from Sec. 4.3, here we also include a certain number of tt background events: 5000 in the left panel and 6000 in the right panel. The fact that theΣ peak is again obtained in the correct location indicates that our method can be viable in the presence of combinatorial background due to pair production of particles D.

Conclusions
In this paper, we reconsidered the classic endpoint method for particle mass determination in SUSY-like decay chains like the one shown in Fig. 1. Our main points are: • We have identified a "flat direction" in the mass parameter space (1.1), along which mass differences can be measured relatively well, but the overall mass scale remains poorly constrained. (The analytical formulas parametrizing this flat direction can be found in Appendix A.) We quantified the problem with examples of specific study points, P 31 and P 32 , considered in Sections 3 and 4, respectively.
• We then proposed a new method for mass measurements in general, and for extracting the mass scale along the flat direction, in particular. The method takes advantage of the changes in the shape of the two-dimensional kinematic boundary surface within the fully differential three-dimensional space of observables, as one moves along the flat direction. We have tested our Voronoi-based algorithm [93,95] for detecting the boundary surface and demonstrated that it can be usefully applied in order to lift the degeneracy along the flat direction. This approach represents the natural extension of the one-dimensional kinematic endpoint method to the relevant three dimensions of invariant mass phase space.
• We introduced a new variable,Σ, which is the average RSD per unit area, calculated over the hypothesized kinematic boundary. We showed that the location of theΣ maximum correlates very well with the true values of the new particle masses, see Figs. 11, 13, 23, 25, and 27.
The work reported here can be extended in several directions. First of all, the method can be readily generalized to longer decay chains with more visible particles, where the boundary enhancement is even more pronounced [94], and therefore, the detection of the boundary surface should be in principle easier. One could also try to apply Voronoi-based boundary detection algorithms for the discovery of new physics. It is also interesting to develop a general and universal method for estimating the statistical significance of the local peaks found in Figs. 11, 13, 23 and 25, and hence the statistical precision of our mass measurement. These, along with many other interesting questions, will be investigated in future studies [109].

A Inverse formulas
In this appendix we derive the inverse relations which define m B , m C and m D in terms of the three measured endpoints The case of region (3,2) Region (3,2) is defined by the following two conditions: The kinematic endpoints are given by the following formulas:  The case of region (2, 3).
Region (2,3) is defined by the following condition: The kinematic endpoints are given by the following formulas: The masses of B, C and D are given by , and while the term (2c − b + a) can have either sign, the discriminant (i.e., the term inside the square root in (A.26)) is always larger than a 2 (2c − b + a) 2 , which guarantees a nonnegative result for m 2 B . In this region, the relation (2.12) is again satisfied, so that d is again given by (A.11), which can be used to cross-check the result (A.26-A.28).
Region (4, 1) is defined by the following conditions: For the case (4,1) the endpoints are given by the following formulas: The case of region (4, 2).
Region (4,3) is defined by the following conditions: For the case (4,3) the endpoints are given by the following formulas: where again the masses are calculated in the order m D , m B and then m C . The d endpoint is given by