11institutetext: Centre for Cosmology, Particle Physics and Phenomenology
(CP3),
Université catholique de Louvain, Chemin du Cyclotron 2, B-1348 Louvain-la-
Neuve, Belgium

# MoMEMta, a modular toolkit for the Matrix Element Method at the LHC

Sébastien Brochet  Christophe Delaere  Brieuc François  Vincent Lemaître
Alexandre Mertens  Alessia Saggio  Miguel Vidal Marono  Sébastien Wertz

(January 2019)

###### Abstract

The Matrix Element Method has proven to be a powerful method to optimally
exploit the information available in detector data. Its widespread use is
nevertheless impeded by its complexity and the associated computing time.
MoMEMta, a C++ software package to compute the integrals at the core of the
method, provides a versatile implementation of the Matrix Element Method to
both the theory and experiment communities. Its modular structure covers the
needs of experimental analysis workflows at the LHC without compromising ease
of use on simpler and smaller simulated samples used for phenomenological
studies. With respect to existing tools, MoMEMta improves on usability and
flexibility. In this paper, we present version 1.0 of MoMEMta, together with
examples illustrating the wide range of applications at the LHC accessible for
the first time with a single tool.

††journal: Eur. Phys. J. C

##  1 Introduction

The discovery of the Higgs boson by the ATLAS and CMS experiments in 2012
HiggsAtlas ; HiggsCms  opened a new era in particle physics. More than just a
new particle, a new set of interactions needs to be characterised. The LHC
physics program therefore includes precision measurements of standard model
(SM) processes (in particular in the top-quark and Higgs sectors) and the
search for rare production mechanisms or rare decay channels. The absence so
far of any obvious sign of physics beyond the SM further increases the need to
look in places where the backgrounds are large and the effect of new physics
subtle.

In all these studies, it is of the uttermost importance to fully exploit the
potential of the large data set collected. For most of the analyses performed
in high energy physics (HEP), obtaining an optimal result implies the
treatment of multiple correlated quantities in a multivariate setting. The
most popular methods for multivariate analysis in HEP are machine learning
techniques, such as boosted decision trees and neural networks. These
approaches require large training data sets (usually obtained by Monte Carlo
techniques) in order to learn the structure of the data. On the contrary, the
Matrix Element Method (MEM) uses directly our theoretical knowledge of a
process to assign to each event a probability that measures the compatibility
of experimental data with a given hypothesis. There is no training, since the
underlying Lagrangian, from which the matrix element of the partonic process
is derived, is known.

The MEM, originally designed at the Tevatron experiments DØ and CDF for top
quark mass measurements in $\mathup{{{t}}}\mathup{{\overline{{t}}}}$
production Dalitz:1998zn ; Aaltonen:2008bg ; Aaltonen:2008bd ; Aaltonen:2011di
; Collaboration:2012hz ; Abazov:2015spa ; D0:2016ull , is nowadays a common
technique in particle physics. Recent examples of its use at the LHC are
searches for $\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$
Aad:2015gra ; Aaboud:2017rss ; Khachatryan:2015ila ; Sirunyan:2018mvw ;
Sirunyan:2018ygk ; Aaboud:2017vzb ; Aaboud:2017jvq ; Sirunyan:2018shy  and
single top quark production Aad:2015upn , and a measurement of spin
correlations in $\mathup{{{t}}}\mathup{{\overline{{t}}}}$ production
Khachatryan:2015tzo . Nevertheless, while it can be used for a wide variety of
studies, the practical application of the MEM has been impeded by its
complexity and by the associated computing time. In order to evaluate the
probability under a given theoretical hypothesis of a given experimental
event, a difficult convolution of the theoretical information on the hard
scattering (i.e. the matrix element squared) with the experimentally available
information on the final state (encoded in the so-called transfer functions)
has to be performed. The corresponding integrand varies by several orders of
magnitudes in different regions of the phase space, which requires the use of
adaptive numerical integration techniques together with a smart choice of
integration variables. A general algorithm has been proposed in Ref. Madweight
, which involves optimised phase-space mappings designed to remove as much as
possible the peaks in the integrand. However, the corresponding implementation
(MadWeight) is not supported anymore and suffers from a lack of flexibility
that prevents—or significantly limits—its use in large scale analyses of LHC
data by the collaborations, and does not allow the user to implement
simplifying assumptions.

In this paper, we present MoMEMta, a modular C++ software package to compute
the convolution integrals at the core of the method. Its modular structure
covers the needs of experimental analysis workflows at the LHC without
compromising the ease of use on simpler and smaller simulated samples used for
phenomenological studies. It relies on the same approach as MadWeight to
address the parameterisation of the phase space but leaves more freedom to the
user. Since it follows the same approach, MoMEMta’s performance in terms of
accuracy and CPU time is similar to that of MadWeight. But contrarily to its
predecessor, it adapts to any process and can be fitted to any C++ or Python
analysis workflow. Modularity and flexibility also open the door to specific
optimisations either when designing the integration structure, or when
choosing the integration engine. It is also possible to provide a custom
(optimised) matrix element when appropriate, without loosing all the
advantages of the block decomposition described in Sec. 3.

In the following, we will first briefly review the MEM, with an emphasis on
the assumptions made in MoMEMta, before presenting shortly the philosophy of
the implementation. We will then concentrate on a few concrete use cases that
illustrate the variety of problems that can be tackled using MoMEMta, and how
the modularity can best be exploited to adapt to these problems.

##  2 The matrix element method

The MEM is a technique to calculate the conditional probability density
$P(x|\alpha)$ to observe an experimental event $x$, given a specific
theoretical hypothesis $\alpha$. Details about the method can be found for
example in Ref. Fiedler:2010sg . We will here concentrate on the main aspects.

The likelihood for a partonic final state $y$ to be produced in the hard-
scattering process is proportional to the differential cross section ${\rm
d}\sigma_{\alpha}$ of the corresponding process, given by

$\displaystyle{\rm
d}\sigma_{\alpha}(q_{1},q_{2},y)=\frac{(2\pi)^{4}\\!\left|\mathscr{M}_{\alpha}\left(q_{1},q_{2},y\right)\right|^{2}}{q_{1}q_{2}s}{\rm
d}\Phi(y)\ ,$ (1)

where $q_{1}$ and $q_{2}$ stand for the initial state parton momentum
fractions, $s$ stands for the hadronic centre-of-mass energy and $y$ stands
for the kinematics of the final state.

The central element in that expression is the squared matrix element for
process $\alpha$, denoted
$|\mathscr{M}_{\alpha}\left(q_{1},q_{2},y\right)|^{2}$, where the summation
over spin and colour states is understood. It can be obtained either
analytically or numerically through packages like MG5_aMC@NLO madgraph5  or
MCFM MCFM . Because of the intrinsic theoretical difficulty to identify final
state particles with partons at NLO, the leading order matrix element is used
in most applications. The n-body phase space ${\rm d}\Phi(y)$ must also be
considered in the calculation, as it plays an important role in any change of
variable needed for the integration of the differential cross section.

To obtain the differential cross section ${\rm d}\sigma_{\alpha}(y)$ in hadron
collisions, (1) is convoluted with the parton density functions (PDF) and
summed over all possible flavour compositions of the colliding partons,

${\rm d}\sigma_{\alpha}(y)=\int\limits_{q_{1},q_{2}}\\!\sum_{a_{1},a_{2}}{\rm
d}q_{1}\,{\rm d}q_{2}\,f_{a_{1}}(q_{1})\,f_{a_{2}}(q_{2})\ {\rm
d}\sigma_{\alpha}(q_{1},q_{2},y),$ (2)

where $f_{a_{1}}(q_{1})$ and $f_{a_{2}}(q_{2})$ are the PDFs for a given
flavour $a_{i}$ and momentum fraction $q_{i}$.

The evolution of the parton-level configuration $y$ into a reconstructed event
$x$ in the detector is modelled by a transfer function $T(x|y)$, normalised as
a probability density over $x$, that describes how the partonic final state
$y$ is reconstructed as $x$ in the detector. This includes the effects from
the parton shower, hadronisation, and the limited detector resolution. The
efficiency $\epsilon(y)$, i.e. the probability to reconstruct and select a
specific partonic configuration $y$, also needs to be taken into account. This
includes geometrical acceptance effects.

The transfer function and efficiency are assumed to factorise into
contributions from each measured final-state particle, and each of these
contributions are often assumed to further factorise into simple direction-
and momentum-dependent terms. For most applications, it is realistic to then
assume that particle directions are perfectly reconstructed. The transfer
function in that case is a Dirac delta function on the angular variables, and
takes a non-trivial form only for the energy (or transverse momentum) degree
of freedom. All these assumptions may not be valid in cases where objects,
especially jets, are close to each other or reconstructed together (e.g. in
boosted topologies), which will then result in a less accurate result.

Once the constraints related to the kinematics of observed particles have been
taken into account, there may remain unobserved degrees of freedom, such as
those pertaining to neutrinos (or any other invisible particles) or to
unreconstructed objects outside of detector acceptance, as well as to the
initial-state partons. Some of these degrees of freedom can be removed by
enforcing the conversation of total 4-momentum in the initial and final
states; the remaining ones need to be marginalised, resulting in a potentially
large volume of phase space over which to integrate. Additional constraints
may then be used to reduce this volume, such as mass constraints from
intermediate, narrow resonances in the considered process, assumed to be on
their mass shell, or the experimentally measured total transverse momentum
$\vec{p}^{\textrm{miss}}_{\textrm{T}}$ of the missing particles in the event.
Note that while a resolution function can be built on
$\vec{p}^{\textrm{miss}}_{\textrm{T}}$, it is not accurate to trivially
factorise the transfer function on the two components of
$\vec{p}^{\textrm{miss}}_{\textrm{T}}$ from the terms relative to visible
particles in the final state, since the experimental error in measuring
$\vec{p}^{\textrm{miss}}_{\textrm{T}}$ is correlated with the error made in
measuring all the other particles in the event.

Aspects to be considered in the transfer function and efficiency are the
measurement of the momentum of a particle as well as its (mis-)identification.
This latter point might be relevant for b quarks and $\tau$ leptons, and
allows in principle to combine different partonic final state hypotheses for
the same event.

After convolution, the full expression reads

$\displaystyle P(x|\alpha)$
$\displaystyle=\frac{1}{\sigma_{\alpha}^{\text{vis}}}\int\\!{\rm
d}\sigma_{\alpha}(y)\,T(x|y)\,\epsilon(y)$ (3a)
$\displaystyle=\frac{1}{\sigma_{\alpha}^{\text{vis}}}\int\limits_{q_{1},q_{2}}\\!\sum_{a_{1},a_{2}}\int\limits_{y}\\!{\rm
d}\Phi(y)\,{\rm d}q_{1}{\rm
d}q_{2}f_{a_{1}}(q_{1})\,f_{a_{2}}(q_{2})\,|\mathcal{M}_{\alpha}(q_{1},q_{2},y)|^{2}\,T(x|y)\,\epsilon(y),$
(3b)

where $\sigma_{\alpha}^{\text{vis}}$ is a normalisation factor that ensures
$P(x|\alpha)$ is a probability density over $x$. While that factor can be
computed by explicitly integrating $P(x|\alpha)$ over the phase space of
reconstructed events $x$, it is often more practical to estimate it using a
sample of events simulated under hypothesis $\alpha$, in which case
$\sigma_{\alpha}^{\text{vis}}=\sigma_{\alpha}\cdot\expectationvalue{\epsilon}_{\alpha}$,
where $\sigma_{\alpha}=\int{\rm d}\sigma_{\alpha}(y)$, and
$\expectationvalue{\epsilon}_{\alpha}$ is the average reconstruction and
selection efficiency of the simulated events. Finally, one has also to take
into account the fact that some of the particles measured in the detector
cannot be assigned unambiguously to specific final-state partons. Generally,
all possible combinations have then to be considered and the resulting values
for $P(x|\alpha)$ averaged.

The information contained in (3b) can be exploited in different ways, from the
extraction of the most probable value of theory parameters through a
likelihood maximisation method (see e.g. FerreiradeLima:2017iwx ), for which
the dependence of the normalisation constant $\sigma_{\alpha}^{\text{vis}}$ on
the considered hypothesis has to be properly taken into account, to the bare
use of the integral result without normalisation
$\sigma_{\alpha}^{\text{vis}}$, referred to in the literature as matrix
element weight, $W(x|\alpha)$.

The integral defined in (3b) is typically a small number that varies over
several orders of magnitudes from event to event. It is therefore common to
use instead the _event information_ defined by $I_{\alpha}\equiv-\log
P(x|\alpha)$. When computed from the weight instead of the probability, the
information is only modified by an additive constant, with no consequence in
many applications. We will denote this quantity
$I^{\prime}_{\alpha}\equiv-\log W(x|\alpha)$.

In the limit where all the quantities and functions in (3b) are known with
perfect accuracy, $P(x|\alpha)$ is a likelihood. By the Neyman–Pearson lemma,
the ratio between the likelihoods obtained under two different hypotheses
$\alpha$ and $\alpha^{\prime}$ is the most powerful test statistic to
discriminate one from the other NeymanPearson . Hence, if it can be
implemented, the MEM should provide optimal experimental sensitivity. In
practice, we are limited by the use of leading-order matrix elements, or by
assumptions made in constructing the transfer function and efficiency term.
The quantity (3b) is therefore not a true likelihood, and the Neyman–Pearson
lemma does not strictly apply. For discrimination purposes, it is then common
to use the event information as input of another multivariate method
(typically a boosted decision tree or a neural network).

##  3 Implementation

The MEM is used in HEP by both theoretical and experimental communities with
different purposes and levels of complexity ranging from the evaluation of a
matrix element on a reconstructed event to the precise evaluation of model
parameters (e.g. the top quark mass) through the use of the properly
normalised likelihood derived from (3b). Note that in the former case no
integration process is required and effects related to parton showering,
hadronisation, and finite detector resolution are explicitly neglected.

In order to adapt to these very different use cases, a novel modular design
has been adopted for MoMEMta. The core library is written in C++ and provides
modules for various purposes: to represent and evaluate the matrix element and
parton density functions, to represent and evaluate transfer functions, to
perform changes of variables, to handle the combinatorics of the final state,
etc. That way, every term of (3b) is treated as a module that can be
configured by the user. Weights are computed for a given process by calling
and linking the proper set of modules in a configuration file written in the
Lua scripting language lua . Thanks to this modularity, the user is free to
substitute any module provided with a custom implementation without loosing
the benefits of other parts of the tool. The resulting object can be called
from any C++ or Python code, which means that it seamlessly integrates into
the complex analysis environment of the large experimental collaborations but
can also be used within small programs reading events from files in any format
(e.g. a custom text file, or a file in the Root ROOT , HepMC hepmc , Lhco lhco
or StdHEP stdhep  format). In case the modules shipped out-of-the-box are not
sufficient for a particular application, it is straightforward for the user to
extend MoMEMta’s functionalities by adding new modules handling a specific
task, still profiting from the existing infrastructure provided by the tool.

The computation of the weights requires, in most cases, the evaluation of
multidimensional integrals via adaptive Monte Carlo techniques. The efficiency
in computing these integrals depends on the parameterisation of the phase-
space measure used in the integration. In order to map in an efficient way all
the structures in the integrand, MoMEMta follows the philosophy introduced by
MadWeight Madweight . In this approach, starting from a standard
parameterisation, the phase-space measure is optimised by using a finite
number of _analytic_ transformations over subsets of the integration
variables, called “blocks”. A list of the blocks available in MoMEMta along
with the addressed event topologies, and the integration variables removed and
introduced by the changes of variables, is shown in Tabs. 1 and 2. For
consistency, we have adopted the same terminology as in Ref. Madweight . Since
the considered transformations are nonlinear, specifying the value of the new
integration variables typically yields several solutions for the canonical
variables. The matrix element, transfer function, efficiencies and PDFs all
need to be evaluated on each of those solutions, and summed to define the
final integrand. The first table lists “Main Blocks”, i.e. changes of
variables that allow to integrate out the four-dimensional Dirac delta
function present in the phase-space density term ${\rm d}\Phi(y)$, that
enforces conservation of momentum between the initial and final states. The
second table lists “Secondary Blocks”, i.e. simple changes of variables that
do not remove any degree of freedom. Main and secondary blocks are implemented
as dedicated MoMEMta modules, which can be chained to perform the change of
integration variables that is best suited for the problem at hand. These
modules also take care of computing the jacobian factors required by the
changes of integration variables, to be multiplied with the considered
integrand. Examples of using these blocks are given in Sec. 4.

Table 1:  Set of MoMEMta Main Blocks. Each block performs a specific change of integration variables, and removes four degrees of freedom by enforcing momentum conservation between the initial and final states. The third and fourth columns show the integration variables respectively removed and introduced in each block definition. Block G is introduced in this work, in addition to those originally defined in MadWeight. We denote by $q_{i}$ the Bjorken fractions of the initial-state partons, and by $p_{i}$ the 4-momentum of a final-state particle, parameterised in polar coordinates by $|p_{i}|$, $\theta_{i}$ and $\phi_{i}$. An off-shell particle decaying to a set of on-shell final-state particles is written $s_{i\dots j}(\to p_{i}\dots p_{j})$, and the quantity $s_{i\dots j}$ is defined as $(p_{i}+\dots+p_{j})^{2}$. In block E, $y$ denotes the rapidity of the total partonic system. Removing a particle, $p_{i}$, means removing all three degrees of freedom associated with that particle. Variables that are not explicitly removed are understood to remain present as in the standard polar phase-space parameterisation. Similarly, additional final-state particles not mentioned in the block topology are allowed. Secondary blocks may be used to transform the corresponding phase-space variables.  Main | Topology | Removes… | For
---|---|---|---
block |  |  |
A | $(q_{1},q_{2})\to p_{1}+p_{2}$ |  $q_{1}$, $q_{2}$, $|p_{1}|$, $|p_{2}|$ |
B | $(q_{1},q_{2})\to s_{12}(\to\,p_{1}+p_{2})$ |  $q_{1}$, $q_{2}$, $p_{1}$ | $s_{12}$
C | $(q_{1},q_{2})\to s_{123}\to p_{3}+s_{12}(\to\,p_{1}+p_{2})$ |  $q_{1}$, $q_{2}$, $p_{1}$, $|p_{3}|$ |  $s_{12}$, $s_{123}$
D | $(q_{1},q_{2})\to s_{134}(\to p_{4}+s_{13}(\to\,p_{1}+p_{3}))+$ |  $q_{1}$,$q_{2}$, $p_{1}$, $p_{2}$ |  $s_{13}$, $s_{134}$, $s_{25}$, $s_{256}$
| $s_{256}(\to p_{6}+s_{25}(\to\,p_{2}+p_{5}))$ |  |
E | $(q_{1},q_{2})\to(s_{1234},y)\to s_{13}(\to\,p_{1}+p_{3})+$ |  $q_{1}$, $q_{2}$, $p_{1}$, $p_{2}$ |  $s_{1234}$, $y$, $s_{13}$, $s_{24}$
| $s_{24}(\to\,p_{2}+p_{4})$ |  |
F | $(q_{1},q_{2})\to s_{13}(\to\,p_{1}+p_{3})+s_{24}(\to\,p_{2}+p_{4})$ |  $p_{1}$, $p_{2}$ |  $q_{1}$, $q_{2}$, $s_{13}$, $s_{24}$
G | $(q_{1},q_{2})\to s_{12}(\to p_{1}+p_{2})+s_{34}(\to p_{3}+p_{4})$ |  $q_{1}$, $q_{2}$, $|p_{1}|$, $|p_{2}|$, $|p_{3}|$, $|p_{4}|$ |  $s_{12}$, $s_{34}$
Table 2: Set of MoMEMta Secondary Blocks. Each block performs a specific change of integration variables, acting exclusively on final-state particles. The third and fourth columns show the integration variables respectively removed and introduced by each change of variables. Removing a particle, $p_{i}$, means removing all three degrees of freedom associated with that particle. Variables that are not explicitly removed are understood to remain present as in the standard polar phase-space parameterisation. Secondary | Topology | Removes… | For
---|---|---|---
block |  |  |
A | $s_{1234}\to(s_{123}\to s_{12}(\to\,p_{1}+p_{2})+p_{3})+p_{4}$ | $p_{1}$ |  $s_{1234}$, $s_{123}$, $s_{12}$
B | $s_{123}\to s_{12}(\to\,p_{1}+p_{2})+p_{3}$ |  $|p_{1}|$, $\theta_{1}$ |  $s_{12}$, $s_{123}$
C/D | $s_{12}\to\,p_{1}+p_{2}$ | $|p_{1}|$ | $s_{12}$
E | $s_{123}\to s_{12}(\to\,p_{1}+p_{2})+p_{3}$ |  $|p_{1}|$, $|p_{2}|$ |  $s_{12}$, $s_{123}$

In order to efficiently handle the potentially large combinatorial ambiguity
in the assignment between reconstructed final-state objects and partons in the
matrix element, we have included a dedicated module to average over all
permutations between a given set of particles. This module requires an
additional dimension for the integrated phase space, so that the associated
variable governs which assignment should be used for the computation of the
integrand, and the resulting integral corresponds to the average weight over
the considered permutations. However, compared to a naive averaging of the
possible assignments, this scheme allows adaptive integration algorithms to
concentrate on those yielding the largest contribution to the final result.
For a fixed number of evaluations of the computationally expensive parts of
the integrand, such as the squared matrix element, the precision on the result
is thereby increased.

MoMEMta ships with matrix elements for a few processes, but any leading-order
process handled by MG5_aMC@NLO madgraph5  can be added using a Matrix Element
Exporter plugin provided MoMEMta-MaGMEE . Native support for other matrix
element generators is planned for future releases, but the modular structure
already enables the user to wrap any C++ code that computes a matrix element
to be used with MoMEMta. This novel feature can potentially speed up the
computation by a substantial amount, since the evaluation of the matrix
element largely dominates the computation time.

Parton density functions are obtained from Lhapdf6 Buckley:2014ana  and the
integration is done using the Cuba library Hahn:2004fe , that offers a choice
of four independent routines for multidimensional numerical integration: Vegas
vegas , Suave Hahn:2004fe , Divonne divonne , and Cuhre cuhre1 ; cuhre2 .

The MoMEMta implementation MoMEMtaImpl  is publicly available, is licensed
under the GLPv3, and comes together with an online documentation
momemta.github.io  and tutorials tutorials .

##  4 MEM use cases

One common application of the MEM is parameter estimation, through which one
can extract a parameter of interest by means of likelihood maximisation.
Nowadays, applications of the MEM in high-energy physics are most often
restricted to computing weights $W(x|\alpha)$ under several hypotheses, to
discriminate a signal from one or several backgrounds. MoMEMta fits the needs
for either purpose, since it is specifically designed to efficiently compute
integrals as defined in (3b). In complex situations with several reconstructed
objects and unconstrained degrees of freedom, there is no unique solution to
the problem of efficiently and precisely computing $W(x|\alpha)$, and the user
has to play an active role in defining how to evaluate $W(x|\alpha)$ and at
which accuracy.

In this section we describe a few use cases of the MEM for signal extraction
in LHC analyses. The examples illustrate various levels of complexity, from
the simplest case with a precisely reconstructed final state for which no
integration is needed, to complex final states including six reconstructed and
two unobserved objects. The Lua configurations for each of the examples can be
found together with the MoMEMta tutorials tutorials . For all the examples,
simulated events are generated using MG5_aMC@NLO madgraph5 , Pythia pythia82
and Delphes delphes3 .

The computation times vary by several orders of magnitude among the different
use cases, and strongly depend on the choice of parameters governing the
integration procedure. Indicative performance figures are given in Sec. 5.

###  4.1 Discovery and characterisation of the Higgs boson

The MEM was instrumental for the CMS collaboration in the discovery of the
Higgs boson in the $\rm{H}\to\mathup{{{Z}}}\mathup{{{Z}}}^{*}\to
4\mathit{{\ell}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
channel HiggsCms . Likewise, the characterisation by ATLAS Aad:2014eva ;
Aad:2015mxa  and CMS Sirunyan:2017exp ; Sirunyan:2017tqd  of the discovered
resonance in terms of coupling structure, spin and parity, has relied on
matrix-element techniques as suggested in Refs. Bolognesi:2012mm ;
Anderson:2013afp . In this channel, all final-state particles can be detected
and there are no unobserved degrees of freedom over which to integrate. Given
the good experimental resolution on muon and electron direction and momentum,
it is reasonable to approximate the transfer function in (3b) by
$T(x,y)=\delta(x,y)$. Hence, the integral reduces to a simple evaluation of
the matrix element squared and the PDFs using the measured momenta in the
event. In this framework, dubbed matrix element likelihood analysis (MELA), it
is straightforward to build a discriminating variable between the signal and
the
$\mathit{{q}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{\mkern
4.0mu\overline{\mkern-4.0mu{q}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}\to\mathup{{{Z}}}\mathup{{{Z}}}/\mathup{{{Z}}}\mathup{{{\gamma}}}^{*}\to
4\mathit{{\ell}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
background by considering the matrix elements of these two hypotheses:

$\mathcal{D}_{\text{bkg}}(x)=\left(1+\frac{P(x|\text{bkg})}{P(x|\text{sig})}\right)^{-1}.$
(4)

Similar variables can be constructed to discriminate, for instance, a SM Higgs
boson ($J^{P}=0^{+}$) from a resonance of the same mass but different spin
and/or opposite parity:

$\mathcal{D}_{J^{P}}(x)=\left(1+\frac{P(x|J^{P})}{P(x|0^{+})}\right)^{-1}.$
(5)

Although MoMEMta was designed to handle more complex final states, its
flexibility allows the user to easily implement a MELA-like analysis. To
illustrate this fact, we have simulated events for the
$\rm{gg}\to\rm{H}\to\mathup{{{Z}}}\mathup{{{Z}}}^{*}\to 4\rm{\mu}$ and
$\mathit{{q}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{\mkern
4.0mu\overline{\mkern-4.0mu{q}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}\to\mathup{{{Z}}}\mathup{{{Z}}}/\mathup{{{Z}}}\mathup{{{\gamma}}}^{*}\to
4\rm{\mu}$ processes. The SM Higgs sample, as well as the production and decay
of a resonance of spin/parity $J^{P}=0^{-}$ were generated using the Higgs
characterisation framework higgs_character . With MoMEMta’s plugin for
MG5_aMC@NLO, the corresponding matrix elements can be exported in a format
suitable for MoMEMta. The configuration of MoMEMta in this use case only
requires a single module, which returns the product of the matrix element and
the PDFs evaluated on a given event. The phase-space density term present in
(3b) does not need to be included, since it cancels in the ratios in (4) and
(5). With $P(x|\text{bkg})$, $P(x|\text{sig})$ and $P(x|0^{-})$ computed by
MoMEMta, the discriminant variables $\mathcal{D}_{\text{bkg}}$ and
$\mathcal{D}_{0^{-}}$ can be built. The distributions of these variables, for
the different processes considered, are shown on Fig. 1 after an event
selection closely following the analysis in Ref. Sirunyan:2017tqd . The
discrimination power between the competing hypotheses is comparable to what is
obtained in Refs. Sirunyan:2017exp ; Sirunyan:2017tqd .

Figure 1: Distribution of the $\mathcal{D}_{\text{bkg}}$ (left) and
$\mathcal{D}_{0^{-}}$ (right) variables for the
$\rm{gg}\to\rm{H}\to\rm{ZZ}^{*}\to 4\rm{\mu}$ and
$\rm{q}\bar{\rm{q}}\to\rm{ZZ}/\rm{Z\gamma}^{*}\to 4\rm{\mu}$ processes, where
the resonance $H$ is taken to be the SM Higgs or a pseudoscalar of the same
mass ($0^{-}$). For the right-hand figure, we require
$\mathcal{D}_{\text{bkg}}>0.5$. Note that both the background and the SM Higgs
processes have similar distributions of the $\mathcal{D}_{0^{-}}$
discriminant. All distributions are normalised to unit area.

###  4.2 Charge identification in $\mathup{{{t}}}\mathup{{{W}}}$ production

The MEM has been extensively used in the study of single top quark production
processes at the Tevatron Aaltonen:2009jj ; Abazov:2009ii , and was
instrumental in the most sensitive search for $s$-channel single top
production at the LHC Aad:2015upn . Incidentally, single top and
$\mathup{{{W}}}$ boson associated production ($\mathup{{{t}}}\mathup{{{W}}}$)
provides a good opportunity to showcase MoMEMta’s abilities. This process
features three propagators in the matrix element and, in the case where both
the top quark and $\mathup{{{W}}}$ boson decay leptonically (dilepton
channel), missing information due to the presence of two neutrinos in the
final state.

We consider the charge-conjugate processes,
$\mathup{{{t}}}\mathup{{{W}}^{\scriptstyle{-}}}$ and
$\mathup{{\overline{{t}}}}\mathup{{{W}}^{\scriptstyle{+}}}$, which yield the
same visible final state and have practically the same rate at the LHC. It has
been suggested to measure the CKM matrix element
$|V_{\mathup{{{t}}}\mathup{{{d}}}}|$ at the LHC using the charge asymmetry
between these two processes Alvarez:2017ybk . This requires the ability to
efficiently disentangle them, a task made difficult by the system not being
entirely reconstructible in the dilepton channel.

We thus suggest to construct a MEM-based observable as:

$\mathcal{D}_{\pm}(x)=\frac{W(x|\mathup{{\overline{{t}}}}\mathup{{{W}}^{\scriptstyle{+}}})-W(x|\mathup{{{t}}}\mathup{{{W}}^{\scriptstyle{-}}})}{W(x|\mathup{{\overline{{t}}}}\mathup{{{W}}^{\scriptstyle{+}}})+W(x|\mathup{{{t}}}\mathup{{{W}}^{\scriptstyle{-}}})}.$
(6)

The charge asymmetry can then be defined by counting the number of events for
which either $\mathcal{D}_{\pm}(x)<0$ or $\mathcal{D}_{\pm}(x)>0$.

Figure 2: Normalised distribution of the $\mathcal{D}_{\pm}$ asymmetry
observable, as defined in (6), for the charge-conjugate processes
$\mathup{{{t}}}\mathup{{{W}}^{\scriptstyle{-}}}$ and
$\mathup{{\overline{{t}}}}\mathup{{{W}}^{\scriptstyle{+}}}$. Backgrounds such
as $\mathup{{{t}}}\mathup{{\overline{{t}}}}$ production are expected to be
distributed symmetrically around zero.

Computing the weights
$W(x|\mathup{{\overline{{t}}}}\mathup{{{W}}^{\scriptstyle{+}}})$ and
$W(x|\mathup{{{t}}}\mathup{{{W}}^{\scriptstyle{-}}})$ requires a careful
consideration of the constraints and degrees of freedom at hand. We start by
assuming that the directions of all “visible” objects ($\mathup{{{b}}}$ quark,
leptons) are perfectly reconstructed, so that the transfer function reduces to
factorised parameterisations of the resolution on their energies. Thus, 11
degrees of freedom are present in the system: the longitudinal momentum of the
initial-state partons (2), the energies of the visible particles in the final
state (3), as well as the directions and energies of the two neutrinos (6).
Enforcing conservation of energy and momentum between the initial and final
state will remove four of these, so that we end up with seven dimensions over
which to integrate.

The numerical integration will be most efficient if the integration variables
are mapped to the six peaks in the integrand generated by the top quark and
$\mathup{{{W}}}$ boson propagators and by the transfer functions on the
energies of the visible particles; the remaining degree of freedom can be
chosen freely. This can be easily achieved in MoMEMta by pairing the
“Secondary Block B” with the “Main Block B”.

We identify the decay chain
$\mathup{{{t}}}\to\mathup{{{W}}}_{\mathup{{{t}}}}(\to\mathit{{\nu}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}_{\mathup{{{t}}}}+\mathit{{\ell}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}_{\mathup{{{t}}}})+\mathup{{{b}}}$
with the notation $s_{123}\to s_{12}(\to p_{1}+p_{2})+p_{3}$ for the chosen
secondary block in Tab. 2. The secondary block does not remove any degree of
freedom and simply exchanges the energy and polar direction of
$\mathit{{\nu}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}_{\mathup{{{t}}}}$
for the squared invariant masses of $\mathup{{{t}}}$ and
$\mathup{{{W}}}_{\mathup{{{t}}}}$, which are taken as integration variables.
Since the energies of
$\mathit{{\ell}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}_{\mathup{{{t}}}}$
and $\mathup{{{b}}}$ are associated with peaks in the transfer function, they
should be retained as integration variables, which is straightforward since
the chosen block does not affect these quantities. There remains a free
variable over which to integrate, the azimuthal direction of
$\mathit{{\nu}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}_{\mathup{{{t}}}}$,
which is not directly associated with any peak in the integrand and can be
kept as is. Given fixed values for $s_{123}$, $s_{12}$ and $\phi_{1}$, as well
as $p_{2}$ and $p_{3}$, the block solves the following nonlinear system:

$\displaystyle s_{12}$ $\displaystyle=(p_{1}+p_{2})^{2}$ (7) $\displaystyle
s_{123}$ $\displaystyle=(p_{1}+p_{2}+p_{3})^{2}$ (8) $\displaystyle p_{1}^{2}$
$\displaystyle=m_{1}^{2}=0.$ (9)

Thus, the block returns up to two solutions for the full kinematics of
$p_{1}$, as well as the jacobian factor associated with this change of
variables. Each of these solutions is used as input for the rest of the
computation, described below.

The remaining phase-space variables to be considered are related to the other
$\mathup{{{W}}}(\to\mathit{{\nu}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}+\mathit{{\ell}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}})$
decay as well as to the initial-state partons. This system can be identified
with the topology $(q_{1},q_{2})\to s^{\prime}_{12}(\to
p^{\prime}_{1}+p^{\prime}_{2})$ in Tab. 1. Again, the energy of the charged
lepton
$\mathit{{\ell}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
is chosen as integration variable. The main block removes four degrees of
freedom by enforcing the conservation of total 4-momentum in the initial and
final states, which results in a single variable left to integrate over,
chosen as the squared invariant mass $s^{\prime}_{12}$ of the $\mathup{{{W}}}$
boson, conveniently aligned with the last peak in the integrand. The system
solved by the main block is:

$\displaystyle s^{\prime}_{12}$
$\displaystyle=(p^{\prime}_{1}+p^{\prime}_{2})^{2}$ (10)
$\displaystyle(p^{\prime}_{1})^{2}$ $\displaystyle=(m^{\prime}_{1})^{2}=0$
(11) $\displaystyle p^{\prime}_{1x}$
$\displaystyle=\vec{p}^{\textrm{miss}}_{x}-p_{1x}$ (12) $\displaystyle
p^{\prime}_{1y}$ $\displaystyle=\vec{p}^{\textrm{miss}}_{y}-p_{1y},$ (13)

where $\vec{p}^{\textrm{miss}}_{x,y}$ are the measured components of the
missing transverse momentum along the $x$ and $y$ axes. The block yields up to
two solutions for the neutrino $p^{\prime}_{1}$ and evaluates the jacobian
factor associated with the integration of a four-dimensional Dirac delta as
well as with the change of variables.

Using the solutions obtained for the final-state partonic systems, the
kinematics of the initial-state partons can now be computed using a dedicated
module. Since the measured $\vec{p}^{\textrm{miss}}_{\textrm{T}}$ has been
used to constrain the $\vec{p}_{\textrm{T}}$ of the pair of neutrinos, the
total $\vec{p}_{\textrm{T}}$ of the system is not guaranteed to vanish. As
suggested in Ref. Alwall:2010cq , we thus apply a transverse boost on the
final-state system to a frame of reference where its total
$\vec{p}_{\textrm{T}}$ is zero, compute the longitudinal components of the two
initial-state partons in that frame, and boost the complete system back to the
original laboratory frame. This procedure can be understood as a way to
correct the effect of initial-state radiation in the observed event. Finally,
using the (up to) four obtained solutions for the full partonic system, the
transfer functions, jacobians, PDFs and squared matrix element are evaluated
and the results are summed to define the desired integrand function. Note that
the enhancements in the matrix elements due to the top quark and
$\mathup{{{W}}}$ boson propagators can easily be removed by further well-known
transformations applied to $s_{123}$, $s_{12}$ and $s^{\prime}_{12}$. These
transformations are handled by specialised modules.

The distribution of the $\mathcal{D}_{\pm}$ discriminant is shown on Fig. 2.
About 75% of events from either process can be retained on each side of
$\mathcal{D}_{\pm}=0$, which by symmetry leads to a corresponding mistag rate
of 25%. Depending on the analysis needs, the purity can be further improved at
the cost of efficiency (e.g. we obtain 1.5% mistag rate for 25% efficiency).

The strategy adopted above for the phase-space integration is by no means
unique. We stress that thanks to the modularity of MoMEMta, it is easy for the
user to quickly test alternate approaches. For instance, working in the
narrow-width approximation (NWA) is simply achieved by removing the modules
handling the integration over the propagator invariant masses and fixing these
to chosen pole masses. Thanks to the changes of variables applied, the
kinematic constraints in the system are automatically satisfied. Modifying the
assumptions underlying the transfer functions is equally easy, by configuring,
adding or removing modules representing the finite resolution on the
kinematics of final-state partons. The user might also choose to enforce that
the total $\vec{p}_{\textrm{T}}$ of the partonic system be zero in the
laboratory frame, which is achieved by configuring the main block so that the
$\vec{p}_{\textrm{T}}$ of the neutrino represented by $p^{\prime}_{1}$
balances that of all other final-state particles.

###  4.3 $\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$ production

One of the most successful uses of the MEM at the LHC can be found in the
searches for $\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$
production. The ATLAS and CMS collaborations have applied the MEM in final
states with $\mathup{{{H}}}\to\mathup{{{b}}}\mathup{{\overline{{b}}}}$
Aad:2015gra ; Aaboud:2017rss ; Khachatryan:2015ila ; Sirunyan:2018mvw ;
Sirunyan:2018ygk , and multi-lepton final states with either
$\mathup{{{H}}}\to\mathrm{V}\mathrm{V}^{*}$, where $\mathrm{V}=\mathup{{{W}}}$
or $\mathup{{{Z}}}$, or $\mathup{{{H}}}\to\mathup{{{\tau}}}\mathup{{{\tau}}}$
Aaboud:2017vzb ; Aaboud:2017jvq ; Sirunyan:2018shy .

Here we demonstrate MoMEMta’s ability to efficiently handle processes as
complex as $\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$, featuring
a large final-state multiplicity, several propagator enhancements in the
matrix element, missing information due to neutrinos, and many possible jet-
parton assignments. We consider the channel where the Higgs boson decays to
$\mathup{{{b}}}\mathup{{\overline{{b}}}}$ and both top quarks decay
leptonically, for which the main irreducible background consists of
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{b}}}\mathup{{\overline{{b}}}}$
associated production. The relevance of the MEM in this channel was first
demonstrated in Ref. Artoisenet:2013vfa . We generate samples for signal and
background processes and select events with two opposite-charge leptons and at
least four $\mathup{{{b}}}$-tagged jets.

Weights are computed with MoMEMta under two hypotheses,
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}(\mathup{{{b}}}\mathup{{\overline{{b}}}})$
and
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{b}}}\mathup{{\overline{{b}}}}$.
The strategy adopted to parameterise the phase space is in many ways similar
to that described in Sec. 4.2, and will only be briefly summarised here. The
assumptions related to the transfer function are the same as those considered
for the previous example, with the exception of what concerns the energy of
the two charged leptons, assumed to be perfectly measured. For both
hypotheses, the energies of the $\mathup{{{b}}}$ quarks coming from the decays
of the top quarks are retained as integration variables. In order to reduce
the number of dimensions over which to integrate, we work in the narrow-width
approximation (NWA), by which the $\mathup{{{W}}}$ boson and top quark
propagators are approximated by Dirac delta functions.

The momenta of the two unobserved neutrinos (six degrees of freedom) can be
fixed using four constraints corresponding to the top quark and W boson
invariant masses, as well as by the requirement that their combined transverse
momentum equals the observed transverse missing momentum in the event. This
solving strategy can be implemented by applying the change of variable “Main
Block D” on the standard phase-space parameterisation for the decay products
of the top quarks, and fixing the invariants associated with the top quark and
W boson propagators ($s_{134}$, $s_{256}$, $s_{13}$ and $s_{25}$ in Tab. 1) to
their respective pole masses. The remaining degrees of freedom are handled
differently, depending on the hypothesis:

  * •

$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{b}}}\mathup{{\overline{{b}}}}$:
The standard polar phase-space parameterisation for the two extra
$\mathup{{{b}}}$ quarks is retained, i.e. we integrate over both their
energies.

  * •

$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$: We integrate over the
energy of one of the $\mathup{{{b}}}$ quarks from the Higgs boson decay. The
other quark’s energy is fixed by the requirement that the pair’s invariant
mass be equal to the true Higgs boson mass (NWA). In MoMEMta, this is achieved
by applying the transformation of the “Secondary Block C/D”, and fixing
$s_{12}$ to $m^{2}_{\mathup{{{H}}}}$.

Using the above parameterisation, the peaks in the integrand remain mapped to
the integration variables, and the unobserved degrees of freedom due to the
two neutrinos in the final state are effectively removed. Finally, the
integrand needs to be averaged over every one of the $4!=24$ possible
assignments between jets and partons. This task is efficiently handled by a
dedicated module that concentrates on the assignments dominating the average,
as described in Sec. 3.

Figure 3:  Left: signal information under the signal hypothesis
($\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$). Right:
discriminating variable built from the weights in the signal and background
hypotheses. All distributions are normalised to unity.

Figure 3 shows the normalised distributions of the event information
$I^{\prime}_{\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}}$ and a
discriminating variable defined as

$\mathcal{D}_{\text{sig}}(x)=\left(1+\frac{P(x|\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{b}}}\mathup{{\overline{{b}}}})}{P(x|\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}})}\right)^{-1}.$
(14)

By applying a requirement on $\mathcal{D}_{\text{sig}}$ such that 50% of the
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$ signal is retained,
83% of the $\mathup{{{t}}}\mathup{{\overline{{t}}}}$ \+ jets background can be
rejected. As a comparison, using the invariant masses of pairs of b-tagged
jets in the events, by choosing the pair of jets with mass closest to the true
Higgs boson mass, would only reject 65% of the background for the same signal
efficiency. The distributions shown in Fig. 3 (right) can be compared to those
of Refs. Khachatryan:2015ila ; Sirunyan:2018mvw ; Artoisenet:2013vfa .

##  5 Indicative performance figures

We give approximate performance figures observed when computing weights for
the different use cases presented above. It should be clear that those numbers
are indicative only, as the computation time strongly depends on the
considered hypothesis and the parameters of the numerical integration
algorithm. Furthermore, these results were obtained using the functionalities
available out-of-the-box in MoMEMta, and with matrix elements generated by our
plugin for MG5_aMC@NLO, which means no attempt whatsoever was made towards
optimising the computation for these particular cases. MoMEMta has been
designed with the aim of being flexible, enabling users to implement
simplifications or optimisations fit for their needs. Note that tuning the
parameters of the algorithms used for the numerical integration of the weights
can have a strong impact on both the precision $p$ of the resulting integrals
(which in turns impacts the power of the discriminant built from the weights),
and the overall computation time $T$. Generally, all other things being equal,
the evaluation time $T$ scales roughly as $T\propto p^{-2}$. The computation
of the weights used in Sec. 4.3 was carried out using two different
integration algorithms available in the Cuba library: Vegas and Divonne. The
latter was found to yield substantially shorter completion times, without
compromising the discrimination between signal and background with respect to
the former.

In Tab. 3 we give the average per-event computation times for the weights used
in the examples of Sec. 4, along with the average relative precision on these
weights reported by the integration algorithm. These results were obtained on
a computer cluster with an average per-core HS06
score111https://w3.hepix.org/benchmarking.html of 9.1. Table 4 shows how much
time is spent on the main elements of the computation. These fractions are
indicative and vary from event to event, but show that for complex hypotheses
such as those considered in Sec. 4.3, the bottleneck in the computation is due
to the evaluation of the matrix element.

The memory consumption of MoMEMta is strongly linked to the way the
integration algorithm is configured. In practice, for the examples shown here,
memory consumption was observed never to exceed 200 MB.

Table 3: Average computation time of weights under the different hypotheses used in Sec. 4. The average relative precision on the resulting weights are also given. Hypothesis | Avg. time | Avg. precision
---|---|---
$\mathup{{{g}}}\mathup{{{g}}}\to\rm{H}\to\mathup{{{Z}}}\mathup{{{Z}}}^{*}\to 4\mathup{{{\mu}}}$, $\mathit{{q}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}\mathit{{\mkern 4.0mu\overline{\mkern-4.0mu{q}}}{}_{\mspace{-2.0mu}\scriptstyle{}}^{\mspace{0.0mu}\scriptstyle{}}}\to\mathup{{{Z}}}\mathup{{{Z}}}/\mathup{{{Z}}}\mathup{{{\gamma}}}^{*}\to 4\mathup{{{\mu}}}$ | 0.6 ms | /
$\mathup{{{t}}}\mathup{{{W}}^{\scriptstyle{-}}}$, $\mathup{{\overline{{t}}}}\mathup{{{W}}^{\scriptstyle{+}}}$ | 4.5 s | 2.3%
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}(\mathup{{{b}}}\mathup{{\overline{{b}}}})$ (Vegas) | 140 s | 0.9%
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{b}}}\mathup{{\overline{{b}}}}$ (Vegas) | 700 s | 1.0%
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}(\mathup{{{b}}}\mathup{{\overline{{b}}}})$ (Divonne) | 90 s | 0.5%
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{b}}}\mathup{{\overline{{b}}}}$ (Divonne) | 600 s | 0.6%
Table 4: Indicative shares of computation time due to the various elements entering the evaluation of weights under different hypotheses. The elements shown are the evaluation of the matrix element, the parton distribution functions, the transfer functions, and the generation of the phase space (including the changes of variables introduced in Sec. 3).  Hypothesis | Matrix element | PDF | Transfer functions | Phase space
---|---|---|---|---
$\mathup{{{t}}}\mathup{{{W}}^{\scriptstyle{-}}}$, $\mathup{{\overline{{t}}}}\mathup{{{W}}^{\scriptstyle{+}}}$ | 28% | 14% | 30% | 28%
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}(\mathup{{{b}}}\mathup{{\overline{{b}}}})$ | 93% | 1.3% | 2.7% | 2.6%
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{b}}}\mathup{{\overline{{b}}}}$ | 98% | 0.6% | 0.6% | 0.5%

##  6 Summary

We have presented MoMEMta, a modular software package to compute the
convolution integrals at the core of the MEM. Its modular structure covers the
needs of experimental analysis workflows at the LHC without compromising the
ease of use on simpler and smaller simulated samples used for phenomenological
studies.

The MEM has been used in HEP by both theoretical and experimental communities
with different purposes and levels of complexity ranging from the evaluation
of a matrix element on a reconstructed event to the precise evaluation of
model parameters through the use of the properly normalised likelihood. We
have described a few use cases of the MEM for signal extraction in LHC
analyses showcasing different levels of complexity. From the most simple
implementation, for Higgs boson characterisation in the
$\rm{H}\to\rm{ZZ}^{*}\to
4\mathit{{\ell}{}_{\mspace{-3.0mu}\scriptstyle{}}^{\mspace{-1.0mu}\scriptstyle{}}}$
channel, to complex final states such as
$\mathup{{{t}}}\mathup{{\overline{{t}}}}\mathup{{{H}}}$ production, MoMEMta
has proven to be sufficiently flexible to properly handle these different
situations.

The main advantage of MoMEMta over past and existing tools comes from its
modular design, that greatly improves on usability and flexibility. MoMEMta is
designed to offer a versatile and reusable framework for a wide range of
applications of the MEM. While it is able to cover numerous use cases out of
the box, the modular architecture of MoMEMta also enables users to easily
extend its functionalities to handle situations we have not considered, while
keeping the benefits of the I/O, configuration, and integration framework (a
typical example would be the use of an optimised or simplified matrix element
implementation). As possible future developments, we are considering adding an
interface to other matrix element libraries such as MCFM MCFM  or Sherpa
Gleisberg:2008ta , or enhancing the performance of the integration itself
through the use of vector integrand with modified transfer function to
evaluate the effect of systematic uncertainties, or through the use of
machine-learning inspired integration algorithms Bendavid:2017zhk .

##  7 Acknowledgments

We warmly thank Andrea Giammanco and Olivier Mattelaer for their valued
feedback. This project is funded by FRS-FNRS (Belgian National Scientific
Research Fund) IISN projects 4.4503.17 and 4.4503.16. MoMEMta is part of
AMVA4NP, a project that has received funding from the European Horizon 2020
research and innovation programme under grant agreement №675440. SW is
supported through a FRIA grant by the F.R.S.-FNRS. Computational resources
have been provided by the supercomputing facilities of the Université
catholique de Louvain (CISM/UCL) and the Consortium des Équipements de Calcul
Intensif en Fédération Wallonie Bruxelles (CÉCI) funded by the Fond de la
Recherche Scientifique de Belgique (F.R.S.-FNRS) under convention 2.5020.11.
This work would not have been possible without the help of the MadWeight team.
Special thanks to Matthias Komm who designed our logo.

This is a post-peer-review, pre-copyedit version of an article published in
Eur. Phys. J. C. The final authenticated version is available online at:
http://dx.doi.org/10.1140/epjc/s10052-019-6635-5.

## References

  * (1) G. Aad et al. [ATLAS Collaboration], “Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC,” Phys. Lett. B 716 (2012) 1 doi:10.1016/j.physletb.2012.08.020 [arXiv:1207.7214 [hep-ex]]
  * (2) S. Chatrchyan et al. [CMS Collaboration], “Observation of a new boson at a mass of 125 GeV with the CMS experiment at the LHC,” Phys. Lett. B 716 (2012) 30 doi:10.1016/j.physletb.2012.08.021 [arXiv:1207.7235 [hep-ex]]
  * (3) R. H. Dalitz and G. R. Goldstein, “Test of analysis method for top-antitop production and decay events,” Proc. Roy. Soc. Lond. A 455 (1999) 2803 doi:10.1098/rspa.1999.0428 [arXiv:9802249 [hep-ph]] and references therein
  * (4) T. Aaltonen et al. [CDF Collaboration], “Top Quark Mass Measurement in the $\mathrm{t\bar{t}}$ All Hadronic Channel using a Matrix Element Technique in $\mathrm{p\bar{p}}$ Collisions at $\sqrt{s}=1.96$ TeV,” Phys. Rev. D 79 (2009) 072010 doi:10.1103/PhysRevD.79.072010 [arXiv:0811.1062 [hep-ex]]
  * (5) T. Aaltonen et al. [CDF Collaboration], “Measurement of the top quark mass with dilepton events selected using neuroevolution at CDF,” Phys. Rev. Lett. 102 (2009) 152001 doi:10.1103/PhysRevLett.102.152001 [arXiv:0807.4652 [hep-ex]]
  * (6) T. Aaltonen et al. [CDF Collaboration], “Measurement of the top-quark mass in the lepton+jets channel using a matrix element technique with the CDF II detector,” Phys. Rev. D 84 (2011) 071105 doi:10.1103/PhysRevD.84.071105 [arXiv:1108.1601 [hep-ex]]
  * (7) T. Aaltonen et al. [CDF Collaboration], “Measurements of the Top-quark Mass and the $\mathrm{t\bar{t}}$ Cross Section in the Hadronic $\mathup{{{\tau}}}+$ Jets Decay Channel at $\sqrt{s}=1.96$ TeV,” Phys. Rev. Lett. 109 (2012) 192001 doi:10.1103/PhysRevLett.109.192001 [arXiv:1208.5720 [hep-ex]]
  * (8) V. M. Abazov et al. [D0 Collaboration], “Precision measurement of the top-quark mass in lepton+jets final states,” Phys. Rev. D 91 (2015) 112003 doi:10.1103/PhysRevD.91.112003 [arXiv:1501.07912 [hep-ex]]
  * (9) V. M. Abazov et al. [D0 Collaboration], “Measurement of the Top Quark Mass Using the Matrix Element Technique in Dilepton Final States,” Phys. Rev. D 94 (2016) 032004 doi:10.1103/PhysRevD.94.032004 [arXiv:1606.02814 [hep-ex]]
  * (10) G. Aad et al. [ATLAS Collaboration], “Search for the Standard Model Higgs boson produced in association with top quarks and decaying into $\mathrm{b\bar{b}}$ in pp collisions at $\sqrt{s}=8\leavevmode\nobreak\ \mathrm{TeV}$ with the ATLAS detector,” Eur. Phys. J. C 75 (2015) no. 7, 349 doi:10.1140/epjc/s10052-015-3543-1 [arXiv:1503.05066 [hep-ex]]
  * (11) M. Aaboud et al. [ATLAS Collaboration], “Search for the standard model Higgs boson produced in association with top quarks and decaying into a $\mathrm{b\overline{b}}$ pair in pp collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector,” Phys. Rev. D 97 (2018) 072016 doi:10.1103/PhysRevD.97.072016 [arXiv:1712.08895 [hep-ex]]
  * (12) V. Khachatryan et al. [CMS Collaboration], “Search for a Standard Model Higgs Boson Produced in Association with a Top-Quark Pair and Decaying to Bottom Quarks Using a Matrix Element Method,” Eur. Phys. J. C 75 (2015) no.6, 251 doi:10.1140/epjc/s10052-015-3454-1 [arXiv:1502.02485 [hep-ex]]
  * (13) A. M. Sirunyan et al. [CMS Collaboration], “Search for $\mathrm{t\overline{t}}$H production in the H$\to\mathrm{b\overline{b}}$ decay channel with leptonic $\mathrm{t\overline{t}}$ decays in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” Submitted to JHEP. [arXiv:1804.03682 [hep-ex]]
  * (14) A. M. Sirunyan et al. [CMS Collaboration], “Search for $\mathrm{t}\overline{\mathrm{t}}$H production in the all-jet final state in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” JHEP 1806 (2018) 101 doi:10.1007/JHEP06(2018)101 [arXiv:1803.06986 [hep-ex]]
  * (15) M. Aaboud et al. [ATLAS Collaboration], “Measurement of the Higgs boson coupling properties in the H$\rightarrow\mathrm{ZZ}^{*}\rightarrow 4\ell$ decay channel at $\sqrt{s}$ = 13 TeV with the ATLAS detector,” JHEP 1803 (2018) 095 doi:10.1007/JHEP03(2018)095 [arXiv:1712.02304 [hep-ex]]
  * (16) M. Aaboud et al. [ATLAS Collaboration], “Evidence for the associated production of the Higgs boson and a top quark pair with the ATLAS detector,” Phys. Rev. D 97 (2018) 072003 doi:10.1103/PhysRevD.97.072003 [arXiv:1712.08891 [hep-ex]].
  * (17) A. M. Sirunyan et al. [CMS Collaboration], “Evidence for associated production of a Higgs boson with a top quark pair in final states with electrons, muons, and hadronically decaying $\tau$ leptons at $\sqrt{s}=$ 13 TeV,” JHEP 1808 (2018) 066 doi:10.1007/JHEP08(2018)066 [arXiv:1803.05485 [hep-ex]]
  * (18) G. Aad et al. [ATLAS Collaboration], “Evidence for single top-quark production in the $s$-channel in proton-proton collisions at $\sqrt{s}=$8 TeV with the ATLAS detector using the Matrix Element Method,” Phys. Lett. B 756 (2016) 228–246 doi:10.1016/j.physletb.2016.03.017 [arXiv:1511.05980 [hep-ex]]
  * (19) V. Khachatryan et al. [CMS Collaboration], “Measurement of spin correlations in $\mathrm{t\bar{t}}$ production using the matrix element method in the muon+jets final state in pp collisions at $\sqrt{s}=$ 8 TeV,” Phys. Lett. B 758 (2016) 321 doi:10.1016/j.physletb.2016.05.005 [arXiv:1511.06170 [hep-ex]].
  * (20) P. Artoisenet, V. Lemaitre, F. Maltoni and O. Mattelaer, “Automation of the matrix element reweighting method,” JHEP 1012 (2010) 068 doi:10.1007/JHEP12(2010)068 [arXiv:1007.3300 [hep-ph]]
  * (21) F. Fiedler, A. Grohsjean, P. Haefner and P. Schieferdecker, “The Matrix Element Method and its Application in Measurements of the Top Quark Mass,” Nucl. Instrum. Meth. A 624 (2010) 203 doi:10.1016/j.nima.2010.09.024 [arXiv:1003.1316 [hep-ex]]
  * (22) J. Alwall et al., “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,” JHEP 1407 (2014) 079 doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]]
  * (23) J. M. Campbell, R. K. Ellis and W. T. Giele, “A Multi-Threaded Version of MCFM,” Eur. Phys. J. C 75 (2015) no.6, 246 doi:10.1140/epjc/s10052-015-3461-2 [arXiv:1503.06182 [physics.comp-ph]]
  * (24) D. E. Ferreira de Lima, O. Mattelaer and M. Spannowsky, “Searching for processes with invisible particles using a matrix element-based method,” Phys. Lett. B 787 (2018) 100 doi:10.1016/j.physletb.2018.10.044 [arXiv:1712.03266 [hep-ph]]
  * (25) J. Neyman, E. S. Pearson, “On the problem of the most efficient tests of statistical hypotheses,” Phil. Trans. R. Soc. Lond. A (1933) 231 289–337 doi: 10.1098/rsta.1933.0009
  * (26) R. Ierusalimschy, L. Henrique de Figueiredo,and W. Celes Filho, “Lua — an extensible extension language,” Software: Practice & Experience 26 no. 6 (1996) 635–652 doi:10.1002/(SICI)1097-024X(199606)26:6<635::AID-SPE26>3.0.CO;2-P
  * (27) R. Brun and F. Rademakers, “ROOT - An Object Oriented Data Analysis Framework”, Proceedings AIHENP’96 Workshop, Lausanne, Sep. 1996, Nucl. Inst. & Meth. in Phys. Res. A 389 (1997) 81-86. See also http://root.cern.ch/
  * (28) M. Dobbs and J. B. Hansen, “The HepMC C++ Monte Carlo event record for High Energy”, Physics Comput. Phys. Commun. 134, 41 (2001).
  * (29) J. Thaler, “How to Read LHC Olympics Data Files”, http://madgraph.phys.ucl.ac.be/Manual/lhco.html (2006), Accessed: February 2018.
  * (30) L. Garren, P. Lebrun, “StdHep User Manual”, http://cepa.fnal.gov/psm/stdhep/, Accessed: February 2018.
  * (31) S. Brochet, S. Wertz, J. de Favereau “MoMEMta - MadGraph Matrix Element Exporter”, May 2018, doi:10.5281/zenodo.1250685
  * (32) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr and G. Watt, “LHAPDF6: parton density access in the LHC precision era,” Eur. Phys. J. C 75 (2015) 132 doi:10.1140/epjc/s10052-015-3318-8 [arXiv:1412.7420 [hep-ph]]
  * (33) T. Hahn, “CUBA: A Library for multidimensional numerical integration,” Comput. Phys. Commun. 168 (2005) 78 doi:10.1016/j.cpc.2005.01.010 [arXiv:0404043 [hep-ph]]
  * (34) G.P. Lepage, “A New Algorithm for Adaptive Multidimensional Integration,” J. Comput. Phys. 27 (1978) 192 doi:10.1016/0021-9991(78)90004-9
  * (35) J.H. Friedman, M.H. Wright, “A Nested Partitioning Procedure for Numerical Multiple Integration and Adaptive Importance Sampling,” ACM Trans. Math. Software 7 (1981) 76 doi:10.1145/355934.355939
  * (36) J. Berntsen, T.O. Espelid and A. Genz, “An adaptive algorithm for the approximate calculation of multiple integrals,” ACM Trans. Math. Software 17 (1991) 437–451 doi:10.1145/210232.210233
  * (37) J. Berntsen, T.O. Espelid and A. Genz, “Algorithm 698: DCUHRE: an adaptive multidemensional integration routine for a vector of integrals,” ACM Trans. Math. Software 17 (1991) 452–456 doi:10.1145/210232.210234
  * (38) S. Brochet, S. Wertz, M. Vidal, B. François, A. Saggio, C. Delaere, and V. Lemaître, “Momemta/momemta: 1.0.0”, May 2018, doi:10.5281/zenodo.1250697
  * (39) S. Brochet, S. Wertz, M. Vidal, B. François, A. Saggio, C. Delaere, and V. Lemaître, “The MoMEMta project website”, May 2018, doi:10.5281/zenodo.1250743
  * (40) S. Brochet, S. Wertz, A. Saggio “MoMEMta/Tutorials (Version v1.0.0)”, May 2018, doi:10.5281/zenodo.1250682
  * (41) T. Sjöstrand et al., “An Introduction to PYTHIA 8.2,” Comput. Phys. Commun. 191 (2015) 159 doi:10.1016/j.cpc.2015.01.024 [arXiv:1410.3012 [hep-ph]]
  * (42) J. de Favereau et al., “DELPHES 3, A modular framework for fast simulation of a generic collider experiment,” JHEP 1402 (2014) 057 doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346 [hep-ex]]
  * (43) G. Aad et al. [ATLAS Collaboration], “Measurements of Higgs boson production and couplings in the four-lepton channel in pp collisions at center-of-mass energies of 7 and 8 TeV with the ATLAS detector,” Phys. Rev. D 91 (2015) no. 1, 012006 doi:10.1103/PhysRevD.91.012006 [arXiv:1408.5191 [hep-ex]]
  * (44) G. Aad et al. [ATLAS Collaboration], “Study of the spin and parity of the Higgs boson in diboson decays with the ATLAS detector,” Eur. Phys. J. C 75 (2015) no. 10, 476 Erratum: [Eur. Phys. J. C 76 (2016) no. 3, 152] doi:10.1140/epjc/s10052-015-3685-1, 10.1140/epjc/s10052-016-3934-y [arXiv:1506.05669 [hep-ex]]
  * (45) A. M. Sirunyan et al. [CMS Collaboration], “Measurements of properties of the Higgs boson decaying into the four-lepton final state in pp collisions at $\sqrt{s}=13$ TeV,” JHEP 1711 (2017) 047 doi:10.1007/JHEP11(2017)047 [arXiv:1706.09936 [hep-ex]]
  * (46) A. M. Sirunyan et al. [CMS Collaboration], “Constraints on anomalous Higgs boson couplings using production and decay information in the four-lepton final state,” Phys. Lett. B 775 (2017) 1 doi:10.1016/j.physletb.2017.10.021 [arXiv:1707.00541 [hep-ex]]
  * (47) S. Bolognesi, Y. Gao, A. V. Gritsan, K. Melnikov, M. Schulze, N. V. Tran and A. Whitbeck, “On the spin and parity of a single-produced resonance at the LHC,” Phys. Rev. D 86 (2012) 095031 doi:10.1103/PhysRevD.86.095031 [arXiv:1208.4018 [hep-ph]]
  * (48) I. Anderson et al., “Constraining anomalous HVV interactions at proton and lepton colliders,” Phys. Rev. D 89 (2014) no. 3, 035007 doi:10.1103/PhysRevD.89.035007 [arXiv:1309.4819 [hep-ph]]
  * (49) P. Artoisenet et al., “A framework for Higgs characterisation,” JHEP 1311 (2013) 043 doi:10.1007/JHEP11(2013)043 [arXiv:1306.6464 [hep-ph]]
  * (50) T. Aaltonen et al. [CDF Collaboration], “First Observation of Electroweak Single Top Quark Production,” Phys. Rev. Lett. 103 (2009) 092002 doi:10.1103/PhysRevLett.103.092002 [arXiv:0903.0885 [hep-ex]]
  * (51) V. M. Abazov et al. [D0 Collaboration], “Observation of Single Top Quark Production,” Phys. Rev. Lett. 103 (2009) 092001 doi:10.1103/PhysRevLett.103.092001 [arXiv:0903.0850 [hep-ex]]
  * (52) E. Alvarez, L. Da Rold, M. Estevez and J. F. Kamenik, “Measuring $|V_{\mathrm{td}}|$ at the LHC,” Phys. Rev. D 97 (2018) 033002 doi:10.1103/PhysRevD.97.033002 [arXiv:1709.07887 [hep-ph]]
  * (53) J. Alwall, A. Freitas and O. Mattelaer, “The Matrix Element Method and QCD Radiation,” Phys. Rev. D 83 (2011) 074010 doi:10.1103/PhysRevD.83.074010 [arXiv:1010.2263 [hep-ph]]
  * (54) P. Artoisenet, P. de Aquino, F. Maltoni and O. Mattelaer, “Unravelling $\mathrm{t\overline{t}H}$ via the Matrix Element Method,” Phys. Rev. Lett. 111 (2013) no.9, 091802 doi:10.1103/PhysRevLett.111.091802 [arXiv:1304.6414 [hep-ph]].
  * (55) T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert and J. Winter, “Event generation with SHERPA 1.1,” JHEP 0902 (2009) 007 doi:10.1088/1126-6708/2009/02/007 [arXiv:0811.4622 [hep-ph]]
  * (56) J. Bendavid, “Efficient Monte Carlo Integration Using Boosted Decision Trees and Generative Deep Neural Networks,” [arXiv:1707.00028 [hep-ph]]


# Integrated Sensing and Communication with Multi-Domain Cooperation

Jie Yang, Xi Yang, Chao-Kai Wen, and Shi Jin  Jie Yang, Xi Yang, and Shi Jin
are with the National Mobile Communications Research Laboratory, Southeast
University, Nanjing, China (e-mail: {yangjie;yangxi;jinshi}@seu.edu.cn). Chao-
Kai Wen is with the Institute of Communications Engineering, National Sun Yat-
sen University, Kaohsiung, 804, Taiwan (e-mail:
chaokai.wen@mail.nsysu.edu.tw).

###### Abstract

With the continuous increase of the spectrum and antennas, endogenous sensing
is now possible in the fifth generation and future wireless communication
systems. However, sensing is a highly complex task for a heterogeneous
communication network with massive connections. Seeking multi-domain
cooperation is necessary. In this article, we present an integrated sensing
and communication (ISAC) system that performs active, passive, and interactive
sensing in different stages of communication through hardware and software. We
also propose different methods about how multi-user and multi-frequency band
cooperate to further enhance the ISAC system’s performance. Finally, we
elaborate on the advantages of multi-domain cooperation from the physical
layer to the network layer for the ISAC system.

##  I Introduction

The future wireless communication system is expected to provide ubiquitous
connectivity with ultra-high throughput, ultra-low latency, and ultra-high
reliability. It will also go far beyond communication to realize the ability
to sense, control, and even optimize wireless environments [1]. Because, with
the rapid advent of the intelligent age, more and more new applications emerge
out, such as extended reality, holographic communication, autonomous driving,
smart medical, intelligent industry, _etc._ , which require mass data
transmission, centimeter-level localization, and highly fine-grained
environmental information. To bring the applications above into reality and
fulfill the vision that everything will be sensed, connected, and intelligent,
integrated sensing and communication (ISAC) is anticipated to play a pivotal
role.

The continuous evolution of wireless communication is conspicuously witnessed.
On the one hand, millimeter wave (mmWave) and even Tera-Hertz (THz)
frequencies are further exploited [2], which coincides with the spectrum of
mmWave radar and high-resolution THz imaging radar. Large bandwidth brings
high time resolution. On the other hand, the antenna deployment is further
scaled up, which enables unprecedented angular resolution [3]. Coupled with
ultra-dense network deployment, the future wireless communication system is
enabled a paradigm shift in sensing capabilities. Moreover, reconfigurable
intelligent surfaces [4, 5], artificial intelligence [6, 7], and integrated
terrestrial-aerial-satellite networks [8] are regarded as promising techniques
to boost both sensing and communication performance. It is feasible to
integrate the advanced technologies by sharing the infrastructure and time-
frequency-space resources of wireless communications. Integrating sensing and
communication functions into one hardware system contributes to reducing cost,
power consumption and deployment complexity. The cooperation and mutual
assistance of the two functions can be conveniently achieved by utilizing the
high-throughput and low-latency information sharing ability of the wireless
communication systems [9]. Therefore, ISAC is the future trend of technology
development, which breaks the traditional pattern of the isolated design of
communication and sensing. ISAC can build a bridge between the physical and
digital worlds and lay a solid foundation for the intelligence of everything.

Although considerable advantages of ISAC have been foreseen, deep integration
of sensing and communication still has a long way to go. Start from the
physical layer, frame structure should be designed, and the time-frequency-
space resources should be allocated to make communication and sensing share
the same physical layer processing. Moreover, it is difficult for a single
network, not to mention a single user, to complete such a complex task of
realizing high-quality communication and high-accuracy sensing [10, 11].
Multi-domain cooperation can provide different observations of common targets,
resulting in more accurate sensing decisions. Local sensing results can be
combined to synthesize more comprehensive global results through multi-domain
cooperation. Multi-domain cooperation can also provide an effective way to
overcome the weaknesses of a single network or user. For example, although the
mmWave massive multi-input multi-output network has high sensing resolution,
its coverage is limited, sub-6 GHz can ensure reliable links to assist mmWave,
wireless local area networks (WLANs) and sensor networks can also provide side
information to mmWave networks [12, 13]. Therefore, it is of paramount
significance to develop flexible and efficient cooperation among multi-domain.
Correspondingly, a protocol interface and information fusion scheme should be
developed to enable information sharing across users, access points, and even
heterogeneous networks.

In this article, we propose a series of methods to enhance the ISAC in the
perspective of multi-domain cooperation, including various sensing types
cooperation with the corresponding signal processing and physical layer
design, crowdsourcing-based multi-user cooperation with the corresponding
information fusion and sharing schemes, and different frequencies cooperation
in heterogeneous networks with the corresponding network topology design and
protocol interface development. The proposed methods also concerning how the
sensing result can be used to speed up the establishment of communication
links and reduce the overhead of beam training and channel sounding.

##  II Fundamentals of the ISAC System

In this section, we introduce the fundamental concepts of the ISAC system. The
sensing function of the ISAC system means obtaining the direction, distance,
velocity, position, and even image of the target and map of the propagation
environment. Simultaneous localization and mapping (SLAM) is one of the
promising techniques to achieve the sensing function in the ISAC system. SLAM
mainly exploits the unchanging nature of the position and state of radio
features in the propagation environment. Through accumulating measurements in
consecutive time slots, each radio feature can be mapped to a geometric place
(i.e., sensing); then, the mapped results assist the localization of the
moving user, thereby extending a single-shot localization solution to the
multiple-shot SLAM [14].

Figure 1: Illustration of ISAC systems in the heterogeneous network. Active,
passive, and interactive sensing correspond to beam scanning, channel
sounding, and data transmission in communication processes, respectively.

In general, we can classify the sensing in the ISAC system into three types,
namely, active, passive, and interactive types, as shown in Fig. 1. For active
type, the system sends radio frequency (RF) signals and then captures targets
from the echo signals. The system serves as a conventional radar with sensing
ability in the active type. For passive-type, the system does not send sensing
signals while capturing targets through the received RF signals sent or
reflected from the targets. The system functions as SLAM in the passive type.
For interactive-type, position and environment information is transmitted
among transceivers during the data transmission stage. The system acts as a
traditional communication system to share information in the interactive type.
Existing research focuses on a single type. In fact, the three sensing types
can be simultaneously achieved. Specifically, as shown in Fig. 1, active
sensing, passive sensing, and interactive sensing may correspond to beam
scanning, channel sounding, and data transmission in the communication
processes, respectively. One of the fundamental problems of ISAC is how to
allocate the resource between sensing and communication. Therefore,
theoretical guidance of the trade-off between the performance and the resource
should be established.

In the era of the internet of everything, the ISAC system should provide
services for a large number of users, including vehicles, mobile phones,
internet of things devices, etc. These users can also act as sensing result
providers, thus forming a framework for crowdsourcing. Through multi-user
cooperation, the ISAC system will be cost-effective. In addition, the ISAC
system is expected to involve satellite, unmanned aerial vehicle (UAV), radar,
cellular, WLAN, and sensor networks, which together form a vast heterogeneous
network (Fig. 1). Each kind of network has its advantages and disadvantages.
For example, satellite-based global positioning system (GPS) can provide
position information outdoors, but it is ineffective indoors; WLAN can be
economically deployed indoors, whereas its sensing resolution is limited.
Therefore, cooperation among heterogeneous networks is imperative.

The ultimate goal of the ISAC system is in two aspects. On the one hand, the
cellular system obtains new functions, which can track devices, localize
targets, identify objects, and map radio environments. On the other hand,
communication performance is enhanced by cellular sensing. In particular, we
aim to improve the sensing performance with limited resources, and then
accurate sensing results help reduce the communication overhead in beam
training and channel sounding. For example, when the transceiver positions and
the multi-path propagation environment are known, the effective beam
directions can be predicted, which speeds up establishing the mmWave
communication links. Therefore, we put forward new ideas for the physical
layer and network layer design in ISAC systems, the benefits and challenges
are explained in detail in the following section.

##  III ISAC System with Multi-domain Cooperation

We present a novel ISAC system with multi-domain cooperation, such as
cooperation through active and passive sensing, multi-user, and multi-
frequency band networks. With the proposed ISAC system, diverse sensing types,
numerous terminals, and heterogeneous networks can improve each other.

###  III-A Multi-Sensing Type Cooperation

Figure 2: The active and passive combined ISAC system with the corresponding
frame structure. (a) illustrates the radio feature map, where the gray lines,
blue dots, and green dots represent the scatterers, AP and the corresponding
VAs, and UE positions, respectively. (b) presents the process of active
sensing, where red and yellow dots denote the mapped scatterers. (c) shows an
example of active and passive sensing cooperation in the diffuse reflection
scenario, where the blue line is the scatterer shape obtained from
interpolation, the green dot denotes the unique solution of UE position, and
the red dot is a wrong UE position estimate. (d) gives an example of active
and passive sensing cooperation with specular reflection. The gray ellipse and
blue dot denote the VA confidence region of active and passive sensing,
respectively.

Fig. 2 shows the frame structure of the proposed ISAC system. In the beam
scanning stage, the access point (AP) sweeps a bank of discrete beam
directions with RF signal and the user equipment (UE) detects the amplitudes
of each beam to determine the proper transmission and receive beams.
Meanwhile, the AP can listen to the echo signals to accomplish active sensing.
For example, the time-of-arrival (TOA) of each beam direction can be estimated
from the echo signals. Thus, the surrounding scatterers around the AP can be
built. In the channel sounding stage, the UE may transmit the pilot signals
from several beam directions. The AP can extract the propagation parameters
such as angle-of-arrival (AOA), angle-of-departure (AOD), and TOA of
propagation paths from the received signal. With the estimated propagation
parameters, the AP can perform passive sensing by localizing the UE and
mapping the scatterers between UE and AP. Note that the UE can also perform
beam scanning to sense the surrounding radio environment actively. With the
cooperation between active sensing and passive sensing, a radio feature map
with UE position can be charted. In particular, the legacy scatterers on the
map will be refined, and the new scatterers will be added to the map over
time. The radio feature map not only enables various UE position-based
applications but also improves communication. For example, the feasible UE
position region of the next time slot can be predicted through that of the
current time slot with the velocity and acceleration. Therefore, the overhead
of beam scanning can be reduced by limiting the scanning range to a few
possible beam directions.

An instance of active sensing and passive sensing is depicted in Fig. 2(a)
with line-of-sight (LOS) path and diffuse and specular non-line-of-sight
(NLOS) reflected paths, where the active sensing is performed at AP. In Fig.
2(a), the virtual anchor (VA) represents the mirror of AP with respect to the
scatterer [14]. Because the impact of high-order reflections can be neglected
due to the severe attenuation, only the LOS and the first-order NLOS reflected
paths are considered. For active sensing, the scatterers that surround the AP
are mapped through multiple scanning beams, which in total cover 360 degrees
(for 2-dimensional case). In Fig. 2(b), the dots represent the mapping
results, and the dot with darker color corresponds to the echo signal with a
larger amplitude which implies the better scatterer estimate. Clearly, active
sensing shows the panorama of the radio environment around AP. In contrast,
passive sensing only obtains the radio environment of a few propagation paths,
as shown in Fig. 2(a). Passive sensing is capable of providing gridless AOAs
and AODs of propagation paths, which implies better UE localization
performance through advanced spectral estimation algorithms [15], whereas the
UE localization performance of active sensing is limited because the beam
directions are restricted to a set of grid points. The SLAM can be effectively
realized by the cooperation between active and passive sensing. Next, we
present two cooperation methods and the corresponding advantages for scenarios
with diffuse and specular reflections.

Fig. 2(c) illustrates a general cooperation method regardless of diffuse and
specular reflections, where the active sensing provides the scatterer position
to assist the UE localization in passive sensing. During the beam scanning
stage, a series of discrete reflection points are provided by active sensing,
and the continuous scatterer shape can be obtained through interpolation.
Furthermore, the association of the reflection points of the current time slot
and that of the previous time slots on the radio feature map yields more
accurate scatterer shape. Given the scatterers and the estimates of
propagation parameters (AOD $\varphi$, AOA $\theta$, and TOA $\tau$) of
diffuse NLOS path obtained from the channel sounding stage, the unique
solution of UE position can be estimated. However, the unique solution does
not exist if the scatterer position is unknown, where a wrong UE estimation is
obtained in Fig. 2(c). According to the UE and AP positions, the reflection
point at the scatterer of each propagation path can be obtained and saved in
the radio feature map.

Another cooperation method can be used to enhance the performance of SLAM by
using specular reflections. As shown in Fig. 2(d), the UE, the VA, and the
reflection point are collinear according to the law of reflection. The AP and
VA remain constant even when the UE moves. Given VA, the reflection point at
scatterer can be obtained accordingly using the law of reflection. Therefore,
the environment sensing is equivalent to find the VA. The scatterer shape is
interpolate by the discrete reflection points obtained by active sensing, and
the corresponding VA can be obtained according to the law of reflection.
Specifically, those reflection points that belong to any legacy VA are used to
refine the estimate of this VA through data fusion. New VAs are obtained from
the remaining reflection points. Next, in the channel sounding stage, the
multiple paths with extracted propagation parameters are associated with the
legacy VAs. The SLAM performs the following steps sequentially: 1) UE
position: the LOS path and those reflected paths with associated VA are fused
to estimate the UE position. 2) Legacy VA refinement: after obtaining the UE
position, the legacy VA of each reflected path is refined. 3) New VA: the new
VAs are estimated using the remaining (unused) reflected paths.

Many topics are wide open for the cooperation between active sensing and
passive sensing. First, the diffuse and specular reflections should be
distinguished in order to use the appropriate cooperation method. Second, the
design of high-performance and low-complexity algorithms, such as the
association algorithm between the current scatterers/paths and the legacy
scatterers/VAs, the fusion algorithm of multiple paths for UE localization,
and the fusion algorithm of current and legacy scatterers/VAs for mapping, are
fundamental and critical for SLAM.

###  III-B Multi-User Cooperation

The connectivity of humans and machine continues to increase in the age of the
internet of everything. Thus, the terminals are expected to share the
information for cooperation. However, the potential benefits provided by
multiple terminals are not fully utilized in the conventional system. As
discussed in Section III-A, each UE can sense the local radio environment of
the small surrounding area in the communication process and charts the local
radio feature map in the ISAC system by active and passive sensing. Because
each UE can provide position and trajectory information, the global radio
feature map of a large area can be quickly established and refined using the
data collected by multiple UEs. For the UEs that have overlapped sensing
areas, the cooperation helps each other to refine the shared part of the local
radio feature map. Note that the sensing performance is closely related to the
relative positions among the UE, the scatterer, and the AP. If some UEs detect
a scatterer with high confidence level, other UEs that observe this scatterer
can directly inherit the reliable scatterer estimate.

Figure 3: The crowdsourcing-based ISAC system and the corresponding frame
structure. (a) shows the global radio feature map (RF-Map), which is fused by
four local RF-Map provided by users with different positions and trajectories.
(b) presents an example of cooperation in the specular reflection scenario,
where different users share a common VA. (c) expresses an example of
cooperation in the diffuse reflection scenario, where different users share
common scatterers.

As illustrated in Fig. 3(a), the local radio feature maps obtained from four
UEs catch a subset of radio features (VAs or scatterers) in the region of
interest (ROI). By uploading the local radio feature maps to the cloud
database, a global radio feature map can be built by fusing the local maps.
Because multiple UEs can sense a common ROI, the radio feature that only
exists in a short time is more likely to be caught and recorded from the
cooperation, and the labor-intensive data collection is circumvented in the
multiuser communication. Although the radio features of each local map may
have a low confidence level, the fusion of the crowdsourcing data of a shared
radio feature (recorded by different maps) can obtain a refined radio feature
with high reliability in the global map. Then, the refined radio features of
the global map can be downloaded by the UEs. Each UE continuously refines the
legacy radio features, adds the new radio features, and reports the local map
to the cloud database, while the cloud refines and distributes the shared
features to the corresponding UEs. In addition, each UE can actively detect a
radio feature to check whether it still exists or not. If a radio feature
vanishes, UE can report to the cloud to delete this radio feature in the
global map. For a newly accessed UE, its local map can be initialized by
downloading the selected part of the global map to reduce sensing overhead. We
present the cooperation methods for a common feature shared by different UEs
in the following discussions.

First, we consider the case with common VAs. According to the discussions in
Section III-A, the VA describes the characteristics of scatterer for specular
reflection. In Fig. 3(b), UE 1 and UE 2 share a common VA (scatterer). UE 1
and UE 2 provide different observations to the shared VA from different
positions. Assuming that the AOD estimates ${\hat{\varphi}}_{1}$ and
${\hat{\varphi}}_{2}$ and the TOA estimates ${\hat{\tau}}_{1}$ and
${\hat{\tau}}_{2}$ of paths 1 and 2 are extracted from the received signals.
According to the studies about the multipath channel sounding [7], the
extracted parameter is generally a Gaussian-corrupted observation of true
value such as ${\hat{\varphi}}_{1}=\varphi_{1}+w$, where $w$ is Gaussian noise
with zero mean and certain variance. Therefore, ${\hat{\varphi}}_{1}$ and
${\hat{\tau}}_{1}$ form a 2-dimensional Gaussian-corrupted observation of VA
position. According to the properties of Gaussian distribution, the confidence
region of VA is represented by an ellipse. In practice, the resolutions of AOD
and TOA are usually different. Fig. 3(b) describes the case that the TOA
measurements are more accurate than the AOD measurements, which means the
accuracy in $x^{\prime}$-axis is superior to that of $y^{\prime}$-axis for UE
1. Moreover, ${\hat{\varphi}}_{2}$ and ${\hat{\tau}}_{2}$ form another
2-dimensional Gaussian-corrupted observation from the perspective of UE 2,
where the accuracy of $y^{\prime}$-axis is superior to that of
$x^{\prime}$-axis. Hence, the data fusion of the two confidence regions of the
VA takes the advantages of both sides and yields a much more accurate VA
estimate.

Next, we consider the case with common scatterer. Fig. 3(c) shows an instance
of the propagation environment, where UE 3 and UE 4 perform beam scanning to
sense the local radio environment. As shown in Fig. 3(c), UE 3 and UE 4 share
two common scatterers A and B. Each UE observes a series of reflection points,
and the cooperation between UE 3 and UE 4 can improve the sensing performance
of both UE. Note that the distance between UE 3 and scatterer A is much
smaller than that between UE 4 and scatterer A. This implies that the sensing
performance of scatterer A provided by UE 3 is much better than that provided
by UE 4 because a smaller distance means larger echo signal amplitude.
Therefore, the scatterer A on the local radio feature map of UE 4 can be
improved using the local map of UE 3. Similarly, the information of scatterer
B provided by UE 4 can help UE 3 to refine its local estimate about scatterer
B. Although the local information is incomplete and biased, the cooperation
helps each UE obtain comprehensive global information. The abovementioned
examples show that the cooperation refines the shared information of each UE,
which significantly speeds up the mapping of the global radio feature map.

One of the main challenges for the multi-user cooperation is to measure the
confidence level (reliability) of detected radio features. According to the
confidence level, unreliable radio features can be filtered out to reduce data
processing overhead. In addition, a good pruning algorithm should be designed
to reduce the amount of radio map data transmitted between UE and cloud to
save the overhead of wireless fronthaul and backhaul.

###  III-C Multi-Frequency Band Cooperation

In Sections III-A and III-B, we present the fundamental localization and
mapping methods and the corresponding cooperation strategies in the physical
layer. As discussed above, the different perspectives provide distinct but
incomplete and biased side information, and the cooperation of side
information shows the possibility of revealing the panorama of global
information. Another strategy is to fuse the heterogeneous networks, such as
sub-6 GHz and mmWave cellular systems, WLANs, sensor networks, UAVs,
satellites, radar systems, etc., in the network layer to realize the multi-
frequency band cooperation. The performance of each network can be improved by
taking advantage of multi-domain information. In the following discussions, we
present the cooperation method for sub-6 GHz and mmWave systems.

Figure 4: The cooperation of sub-6 GHz and mmWave in the ISAC systems. (a)
presents an example of the collocated deployment. (b) provides an instance of
the distributed deployment.

The mmWave system offers a high data rate using large bandwidth in high
frequency. However, the severe attenuation and shadowing imply poor quality of
service. Compared with the mmWave communication, the sub-6 GHz system provides
broad coverage and reliable link quality. The future cellular system is
expected to contain both sub-6 GHz and mmWave APs, where the sub-6 GHz APs
ensure seamless communication, and the mmWave APs provide high-quality access
in hot spot areas. The latest release of 5G New Radio (NR) supports dual-
connectivity, where the mobile terminal is able to connect multiple nodes such
as eNB for Long Term Evolution (LTE) and gNB for 5G NR. There are two types of
deployment strategies for sub-6 GHz and mmWave APs: 1) collocated deployment
and 2) distributed deployment. For collocated deployment, the sub-6 GHz and
mmWave APs are located in the same position, where the former mainly covers
the cell-edge terminals, and the latter supports high data rate transmission
for cell-center terminals. The sub-6 GHz and mmWave APs are located in
different positions in a distributed deployment. The mmWave AP serves the
local terminals in the hot spot area as micro AP and the sub-6 GHz AP provides
reliable access to the terminal as macro AP. Considering the propagation
characteristics, different cooperation methods should be used for different
deployments.

For collocated deployment, the sub-6 GHz and mmWave APs share some common
scatterers. However, the propagation phenomenon via shared scatterer is
different in sub-6 GHz and mmWave communications. In practice, the sub-6 GHz
usually observes a bundle of clustering paths, whereas the multipath in mmWave
is much sparser, and the specular reflection can be observed. As shown in Fig.
4(a), for the common scatterer, mmWave transceivers only observe a single
specular reflection path that associates to a VA, whereas sub-6 GHz
transceivers obtain a batch of clustering paths. The means and the spreads of
the AOA, AOD, and TOA of clustering paths with respect to the shared scatterer
can be extracted from the received signal at sub-6 GHz transceivers. In Fig.
4(a), a feasible VA region is obtained using the means and the spreads of
extracted parameters. According to the feasible VA region, the overhead of
beam scanning in mmWave communication can be considerably reduced. This
approach is especially suitable for a newly accessed terminal, where the local
radio feature map is not yet charted, and the reliable radio environment
information is unavailable. After the initial access, the accurate VA estimate
is obtained through mmWave communication and recorded on the local radio
feature map. Given the VA, if the link between the terminal and mmWave AP
undergoes severe attenuation or shadowing, the terminal can switch to the
sub-6 GHz communication and use the extracted parameters to obtain a feasible
terminal region. Once the link quality between the terminal and mmWave AP is
available, the terminal can switch to the mmWave communication and use the
feasible terminal region to reduce the overhead of beam scanning.

In the distributed deployment, each mmWave AP serves the terminals in the
local hot spot area and establishes the local radio feature map. The sub-6 GHz
AP has broad coverage and reliable link quality, taking charge of the up-
loading/down-loading of the radio feature map and the network control. As
shown in Fig. 4(b), the local radio feature maps of mmWave transceivers can be
uploaded to the sub-6 GHz AP, the accurate position and trajectory of the
terminal obtained by mmWave transceivers help the localization of sub-6 GHz AP
and VAs. The sub-6 GHz AP can also generate the global radio feature map by
fusing the local maps and distribute the selected radio features to each
terminal using the wireless fronthaul. In particular, when a terminal enters
the service area of the heterogeneous network, sub-6 GHz AP can quickly
establish the communication link and distribute the local radio feature map to
the terminal. The terminal can quickly access the local mmWave AP by matching
its rough position obtained from GPS with the downloaded radio feature map. In
addition, only the selected radio features of the local map with respect to
the terminal position are transmitted to avoid the overhead of wireless
fronthaul. If the terminal is successfully access to the mmWave AP, the
subsequent map downloading can be completed through mmWave communication.

In the above discussions, we only present two possible cooperation methods for
sub-6 GHz and mmWave systems. It is critical to analyze the specific
correlations between various communication and sensing systems and design a
general data fusion algorithm framework to fuse the correlated information
from multi-domain effectively. In addition, the radio features which may
contain the UE position information are inevitably uploading and downloading
through the wireless fronthaul and backhaul. Hence, how to protect privacy is
a crucial problem for the ISAC system.

##  IV Conclusions and Future Directions

In this article, we have presented a novel ISAC system that performs active,
passive, and interactive sensing in different stages of communication through
hardware and software. We also elaborated on how multi-domain cooperation can
be achieved to enhance the performance of the ISAC system through active and
passive sensing, multi-user, and multi-frequency band networks.

The ISAC system with RIS and AI are potential research directions. On the one
hand, RIS can manipulate electromagnetic waves in a cost-effective and energy-
efficient manner. Therefore, RIS promises to extend the wireless communication
range, facilitate NLOS communications, and provide low-cost cooperative
localization and sensing opportunities. The ISAC system with RIS assistance
can open another dimension in environment sensing. On the other hand, AI has
good learning ability, which can extract useful features from massive data and
alleviate modeling issues of ISAC systems in complex environments. Hybrid data
and model-based AI methods can reduce the training overhead and enhance the
robustness of the ISAC systems.

## References

  * [1] M. Latva-aho and K. Leppänen, Key drivers and research challenges for 6G ubiquitous wireless intelligence. University of Oulu, 2019. [Online]. Available: http://urn.fi/urn:isbn:9789526223544
  * [2] O. Kanhere and T. S. Rappaport, “Position location for futuristic cellular communications: 5G and beyond,” IEEE Commun. Mag., vol. 59, no. 1, pp. 70-75, Jan. 2021.
  * [3] J. Yang, Y. Zeng, S. Jin, C. -K. Wen and P. Xu, “Communication and localization with extremely large lens antenna array,” _IEEE Trans. Wireless Commun._ , Early access.
  * [4] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” IEEE Commun. Mag., vol. 58, no. 1, pp. 106-112, Jan. 2020.
  * [5] W. Tang et al., “Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 1, pp. 421-439, Jan. 2021.
  * [6] A. Zappone, M. Di Renzo, and M. Debbah, “Wireless networks design in the era of deep learning: Model-based, AI-based, or both?” IEEE Trans. Commun., vol. 67, no. 10, pp. 7331-7376, Oct. 2019.
  * [7] J. Yang, S. Jin, C. -K. Wen, J. Guo, M. Matthaiou and B. Gao, “Model-based learning network for 3-D localization in mmWave communications,” _IEEE Trans. Wireless Commun._ , Early access.
  * [8] Y. Zeng and X. Xu, “Toward environment-aware 6G communications via channel knowledge map,” _IEEE Wireless Commun._ , Early access.
  * [9] R. D. Taranto, S. Muppirisetty, R. Raulefs, D. Slock, T. Svensson, and H. Wymeersch, “Location-aware communications for 5G networks: How location information can improve scalability, latency, and robustness of 5G,” _IEEE Signal Process. Mag._ , vol. 31, no. 6, pp. 102-112, Nov. 2014.
  * [10] M. Z. Win, F. Meyer, Z. Liu, W. Dai, S. Bartoletti, and A. Conti, “Efficient multisensor localization for the internet of things: Exploring a new class of scalable localization algorithms,” IEEE Signal Process. Mag. vol. 35, no. 5, pp. 153-167, Sept. 2018.
  * [11] X. Chu, Z. Lu, D. Gesbert, L. Wang, and X. Wen, “Vehicle localization via cooperative channel mapping,” Feb. 2021. [Online]. Available: https://arxiv.org/abs/2102.04852
  * [12] N. Gonzalez-Prelcic, A. Ali, V. Va, and R. W. Heath, “Millimeter-wave communication with out-of-band information,” IEEE Commun. Mag., vol. 55, no. 12, pp. 140-146, Dec. 2017.
  * [13] R. Liu, S. H. Marakkalage, M. Padmal, T. Shaganan, C. Yuen, Y. L. Guan, and U. Tan, “Collaborative SLAM based on WiFi fingerprint similarity and motion information,” _IEEE Internet Things J._ , vol. 7, no. 3, pp. 1826-1840, Mar. 2020.
  * [14] E. Leitinger, F. Meyer, F. Hlawatsch, K. Witrisal, F. Tufvesson, and M. Z. Win, “A belief propagation algorithm for multipath-based SLAM,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 12, pp. 5613-5629, Sept. 2019.
  * [15] Y. Han, T. Hsu, C. Wen, K. Wong and S. Jin, “Efficient downlink channel reconstruction for FDD multi-antenna systems,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 6, pp. 3161-3176, Jun. 2019.