Using sorted invariant mass variables to evade combinatorial ambiguities in cascade decays

The classic method for mass determination in a SUSY-like cascade decay chain relies on measurements of the kinematic endpoints in the invariant mass distributions of suitable collections of visible decay products. However, the procedure is complicated by combinatorial ambiguities: e.g., the visible final state particles may be indistinguishable (as in the case of QCD jets), or one may not know the exact order in which they are emitted along the decay chain. In order to avoid such combinatorial ambiguities, we propose to treat the final state particles fully democratically and consider the sorted set of the invariant masses of all possible partitions of the visible particles in the decay chain. In particular, for a decay to N visible particles, one considers the sorted sets of all possible n-body invariant mass combinations (2<= n<= N) and determines the kinematic endpoint m_(n,r)^max of the distribution of the r-th largest n-body invariant mass m_(n,r) for each possible value of n and r. For the classic example of a squark decay in supersymmetry, we provide analytical formulas for the interpretation of these endpoints in terms of the underlying physical masses. We point out that these measurements can be used to determine the structure of the decay topology, e.g., the number and position of intermediate on-shell resonances.


Introduction
Now that the long awaited Higgs boson of the Standard Model (SM) appears to have been discovered [1][2][3], the best evidence for new physics Beyond the Standard Model (BSM) is provided by the dark matter problem [4]. Dark matter particles can be produced directly at high energy colliders like the Large Hadron Collider (LHC) at CERN [5][6][7]. However, the expected rates are relatively low, since dark matter is (super)weakly interacting. In generic models, therefore, the indirect dark matter production at colliders is much more copious, with dark matter particles appearing in the decay chains of heavier (perhaps colored or charged) new particles. One such possible decay chain is shown in Fig. 1, where a heavy new particle D decays successively to lighter particles C, B, and A, the latter perhaps being a dark matter candidate. This particular decay chain is rather ubiquitous in models of low energy supersymmetry 1 (SUSY), where particle D is a squark, C is a heavy neutralino, B is a slepton, and A is the lightest neutralino (the typical dark matter candidate in SUSY). 1 An analogous decay chain can be present in many other BSM models, e.g. in Universal Extra Dimensions (UED) [8], where D is a Kaluza-Klein (KK) quark, C is a KK Z-boson, B is a KK lepton, and A is a KK photon [9]. 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 [50-76] 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 - Figure 1. The typical cascade decay chain used in mass measurement studies. 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. Results for the masses of the new particles are often quoted in terms of the D mass m D and the three dimensionless squared mass ratios R CD , R BC and R AB .
The main difficulty in the analysis of the decay chain of Fig. 1 stems from the fact that particle A, being a dark matter candidate, is invisible in the detector, hence its energy and momentum are not directly measured. This makes the problem of determining the masses and spins of the new particles A through D rather challenging. Over the last 20 years, a fairly large body of literature 2 has been devoted to this problem. Among the different approaches which have been proposed, the classic method of kinematic endpoints is arguably the most popular and robust technique for mass determination. With this method, one studies the invariant mass distributions of different combinations of the visible decay products and attempts to locate their kinematic endpoints (generally in the presence of some background continuum). In the example of the decay chain in Fig. 1, one can form three 2-body invariant mass variables (m j n , m j f , and m ) and one 3-body invariant mass variable, m j , so that the basic set of invariant mass variables is m j n , m j f , m , m j . (1.1) Ideally, one would like to study each of the variables (1.1) individually and obtain the corresponding kinematic endpoints m max j n , m max j f , m max , and m max j , which can be simply related to the unknown masses {m A , m B , m C , m D } by analytic expressions available in the literature (see, e.g., [12][13][14]). Unfortunately, the situation is not that simple, as it becomes muddled by various combinatorial ambiguities: • Partitioning ambiguity. In general, in addition to the three visible objects from the decay chain in Fig. 1, there will be a number of additional objects in the event -the decay products from the other 3 decay chain in the event, jets from initial state radiation, pile-up, etc. The question then becomes, how to select the correct objects j, n and f for the analysis variables (1.1). This is a really challenging problem, for which no universal solution exists, although several ideas have been tried, including the "hemisphere method" [15,16], (a combination of) invariant mass and M T 2 cuts [17][18][19][20], and neural networks [21]. In this paper, we will ignore the partitioning ambiguity and instead focus on the ordering ambiguity described next.
Our assumption is justified in the case when particle D is produced singly, or when D is produced in association 4 with the stable particle A, so that a second decay chain simply does not exist.
• Ordering ambiguity. Having selected the correct visible objects arising in a given decay chain, we still have to decide on the order in which they are emitted along the chain. For example, motivated by SUSY, in Fig. 1 one makes the specific assumption that the jet comes first, followed by the two leptons. However, even with this extra assumption, the ambiguity is not completely resolved, as we still do not know the exact ordering of the two leptons. In other words, we are not justified in using the labels "near" and "far" to denote the two leptons, which makes it impossible to study separately the distributions of m j n and m j f in the real experiment.
Two possible ways out of this conundrum have been suggested. The standard approach [12] is to trade the variables m j n and m j f for their ordered cousins Then, instead of the set (1.1), one can consider the alternative set of variables measure their respective endpoints, and from those extract the physical mass spectrum [12][13][14]. A more recent alternative approach [24] introduces 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.
However, while both of these approaches are designed to address the ordering ambiguity problem, it is our view that they do not go far enough -in the sense that the assumption of the jet being the first emitted particle is still hardwired in the analysis from the very beginning. From the point of view of an experimenter whose duty is to perform unbiased measurements without theoretical prejudice, there is no compelling reason to make that assumption. One can easily construct theory models 5 in which the jet is the second (or the third) visible particle in the diagram of Fig. 1. In order to account for such scenarios, one needs to further generalize the reordering in Eqs. (1.2,1.3) to include swapping the jet with one of the leptons. To be concrete, for the example of Fig. 1, we propose to further 4 A well-known such example in SUSY is provided by the associated squark-neutralino production [22,23]. 5 Granted, such models will contain intermediate particles with unusual, "leptoquark", quantum numbers.
The generic decay of a heavy resonance D to N indistinguishable visible particles (solid black lines) and one invisible particle A (red dashed line).
replace the two-body invariant mass variables m j (lo) , m j (hi) and m from (1.4) with the three ordered variables where we have used the function max r to denote the r-th largest among a given set of elements 6 . The variable set (1.4) is then replaced by m (2,1) , m (2,2) , m (2,3) , m j , (1.8) whose kinematic endpoints can then be used to extract the physical spectrum. To this end, however, one would first need to derive the analytic formulas for these new endpoints in terms of the physical masses, and this will be one of the main goals of this paper. In summary, in this paper we propose to evade the ordering ambiguity problem by considering all possible partitions of the visible particles resulting from a given cascade decay chain, and then studying the kinematic endpoints of the corresponding sorted invariant mass variables. We shall try to keep our discussion as general as possible -for example, in defining the sorted invariant mass variables, we shall have in mind the generic diagram in Fig. 2 instead of the specific SUSY-inspired example of Fig. 1. In the case of Fig. 2, a heavy resonance D decays to N indistinguishable visible particles (denoted by solid black lines) and one invisible particle A (denoted by a red dashed line). We first form the set of all possible n-body invariant mass combinations m v i 1 v i 2 ...v in for a given n ∈ [2, N ]. The total number C N n of elements in the set S N n depends on the choice of n and is given by the binomial coefficient (1.10) We can now uniquely and unambiguously label the members of the set S N n by defining sorted 7 invariant mass variables in analogy to Eqs. (1.5-1.7): (1.11) Using Eq. (1.10), it is easy to see that for a given N , there are a total of N n=2 such variables. The physical meaning of the variable m (n,r) is that it is the r-th largest among all possible n-body invariant mass combinations of visible particles in Fig. 2. From the general definition (1.11) it is easy to make contact with the previously considered N = 3 case of Fig. 1: in the notation of (1.11), the variable set (1.8) corresponds to m (2,1) , m (2,2) , m (2,3) , m (3,1) .
In introducing the more general variables (1.11) we are motivated by several factors: 1. Often the visible particles in the cascade decay are indistinguishable, e.g. they are all QCD jets. This can easily be the case even with the SUSY example of Fig. 1, whenever the second-lightest neutralino (particle C) decays predominantly hadronically to 2 jets and the lightest neutralino (particle A). Another well motivated SUSY example is a squark-gluino-neutralino decay chain where again all three visible particles are jets. In such scenarios, a priori there is no way to single out any particular jet, and the set (1.13) is the only one which makes physics sense.
2. In a purely off-shell scenario, where particle D decays directly to particle A plus N visible particles, it is not possible to assign any specific order to the decay products, even if they are distinguishable experimentally.
3. Even when the decay of particle D proceeds through intermediate narrow resonances, so that the visible decay products are emitted in some well-defined order, this true order is unknown to the experimentalist, and can only be hypothesized. In general, alternative theory models, with alternative orderings of the same final state particles, are also possible. Therefore, assuming a specific ordering throughout the analysis is dangerous and may lead to wrong conclusions.
4. Finally, by staying clear of any theory bias, and considering the most general case of Fig. 2, we will be able to derive (see Section 2 below) the necessary relations which must be obeyed by the kinematic endpoints of the variables (1.11) in the case of a pure off-shell decay (i.e., with no intermediate resonances). Any observed deviations from those predictions will signal the presence of new particles in addition to the mother D and daughter A. Furthermore, as illustrated in Fig. 26 below, the measured relations among the kinematic endpoints are indicative of the particular on-shell event topology at hand.
In this paper we begin the investigation of the mathematical properties of the variables (1.11), and, in particular, their kinematic endpoints. In Section 2 we focus on the general case of an off-shell 1 → N + 1 decay as depicted in Fig. 2. We derive the formulas for the kinematic endpoints of the sorted invariant mass variables (1.11) in terms of the relevant physical mass parameter, the mass difference m D − m A . For N > 2, the number of variables given by (1.12) already exceeds the number of input parameters, which implies certain specific relations among the measured kinematic endpoints. Among the main results of Section 2 is the derivation of these relations, as they provide a stringent test of the offshellness of the decay topology. 8 Having dealt with the general case of arbitrary N in Section 2, in the following sections we return to the SUSY-motivated case of N = 3. We shall similarly study the dependence of the sorted invariant mass endpoints on the physical mass parameters, in the presence of intermediate on-shell resonances. The relevant special cases are discussed in Sections 3, 4 and 5. Section 6 is reserved for our summary and conclusions.
2 The pure off-shell case of (N + 1)-body decay Consider the decay of a heavy resonance D into N massless visible particles and one massive invisible particle A, as shown in Fig. 2: (2.1) In this section we shall assume that the decay (2.1) proceeds in one step, i.e., exactly as depicted in Fig. 2. In other words, any virtual particles hiding behind the circular blob in Fig. 2 are sufficiently heavy and can be integrated out to give rise to an effective contact interaction as shown in the figure.
As mentioned in the Introduction, we treat all visible particles in the final state as indistinguishable, so that we do not know a priori the sequence in which the visible particles are emitted. This motivates us to consider the sorted invariant mass variables m (n,r) defined in Eq. (1.11). In what follows, we shall refer to the first index n as the "order" of the variable, while the second index r will denote its "rank". Obviously, the order n can be chosen to be any integer from 2 to N ; for completeness we shall consider all possible values of n, i.e., we shall construct two-body, three-body, etc. invariant masses of visible particles. For a given order n, the rank r then takes values from 1 to C N n .
Our main goal in this section is to provide the analytic expressions for the kinematic endpoints m max (n,r) in terms of the physical masses m D and m A . In all results below, we shall always factor out the parent mass m D and write the formulas in terms of the dimensionless squared mass ratios 9 where by assumption the particle masses obey the hierarchy In the case of the pure off-shell process (2.1), the only two masses entering the problem are m D and m A , thus our results will be functions of R AD < 1.

Invariant masses of order n = 2
We first discuss the sorted invariant masses of order 2 (i.e., the two-body invariant masses), for which a useful sum rule can be derived as follows. Using 4-momentum conservation for the reaction (2.1) 10 we can write Since the visible particles are assumed massless, p 2 i = 0, and furthermore, 2p i · p j is simply (p i + p j ) 2 = m 2 ij so that the above relation can be rewritten as The left-hand side in this equation may be interpreted as the "total available invariant mass squared" which is allocated to m 2 (2,r) 's in a given event. Since Eq. (2.6) is Lorentzinvariant, we can evaluate its left-hand side in any frame. It is convenient to choose the rest frame of particle D, where We are interested in kinematic endpoints, i.e., the cases in which a particular variable m (2,r) is maximized. Eq. (2.6) suggests that in order to maximize an individual variable 9 Note the following transitive and inversion properties  m (2,r) , we must necessarily maximize the total quantity (2.7) as well. It is easy to see that (2.7) is maximized when A is produced at rest and its energy E A = m A . We therefore conclude that for an event which yields a kinematic endpoint of m (2,r) for some value of r, the following sum rule holds We are now in position to derive the formulas for the kinematic endpoints m max (2,r) for various r.
Rank r = 1. The largest possible value of m (2,1) is obtained when all other invariant mass combinations are vanishing, i.e., for events in which m (2,1) = m max (2,1) and m (2,r) = 0, for r = 2, . . . , C N 2 . (2.9) The momentum configuration of such an event is shown in Fig. 3 -two of the visible particles are exactly back-to-back, while the remaining visible particles, together with the invisible particle A, are all at rest. The m max (2,1) endpoint is now obtained by substituting (2.9) into (2.8): Rank r = 2. According to Eq. (2.8), in order to maximize m (2,2) , we need to minimize both m (2,1) and m (2,r) for r ≥ 3. However, by definition m (2,1) cannot be less than m (2,2) , thus for events giving the largest possible value of m (2,2) we expect to have m (2,1) = m (2,2) = m max (2,2) and m (2,r) = 0, for r = 3, . . . , C N 2 . (2.11) The momentum configurations of such events are also collinear, as shown in Fig. 3. Now there are three visible particles with non-zero momenta, two of them having equal momenta and recoiling against the third. From (2.11) and (2.8) we obtain Higher ranks (r > 2). Proceeding analogously, one might naïvely expect that for an arbitrary rank r, Eqs. (2.10) and (2.12) generalize to m max (2,r) However, one has to be careful to check whether physical momentum configurations exist such that This check is not aways trivial. As a concrete example, consider the rank 10 variable m (2,10) in the case of N = 5 total visible particles. It is not difficult to see that in three spatial dimensions, there are no allowed momentum configurations which would give the required case with m (2,1) = m (2,2) = · · · = m (2,10) . As a result, in this case the conjectured endpoint m max (2,10) will not be saturated, and the true endpoint will appear at slightly lower values of m (2,10) . However, in such situations where the general formula (2.13) happens to overestimate the kinematic endpoint, it is nevertheless possible to obtain the correct answer by inspecting the candidate momentum configurations for the visible particles in the rest frame of the mother particle D.
Let us illustrate the procedure with the above example of N = 5 visible particles, and try to obtain the exact upper bound for m (2,10) . For this purpose, we look for momentum configurations, in which the visible particles v i are maximally "spread out" in the rest frame of the mother particle D, while the massive invisible particle A is still at rest. The main idea is that the desired configurations will exhibit a certain level of symmetry, as demonstrated in Figs. 4 and 5. The left panels in the figures depict the momentum configurations of the visible particles in the D rest frame, whereby momenta in blue (red) have the same magnitude a (b). Energy and momentum conservation imply certain relations among the magnitudes and the relative orientation of the momentum vectors as shown. In each case, the remaining sole degree of freedom can be varied in order to find the maximum of the ranked variable m (2,10) as As expected, this bound is tighter than the naïve expectation (2.15).  Figure 5. The same as Fig. 4, but for an alternative, pyramid-like, momentum configuration. Energy conservation implies 4b + a = m D − m A , while momentum conservation along z requires 4b cos θ = a, leaving only one independent degree of freedom, which in the right panel is chosen as cos θ. In this configuration, we have One can similarly analyze the case of N = 6, where we find three competing momentum configurations: a pentagonal pyramid, a square bipyramid, and a triangular prism wedge. The latter provides the true maximum value of the m (2,15) ranked variable: while the bipyramid-like configuration provides the maximum of m (2,14) : In summary, for large enough ranks r (as in the examples just considered), the upper bound provided by Eq. (2.13) will not be saturated, and the kinematic endpoint m max (2,r) will be found at somewhat lower values. The lowest value of r at which the prediction (2.13) begins to deviate from the true answer, in general depends on the number of visible particles N . We have checked that Eq. (2.13) can be trusted up to the following rank where k is a non-zero integer. For ranks r higher than the rank r given by Eq. (2.19), the expression (2.13) provides simply an upper bound on the kinematic endpoint m max (2,r) .

Invariant masses of order n > 2
Once we consider more than two particles at a time, the situation becomes more complicated. As a concrete example, let us take N = 4 visible particles and investigate the third order (n = 3) variables m 123 , m 124 , m 134 , and m 234 . Momentum conservation (2.4) now leads to the following relation (compare to Eqs. (2.5) and (2.6)) As before, the kinematic endpoints are attained when particle A is produced at rest in the rest frame of particle D, so that the above equation reduces to the following analogue of Eq. (2.8) Rank r = 1. There are two types of events which determine the endpoint of m (3,1) . The first type of events have two visible particles moving back-to-back in the D rest frame, while the other two visible particles are at rest. In this configuration, we have In the other configuration, three visible particles with equal energies are moving in a plane at 120 • with respect to each other, while the fourth visible particle is at rest. This in turn implies that Rank r = 2. The maximal value for m (3,2) is obtained for the momentum configuration given by (2.22), thus the endpoint is the same as (2.24): Ranks r = 3 and r = 4. For the third-and fourth-ranked three-body invariant masses, we apply the reasoning from Section 2.1 to obtain similar expressions. Our final answer for the case of N = 4 visible particles is thus (2.26) Having worked out this simple example, we can now generalize (2.26) to higher orders n > 3. Suppose that there are N visible particles as usual, and we consider invariant mass combinations of order n. Fixing two visible particles, say i and j, the term 2p i · p j appears in C N −2 n−2 invariant mass variables out of the total possible number C N n . We have already seen that summing 2p i · p j over all possible pairs of indices i and j is related to a sum over all possible n-body invariant mass combinations 11 : The kinematic endpoints that we are interested in are obtained when the right-hand-side of this relation is maximized (by virtue of particle A being at rest in the D rest frame): Retracing the steps which led to Eq. (2.26), we get (2.29) As already discussed at the end of Section 2.1, this formula provides the exact maximum only up to some rank, i.e., for sufficiently high ranks, it only gives an upper bound. However, even with such high ranks, the true endpoint will still be proportional to the , only with a pre-factor which is somewhat smaller than 11 See, e.g., the special case in Eq. (2.20).

Testing for off-shellness
Armed with the general result (2.29), one can now design a specific test to verify that the decay topology is indeed a purely off-shell one as hypothesized in Fig. 2. The main observation is that in a purely off-shell topology the kinematic endpoints of all invariant mass variables are functions of a single degree of freedom, m D − m A . This implies certain relationships, or "sum rules" for short, among the kinematic endpoints. These sum rules are quantitatively predicted by Eq. (2.29). Note that by introducing the sorted variables (1.11), we are considering the largest possible number of invariant mass variables, and therefore we obtain the largest possible number of sum rules. 12 For illustration, let us consider the simplest case of N = 3 visible particles, as in the SUSY-like decay chain of Fig. 1. There are 2 3 − 3 − 1 = 4 sorted variables given by (1.13), and one unknown degree of freedom, m D − m A , which leaves us with three sum rules. Therefore, if this were a purely off-shell process, the following relations must hold The violation of one or more of these relations would indicate one of two things -either the presence of intermediate on-shell resonances, or some sort of momentum-dependent couplings (form-factors).
The relationships (2.30-2.32) are illustrated in Fig. 6. We consider a four-body pure off-shell decay, i.e., the reaction (2.1) with N = 3, and plot the distributions of the four relevant sorted invariant mass variables: m (2,1) (red dashed line), m (2,2) (green dotted line), m (2,3) (blue dot-dashed line) and m (3,1) (black solid line). Eq. (2.29) predicts that their respective endpoints will be located at 450 GeV, 318 GeV, 260 GeV and 450 GeV. Fig. 6 shows that the endpoint structure for m (3,1) can be identified very well and the value of the endpoint clearly agrees with the theoretical prediction. The two-body invariant mass distributions for m (2,i) also saturate the theoretical bounds (2.29). However, we note that those distributions are relatively shallow near their upper kinematic endpoints [25][26][27], which might make the experimental extraction of those endpoints rather challenging [28]. 13 3 The pure on-shell decay topology with N = 3 visible particles Having considered the purely off-shell case in complete generality in the previous section, we now turn our attention to decay topologies with intermediate on-shell resonances. Un- 12 Recall from Eq. (1.12) that the number of sorted variables is 2 N − N − 1. Therefore the total number of sum rules in the purely off-shell case is 2 N − N − 2.
13 In principle, one should be able to benefit from the knowledge of the asymptotic behavior of the distribution near the endpoint, which, however, is only known for the cases of m (2,1) and m (3,1) [26].  fortunately, the general analysis for an arbitrary number of visible particles N gets quite complicated, which is why in this and the subsequent sections we shall limit our discussion to the most interesting case of N = 3, as in the SUSY-like decay chain from Fig. 1. In particular, we shall focus on the four possibilities depicted in Fig. 7. The event topology of Fig. 7(a) is simply a special case of a purely off-shell decay already considered in the previous section. The event topology of Fig. 7(b) is the typical SUSY scenario from Fig. 1 and will be the main subject of this section. It involves a sequence of three 2-body decays, where each decay produces one visible particle. We shall sometimes refer to the diagram of Fig. 7(b) as a (1, 1, 1 Figure 8. (a) The (1, 1, 1) decay topology of a heavy particle D into a lighter particle A and three massless visible particles v 1 , v 2 and v 3 . (b) The relevant kinematic variables: α is the angle between the momenta of v 1 and v 2 in the rest frame of particle C, θ is the polar angle of the v 3 momentum with respect to the direction of C in the B rest frame, and φ is the azimuthal angle between the planes defined by the momenta of a) v 1 and v 2 and b) v 2 and v 3 , in the C rest frame.
"hybrid" event topologies in the sense that they involve both a two-body and a three-body decay. Correspondingly, the diagram of Fig. 7(c) will be referred to as a (2, 1) topology and will be studied in the next Section 4, while the diagram of Fig. 7(d) will be labelled as a (1, 2) topology and will be considered in Section 5.

3.1
The phase space structure of the (1, 1, 1) decay topology Before discussing the properties of the ranked variables (1.11), it is instructive to review the properties of the allowed phase space in terms of the original variables (1.1) [29,30]. The relevant kinematic variables for the (1, 1, 1) decay of a heavy particle D into an invisible particle A and visible particles v 1 , v 2 and v 3 are depicted in Fig. 8. We note that the (1, 1, 1) decay is most conveniently described in the rest frame of particle C as shown in Fig. 8(b). Naïvely, the total number of degrees of freedom is four, but one of them (here the overall azimuthal angle β) can be neglected taking into account the azimuthal symmetry of this phase space. The remaining three degrees of freedom are: the polar angle α of the momentum of v 2 with respect to the direction of v 1 , the polar angle θ of the momentum of v 3 with respect to the direction of C in the rest frame of particle B, and the azimuthal angle φ between the planes defined by (v 1 , v 2 ) and (v 2 , v 3 ) in the C rest frame. This phase space can be equivalently described in terms of the invariant mass variables However, to simplify notation, from here on we shall work with the dimensionless variables instead of the dimensionful set (3.1). In terms of the SUSY-like decay chain of Fig. 1, the variable x corresponds to m 2 , the variable y is the analogue of m 2 j n , while z represents m 2 j f . The three angular degrees of freedom θ, α and φ from Fig. 8(b) can be mapped to the three dimensionless invariant mass variables (3.2-3.4) as follows: and The allowed phase space spanned by (3.5-3.7) is shown in Fig. 9. Its shape has been likened to that of a samosa [30] and can be parametrized by two functions, z + (x, y) for the top surface (colored in red) and z − (x, y) for the bottom surface (colored in blue). In order to find the explicit form of z ± (x, y), we note that the angle φ enters only the definition of z in Eq. (3.4). Thus the extreme values of z (for a fixed x and y) are found for the extreme values of φ, namely, φ = 0 and φ = π [29]: The exact shape of the "samosa" is determined by the location of the four "corner points" which in Fig. 9 are denoted as {O, P 1 , P 2 , P 3 }

Computer tomography of the allowed phase space
In order to analyze the allowed phase space from Fig. 9 in terms of the sorted variables (1.11), we need to rank the variables x, y and z among themselves. We do this by performing a computerized tomography (CT) scan in which the relevant cross sectional CT images are obtained at a fixed value z = z s along the z-axis (see the green plane in Fig. 9). In this CT scan process the point P 1 plays a special role because it divides the obtained images into two groups. Whenever the scan image is taken "below" P 1 , i.e., at a value of z s smaller than the z component of P 1 the boundary of the image consists of two segments obtained from setting z + (x, y) = z s , interspersed with another two segments given by z − (x, y) = z s (see Fig. 9). On the other hand, when the image is taken "above" P 1 , i.e., when z s > z c , the corresponding image boundary is made up of only one segment from each surface (top and bottom). The basic procedure of ranking x, y and z with the CT scan method is the following. For a fixed z = z s , the intersection of the green plane shown in Fig. 9 with the interior of the samosa determines the allowed values x(z s ) and y(z s ) at this particular value of z = z s . We then sort x(z s ) and y(z s ) and find their respective maxima: Once this is done, we need to compare the thus obtained values of r 1 and r 2 to the value of z, so we sort again {r 1 , r 2 , z s } by magnitude for all possible r 1 and r 2 to obtain the "local" In evaluating r 2 , it is convenient to fold the (x, y) plane along the x = y line [14]. This motivates us to treat separately the cases of x 0 < y 0 and x 0 ≥ y 0 .

CT scans for the case of x 0 < y 0
We first discuss the case of x 0 < y 0 , i.e., when the range of possible y values is larger than the range of possible x values. According to Eqs. (3.8) and (3.9), this occurs for mass spectra obeying the relation At those values of z s , one of the r i variables, either r 1 or r 2 , when considered as a function of z s , exhibits some interesting behavior. This is depicted more clearly in Fig. 11(a), where we track the functional dependence of r 1 (z s ) and r 2 (z s ). At first, for very small values of z s (panels (a) and (b) in Fig. 10), both r 1 (z s ) and r 2 (z s ) stay constant at r 1 (z s ) = y 0 and r 2 (z s ) = x 0 15 . As the value of z s is being increased, the blue point marking the location of r 2 is lowered until it eventually reaches the diagonal line of x = y. This occurs at a special value of z s ≡ z 1 such that As we continue to increase z s beyond z 1 and up to z c (panels (c) and (d) in Fig. 10), two effects take place. First, the value of r 2 is not given by x 0 any more, but is obtained from 14 We remind the reader that we are using the notation maxr to indicate the r-th largest among a given set of elements. The index r in Eq. (3.20) is thus the rank index which in this case takes values r = 1, 2, 3. 15 Recall that in this subsection we have assumed that y0 > x0 Figure 10. For the case of x 0 < y 0 , we show the eight characteristic CT images obtained at different fixed values of z s . Each CT image typically consists of a green region, in which x < y, and a grey region, in which x > y. In order to rank x and y, we fold each CT image along the diagonal line x = y, thus mapping each grey region onto a corresponding green hatched region. The extremal values of r 1 and r 2 are then found by examining the two green regions (with and without a grid hatch). In each panel, the red (blue) dot indicates the location of the point giving the value of r 1 (r 2 ). the folding along the x = y line. The functional dependence r 2 (z s ) is thus given implicitly by the equation z + (r 2 , r 2 ) = z s .

(3.24)
Second, the red point in Fig. 10(a-d) indicating the value of r 1 moves to the left, until it eventually reaches the point (x, y) = (0, y 0 ), where the z coordinate is given by Comparing to (3.15), we see that the scan depicted in Fig. 10(d) is taken through point P 1 in Fig. 9, and from that point on, the value of r 1 will begin to decrease with z s . Indeed, Figure 11. as the value of z s increases further from z c to z 2 (panels (e) and (f) in Fig. 10), the green region shrinks and r 1 (z s ) decreases linearly with z s as until r 1 (z s ) reaches the value x 0 . This occurs at a value of z s = z 2 such that r 1 (z 2 ) = x 0 . Inverting (3.26) and solving for z 2 , we obtain Finally, for z s > z 2 , r 1 (z s ) again becomes constant at r 1 = x 0 . Thus the complete functional behavior of r 1 (z s ) is given by x 0 for z 2 < z s ≤ z 0 , as illustrated by the red solid line in Fig. 11(a). In order to complete the discussion of Figs. 10 and 11, we again turn our focus to r 2 (z s ), which was given implicitly by (3.24). However, this is true only as long as the CT images are crossed by the diagonal line x = y. Eventually, as z s approaches its maximal value z 0 , the CT image is confined to the region 16 with x ∼ x 0 and y ∼ 0 and may not 16 Recall that the "tip" of the samosa is located at point P2 in Fig. 9, whose coordinates are given by Eq. (3.16). extend all the way up to the x = y line, as depicted in Fig. 10(h). As shown in Fig. 10(g), the image becomes "disconnected" from the diagonal line x = y at the point From that point on, for z s > z t , the value of r 2 is determined by the rightmost point of the hatched green CT image, i.e., the blue point in Fig. 10(h). Thus we find that for z s > z t , the functional dependence of r 2 (z s ) is given by Interestingly, this is the same function as (3.26). This can also be seen by inspecting Fig. 11(a), where the blue and red slanted straight segments are lined up.
In conclusion, we comment on the relative ordering of the special points {z 1 , z c , z 2 , z t }. First, the definition (3.29) can be rewritten as which makes it evident that z t > z 2 . Furthermore, the definition (3.27) can be rewritten as Since e −2η = R BC ≤ 1 by definition and x 0 < y 0 by the assumption of this subsection, Eqs. (3.32) imply that z 2 is larger than both z 1 and z c . Finally, the relative size of z c and z 1 is not predetermined, but depends on the mass spectrum. Therefore, the ordering is We are now in position to derive the endpoints of the ranked invariant mass variables. We first add the straight line z = z s as the green line in Fig. 11(a), and proceed to order {r 1 (z s ), r 2 (z s ), z s } for each value of z s as shown in Fig. 11(b), where the red, blue, and green colored curves track the values of the largest, the second largest, and the smallest invariant mass combination for each z s . Therefore, the endpoint of each ranked invariant mass arises at the maximum of each colored curve. Since r 1 and r 2 have already been ordered among themselves in accordance with (3.19), all that is left to do now is to compare r 1 (z s ) and r 2 (z s ) to the value of z s itself. This leads to two interesting intersection points seen in Fig. 11(a): first, the intersection point between r 1 (z s ) and the straight line z s Figure 12. The same as Fig. 10, but for the case of x 0 > y 0 . and the corresponding intersection point between r 2 (z s ) and the straight line z s Armed with these results, it is now straightforward to determine the kinematic endpoints of the sorted invariant mass variables on a case-by-case basis. We postpone the presentation of the relevant results until Section 3.3, after we have had the chance to also discuss the case of x 0 > y 0 , which is the subject of the next subsection.

CT scans for the case of x 0 > y 0
We now repeat the analysis of the previous subsection for the case of x 0 > y 0 . Note that the vertices P 2 and P 3 from Fig. 9 have a common x coordinate x = x 0 (see also Eqs. (3.16) and (3.17)). Therefore, for any value of z s , we will have the simple relation (contrast this to (3.28)) This fact can be clearly observed in Figs. 12 and 13, which are the x 0 > y 0 analogues of Figs. 10 and 11, respectively. On the other hand, the behavior of the function r 2 (z s ) is similar to the x 0 < y 0 case from the previous subsection, only now the roles of x and y are reversed. Again, there are two special points, z s = z 3 and z s = z t , where the functional form of r 2 (z s ) changes. Fig. 12(a,b) reveals that as we start increasing z s from 0, the value of r 2 stays constant at r 2 (z s ) = y 0 . Eventually the blue point determining the value of r 2 reaches the diagonal line x = y. This occurs at a special value of z s = z 3 given by Once z s exceeds z 3 , the blue point begins to track the straight line of x = y and r 2 (z s ) is again found from (3.24). As illustrated in Fig. 12(c-f), this trend continues until the contour of z + (x, y) = z s is detached from the x = y line at the point z s = z t , where z t z max (a) Figure 13. The same as Fig. 11, but for the case x 0 > y 0 shown in Fig. 12. is again given by (3.29). Finally, for z s > z t , r 2 (z s ) is the maximal y-coordinate of the contour line z + (x, y) = z s which was already computed earlier in Eq. (3.26).
The rest of the steps are identical to those for the previous case of x 0 < y 0 . Having ranked x and y among themselves in terms of r 1 and r 2 , it remains to rank z relative to r 1 and r 2 , as illustrated in Fig. 13(b). The relevant results are summarized in the next subsection.

Results summary for the (1, 1, 1) decay topology
In this subsection, we collect all relevant results derived from the arguments in the previous two subsections. We have seen that the endpoints of the ranked invariant mass variables are given in terms of the endpoints x 0 , y 0 and z 0 of the original unranked variables (3.1), as well as the two intersection points (3.34) and (3.35). The actual mass spectrum {m D , m C , m B , m A } then determines which of these five expressions applies to which ranked variable. 17 Therefore, in presenting the results, we must first describe how the mass parameter space {m D , m C , m B , m A } is partitioned into domains, and then specify the relevant endpoint formulas for each individual domain.
In order to define the domains, we can factor out the overall scale, say m D , and then trade the remaining three mass variables {m C , m B , m A } for the dimensionless ratios (2.2). Even though the resulting parameter space {R AB , R BC , R CD } is three-dimensional, it turns out that all relevant domains can be exhibited on a suitably chosen 2-dimensional slice at constant R AB , as illustrated in Figs. 14 and 15.
The color coding in Figs. 14 and 15 shows that there are 18 different regions which are needed in order to define the endpoints of the sorted two-body invariant mass variables m (2,r) . Sixteen of those regions are always visible, for any value of the fixed parameter R AB . On the other hand, Region X only appears at low values of R AB , R AB < (3 − (see Fig. 14), while Region V only appears at high values of R AB , R AB > (3 − √ 5)/2 (see Fig. 15). In what follows, each region will be defined by specifying a range of R BC values and a corresponding range for R CD values. Given the geometry of Figs. 14 and 15, it is convenient to define the region boundaries at a fixed R AB by treating R BC as the independent variable, and R CD as the dependent variable, as follows: (3.43) The functions f (R BC ) defining those boundaries are labelled by the replacement which needs to be done in the endpoint formulas (3.46) below as one crosses the corresponding boundary and moves from one region to the next. In addition, Figs. 14 and 15 also display two vertical boundaries Using the definitions of the boundaries (3.38-3.45), the 18 colored regions 18 seen in Figs. 14 and 15 can be explicitly defined as in Table 1. Then, the kinematic endpoints for 18 A close inspection of the endpoint formulas (3.46) given below reveals that in Regions IX and XV the endpoints are given by the same expressions, so those two regions can be effectively merged. The same observation applies to Regions XI and XV II. Thus, strictly speaking, for the purposes of defining the endpoints of the ranked m (2,r) variables, one only needs to consider 16 cases.

Region
Color min max min max the sorted two-body invariant mass variables m (2,r) are given by For completeness, we also list the well known endpoint formulas for the three-body invariant mass m (3,1) [12,13] (3.47) 4 Type (2, 1) cascade decay chain In this section, we proceed to analyze one of the hybrid decay topologies, namely type (2, 1). The relevant decay chain is depicted in Fig. 16(a): a massive particle D decays via a threebody decay into two massless visible particles v 1 and v 2 and an on-shell intermediate particle C. In turn, particle C decays into a massless visible particle v 3 and an invisible particle B. The relevant decay configuration seen in the rest frame of particle D is illustrated in Fig. 16(b). Again, naïvely there exist four degrees of freedom, however, out of the two azimuthal angles, β and φ, only their difference is relevant -we can then safely parametrize it with φ, and set β to zero. The allowed {x, y, z} phase space is shown in Fig. 17. In order to obtain the equation for the surface boundary of the allowed region, we start from the kinematic relation where the superscript on E B implies that this energy is measured in the rest frame of particle D. Notice that the symmetry v 1 ↔ v 2 implies that the variables x and z always enter in the combination x + z. Then, the equation for the boundary surface is where the functions Σ ± (y), which are the analogues of (3.12,3.13), are obtained from (4.1) for the extreme values of cos θ: In order to derive the maximal values of the sorted invariant masses of order 2, we repeat the scanning procedure from the previous section, only now we scan along the y-axis. For a given y = y s , the CT image in the (x, z) plane is an isosceles trapezoid, as shown in Fig. 18. We again rank x(y s ) and z(y s ) for all possible pairs of (x, z) to obtain the ranked variables (4.5) Due to the simple geometry of Fig. 18(a,b), r 1 and r 2 can be readily computed as We then compare the variables r 1 (y s ), r 2 (y s ), and y = y s as before, see Fig. 18(c,d). There are two special points, denoted by y 1 and y 2 , which arise from the intersection of r 2 (y s ) and r 1 (y s ) with the straight line y = y s : The kinematic endpoints for the sorted two-body invariant masses m (2,r) will be given in terms of the endpoints x 0 = z 0 and y 0 of the original unsorted variables and the two special points y 1 and y 2 given by (4.8) and (4.9). Fig. 18(d) shows an example where m max (2,1) = m D √ y 0 , m max (2,2) = m D √ y 2 and m max (2,3) = m D √ y 1 . The relevant formulas for the general case are collected in Section 4.1.
In conclusion of this subsection, we discuss the derivation of the endpoint of the m (3,1) variable for the case of a type (2, 1) decay topology. Since in supersymmetry models Figure 19. The allowed phase space (yellow-shaded region) of a Type (2, 1) decay chain in the (y, x + z) plane. The kinematic boundary is given by Σ ± (x,z) (y s ) for a given y s . The two panels show the two distinct cases for m max one typically gets either (1, 1, 1) or (1, 2) decay topologies, this case has not drawn much attention in the existing literature. It is most convenient to analyze the decay kinematics in the plane of (x + z) versus y, as shown in Fig. 19. Since all visible particles are assumed massless, we have the identity which describes a straight line with a slope −1 and intercept k: x + z = −y + k. Given Eq. (4.10), in order to maximize m (3,1) , we need to maximize the intercept, k, with respect to the full allowed phase space of (y, x + z), keeping the slope fixed at −1. As shown in Fig. 19, the boundary of the allowed phase space is given by Σ ± (x,z) (y), and we need to consider two distinct cases, illustrated in panels (a) and (b). Note that Σ + (x,z) (y) is a monotonically decreasing function of y, whose (negative) slope increases in magnitude and approaches −∞ at y = y max . Therefore, if the slope at y = 0 is larger than −1, then Eq. (4.11) will appear as a tangential line on the curve of Σ + (x,z) (y) as shown in Fig. 19(a), and the resulting k will give rise to the maximum m (3,1) . If, however, the slope at y = 0 is already smaller than −1, then no such tangential line is possible, and the line giving the maximum will pass through (y, x + z) = (0, Σ + (x,z) (0)), as shown in Fig. 19(b). A simple algebra results in the following kinematic endpoints: (1 − R BC )(1 − R CD ), otherwise. (4.12) Figure 20. The same as Figs. 14 and 15, but for the Type (2,1) decay topology discussed in Section 4.

Results summary for the (2, 1) decay topology
The endpoints of the m (2,r) and m (3,1) sorted invariant mass variables are given in terms of the nominal endpoints 14) and the intersection points y 1 and y 2 given by (4.8) and (4.9). The endpoint formulas are again piecewise-defined functions, and the boundaries of the defining regions are given by functions R CD = f (R BC ) in analogy to (3.38-3.43): , (4.16) (4.17) The boundaries (4.15-4.17) divide the (R BC , R CD ) parameter space of the (2, 1) decay topology into the 4 color-coded regions shown in Fig. 20. Then, the kinematic endpoints (a) Figure 21. of the ranked two-body invariant mass distributions are given by while the endpoint of the three-body invariant mass is (1 − R BC )(1 − R CD ) otherwise. (4.19) 5 Type (1, 2) cascade decay chain In this section, we analyze the decay topology of type (1, 2), which is depicted in Fig. 21(a). First, a massive particle D decays into a visible particle v 1 along with an on-shell intermediate particle C, and in turn, particle C decays into two visible particles v 2 and v 3 and an invisible particle B via a three-body decay. It is convenient to analyze the kinematics in the rest frame of particle C, as illustrated in Fig. 21(b). As one might expect, most of the analysis leading to the final formulae will be similar to that of the preceding section. The allowed {x, y, z} phase space is illustrated in Fig. 22. Notice the y ↔ z symmetry which follows from the v 2 ↔ v 3 symmetry. The boundary of the allowed region can be derived from a kinematic relation analogous to (4.1): The boundary equation is obtained from here by taking cos θ = ±1. One finds where In order to find the largest two-body invariant masses, we again scan the allowed phase space shown in Fig. 22, this time at fixed values for x = x s . The obtained images are again isosceles trapezoids with bases of length Σ + (y,z) (x s ) and Σ − (y,z) (x s ), as shown in Fig. 23(a,b). We then rank y(x s ) and z(x s ) in analogy to (3.19) and (4.7) Just like the previous case of type (2,1) topology, due to the symmetric structure of the phase space, the corresponding r 1 and r 2 are simple to evaluate -they are given by Σ + (y,z) (x s ) and Σ + (y,z) (x s )/2, respectively. Finally, the ranking procedure among r 1 (x s ), r 2 (x s ) and x s again introduces two intersection points, x 1 and x 2 , which arise from the crossing of r 2 (x s ) with x = x s , and r 1 (x s ) with x = x s , respectively (see Fig. 23(c,d)): The endpoints of the ranked two-body invariant masses m (2,r) will be given in terms of y 0 (= z 0 ), x 0 , y 0 /2, x 1 , or x 2 , depending on the mass spectrum. Fig. 23(c, d) illustrates y y x s x s Finally, we discuss the three-body invariant mass m (3,1) . The decay topology (1, 2) is very common in supersymmetry, e.g., in the decayq → qχ 0 2 of a squarkq to the second-tolightest neutralinoχ 0 2 , which in turn decays by a three-body decay to the lightest neutralinõ χ 0 1 and a couple of jets or leptons:χ 0 2 → qqχ 0 1 orχ 0 2 → + −χ0 1 . The expression for the endpoint m max (3,1) is already known (see, e.g. [31]) and here we shall simply re-derive it using the method from Section 4.
Again, it is convenient to study the variable of interest m 2 (3,1) = x + y + z ≡ k (5.8) in the plane of (y + z) versus x, as depicted in Fig. 24, where the allowed region is shaded x max x s ⌃ + (y,z) (x s ) ⌃ (y,z) (x s ) x + y + z ⌃ + (y,z) (x s ) ⌃ (y,z) (x s ) x + y + z x max x s y + z y + z Figure 24. The same as Fig. 19, but for a type (1, 2) decay topology, where the relevant phase space is best viewed in the (x, y + z) plane.
in yellow. In this plane, the relation (5.8) again represents a straight line with constant negative slope −1 and intercept k. Just like in the previous section, the task is to find the point on the phase space boundary which maximizes the intercept, k, for a fixed slope −1. As the two panels of Fig. 24 show, one again has to consider two cases, depending on whether the slope of the boundary curve Σ + (y,z) (x) at x = 0 is larger or smaller than −1. In the former case, shown in Fig. 24(a), the endpoint m max (3,1) is obtained from a line tangential to the boundary, while in the latter case the endpoint is simply given by y 0 : (1 − R CD )(1 − R BC ) otherwise. (5.10)

Conclusions and outlook
The dark matter problem greatly motivates the search for semi-invisibly decaying resonances in Run II of the LHC. After the discovery of such particles, their masses will most likely have to be measured using the classic kinematic endpoint techniques. In fact, such techniques can already be usefully applied in the current data -for example, following the procedure outlined in [32], the CMS collaboration has published an analysis of simultaneous extraction of the top, W and neutrino masses from the measurement of kinematic endpoints in the tt dilepton system [33].
In this paper, we revisited the classic method for mass determination via kinematic endpoints, where one studies the invariant mass distributions of suitable collections of visible decay products, and extracts their upper kinematic endpoints. We generalized the existing studies on the subject in several ways: • We shied away from making any assumptions about the structure of the decay topology, and considered the invariant masses of all possible sets of visible decay products. This led us to the introduction in Eq. (1.11) of the sorted invariant mass variables m (n,r) , where we consider all possible partitions of n visible particles, and then rank the resulting invariant masses in order. The so defined sorted invariant mass variables allow us to study SUSY-like decay chains in a fully model-independent way.
• In Section 2 we considered a completely general semi-invisible decay with no intermediate resonances, where a heavy particle D decays directly to an arbitrary number N of massless SM particles and a single massive NP particle A. For this very general case, we derived the corresponding formulas for the endpoints of the sorted invariant mass variables, Eq. (2.29). The importance of those results lies in the fact that they allow the experimenter to test for the presence of intermediate on-shell resonances between particles D and A -in the absence of such resonances, the ratios of all endpoints are uniquely predicted by Eq. (2.29). Any measured deviation from those ratios will signal the presence of some other new intermediate particles.