###### Abstract

Various physics observables can be determined from the localisation of
distinct edge-like features in distributions of measurement values. In this
paper, we address the observation that neither differentiating nor fitting the
measured distributions is robust against significant fluctuations in the
experimental data. We propose the application of Finite Impulse Response (FIR)
filters instead. To demonstrate the method, we consider the typical case in
particle physics in which the precise localisation of kinematic edges, often
blurred by e.g. background contributions and detector effects, is crucial for
determining particle masses. We show that even for binned data, typical for
high energy physics, the optimal FIR filter kernel can be approximated by the
first derivative of a Gaussian (FDOG). We study two highly complementary
supersymmetric scenarios that, if realised in nature, could be observed at a
future high-energy $e^{+}e^{-}$ collider such as the International Linear
Collider (ILC) or the Compact Linear Collider (CLIC). The first scenario
considers the production of $\tilde{e}^{\pm}-$pairs while the second focuses
on the $\tilde{\chi}^{\pm}_{1}$ and $\tilde{\chi}^{0}_{2}$-pair production. We
demonstrate that the FIR filter method for edge extraction is superior to
previously employed methods in terms of robustness and precision.

##  1 Introduction

Physics observables are commonly extracted by detecting the position of
differential features in distributions of experimental data. An important
example in particle physics is the measurement of edge positions in kinematic
distributions, which can reveal the kinematic end-points of a process. In an
actual experiment, such an edge position will be subject to noise from
statistical fluctuations of the signal process as well as from remaining
background contributions, and will be blurred e.g. by the finite resolution of
the detector.

Various techniques can be used to extract the sought after kinematic feature
from such distributions. The most commonly used technique is to fit the
distribution with a differentiable function that can model the feature to be
measured such that the relevant parameter value(s) can be extracted from the
fitted function. This technique implies a model choice, i.e. an assumption
regarding the analytical form of the feature to detect. In many cases,
especially when there are detector effects and background processes affecting
the distribution, there is no unique, model-independent choice for an analytic
description. It might be also unfeasible to find any viable analytical
function that can model the distribution such that it stays robust against
small changes of its shape.

In this article we describe a computational method based on finite impulse
response filters. Related techniques are not uncommon in digital image
processing, and are also used elsewhere, e.g. in the automatic calibration
system of the LIGO experiment [1]. Here, we adjust the method to the
application of typical particle physics distributions. As an example, we focus
on the measurement of kinematic end-points of physical processes which are
particularly relevant for the physics program of future electron-positron
colliders. Two complementary study cases have been chosen. The first one
exhibits a large signal with a sharp edge over a small background. Here, we
will especially study the systematic uncertainties of the filter method. The
second example has a much smaller signal-to-background ratio and a shallow
edge, making a robust determination of its position challenging.

This paper is structured as follows: in Sec. 2 we describe and characterize
our edge detection method. The two study cases are introduced in Sec. 3, while
the application of the algorithm to the two benchmarks is discussed and
compared to alternative approaches in Sec. 4 and 5, respectively. We conclude
by summarizing our results in Sec. 6.

##  2 A computational approach to edge detection

In the general case of a continuous function $\bm{g(x)}$, an edge-like feature
is characterised by a local maximum or minimum in the first derivative
$\bm{g^{\prime}(x)}$ and a zero-crossing in the second derivative
$\bm{g^{\prime\prime}(x)}$. Unfortunately, in the case of experimental data,
the physical signal is disturbed by noise and small variations in the noise
amplitude can have large impact on the results of the differential operator
[2]. One possibility to regularise this problem is to convolute the target
function $\bm{g(x)}$ with an appropriate function $\bm{f(x)}$ _before_ the
differentiation stage.

Assuming $\bm{f}$ is differentiable and absolutely integrable, and $\bm{g}$ is
absolutely integrable, the derivative of the convolution (*) of these two
functions, denoted here as the response $\bm{R}$, can be written as

$\displaystyle R(x)=\frac{d}{dx}(f(x)*g(x))=\frac{df(x)}{dx}*g(x)=h(x)*g(x).$
(1)

The edge-like features in $\bm{g(x)}$ can be identified from the maxima of
$\bm{h(x)*g(x)}$ or the zero-crossings of $\bm{h^{\prime}(x)*g(x)}$. In this
context, the function $\bm{df(x)/dx=h(x)}$ is called _a matched filter_ and
can be characterised by the result of its convolution with a delta function.
The convolution result is called the algorithm’s _Impulse Response_ (IR).
Since the filter unambiguously defines its impulse response we will use the
two terms as synonyms. When the filter values drop to 0 everywhere outside a
finite range in its domain, that function is known as a _Finite Impulse
Response_ (FIR) filter, otherwise it is called an _Infinite Impulse Response_
(IIR) filter.

In the analysis of physics data though, the typical observable distributions
are discrete, therefore the filters and the convolution algorithm should be
defined in the same domain. For this discussion, we assume that our input
distributions are defined on a regular grid, i.e. that the distance between
the samples in the distribution domain is constant, and that the error
associated to each sample is independent from all the others. A histogram with
constant binning is a special case of such a distribution. Consequently, the
filter can also be represented by a discrete function on the same regular
grid. Considering the input function $\bm{g_{d}(i)}$ and a filter represented
by a discrete function $\bm{h_{d}(k)}$, the _Discrete Response Function_
$\bm{R_{d}(i)}$, that is their discrete convolution, can be written as:

$\displaystyle R_{d}(i)=\sum_{k=-N}^{N}h_{d}(k)\cdot g_{d}(i-k)$ (2)

where N is the value after which the filter amplitude drops to 0.

In the case of a continuous distribution, the optimal characteristics of a
differential filter were described in Ref. [3], which showed that the optimal
differential filter for a continuous distribution disturbed by random noise,
called the _Canny filter_ , is very similar to the _First Derivative of a
Gaussian_ (FDOG), which will be defined in Sec. 2.1. In the case of a discrete
distribution the results of the optimization can differ substantially, as
described extensively in [4], and different types of filters may be used to
obtain the optimal results. Therefore we will redefine those optimization
criteria for the specific case of the detection of an edge-like feature in a
one-dimensional distribution and develop an operational procedure to evaluate
the algorithm performance to optimize the parameters of its operation in the
following subsections.

Using those ideas, given a discrete physical distribution defined on a regular
grid, e.g., a histogram with uniform binning, and a filter defined by its
$2N+1$ coefficients $\bm{h_{d}(k)}$, we can define the position of an edge-
like feature as the point where the input function’s first derivative is
maximum. Applying our algorithm we can compute this position calculating the
position of the maximum of the response function $\bm{R_{d}(i)}$. A simple
example of a noisy step function as test function $g_{d}$ and the response of
a FDOG filter $R_{d}$ is shown in Fig. 1.

Figure 1: Example of the application of the edge finding algorithm to a noisy
step function with a signal-to-noise ratio of 1.0. The test function to which
the algorithm is applied is obtained by overlaying a white-noise of $\sigma=1$
to an ideal step function of amplitude $A=1$ at $i=0$. We processed the test
function with a FDOG filter with $\sigma_{f}=5$. Using our algorithm the
maximum positive response which identifies the step position was detected at
$x=-0.4$

This technique is common in data processing. It is employed, e.g., for finding
edges in digital images or for identifying a significant pulse in a noisy
background when performing signal processing. In those contexts, the emphasis
is to find an appropriate solution in large number of independent datasets,
e.g., to identify the edges in millions of different images. In contrast, our
study aims at defining the optimal filter to measure, as precisely as
possible, the position of a specific feature, i.e. the edge, in a single
distribution, taking into account the measurement errors of all samples.
Furthermore, it is also crucial to determine a robust and reliable estimate
for the uncertainty associated with the edge’s position measurement.
Therefore, the most important element of our study is the establishment and
application of an optimization and uncertainty estimation procedure.

###  2.1 Definition of Data Properties and Filter Types

Figure 2: Definition of the important parameters used to characterize the
edges. L represents the edge position, A the amplitude, W the width.

The data distributions to be analysed for the presence of an edge can be
characterized by the following properties, visualized in Fig. 2:

Location of the edge (L):

The location $L$ of the edge is the main target observable. For an analytical
function, this corresponds to the location of the maximum of the first
derivative of the function.

Amplitude of the edge (A):

The amplitude $A$ of the edge is the total height of the edge and a measure of
the signal strength.

Width of the edge (W):

The edge width $W$ which, in this work, we will define as the range in the
function domain, or, in case of a histogram, the number of bins over which the
distribution rises from $10$% to $90$% (or drops from $90$% to $10$%) of the
edge amplitude $A$. For an ideal step function, $W=1$. When the step is
defined by an analytical function, the width can also be calculated
analytically and be defined as a real number.

Signal-to-Noise Ratio (SNR):

In case additional white noise is overlaid to the ideal function, we define
the _Signal-to-Noise Ratio_ as the ratio between the maximum of the derivative
(S) of the ideal function and the noise RMS.

In order to compare the performance of different filters, it is convenient to
define a standard _Scale Factor_ , computed as the RMS of the _Step Response_
, that is the filter response to an ideal, noiseless step function with
amplitude $A=1$, width $W=1$ and the edge located at $k=0$.
In the following, we will compare four different convolution filters:

Difference of Boxes (DOB):

This filter is equivalent to the execution of a moving average in a window of
size $2N+1$ prior to the derivation stage. The characterizing parameter of
this filter is the size $\bm{N}$ of the averaging window and can be
analytically written as:

$\displaystyle h_{d}(k)=1$ for $\displaystyle-N\leq k<0$ $\displaystyle
h_{d}(k)=0$ for $\displaystyle k=0$ $\displaystyle h_{d}(k)=-1$ for
$\displaystyle 0<k\leq N$ (3)

First Derivative of a Gaussian (FDOG):

This filter is equivalent to the execution of a Gaussian smoothing prior to
the derivative. The filter is parametrised by the $\bm{\sigma}$ parameter of
the Gaussian:

$h_{d}(k)=-k\ e^{-\frac{k^{2}}{2\sigma^{2}}}\quad\text{for}\quad-N\leq k\leq
N$ (4)

In case of the FDOG, the scale factor is equal to the width $\sigma$ of the
Gaussian.

Deriche:

This algorithm, described in [5], is characterized by the following impulse
response, parametrised by the exponential decay length $\bm{\alpha}$:

$h_{d}(k)=-k\ e^{-\alpha|k|}\quad\text{for}\quad-N\leq k\leq N$ (5)

Shen:

Similar to the previous one, this algorithm, described in [6], is equivalent
to an exponential smoothing and presents a discontinuity in the derivative at
$k=0$:

$\displaystyle h_{d}(k)=e^{\alpha k}$ for $\displaystyle-N\leq k<0$
$\displaystyle h_{d}(k)=0$ for $\displaystyle k=0$ $\displaystyle
h_{d}(k)=-e^{-\alpha k}$ for $\displaystyle 0<k\leq N$ (6)

Additionally, we are going to compare the results of those filters to the
numerical derivative, labeled in all figures as “derivative”, equivalent, up
to a constant, to a DOB filter with ${N=1}$. A detailed and extensive
description of the characterization of these filters is given in [7] and we
will here summarize the most relevant findings.

Figure 3: Normalized step response function (SR) of the 4 studied filters, as
compared with the response of the numerical derivative. The characterizing
parameter of each filter was chosen so that they will have a scale factor (RMS
of the SR) of about 5 to make a comparison more meaningful.

Those filters are graphically visualized through their _Step Response_ (SR) in
Fig. 3. Thereby the characterizing parameters of all filters have been
adjusted such that the scale factor is about 5 for each filter. The value of
the scale factor has an impact on the amount of noise that the filter picks up
during the convolution with the input distributions as well as on the
precision of the edge localization. Qualitatively, we expected that when that
value is too small, the filter response will be strongly influenced by the
noise. In contrast, if it is too large, the correct edge will be found above
the noise, but its position determination may be less precise and be
influenced by the presence of other nearby features. We will quantify these
aspects in the next sections.

###  2.2 Optimisation criteria for a discrete function edge finder

Using the tools introduced in the previous subsection we define two benchmark
quantities: the filter efficiency and its localization power. Those are
conceptually similar to those defined in [3] and we used them to characterize
and optimize the filtering algorithm. In all benchmarks presented in this
section we use, as a test function, a Gaussian smoothed step of fixed
amplitude with variable width and SNR. Using this function allows to define
the width and the signal amplitude analytically as a function of the width
parameter of the Gaussian smoothing $\bm{\Sigma}$ as: $\bm{W=2\sqrt{2\ln
10}\Sigma}$ and $\bm{S=A/\sqrt{2\pi}\Sigma}$. The Monte Carlo samples for the
optimization criteria evaluation are generated using the aforementioned smooth
step as a template, adding randomly generated white noise to each point in the
function domain. The position of the grid points is randomized at every
sampling to avoid any bias coming from the relative position of the function
edge within a single subdivision of the grid.

#### Filter efficiency:

The filter efficiency, in the case of a rising edge, is defined as _the
probability that the local maximum of the filter’s response closest to the
true edge position of the input function is the global maximum_. In case of a
falling edge it is sufficient to replace minimum with maximum in the previous
definition. This probability can be calculated by processing the input
function with different random noise samples an arbitrary number of times.
This benchmark quantity was evaluated for the four studied filters and, for
comparison, for the simple numerical derivative as a function of the SNR of
the input distribution. As shown in fig. 4, the efficiency of all filters
increases with the SNR and increases monotonically with the filter scale
factor at all SNR. At low SNR in particular, all studied filters perform
significantly better than the simple derivative, which is more sensitive to
noise than any of the filters.

Figure 4: On the left, the plot shows the efficiency of four different
filters, compared to that of the simple derivative, as function of the Signal-
to-Noise-Ratio (SNR) of a noisy, smooth step function with a width $W=1$. Each
point is obtained using 10000 samples of the same function with random white
noise. The characterizing parameter of each filter was chosen so that they
have a scale factor of about 3. On the right, a similar plot shows the
efficiency of the same filter with changing scale factor.

#### Localization power:

The second parameter, the localization power, is particularly important for
this work as it defines the systematic and statistical uncertainty that can be
associated to the results we obtain. It comprises the _localization bias_ and
the _localization error_ , which are defined as the mean $\bm{\hat{L}}$ and
standard deviation $\bm{\sigma_{L}}$, respectively, of the detected edge
positions obtained when applying the algorithm to many random variations of
the same underlying function. In general, the localization is defined in unit
of the grid granularity, e.g. the width of an histogram bin. The difference
between $\hat{L}$ and the true edge position defines the systematic
uncertainty of the filter output, while $\sigma_{L}$ corresponds to the
statistical uncertainty, whose evaluation does not require any prior
definition of a _true_ edge position.

In our test cases, where the edge is defined on top of a constant baseline and
the additive noise is uniform and doesn’t cause an asymmetric distortion of
the edge, the localization bias $\hat{L}$ is always compatible with zero in
all conditions. This will change in the physics examples discussed in Sec. 4
and 5, which include non-flat backgrounds and non-uniform bin errors. The
localization error $\sigma_{L}$, on the other hand, has important variations
depending on the SNR, on the edge width $W$ and the scale factor of the
filter, which are illustrated in fig. 5 and 6, respectively.

In relation to the noise level we observe that, at large enough SNR, all
filters of similar scale factor are roughly equivalent, with a localization
error decreasing asymptotically towards a value of about 30% of the grid
granularity. Moreover we can note that all are notably superior to a simple
numerical derivative. We also note, looking at the right panel of fig. 5, that
there is an optimal filter configuration that provides the best localization
and that this configuration can be chosen unambiguously using our optimization
technique.

To transition to the discussion of the case studies, in which the function to
analyse is obtained from a histogram with user-defined binning, the plots in
fig. 6 show the ratio between the localization error and the step width, which
we will call _normalized localization_ , rather than the simple localization
error. Moreover the SNR overlaid on the function is proportional to
$1/\sqrt{W}$, starting from an SNR of 3 for $W=2$. In this way we can emulate
the case of an histogram where changing the bin width leads to a correlated
change of the signal and background entries and thus the SNR. Here we notice
again that all filters show a similar behaviour: the normalized localization
improves rapidly with the increasing step width for small values of that
quantity and tends asymptotically to a value that depends on the choice of the
filter and on its scale factor.

Figure 5: On the left, the plot shows the localization error of the four
different filters, compared to that of the simple derivative, as function of
the Signal-to-Noise-Ratio (SNR) of a noisy, smooth step function with a width
$W=1$. Each point is obtained using 10000 samples of the same function with
random white noise. The characterizing parameter of each filter was chosen so
that they have a scale factor of about 3. Note that this scale factor was
chosen for ease of comparison between all plots rather than to select the best
possible filter parametrization. On the right, a similar plot shows the
localization error of the same filter type when changing the scale factor.

Figure 6: On the left, the plot shows the normalized localization error of the
four different filters, as a function of the width of a smooth step of varying
width and Signal-to-Noise-Ratio (SNR), as described in the text. Each point is
obtained using 10000 samples of the same function with random white noise. The
characterizing parameter of each filter was chosen so that they have a scale
factor of about 3. On the right the same is shown for a set of FDOG filters
with different scale factors. The derivative is not shown in those plots as
its localization is so poor to be out of scale.

###  2.3 Filter choice and optimization

With the results shown in the previous sections we can define a procedure to
perform an accurate and optimal measurement of the position of edge-like
features in discrete one-dimensional distributions like those obtained in
particle physics experiments and a robust way to evaluate the statistical and
systematic error to associate to such measurements.

Previously we have shown that the efficiency and localization power of all
filters drop significantly when the SNR is reduced below a threshold that
depends on the filter scale factor and type. Above that threshold larger
filters tend to perform better but improvements quickly become negligible.

As we have shown in the previous section, the algorithm parameters can be
optimized through a Monte-Carlo technique, by considering the input
distribution a template from which an arbitrary number of samples can ge
generated. It is important to note that this whole optimization can be
performed without any knowledge of the specific shape of the feature to detect
or of its theoretical position, if that could be calculated at all. In fact,
this optimization does not require any prior knowledge of the underlying
distribution.

When the grid resolution can be customized by the user, like in the case of an
histogram generated from a large dataset, the first step is thus to optimize
that resolution, e.g. the bin width of the input histogram, such that the
average difference between the value at neighbouring positions in the step
region is in the order of 25 entries, corresponding to a SNR of about 5, to be
able to use the largest possible filter while keeping an optimal efficiency.

With the grid resolution fixed it is then possible to fine tune the algorithm
and determine the localization errors for different filters of different scale
factors, to find the combination which gives the smallest localization error.
Unless the step to be localized is very narrow even after the adjustment of
the grid resolution, the filter of choice is typically the FDOG, where the
filter scale factor corresponds to the parameter $\sigma$.

###  2.4 Calibration procedures for the detected edges

To define the systematic accuracy of the edge position measurement we need to
introduce some theoretical model of the underlying distribution that would
define the _true edge_ position. By using such a model it is possible to use a
toy Monte-Carlo procedure to calibrate the edge position measurement and to
define its systematic accuracy. In general, this can be done by calculating
the localization bias of the edge position measurement for each element in a
set of model parameters influencing the true edge position so that a
correlation between them can be derived. This correlation is model dependent
and can be studied in each specific case to estimate the systematic error
induced to the measurement for each set of model parameters. In the rest of
this article, two physics benchmarks with very different challenges and
different types of systematic uncertainties are studied, with the goal to
validate the proposed method in actual physics analyses.

##  3 Study Cases

Many models of new physics predict additional particles which promptly decay
into more stable particles, of which some can be seen in the detector, while
others may interact too weakly to leave any trace in the detector. If all
decay products can be detected, the mass of the new particle can simply be
calculated as the invariant mass of all decay products. If there are, however,
also “invisible” particles among the decay products, the detection of
kinematic edges in the energy (momentum) or invariant mass spectra of the
visible decay products is a well established technique to determine the mass
of the mother particle. This technique was widely used, for example, by the
LEP experiments when performing supersymmetry searches [8, 9, 10, 11]. Similar
techniques are used in supersymmetry searches at LHC [12, 13, 14, 15, 16, 17].
Beyond supersymmetry, partially invisible decays are expected e.g. in most
extensions of the Standard Model of particle physics which contain dark matter
candidate particles.

In order to test and evaluate the performance of our proposed edge
localisation method, we consider two different supersymmetry scenarios. Since
at lepton colliders the well-known momenta of the colliding elementary
particles allows a more complete reconstruction of the event kinematics, we
especially study the production of supersymmetric particles at a future
$e^{+}e^{-}$ collider. In both scenarios, the analysed reactions are of the
generic type
$e^{+}e^{-}\rightarrow\widetilde{Y}\bar{\widetilde{Y}}\rightarrow\widetilde{\chi}^{0}_{1}Y\widetilde{\chi}^{0}_{1}\bar{Y}$,
where the $\widetilde{\chi}^{0}_{1}$ is the lightest SUSY particle, assumed to
be stable and not detectable, $\widetilde{Y}$ is a generic SUSY particle, and
$Y$ its SM partner.

The kinematic edges comes about as follows: in the rest-frame of a system with
mass $M_{0}$ decaying to two particles with masses $M_{1}$ and $M_{2}$, one
has

$\displaystyle E_{1,(2)}=$
$\displaystyle\frac{M^{2}_{o}+M^{2}_{1,(2)}-M^{2}_{2,(1)}}{2M_{0}}$
$\displaystyle p_{1}=p_{2}=p=$
$\displaystyle\frac{\sqrt{\lambda\left(M^{2}_{0},M^{2}_{1},M^{2}_{2}\right)}}{2M_{0}}=\frac{\sqrt{\lambda_{0,1,2}}}{2M_{0}}$
(7)

where the Källén function $\lambda$ is defined by
$\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc$.

At production, $M_{0}=E_{cms}$, and in the case of production of a particle-
antiparticle pair $\widetilde{Y}\bar{\widetilde{Y}}$ one has
$M_{1}=M_{2}=M_{\widetilde{Y}}$, so that
$E_{1}=E_{2}=E_{\widetilde{Y}}=E_{cms}/2$ and
$p=(E_{cms}/2)\sqrt{1-\left(M_{\widetilde{Y}}/(E_{cms}/2)\right)^{2}}.$ In the
rest-frame of the decaying particle, one would identify $M_{0},M_{1}$ and
$M_{2}$ in eq. 7 with $M_{\widetilde{Y}}$, $M_{\widetilde{\chi}^{0}_{1}}$ and
$M_{Y}$, respectively. The energy of $Y$ in the laboratory frame,
$E_{Y\leavevmode\nobreak\ lab}$, is then given by

$\displaystyle E_{Y\leavevmode\nobreak\ lab}=$ $\displaystyle\gamma
E_{Y}+\gamma\beta\cos{\theta}p_{Y}$ (8)

where $\gamma$ and $\gamma\beta$ are the parameters of the Lorentz-
transformation between the laboratory frame and the rest frame of
$\widetilde{Y}$. The highest and lowest possible energies are obtained for
$\theta=0$ and $\theta=\pi$, respectively. Therefore, the kinematic edges are

$\displaystyle E_{Y\leavevmode\nobreak\ lab^{max}_{(min)}}=$
$\displaystyle\gamma E_{Y}\begin{smallmatrix}+\\\
(-)\end{smallmatrix}\gamma\beta p_{Y}$ (9)

and are independent of the distribution of $\cos{\theta}$.

The relation between the position of the edges and the masses
$M_{\widetilde{Y}}$ and $M_{\widetilde{\chi}^{0}_{1}}$, can be found by noting
that the Lorentz factors of the transformation from the laboratory frame to
the rest-frame of the produced particle are
$\gamma=E/M=E_{CMS}/(2M_{\widetilde{Y}})$ and $\gamma\beta=p/M$
=$E_{cms}/(2M_{\widetilde{Y}})$
$\sqrt{1-\left(M_{\widetilde{Y}}/(E_{CMS}/2)\right)^{2}},$ so that the edges
of the $E_{Y\leavevmode\nobreak\ lab}$ spectrum are given by

$\displaystyle E_{Y\leavevmode\nobreak\ lab^{max}_{(min)}}=$
$\displaystyle\frac{E_{CMS}}{4}\left(\frac{M^{2}_{\widetilde{Y}}+M^{2}_{Y}-M^{2}_{\widetilde{\chi}^{0}_{1}}}{M^{2}_{\widetilde{Y}}}\begin{smallmatrix}+\\\
(-)\end{smallmatrix}\sqrt{1-\left(\frac{M_{\widetilde{Y}}}{E_{CMS}/2}\right)^{2}}\frac{\sqrt{\lambda_{\widetilde{Y},Y,\widetilde{\chi}^{0}_{1}}}}{M^{2}_{\widetilde{Y}}}\right)$
(10)

In the study-case discussed in Sec. 4, one might set $M_{Y}=0$ , in which case
$E_{Y}=p_{Y}=(M^{2}_{\widetilde{Y}}-M^{2}_{\widetilde{\chi}^{0}_{1}})/(2M_{\widetilde{Y}}),$
so that eq. 10 simplifies to

$\displaystyle E_{Y\leavevmode\nobreak\ lab^{max}_{(min)}}=$
$\displaystyle\frac{E_{CMS}}{4}\left(1-\frac{M^{2}_{\widetilde{\chi}^{0}_{1}}}{M^{2}_{\widetilde{Y}}}\right)\left(1\begin{smallmatrix}+\\\
(-)\end{smallmatrix}\sqrt{1-\left(\frac{M_{\widetilde{Y}}}{E_{CMS}/2}\right)^{2}}\right)$
(11)

In eq. 11, it is easy to see that by calculating the sum of the upper and
lower edge of the spectrum ($\Sigma$), $M^{2}_{\widetilde{\chi}^{0}_{1}}$ can
be expressed in $M^{2}_{\widetilde{Y}}$, which can then be used to eliminate
the former in the expression for the difference between the edges ($\Delta$).
From the obtained expression, $M^{2}_{\widetilde{Y}}$ is obtained, and by
back-substitution also $M^{2}_{\widetilde{\chi}^{0}_{1}}$ is obtained.

$\displaystyle M_{\widetilde{Y}}=$
$\displaystyle\frac{E_{CMS}}{2}\frac{\sqrt{\Sigma^{2}-\Delta^{2}}}{\Sigma}$
$\displaystyle M_{\widetilde{\chi}^{0}_{1}}=$
$\displaystyle\frac{E_{CMS}}{2}\sqrt{\frac{E_{CMS}-2\Sigma}{E_{CMS}}}\frac{\sqrt{\Sigma^{2}-\Delta^{2}}}{\Sigma}$
(12)

Similar, albeit more complicated, expressions can also be derived from eq. 10.

The two example scenarios, which will be introduced in more detail in the
following subsections, were studied in the context of the International Linear
Collider (ILC) [18, 19, 20, 21] and the International Large Detector (ILD)
[22, 23].

The ILC is a planned electron-positron linear collider. The construction of
its first stage with
$\sqrt{s}=$250 GeV is under political consideration in Japan. In this article
we consider an ILC centre-of-mass energy of $\sqrt{s}=$500 GeV, an integrated
luminosity of $\int\mathcal{L}dt$ = 500 fb-1 and a beam polarisation
configuration of 80% right-polarised electrons and 30% left-polarised
positrons, i.e., $\mathcal{P}$($e^{-}$, $e^{+}$) = (+80%, -30%). This
corresponds to only small fraction of the total data expected from the full
run plan of the ILC [24, 25]. All Monte-Carlo events used in the two case
studies have been generated with Whizard 1.95 [26], with PYTHIA 6.4 [27] used
for fragmentation and hadronization, and taking into account the beam energy
spectrum, initial and final state radiation, and the effects of
beamsstrahlung.

ILD is one of the two detector concepts proposed for the ILC, which were
specifically designed according to the Particle Flow principles. ILD will be
equipped with a time projection chamber acting as its main tracker and a very
granular calorimetric system, with a fine longitudinal and transverse
segmentation, optimized for particle flow reconstruction. The ILD design is
implemented in a detailed Geant4-based simulation [28, 29], and the simulated
detector response is reconstructed by the algorithms implemented in the Marlin
framework [30], including tracking as well as the particle flow algorithm
PandoraPFA [31, 32, 33].

###  3.1 The STC4 Scenario

The $\tilde{\tau}$-coannihilation scenario (STC) has been defined and
motivated in [34] in the theoretical framework of the phenomenological minimal
supersymmetric Standard Model (pMSSM) [35, 36]. It is characterised by a light
slepton sector and heavy squarks that may have escaped detection at the LHC.
The $\tilde{\tau}$ is the next-to-lightest supersymmetric particle with a mass
approximately 10 GeV larger than that of the $\tilde{\chi}_{1}^{0}$, i.e., the
stable LSP (_M_ ${}_{\tilde{\chi}_{1}^{0}}$ = 95.6 GeV).

Several different STC model points have been defined [37], which differ mainly
in the value assumed for the right-handed third generation squark mass
parameter at 1 TeV. In case of the STC4 point, considered in our evaluation,
the value of this parameter is assumed to be 400 GeV. The part of the spectrum
relevant to an $e^{+}e^{-}$ collider with $\sqrt{s}=$500 GeV depends only
weakly on this parameter and is for all practical purposes the same for all
model points of the STC series.

Our study concentrates on selectron pair production at the ILC and their
subsequent, almost exclusive decay:
$e^{+}e^{-}\rightarrow\tilde{e}\tilde{e}\rightarrow
e^{\pm}\tilde{\chi}_{1}^{0}e^{\mp}\tilde{\chi}_{1}^{0}$. This takes advantage
of the relatively large cross-section that characterises selectron pair
production in the STC4 scenario, due to the t-channel neutralino exchange:
$\mathcal{O}$($10^{2}$ fb) for the $\tilde{e}_{R}\tilde{e}_{R}$ case and
$\mathcal{O}$(10 fb) for the $\tilde{e}_{L}\tilde{e}_{L}$ case. Furthermore,
the event selection also benefits from the higher efficiency to identify
electrons in contrast to $\tau$ decays. Lastly, if _M_
${}_{\tilde{\chi}_{1}^{0}}$ is extracted from the selectron momentum spectrum
the $\tilde{\tau}$ mass can also be determined, using the same data set, from
the endpoint of the spectrum of the $\tau$ decay products.

Based on the assumed STC4 parameter values, the mass spectrum calculator
SPheno [38, 39] predicts a value of _M_ ${}_{\tilde{e}_{L}}$ = 213.1 GeV for
the left-handed selectron mass and _M_ ${}_{\tilde{e}_{R}}$ = 126.2 GeV for
its right-handed counterpart. The CMS experiment [40, 41] has already excluded
left-handed selectrons in the mass range assumed by the STC4 model. However,
their right-handed partners have not yet been excluded. They are only $\sim$30
GeV heavier than the LSP and, due to this small mass difference, the amount of
missing energy in the LHC events would not be sufficient to allow a clear
separation from the Standard Model background.

While theoretically the isotropic decay of a scalar particle in its rest frame
is expected to lead to a pure box-like distribution of the lepton momentum in
the lab frame, the initial state radiation and the interactions of the final-
state electrons in the detector distort the shape of the decay product
momentum distribution so that it cannot be reasonably described by any
specific analytic function calculated from first principles. Therefore this
scenario is an excellent benchmark for our edge detection algorithm. Moreover,
the same simulated data was used to estimate the same quantity using a
different technique [37], which can serve as a direct reference to estimate
the performance of our method. In that publication, the method relied on the
fact that, for an undistorted rectangular distribution without background,
estimators surpassing the CRB (the Cramér-Rao Bound) for the minimal variance
exist for the high and low edges111The reason that the CRB can be surpassed in
the case of the rectangular distribution is that this distribution does not
fulfil one of the conditions of applicability of the Cramér-Rao bound, namely
that the support of the density does not depend on the fitted parameter. They
are simply the smallest and largest observations in the sample, respectively.
These statistics are only asymptotically unbiased, but the finite sample-size
bias is analytically known. The variance of the estimators decreases as the
reciprocal of the _s_ quare of the sample size. To mitigate the distortions
and the presence of background, the approach was to divide the sample in sub-
samples, and estimate the edges, not by the highest and lowest observations,
but to first exclude a fraction of the extreme values. The edge estimates were
then taken as the average of the truncated sub-sample extremes. The two
parameters - sub-sample size and truncation factor - were optimised to give
the smallest variance using a toy Monte-Carlo procedure. The toy MC was also
used to calibrate the estimators in order to remove the bias222The bias due to
a truncation of a pure rectangular distribution is analytically calculable,
but not the effects of background and distortions.. It was found that _M_
${}_{\tilde{e}_{R}}$ could be determined to 210 MeV and _M_
${}_{\tilde{\chi}_{1}^{0}}$ to 160 MeV based on an integrated luminosity of
500 fb-1 with $\mathcal{P}$($e^{-}$, $e^{+}$) = (+80%, -30%) at a center-of-
mass energy of 500 GeV.

###  3.2 The Point 5 Scenario

The supersymmetry scenario denoted as "_Point 5_ " was first defined [42] in
the context of the minimal Supergravity model [43, 44], mSUGRA. It was
originally considered for benchmarking the Particle Flow performance of the
ILD detector concept.

In the Point 5 model, the first chargino ($\tilde{\chi}_{1}^{\pm}$) and the
second neutralino ($\tilde{\chi}_{2}^{0}$) are almost mass degenerate: _M_
${}_{\tilde{\chi}_{1}^{\pm}}$ = 216.5 GeV and _M_ ${}_{\tilde{\chi}_{2}^{0}}$
= 216.7 GeV. They could be produced in pairs at the ILC with an operating
centre-of-mass energy of $\sqrt{s}$ = 500 GeV. The production cross-section is
of the order of $10^{2}$ fb in the chargino case and of the order of $10$ fb
in the neutralino case.

Both gauginos decay with branching ratios larger than 96% to the corresponding
gauge boson and the first neutralino which is the stable LSP:

$\displaystyle e^{+}e^{-}$ $\displaystyle\rightarrow$
$\displaystyle\tilde{\chi}_{1}^{+}\tilde{\chi}_{1}^{-}\rightarrow
W^{+}\tilde{\chi}_{1}^{0}W^{-}\tilde{\chi}_{1}^{0}$ $\displaystyle e^{+}e^{-}$
$\displaystyle\rightarrow$
$\displaystyle\tilde{\chi}_{2}^{0}\tilde{\chi}_{2}^{0}\rightarrow
Z^{0}\tilde{\chi}_{1}^{0}Z^{0}\tilde{\chi}_{1}^{0}$ (13)

Their mass difference with respect to the LSP mass (_M_
${}_{\tilde{\chi}_{1}^{0}}$ = 115.7 GeV) is approximately 100 GeV, i.e. large
enough that the _W_ and _Z_ bosons can be produced on-shell.

Our study aims to reconstruct and determine the masses of all three gauginos:
$\tilde{\chi}_{1}^{\pm}$, $\tilde{\chi}_{2}^{0}$ and $\tilde{\chi}_{1}^{0}$.
Consequently, each one of the two reactions presented above acts as background
for the other. The separation of the chargino and neutralino candidate events
essentially relies on identifying and distinguishing _W_ +_W_ - and _Z_ _Z_
decays. Only the hadronic decays of the gauge bosons are here considered as
signals. This was a deliberate choice, posing a significant challenge both for
the Particle Flow reconstruction and for the proposed edge detection method.
In practice, the (semi-)leptonic modes would of course be included, although
their contribution is statistically limited.

The relevant observables for determining the gaugino masses are the di-jet
energy spectra, where each di-jet corresponds to a gauge boson. Our proposed
edge localisation method must handle the smearing effects on the edges of the
di-jet energy distributions produced by the beam energy spectrum. In addition,
it must also deal with the contributions coming from the intrinsic gauge
bosons’ widths: $\Gamma_{W}\simeq$ 2.1 GeV and
$\Gamma_{Z}\simeq=2.5\,\mathrm{GeV}$ [45]. The additional smearing effect due
to the finite jet energy resolution is rather small compared to these two.

In the original Point 5 study [46], the positions of the edges were determined
by fitting the di-jet energy with a function in which the edge positions were
free parameters. In order to model the signal, a second order polynomial,
describing the part of the spectra between the edges, was used. This was
convoluted with a Voigt function in order to model both the effect of the
width of the gauge bosons and of the uncertainties related to the measurement
of jet energies. The Standard Model background contributions were assumed to
be well understood by the time of ILC operation, and thus they were described
independently by another heuristic function $f_{SM}$ which was fixed in the
edge position fit:

$f_{total}=f_{SM}+\int_{t_{0}}^{t_{1}}(a_{sig}\cdot t^{2}+b_{sig}\cdot
t+c_{sig})Voigt(t-x,\sigma,\Gamma)\mathit{dt}$ (14)

The convolution limits ($t_{0}$ and $t_{1}$) represent the edge positions and
they are free parameters of the function, together with $a_{sig}$, $b_{sig}$,
$c_{sig}$, $\sigma$, and $\Gamma$.

In order to evaluate the robustness of the fitting method, a toy Monte Carlo
study was performed. For this purpose, 104 di-jet energy distributions (for
each of the two gaugino candidates’ samples) were produced as random
variations of the histograms obtained in the full detector simulation-based
analysis. The fitting function given by Eq. 14 was applied for each one of the
newly generated distributions. From 104 performed fits, 8808 converged for the
chargino case and 8282 for the neutralino case.

Figure 7: Results of the toy Monte Carlo study performed to evaluate the
consistency of the fitting method for kinematic edge extraction. The number of
converged fits is illustrated by the colour code. The black, dashed vertical
line shows the location of the model calculated edge.

In Fig. 7, the $\chi^{2}/$ndf ratio, as a measure of the fit performance, is
plotted with respect to the determined edge position for the neutralino upper
edge as an example. The information concerning the number of fits associated
with a specific $\chi^{2}/$ndf - edge position configuration can be read from
the colour code.

The instability of the fitting method is indicated by the wide range of
obtained edge values. It can be seen from Fig. 7 that the obtained edge values
vary over a range of approximately 40 GeV despite the apparently good quality
of the fit ($\chi^{2}/$ndf < 1.5). This effect is significantly larger than
any expected shift in the edge position that could occur from the random
generation of the distribution. Furthermore, the majority of fits falsely
determined that the higher edge positions coincide with the histrograms’ end
points at 170 GeV.

Thus, as a result of the toy Monte Carlo study it could be assessed that the
fitting method is highly sensitive to small fluctuations in the di-jet energy
spectra. This strongly illustrates the need to develop a more robust edge
detection method. We will show in the following sections that the FIR filter
approach is significantly more stable and accurate in measuring the edge
positions.

##  4 Measurement of the $\tilde{e}$ and $\tilde{\chi}^{0}_{1}$ masses at the
ILC

The first study case is the measurement of the masses of the $\tilde{e}_{R}$
and the $\tilde{\chi}^{0}_{1}$, respectively ${M_{\tilde{e}_{R}}}$,
${M_{\tilde{\chi}^{0}_{1}}}$, from the position of the kinematic end-points of
the momentum distribution of the lepton produced in the final state. For this
analysis, the signal and the full SM background were processed by the fast
simulation program SGV [47] which models the response of the tracking system
of the ILD detector from first principles, complemented by a parametrised
response of the particle-flow calorimeters. Those data are finally processed
by an electron identification algorithm, described in [7]. This tool was made
specifically to operate on the SGV output to identify isolated electrons and
the photons they may have radiated before entering the calorimeters.

In the process we are analysing, the mass of the two visible leptons in the
final state is negligible compared to the masses of the other particles
involved. Therefore we can use eq. 11 to estimate the energy distribution of
those electrons, which would ideally be uniform between two sharp, rectangular
edges whose positions we can calculate using eq. 12. In the STC4 scenario,
with $M_{\tilde{\chi}_{1}^{0}}$ = 95.6 GeV and $M_{\tilde{e}_{R}}$ = 126.2 GeV
and a centre-of-mass energy of 500 GeV, the ideal position of the kinematics
edges of the $e^{+}e^{-}\rightarrow\bar{\tilde{e}}_{R}\tilde{e}_{R}\rightarrow
e^{+}e^{-}+2\tilde{\chi}^{0}_{1}$ are $E_{e}^{max}$ = 99.362 GeV and
$E_{e}^{min}$ = 7.298 GeV.

Compared to these ideal values, when considering the realistic distribution of
the centre-of-mass energies, initial-state radiation and detector
uncertainties, we expect the actual positions of the kinematic end-points to
be shifted in a way that cannot be predicted analytically, therefore making
this case an ideal test to measure the performances of our method.

Figure 8: Momentum distribution of the higher-momentum electron after
application of the high-momentum edge selection. This distribution is used as
a template to calculate the position of the high-momentum edge and its
uncertainty.  Figure 9: Momentum distribution of the lower-momentum electron
after application of the low-momentum edge selection described in the text.
This distribution is used as a template to calculate the position of the low-
momentum edge and its uncertainty.

Because the magnitude and sources of the distortions can depend on the energy
of the final-state particle and because the backgrounds are significantly
different at the location of the two end-points, we analyse them separately.
In both cases the selection is based on cuts which identify events with
negligible hadronic activity, large amounts of missing energy and exactly two
isolated electrons of opposite charge. To reduce the $\gamma\gamma$
background, which dominates at low electron momenta, a cut on the acolinearity
and acoplanarity of the two identified electrons with the higher momentum is
applied. The momentum distribution of those two electrons is finally used to
identify the high- and low-momentum end-points. Those distributions, after the
event selection, are shown in Figs. 8 and 9. Taking the whole data set into
account, for the high-momentum selection we achieved a signal purity of 56.5%
with an efficiency of 82% while for the low-momentum selection we achieved a
signal purity of 87.6% with an efficiency of 30%. Full details of the
selection can be found in [7].

###  4.1 Edge position measurement and calibration

(a) High-momentum edge optimization

(b) Low-momentum edge optimization

Figure 10: The plots show the filter optimization graphs for the FDOG, Deriche
and Shen filters, applied to the histograms shown in fig. 8 and fig. 10(b), on
the left and on the right respectively.

To optimize the parameters of the edge detection algorithm we performed the
procedure described in Sec. 2.3. The distributions shown in Figs. 8 and 9 are
used as a template for a toy Monte Carlo with which we optimized the binning
of the distribution and minimized the filter localization error
$\bm{\sigma_{L}}$. As can be seen from Figs. 8 and 9, bin widths of 50 MeV and
250 MeV for the high-momemtum and the low-momentum edges, respectively,
correspond to an average bin-by-bin difference of 25 entries in the edge
region, i.e. to SNR = 5. The plots in Fig. 10 show the localization error of
the algorithm as a function of the filter size for this binning. In case of
the high-momentum edge, Fig. 10(a), the FDOG filter performs best, and for the
final result we use it with a scale factor of 1.2 GeV. Larger scale factors
could have improved the statistical uncertainty a little further, but, as we
will see in the following, the systematic uncertainty on the localisation bias
is of a similar size, and thus we to not push the statistical precision to the
absolute minimum. Due to its steepness, the low-momentum edge shows a
significantly less stable situation in Fig. 10(b). The best performance is
reached with the Shen filter, which is not surprising since it has been
designed specifically for the case of very narrow edges. However, the FDOG
performs worse only by about 10%, relatively, and shows a more stable
behaviour. Therefore and in view of the above mentioned localisation bias
uncertainty, we chose to use the same filter for all cases and thus employ the
FDOG with a scale factor of 0.25 GeV for the low-momentum edge final results.
Using these parameters, we measure the high-momentum edge position to be
$E^{max}=98.748\pm 0.043\,\mathrm{GeV}$ and the low-momentum edge position
$E^{min}=7.409\pm 0.012\,\mathrm{GeV}$.

The second step is the calibration of the edge position measurements, in other
words the determination of the localisation bias and its uncertainty. This was
done by separating the signal data set from the backgrounds and, while keeping
the latter constant, we applied a random offset to the former generating a new
template distribution that we used to recalculate the edge position. This
procedure is repeated many times for different offsets to obtain the
distribution of the localization bias as a function of the applied offset. The
results of the calibration procedure are shown in Fig. 11 where the measured
bias is displayed as a function of the offset chosen for each measurement.

A striking feature observed in both plots is a modulation with a constant
periodicity. This is caused by the residual effect of the underlying
background and is connected with another feature which, albeit not used for
the optimization procedure, we defined and characterized in [7]. That feature
was already described in [3], namely the multiple response sensitivity or,
otherwise said, the minimum distance of two edges which could be
discriminated. In the case of the FDOG filter this distance is roughly equal
to 4 times the scale factor of the filter, which is consistent with the
periodicity observed in the calibration plots.

(a) Low-momentum edge calibration

(b) High-momentum edge calibration

Figure 11: Calibration graphs of the high-and low-momentum edge measurements.
Each data point corresponds to a measurement of the localization parameters on
a different template distribution obtained by applying an offset to the signal
sample. Superimposed to the measurement set we displayed a first degree
polynomial fit to the data showing a long range independence of the bias from
the offset.

Since in an actual measurement one would not know where the true edge position
sits within a period, the localisation bias should not be determined directly
from the value obtained at an offset of 0, but from the value a linear
function fitted to the whole distribution attains at 0 offset. By comparing
this value to the ideal edge positions calculated from the generator masses of
the supersymmetric particles, the average localization bias is $\bm{\hat{L}}$
-0.601 GeV and 0.109 GeV for the two edges respectively. The amplitude of the
periodic structure adds a systematic uncertainty to the localisation
procedure, which we defined as the standard deviation of all bias values from
the calibration procedure, which yields 42 MeV for the high-momentum edge and
8 MeV for the low-momentum edge. Thus, the systematic uncertainty of the
filter method is of the same order or even smaller than the statistical
uncertainty on the edge position. Since the periodic structure originates from
statistical features in the background, the systematic uncertainty is expected
to shrink at a similar rate as the statistical one when considering larger
data sets.

###  4.2 Edge positions and mass measurements

Summarizing the results of the edge analysis, we can write our final edge
measurements results based on $500\,\mathrm{fb^{-1}}$ of data with $+80\%$
electron polarization and $-30\%$ positron polarization as:

$\displaystyle\bm{E^{max}}$ $\displaystyle=99.349\pm 0.043\text{(stat)}\oplus
0.042\text{(bias)}\,\mathrm{GeV}$
$\displaystyle\qquad\bm{E^{max}_{\textbf{true}}}$
$\displaystyle=99.362\,\mathrm{GeV}$ (15) $\displaystyle\bm{E^{min}}$
$\displaystyle=7.300\pm 0.012\text{(stat)}\oplus
0.008\text{(bias)}\,\mathrm{GeV}$
$\displaystyle\qquad\bm{E^{min}_{\textbf{true}}}$
$\displaystyle=7.298\,\mathrm{GeV}$ (16)

with the _true_ positions calculated assuming a nominal centre-of-mass energy
of 500 GeV as well as ${M_{\tilde{e}_{R}}=126.235\,\mathrm{GeV}}$ and
${M_{\tilde{\chi}^{0}_{1}}=95.586\,\mathrm{GeV}}$ as the _true masses_ of the
supersymmetric particles.

Finally, we calculate the masses of the selectron and of the neutralino from
the measured edge positions, in order to compare them to the results of other
studies. As both masses depend on the lower and upper edge, see eqn. 12, their
calculated values will be correlated. Using a Monte Carlo procedure described
in [7], we can extract the masses of the $\tilde{e}_{R}$ and of the
$\tilde{\chi}_{1}^{0}$ with accuracies of:

$\displaystyle\bm{M_{\tilde{e}_{R}}}$ $\displaystyle=126.20\pm
0.11\,\mathrm{GeV}$
$\displaystyle\qquad\bm{M_{\tilde{e}_{R}}^{\textbf{true}}}$
$\displaystyle=126.235\,\mathrm{GeV}$ (17)
$\displaystyle\bm{M_{\tilde{\chi}_{1}^{0}}}$ $\displaystyle=95.56\pm
0.09\,\mathrm{GeV}$
$\displaystyle\qquad\bm{M_{\tilde{\chi}_{1}^{0}}^{\textbf{true}}}$
$\displaystyle=95.586\,\mathrm{GeV}$ (18)
$\displaystyle\text{cov}(M_{\tilde{e}_{R}},M_{\tilde{\chi}_{1}^{0}})$
$\displaystyle=0.0099\,\mathrm{GeV^{2}}$ $\displaystyle\qquad\rho$
$\displaystyle=0.97$ (19)

Here, the given uncertainties comprise both the statistical uncertainty on the
edge position as well as the uncertainty from the bias correction. They are
about a factor of 2 more precise than the best previous results based on the
same data set, but obtained using a different method for locating the edge
positions as explained in Sec. 3.1. The fact that the initial generator values
for the two sparticle masses fall within the error range of the measurements
is a good indication of the reliability of the measurement algorithm and of
the procedure we used to estimate the measurement errors.

##  5 Edge Measurements with FIR Filters in the Point 5 Study Case

The second study case, introduced in Sec. 3.2, aims to measure the masses of
the three lightest gauginos from the kinematic endpoints in
$\tilde{\chi}_{1}^{\pm}$ and $\tilde{\chi}_{2}^{0}$ pair production and was
carried out with the full ILD simulation [28, 29] and reconstruction [31, 32,
33] chain. In particular, the challenging but frequent hadronic decay modes
were considered, which necessitates the full state-of-the-art modelling of the
detector as well as full particle flow reconstruction, as argued in Sec. 3.

Considering the two signal reactions presented in equation 13, the event
selection searches for topologies characterised by: (i) four hadronic jets in
the final state, all localised in the main ILD calorimeters, (ii) little
leptonic activity, and (iii) a large amount of missing energy. For both the
chargino and the neutralino samples, the selection efficiency is approximately
91%, however the purity is very low at this stage [48]. This is not only due
to the SM background contributions, but also the fact that each signal channel
acts as background for the other. To mitigate this, a $\chi^{2}$ test is
carried out, comparing the reconstructed di-jet invariant masses to both the
_W_ and the _Z_ boson masses. Based on the obtained $\chi^{2}$ values,
$\tilde{\chi}_{1}^{\pm}$ and $\tilde{\chi}_{2}^{0}$ candidates are selected
from the events passing the rest of the selection. This enables a final sample
purity of 63% for the charginos and 39% for the neutralino case at
efficiencies of 53% and 30%, respectively.

The jet pairing is based on the output of a kinematic fit applied with an
equal mass constraint, i.e., the fit does not enforce a specific mass value
for the di-jet pairs, but requires that they have equal masses to reflect the
expected presence of the gauge boson pair in the final state. For this
purpose, the MarlinKinFit [49] package was used. The di-jet configuration
providing the largest fit probability was chosen and used further in the
analysis.

The relevant observable for determining the gaugino masses is the di-jet
energy. The two distributions of this observable, obtained from the sets of
$\tilde{\chi}_{1}^{\pm}$ and the $\tilde{\chi}_{2}^{0}$ candidate events,
after applying the kinematic fit, are shown in Fig. 12.

(a) $\tilde{\chi}^{\pm}_{1}$ candidates

(b) $\tilde{\chi}^{0}_{2}$ candidates

Figure 12: The di-jet energy spectra from (a) the chargino and (b) the
neutralino candidate events used in the kinematic edge extraction.
Distributions based on events simulated with a detailed description of the ILD
detector concept. The black vertical lines indicate the positions of the
analytically computed kinematic edges.

Based on the masses (_M_ ${}_{\tilde{\chi}_{1}^{\pm}}$, _M_
${}_{\tilde{\chi}_{2}^{0}}$ and _M_ ${}_{\tilde{\chi}_{1}^{0}}$) assumed in
the Point 5 model, the ideal positions of the kinematic edges can be
calculated analytically. They are indicated by the vertical lines in Fig. 12.
This calculation however does not take into account important aspects like
e.g. the natural width of the decaying gauge bosons, the shape of the beam
energy spectrum or the finite resolution of the detector and the
reconstruction algorithms, which lead to a much more complicated shape of the
realistic spectra, and can even shift the edge positions. All these effects
are taken into account in our simulation and are present in the distributions
in Fig. 12. Therefore, we expect that a calibration procedure will be
necessary in order to retrieve the dependency of the edge positions on the
particle masses.

As in the STC4 analysis, an FIR filter with the first derivative of a Gaussian
(FDOG) kernel was applied to measure the edge positions in the energy spectra
shown in Fig. 12. Before performing the edge measurement, the parameter values
of the FIR filter were optimised following the receipe introduced in Sec. 2.3.

###  5.1 Filter Parameters Optimisation

For the optimisation, 104 variations of the di-jet energy spectra of the
$\tilde{\chi}^{\pm}_{1}$ and the $\tilde{\chi}^{0}_{2}$ candidates were
randomly generated from the original ones (Fig. 12) using a toy Monte Carlo
approach.

The first step was to choose an appropriate binning of the input data. The di-
jet energy spectra upper edges are the noisiest ones not only due to beam
energy spectrum effects, but also because the Standard Model background
contributions become predominant in that region. Particularly in case of the
$\tilde{\chi}^{0}_{2}$ spectrum, shown in Fig. 12(b), a SNR > 5 cannot be
reached unless the whole upper egde is collapsed into 2-3 bins. With a binning
of 1 GeV per bin, the average count difference between neighbouring bins in
the upper edge region is 7.5, corresponding to an SNR $\approx 2.7$, for which
only a small penalty in terms of the localization error is expected, c.f. Fig.
5. Both parameters, the scale factor of the filter and the binning of the
input distribution were varied to verify that the performance expectations
from Sec. 2 also hold in the realistic physics example.

Figure 13: Filter responses obtained when varying the Gaussian $\sigma$
parameter (left) and the input data binning (right) in the context of the
optimisation procedure. The shown distributions were attained using the
simulated $\tilde{\chi}^{0}_{2}$ candidate events. A very similar behaviour
was observed in the case of the $\tilde{\chi}^{\pm}_{1}$ candidates.

Optimisation of the filter scale factor (Gaussian $\sigma$):
To determine the optimal value of $\sigma$, the binning of the input
distributions were made with 100 bins of width 1 GeV. The $\sigma$ value was
varied within a range of $\sigma\in\\{1\dots 10\\}$, in steps of one. The FDOG
filter was then applied with each $\sigma$ value on the toy Monte Carlo
produced $\tilde{\chi}^{\pm}_{1}$ and $\tilde{\chi}^{0}_{2}$ di-jet energy
spectra. The obtained filter responses are shown in Fig. 13.

The very low $\sigma$ values could be immediately excluded since the
corresponding FIR filter response was too noisy for a meaningful edge
extraction. In contrast, all values above 4 provide almost identical filter
responses. Consequently, a value of $\sigma$ = 4 bins (shown in green) was
chosen as most appropriate for the analysis.

Optimisation of the bin size:
To find the optimal bin size for the input di-jet energy spectra, the $\sigma$
parameter was fixed to its previously determined optimal value: $\sigma$ = 4 .
A total of four different bin sizes were considered and the filter responses
obtained for each of them can be seen in Fig. 13.

It was determined that, with a binning of 0.5 GeV/bin, the filter response
could not provide clearly distinguishable and relevant maxima due to too much
noise. For bin sizes larger than 1 GeV/bin, corresponding to smaller step
widths, the peaks marking the edge positions became wider, corresponding to
larger localization errors, in accordance with the expectation from Fig. 6.
Thus, a bin size of 1 GeV/bin was chosen for the final analysis.

###  5.2 Measured Edge Values

Following the optimisation process, the filter $\sigma$ parameter was fixed to
its determined optimal value, $\sigma$ = 4 , and the input data was set to a
bin size of 1 GeV/bin. The FDOG filter was then applied on the di-jet energy
spectra obtained from the $\tilde{\chi}^{\pm}_{1}$ and $\tilde{\chi}^{0}_{2}$
candidate events (Fig. 12) and the filter responses are presented in Fig. 14.
The upper panels of the figure show the original di-jet energy spectra, i.e.
the same distributions from Fig. 12; the lower panels present the computed
filter response, i.e. the result of applying the filter on the di-jet energy
distributions shown above. The measured edge positions and the corresponding
analytically calculated values are indicated by the vertical lines.

The statistical uncertainty was evaluated from 104 toy Monte Carlo experiments
based on the di-jet energy spectra for each of the $\tilde{\chi}^{\pm}_{1}$
and the $\tilde{\chi}^{0}_{2}$ candidate samples. The FDOG filter was applied
on these spectra and both the low and the high edge positions were determined
by performing a Gaussian fit to the maxima present in the filter response.
Consequently, four distributions of measured edge values were obtained, one
for each of the four edges. The widths (RMS) of these distributions represent
the statistical uncertainties. The central values and uncertainties obtained
for $500\,\mathrm{fb^{-1}}$ at $\sqrt{s}$=500 GeV with $\mathcal{P}$($e^{-}$,
$e^{+}$) = ($+80\%$, $-30\%$) are shown in Table 1, along with the calculated
edge positions. There is a bias of up to 2 GeV in the determination of edge
positions, due to the inevitable effects from beam energy spectrum and gauge
boson width introduced above, which demonstrates the need for a calibration
procedure.

Figure 14: Results of applying the FIR filter method for extracting the kinematic edge position in the "Point 5" scenario. The model calculated values are also shown for comparison. Edge |  $\tilde{\chi}^{\pm}_{1}$ low |  $\tilde{\chi}^{\pm}_{1}$ high |  $\tilde{\chi}^{0}_{2}$ low |  $\tilde{\chi}^{0}_{2}$ high
---|---|---|---|---
Model calculated | 80.17 | 131.53 | 93.24 | 129.06
FDOG filter | 79.49 $\pm$ 0.15 | 129.53 $\pm$ 0.74 | 92.11 $\pm$ 0.31 | 128.38 $\pm$ 0.75
Table 1: Kinematic edge positions measured with the FIR filter method, using
the first derivative of a Gaussian as the filter kernel. The analytically
computed edges are also provided - they represent only an estimation of the
real edge positions, as the calculation does not include beam spectrum and
gauge boson widths effects.

###  5.3 Edge Calibration

The mass calculation presented in Sec. 3 entails some approximations. Firstly,
the effects of the beam energy spectrum on the edge localisation are not taken
into account. In addition, in the Point 5 scenario, the gauge bosons’ natural
widths also impact the edge positions. This effect is also not considered in
the calculation. In order to account for these approximations and also for the
reconstruction effects, a mass calibration procedure was performed. However,
in this study, the observables determined directly from the data are the edge
positions and not the masses, hence, it is the former that must be calibrated.
Five new data samples were simulated and reconstructed for this purpose with
the ILD full simulation.

In producing the new data samples, the mass of the LSP,
$\tilde{\chi}^{0}_{1}$, was left unchanged and set to the model value of
$M_{\tilde{\chi}^{0}_{1}}$ = 115.7 GeV. The masses of the two gauginos,
$\tilde{\chi}^{\pm}_{1}$ and $\tilde{\chi}^{0}_{2}$, which are almost
identical in the Point 5 scenario, were then varied simultaneously from 210
GeV to 225 GeV in steps of 3 GeV. The Standard Model background events
produced with the ILD detector full simulation were left unchanged.

Since the LSP mass is not known a priori, a further calibration procedure
should involve fixing the two $\tilde{\chi}^{\pm}_{1}$ and
$\tilde{\chi}^{0}_{2}$ masses to their Point 5 model values while the
$\tilde{\chi}^{0}_{1}$ mass is varied instead. However, this was not performed
in the present analysis and remains the subject of future studies.

After producing the five new Monte Carlo data samples, the corresponding
_measured_ edge positions were obtained by applying the FIR filter method
separately. A calibration function was then obtained for each edge, by
plotting the corresponding measured values with respect to expected edge
values calculated from the input masses. The results are presented in Fig. 15.
Since the $\tilde{\chi}^{\pm}_{1}$ low edge overlaps with the mass of the $W$
boson, it is not useful in the mass calculation and was therefore neglected in
the calibration procedure.

Figure 15: The edge position on reconstruction level versus the model-based
calculated values. The red continuous line represents the obtained calibration
line. The error projection is also shown (black dashed lines).

The ideal slope of the calibration linear fit has a value equal to one. This
would indicate that the effects considered in the calibration would have a
small impact on the edge positions. While the neutralino upper edge is close
to this value, Fig. 15 shows a deterioration that can be observed especially
in the case of the neutralino lower edge. This is not entirely unexpected
since, for the lowest $\tilde{\chi}^{0}_{2}$ mass,
$M_{\tilde{\chi}^{0}_{2}}$=210.7 GeV, the phase space volume for the reaction
$\tilde{\chi}^{0}_{2}\rightarrow Z\tilde{\chi}^{0}_{1}$ tends to zero
($M_{\tilde{\chi}^{0}_{2}}-M_{\tilde{\chi}^{0}_{1}}\approx$ 95 GeV).
Furthermore, since the phase space is small, the influence of the beam energy
spectrum and of the gauge boson width is much larger.

To determine each calibrated edge position, the corresponding measured value,
presented in section 5.2, was projected onto the linear fit (red continuous
line in Fig. 15). The intersection point with the calibration line was then
further projected onto the $x$-axis. This final projection gives the
calibrated value which can be read directly off the horizontal axis. The
uncertainty on the measured values was also taken into account in the
procedure.

|  $\tilde{\chi}^{\pm}_{1}$ (high) (GeV) |  $\tilde{\chi}^{0}_{2}$ (low) (GeV) |  $\tilde{\chi}^{0}_{2}$ (high) (GeV)
---|---|---|---
Reconstructed | 129.53$\pm$0.74 | 92.11$\pm$0.31 | 128.38$\pm$0.75
Calculated | 132.77 | 93.09 | 129.92
Calibrated | 132.46$\pm$0.44 | 92.52$\pm$1.23 | 130.04$\pm$0.77
Table 2: Results of the edge calibration comparing reconstruction level edge
values to the calculated values.

Table 2 presents the edge calibration results and the comparison between the
reconstructed and the calculated edge positions. The cumulative impact of the
natural boson width plus the beam and reconstruction effects can be evaluated
by comparing the “Reconstructed” and the “Calculated” columns. It can be seen
that, after the calibration, the chargino upper edge decreased by 2%. In the
neutralino case, this effect is of the order of 1% for the lower edge while
the upper edge decreases by 0.7%.

Lastly, Table 3 shows the outcome of the calibration procedure when
calculating the gaugino masses. While the mean $\tilde{\chi}^{\pm}_{1}$ is
shifted lower by 1.2%, the final measured mass is well compatible with the
original model values. The neutralino mass was corrected by 1.6% while the LSP
mass was improved by 2%. This is, however, achieved at the cost of an increase
of the uncertainty. But most importantly, the main goal of establishing a
stable and robust method for the edge localization, which was not available
before, has clearly been reached.

Gaugino |  Mass [GeV] (before calib.) |  Mass [GeV] (after calib.) | Model mass [GeV]
---|---|---|---
$\tilde{\chi}^{\pm}_{1}$ | 216.7$\pm$3.1 | 214.1$\pm$4.8 | 216.5
$\tilde{\chi}^{0}_{2}$ | 220.4$\pm$1.3 | 216.9$\pm$3.4 | 216.7
$\tilde{\chi}^{0}_{1}$ | 118.1$\pm$0.9 | 115.5$\pm$1.8 | 115.7
Table 3: Comparison of mass values for the three gaugino masses relevant in
the Point 5 scenario before and after the calibration. The model masses are
also given for reference.

##  6 Conclusions

In this paper we investigated the application of finite impulse response
filters to a typical pattern recognition problem in particle physics, namely
the localisation of edges in kinematic distributions. We showed how to choose
an appropriate filter kernel and its parameters in presence of statistical
fluctuations and irreducible backgrounds in order to optimise the response in
terms of efficiency and localisation power.

Two examples which are highly complementary in signal-to-noise ratio as well
as in steepness of the edge have been studied: pair production of selectrons,
and the production of charginos and neutralinos, decaying to $W$ and $Z$
bosons, respectively. In case of the selectrons, there is an easy-to-find,
very steep edge over a tiny background. The same benchmark has been studied
before with various other techniques to extract the selectron and LSP masses.
We showed that with the FIR method, the edge positions can be extracted with a
twice better precision than with any other method so far investigated. In case
of the charginos and neutralinos, the situation is a priori much more diffcult
due to a much worse signal-to-noise ratio and an edge which is smeared out by
the detector resolution for hadronic jets – and by the finite natural width of
the intermediate gauge boson. Here, other methods, in particular the fit of an
analytic function to the shape of the signal and background parts of the
distribution, fail to give reliable results at all. With the FIR filter
method, the positions of all four edges could be determined reliably to better
than 1% precision. In both cases, the method was calibrated by varying the
input particle masses. After the edge calibration procedure the reconstructed
particle masses agree well with their input values. In case of the selectron
example, precisions of 110 MeV and 90 MeV are obtained on the selectron and
LSP masses, respectively, assuming a data set of $500\,\mathrm{fb^{-1}}$ of
$e^{+}e^{-}$ collisions at a centre-of-mass energy of $\sqrt{s}$=500 GeV and
with beam polarisations of $\mathcal{P}$($e^{-}$, $e^{+}$) = (+80%, -30%). In
case of the gaugino example, the mass of the $\tilde{\chi}^{\pm}_{1}$ was
determined with a precision of 2.2%, while for the two neutralinos precisions
of 1.6% were obtained under the same conditions. In summary, we showed that
the FIR filter technique can be applied very successfully to a variety of
typical particle physics distributions. Other physics analysis could benefit
from applying this edge detection method.

## Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through
the Collaborative Research Centre SFB 676 “Particles, Strings and the Early
Universe”, project B1.

Heister:2002mn,Heister:2001nk,Heister:2002hp
Abbiendi:2002vz,Abbiendi:2003sc,Abbiendi:2002mp Achard:2003ge,Acciarri:1999km
ALEPH ALEPH ALEPH OPAL OPAL OPAL L3 L3

## References

  * [1] Aaron Viets  “Reconstructing the calibrated strain signal in the Advanced LIGO detectors”  In _Class. Quant. Grav._ 35.9, 2018, pp. 095015  DOI: 10.1088/1361-6382/aab658
  * [2] Vincent Torre and Tomaso A. Poggio  “On Edge Detection”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 8.2 Los Alamitos, CA, USA: IEEE Computer Society, 1986, pp. 147–163  DOI: http://doi.ieeecomputersociety.org/10.1109/TPAMI.1986.4767769
  * [3] J. Canny  “A Computational Approach to Edge Detection”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 8.6, 1986, pp. 679–698  DOI: 10.1109/TPAMI.1986.4767851
  * [4] Didier Demigny and Tawfik Kamlé  “A Discrete Expression of Canny’s Criteria for Step Edge Detector Performances Evaluation”  In _IEEE Trans. Pattern Anal. Mach. Intell._ 19.11 Washington, DC, USA: IEEE Computer Society, 1997, pp. 1199–1211  DOI: 10.1109/34.632980
  * [5] Rachid Deriche  “Using Canny’s criteria to derive a recursively implemented optimal edge detector”  In _International Journal of Computer Vision_ 1.2, pp. 167–187  DOI: 10.1007/BF00123164
  * [6] Jun Shen and Serge Castan  “An Optimal Linear Operator for Step Edge Detection”  In _CVGIP: Graph. Models Image Process._ 54.2 Orlando, FL, USA: Academic Press, Inc., 1992, pp. 112–133  DOI: 10.1016/1049-9652(92)90060-B
  * [7] Sabato Stefano Caiazza  “The GridGEM module: a new GEM based readout module for a large TPC & A new algorithm for the measurement of the kinematic edges positions in SUSY decays”, 2018
  * [8] J. Abdallah  “Searches for supersymmetric particles in e+ e- collisions up to 208-GeV and interpretation of the results within the MSSM”  In _Eur. Phys. J._ C31, 2003, pp. 421–479  DOI: 10.1140/epjc/s2003-01355-5
  * [9]
  * [10]
  * [11]
  * [12] I. Hinchliffe et al.  “Precision SUSY measurements at CERN LHC”  In _Phys. Rev._ D55, 1997, pp. 5520–5540  DOI: 10.1103/PhysRevD.55.5520
  * [13] B.. Gjelsten, D.. Miller and P. Osland  “Measurement of SUSY masses via cascade decays for SPS 1a”  In _JHEP_ 12, 2004, pp. 003  DOI: 10.1088/1126-6708/2004/12/003
  * [14] Michael Burns, Konstantin T. Matchev and Myeonghun Park  “Using kinematic boundary lines for particle mass measurements and disambiguation in SUSY-like events with missing energy”  In _JHEP_ 05, 2009, pp. 094  DOI: 10.1088/1126-6708/2009/05/094
  * [15] Konstantin T. Matchev, Filip Moortgat, Luc Pape and Myeonghun Park  “Precise reconstruction of sparticle masses without ambiguities”  In _JHEP_ 08, 2009, pp. 104  DOI: 10.1088/1126-6708/2009/08/104
  * [16] Iftah Galon and Yael Shadmi  “Kinematic Edges with Flavor Splitting and Mixing”  In _Phys. Rev._ D85, 2012, pp. 015010  DOI: 10.1103/PhysRevD.85.015010
  * [17] Morad Aaboud  “Search for new phenomena using the invariant mass distribution of same-flavour opposite-sign dilepton pairs in events with missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions with the ATLAS detector”, 2018  arXiv:1805.11381 [hep-ex]
  * [18] Philip Bambade  “The International Linear Collider: A Global Project”, 2019  arXiv:1903.01629 [hep-ex]
  * [19] Ties Behnke et al.  “The International Linear Collider Technical Design Report - Volume 1: Executive Summary”, 2013  arXiv:1306.6327 [physics.acc-ph]
  * [20] Chris Adolphsen et al.  “The International Linear Collider Technical Design Report - Volume 3.I: Accelerator & in the Technical Design Phase”, 2013  arXiv:1306.6353 [physics.acc-ph]
  * [21] Chris Adolphsen et al.  “The International Linear Collider Technical Design Report - Volume 3.II: Accelerator Baseline Design”, 2013  arXiv:1306.6328 [physics.acc-ph]
  * [22] Toshinori Abe  “The International Large Detector: Letter of Intent”, 2010  DOI: 10.2172/975166
  * [23] Halina Abramowicz  “The International Linear Collider Technical Design Report - Volume 4: Detectors”, 2013  arXiv:1306.6329 [physics.ins-det]
  * [24] Keisuke Fujii  “Physics Case for the 250 GeV Stage of the International Linear Collider”, 2017  arXiv:1710.07621 [hep-ex]
  * [25] T. Barklow et al.  “ILC Operating Scenarios”, 2015  arXiv:1506.07830 [hep-ex]
  * [26] Wolfgang Kilian, Thorsten Ohl and Jurgen Reuter  “WHIZARD: Simulating Multi-Particle Processes at LHC and ILC”  In _Eur. Phys. J._ C71, 2011, pp. 1742  DOI: 10.1140/epjc/s10052-011-1742-y
  * [27] Torbjorn Sjostrand, Stephen Mrenna and Peter Z. Skands  “PYTHIA 6.4 Physics and Manual”  In _JHEP_ 05, 2006, pp. 026  DOI: 10.1088/1126-6708/2006/05/026
  * [28] P. Mora de Freitas and H. Videau  “Detector simulation with MOKKA / GEANT4: Present and future”  In _Linear colliders. Proceedings, International Workshop on physics and experiments with future electron-positron linear colliders, LCWS 2002, Seogwipo, Jeju Island, Korea, August 26-30, 2002_ , 2002, pp. 623–627  URL: http://www-library.desy.de/cgi-bin/showprep.pl?lc-tool03-010
  * [29] S. Agostinelli  “Geant4 - a simulation toolkit”  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 506.3, 2003, pp. 250–303  DOI: https://doi.org/10.1016/S0168-9002(03)01368-8
  * [30] F. Gaede  “Marlin and LCCD: Software tools for the ILC”  In _Advanced computing and analysis techniques in physics research. Proceedings, 10th International Workshop, ACAT05, Zeuthen, Germany, May 22-27, 2005_ A559, 2006, pp. 177–180  DOI: 10.1016/j.nima.2005.11.138
  * [31] J.. Marshall and M.. Thomson  “The Pandora Software Development Kit for Pattern Recognition”  In _Eur. Phys. J._ C75.9, 2015, pp. 439  DOI: 10.1140/epjc/s10052-015-3659-3
  * [32] J.. Marshall, A. Münnich and M.. Thomson  “Performance of Particle Flow Calorimetry at CLIC”  In _Nucl. Instrum. Meth._ A700, 2013, pp. 153–162  DOI: 10.1016/j.nima.2012.10.038
  * [33] M.. Thomson  “Particle Flow Calorimetry and the PandoraPFA Algorithm”  In _Nucl. Instrum. Meth._ A611, 2009, pp. 25–40  DOI: 10.1016/j.nima.2009.09.009
  * [34] Howard Baer and Jenny List  “Post-LHC8 supersymmetry benchmark points for ILC physics”  In _Phys. Rev. D_ 88 American Physical Society, 2013, pp. 055004  DOI: 10.1103/PhysRevD.88.055004
  * [35] Howard Baer, Frank E. Paige, Serban D. Protopopescu and Xerxes Tata  “Simulating Supersymmetry with ISAJET 7.0 / ISASUSY 1.0”  In _Workshop on Physics at Current Accelerators and the Supercollider Argonne, Illinois, June 2-5, 1993_ , 1993, pp. 0703–720  arXiv: http://lss.fnal.gov/cgi-bin/find_paper.pl?other/ssc/sscl-preprint-441
  * [36] Abdelhak Djouadi, Jean-LoÃ¯c Kneur and Gilbert Moultaka  “SuSpect: A Fortran code for the Supersymmetric and Higgs particle spectrum in the MSSM”  In _Computer Physics Communications_ 176.6, 2007, pp. 426–455  DOI: https://doi.org/10.1016/j.cpc.2006.11.009
  * [37] M. Berggren et al.  “Non-simplified SUSY: $\widetilde{\tau}$ -coannihilation at LHC and ILC”  In _Eur. Phys. J._ C76.4, 2016, pp. 183  DOI: 10.1140/epjc/s10052-016-3914-2
  * [38] Werner Porod  “SPheno, a program for calculating supersymmetric spectra, SUSY particle decays and SUSY particle production at e+ e- colliders”  In _Comput. Phys. Commun._ 153, 2003, pp. 275–315  DOI: 10.1016/S0010-4655(03)00222-4
  * [39] W. Porod and F. Staub  “SPheno 3.1: Extensions including flavour, CP-phases and models beyond the MSSM”  In _Comput. Phys. Commun._ 183, 2012, pp. 2458–2469  DOI: 10.1016/j.cpc.2012.05.021
  * [40] Albert M Sirunyan  “Search for supersymmetric partners of electrons and muons in proton-proton collisions at $\sqrt{s}=$ 13 TeV”  In _Submitted to: Phys. Lett._ , 2018  arXiv:1806.05264 [hep-ex]
  * [41] Morad Aaboud  “Search for electroweak production of supersymmetric particles in final states with two or three leptons at $\sqrt{s}=13\,$TeV with the ATLAS detector”, 2018  arXiv:1803.02762 [hep-ex]
  * [42] M. Battaglia et al.  “Physics benchmarks for the ILC detectors”  In _2005 International Linear Collider Workshop : LCWS 2005 : Stanford, California, USA, 18-22 March, 2005_ C050318, 2005, pp. 1602  arXiv:hep-ex/0603010 [hep-ex]
  * [43] Ali H. Chamseddine, Richard L. Arnowitt and Pran Nath  “Locally Supersymmetric Grand Unification”  In _Phys. Rev. Lett._ 49, 1982, pp. 970  DOI: 10.1103/PhysRevLett.49.970
  * [44] Riccardo Barbieri, S. Ferrara and Carlos A. Savoy  “Gauge Models with Spontaneously Broken Local Supersymmetry”  In _Phys. Lett._ 119B, 1982, pp. 343  DOI: 10.1016/0370-2693(82)90685-2
  * [45] M. Tanabashi  “Review of Particle Physics”  In _Phys. Rev. D_ 98.3, 2018, pp. 030001  DOI: 10.1103/PhysRevD.98.030001
  * [46] Taikan Suehara and Jenny List  “Chargino and Neutralino Separation with the ILD Experiment”, 2009  arXiv:0906.5508 [hep-ex]
  * [47] Mikael Berggren  “SGV 3.0 - a fast detector simulation”  In _International Workshop on Future Linear Colliders (LCWS11) Granada, Spain, September 26-30, 2011_ , 2012  arXiv:1203.0217 [physics.ins-det]
  * [48] Madalina Chera  “Particle Flow: From First Principles to Gaugino Property Determination at the ILC”, 2018  DOI: 10.3204/PUBDB-2018-01897
  * [49] B. List, J. List  “An Object-Oriented Kinematic Fitting Package” LC-TOOL-2009-001, 2009
  * [50] A. Heister  “Search for charginos nearly mass degenerate with the lightest neutralino in e+ e- collisions at center-of-mass energies up to 209-GeV”  In _Phys. Lett._ B533, 2002, pp. 223–236  DOI: 10.1016/S0370-2693(02)01584-8
  * [51] A. Heister  “Search for scalar leptons in e+ e- collisions at center-of-mass energies up to 209-GeV”  In _Phys. Lett._ B526, 2002, pp. 206–220  DOI: 10.1016/S0370-2693(01)01494-0
  * [52] A. Heister  “Search for scalar quarks in $e^{+}e^{-}$ collisions at $\sqrt{s}$ up to 209-GeV”  In _Phys. Lett._ B537, 2002, pp. 5–20  DOI: 10.1016/S0370-2693(02)01827-0
  * [53] G. Abbiendi  “Search for nearly mass degenerate charginos and neutralinos at LEP”  In _Eur. Phys. J._ C29, 2003, pp. 479–489  DOI: 10.1140/epjc/s2003-01237-x
  * [54] G. Abbiendi  “Search for chargino and neutralino production at s**(1/2) = 192-GeV to 209 GeV at LEP”  In _Eur. Phys. J._ C35, 2004, pp. 1–20  DOI: 10.1140/epjc/s2004-01758-8
  * [55] G. Abbiendi  “Search for scalar top and scalar bottom quarks at LEP” [Erratum: Phys. Lett.B548,258(2002)]  In _Phys. Lett._ B545, 2002, pp. 272–284  DOI: 10.1016/S0370-2693(02)02808-3, 10.1016/S0370-2693(02)02593-5
  * [56] P. Achard  “Search for scalar leptons and scalar quarks at LEP”  In _Phys. Lett._ B580, 2004, pp. 37–49  DOI: 10.1016/j.physletb.2003.10.010
  * [57] M. Acciarri  “Search for charginos and neutralinos in $e^{+}e^{-}$ collisions at $\sqrt{S}$ = 189-GeV”  In _Phys. Lett._ B472, 2000, pp. 420–433  DOI: 10.1016/S0370-2693(99)01388-X