# High Performance Numerical Computing for High Energy Physics: A New Challenge for Big Data Science

Authors: Florin Pop
Journal: Advances in High Energy Physics
 (2014)

---

## Abstract

Modern physics is based on both theoretical analysis and experimental validation. Complex scenarios like subatomic dimensions, high energy, and lower absolute temperature are frontiers for many theoretical models. Simulation with stable numerical methods represents an excellent instrument for high accuracy analysis, experimental validation, and visualization. High performance computing support offers possibility to make simulations at large scale, in parallel, but the volume of data generated by these experiments creates a new challenge for Big Data Science. This paper presents existing computational methods for high energy physics (HEP) analyzed from two perspectives: numerical methods and high performance computing. The computational methods presented are Monte Carlo methods and simulations of HEP processes, Markovian Monte Carlo, unfolding methods in particle physics, kernel estimation in HEP, and Random Matrix Theory used in analysis of particles spectrum. All of these methods produce data-intensive applications, which introduce new challenges and requirements for ICT systems architecture, programming paradigms, and storage capabilities.


---

## Body

## 1. Introduction

High Energy Physics (HEP) experiments are probably the main consumers of High Performance Computing (HPC) in the area of e-Science, considering numerical methods in real experiments and assisted analysis using complex simulation. Starting with quarks discovery in the last century to Higgs Boson in 2012 [1], all HEP experiments were modeled using numerical algorithms: numerical integration, interpolation, random number generation, eigenvalues computation, and so forth. Data collection from HEP experiments generates a huge volume, with a high velocity, variety, and variability and passes the common upper bounds to be considered Big Data. The numerical experiments using HPC for HEP represent a new challenge for Big Data Science.Theoretical research in HEP is related to matter (fundamental particles and Standard Model) and Universe formation basic knowledge. Beyond this, the practical research in HEP has led to the development of new analysis tools (synchrotron radiation, medical imaging or hybrid models [2], wavelets-computational aspects [3]), new processes (cancer therapy [4], food preservation, or nuclear waste treatment), or even the birth of a new industry (Internet) [5].This paper analyzes two aspects: the computational methods used in HEP (Monte Carlo methods and simulations, Markovian Monte Carlo, unfolding methods in particle physics, kernel estimation, and Random Matrix Theory) and the challenges and requirements for ICT systems to deal with processing of Big Data generated by HEP experiments and simulations.The motivation of using numerical methods in HEP simulations is based on special problems which can be formulated using integral or differential-integral equations (or systems of such equations), like quantum chromodynamics evolution of parton distributions inside a proton which can be described by the Gribov-Lipatov-Altarelli-Parisi (GLAP) equations [6], estimation of cross section for a typical HEP interaction (numerical integration problem), and data representation using histograms (numerical interpolation problem). Numerical methods used for solving differential equations or integrals are based on classical quadratures and Monte Carlo (MC) techniques. These allow generating events in terms of particle flavors and four-momenta, which is particularly useful for experimental applications. For example, MC techniques for solving the GLAP equations are based on simulated Markov chains (random walks), which have the advantage of filtering and smoothing the state vector for estimating parameters.In practice, several MC event generators and simulation tools are used. For example, HERWIG (http://projects.hepforge.org/herwig/) project considers angular-ordered parton shower, cluster hadronization (the tool is implemented using Fortran), PYTHIA (http://www.thep.lu.se/torbjorn/Pythia.html) project is oriented on dipole-type parton shower and string hadronization (the tool is implemented in Fortran and C++), and SHERPA (http://projects.hepforge.org/sherpa/) considers dipole-type parton shower and cluster hadronization (the tool is implemented in C++). An important tool for MC simulations is GATE (GEANT4 Application for Tomographic Emission), a generic simulation platform based on GEANT4. GATE provides new features for nuclear imaging applications and includes specific modules that have been developed to meet specific requirements encountered in SPECT (Single Photon Emission Tomography) and PET (Positron Emission Tomography).The main contributions of this paper are as follows:(i)
introduction and analysis of most important modeling methods used in High Energy Physics;
(ii)
identifying and describing of the computational numerical methods for High Energy Physics;
(iii)
presentation of the main challenges for Big Data processing.The paper is structured as follows. Section2 introduces the computational methods used in HEP and describes the performance evaluation of parallel numerical algorithms. Section 3 discusses the new challenge for Big Data Science generated by HEP and HPC. Section 4 presents the conclusions and general open issues.


## 2. Computational Methods Used in High Energy Physics

Computational methods are used in HEP in parallel with physical experiments to generate particle interactions that are modeled using vector of events. This section presents general approach of event generation, simulation methods based on Monte Carlo algorithms, Markovian Monte Carlo chains, methods that describe unfolding processes in particle physics, Random Matrix Theory as support for particle spectrum, and kernel estimation that produce continuous estimates of the parent distribution from the empirical probability density function. The section ends with performance analysis of parallel numerical algorithms used in HEP.


### 2.1. General Approach of Event Generation

The most important aspect in simulation for HEP experiments is event generation. This process can be split into multiple steps, according to physical models. For example, structure of LHC (Large Hadron Collider) events:(
1
) hard process; (
2
) parton shower; (
3
) hadronization; (
4
) underlying event. According to official LHC website (http://home.web.cern.ch/about/computing): “approximately 600 million times per second, particles collide within the LHC
…
Experiments at CERN generate colossal amounts of data. The Data Centre stores it, and sends it around the world for analysis.” The analysis must produce valuable data and the simulation results must be correlated with physical experiments.Figure1 presents the general approach of event generation, detection, and reconstruction. The physical model is used to create simulation process that produces different type of events, clustered in vector of events (e.g., the fourth type of events in LHC experiments).Figure 1
General approach of event generation, detection, and reconstruction.In parallel, the real experiments are performed. The detectors identify the most relevant events and, based on reconstruction techniques, vector of events is created. The detectors can be real or simulated (software tools) and the reconstruction phase combines real events with events detected in simulation. At the end, the final result is compared with the simulation model (especially with generated vectors of events). The model can be corrected for further experiments. The goal is to obtain high accuracy and precision of measured and processed data.Software tools for event generation are based on random number generators. There are three types of random numbers: truly random numbers (from physical generators), pseudorandom numbers (from mathematical generators), and quasirandom numbers (special correlated sequences of numbers, used only for integration). For example, numerical integration using quasirandom numbers usually gives faster convergence than the standard integration methods based on quadratures. In event generation pseudorandom numbers are used most often.The most popular HEP application uses Poisson distribution combined with a basic normal distribution. The Poisson distribution can be formulated as(1)
P
[
X
=
k
]
=
μ
k
k
!
exp
⁡
{
-
μ
}
,
k
=
0,1
,
…
,
with E
(
k
)
=
V
(
k
)
=
μ (V is variance and E is expectation value). Having a uniform random number generator called RND() (Random based on Normal Distribution) we can use the following two algorithms for event generation techniques.The result of running Algorithms1 and 2 to generate around 1
0
6 random numbers is presented in Figure 2. In general, the second algorithm has better result for Poisson distribution. General recommendation for HEP experiments indicates the use of popular random number generators like TRNG (True Random Number Generators), RANMAR (Fast Uniform Random Number Generator used in CERN experiments), RANLUX (algorithm developed by Luscher used by Unix random number generators), and Mersenne Twister (the “industry standard”). Random number generators provided with compilers, operating system, and programming language libraries can have serious problem because they are based on system clock and suffer from lack of uniformity of distribution for large amounts of generated numbers and correlation of successive values.Algorithm 1:Random number generation for Poisson distribution using many random generated numbers with normal distribution (RND).
(1)procedureRa
ndom_Generator_Poisson(μ)
(2)n
u
m
b
e
r
←
-
1;
(3)a
c
c
u
m
u
l
a
t
o
r
←
1.0;
(4)q
←
exp
⁡
{
-
μ
};
(5)while
a
c
c
u
m
u
l
a
t
o
r
>
q
do
(6)r
n
d
_
n
u
m
b
e
r
←
R
N
D
(
);
(7)a
c
c
u
m
u
l
a
t
o
r
←
a
c
c
u
m
u
l
a
t
o
r
*
r
n
d
_
n
u
m
b
e
r;
(8)n
u
m
b
e
r
←
n
u
m
b
e
r
+
1;
(9)end while
(10)return number;
(11)end procedureAlgorithm 2:Random number generation for Poisson distribution using one random generated number with normal distribution.
(1)procedure Ra
ndom_Generator_Poisson_RND(
μ
,
r
)
(2)n
u
m
b
e
r
←
0;
(3)q
←
exp
⁡
{
-
μ
};
(4)a
c
c
u
m
u
l
a
t
o
r
←
q;
(5)p
←
q;
(6)while
r
>
a
c
c
u
m
u
l
a
t
o
r
do
(7)n
u
m
b
e
r
←
n
u
m
b
e
r
+
1;
(8)p
←
p
*
μ
/
n
u
m
b
e
r;
(9)a
c
c
u
m
u
l
a
t
o
r
←
a
c
c
u
m
u
l
a
t
o
r
+
p;
(10)end while
(11)return number;
(12)end procedureGeneral approach of event generation, detection, and reconstruction.
(a)
(b)The art of event generation is to use appropriate combinations of various random number generation methods in order to construct an efficient event generation algorithm being solution to a given problem in HEP.


### 2.2.  Monte Carlo Simulation and Markovian Monte Carlo Chains in HEP

In general, a Monte Carlo (MC) method is any simulation technique that uses random numbers to solve a well-defined problem,P. If F is a solution of the problem P (e.g., F
∈
R
n or F has a Boolean value), we define F
^, an estimation of F, as F
^
=
f
(
{
r
1
,
r
2
,
…
,
r
n
}
,
…
), where {
r
i
}
1
≤
i
≤
n is a random variable that can take more than one value and for which any value that will be taken cannot be predicted in advance. If ρ
(
r
) is the probability density function, ρ
(
r
)
d
r
=
P
[
r
<
r
′
<
r
+
d
r
], the cumulative distributed function is
(2)
C
(
r
)
=
∫
-
∞
r
‍
ρ
(
x
)
d
x
⟹
ρ
(
r
)
=
d
C
(
r
)
d
r
.C
(
r
) is a monotonically nondecreasing function with all values in [
0,1
]. The expectation value is
(3)
E
(
f
)
=
∫
‍
f
(
r
)
d
C
(
r
)
=
∫
‍
f
(
r
)
ρ
(
r
)
d
r
.
And the variance is
(4)
V
(
f
)
=
E
[
f
-
E
(
f
)
]
2
=
E
(
f
2
)
-
E
2
(
f
)
.


#### 2.2.1. Monte Carlo Event Generation and Simulation

To define a MC estimator the “Law of Large Numbers (LLN)” is used. LLN can be described as follows: let one choosen numbers r
i randomly, with the probability density function uniform on a specific interval (
a
,
b
), each r
i being used to evaluate f
(
r
i
). For large n (consistent estimator),
(5)
1
n
∑
i
=
1
n
‍
f
(
r
i
)
⟶
E
(
f
)
=
1
b
-
a
∫
a
b
‍
f
(
r
)
d
r
.The properties of a MC estimator are being normally distributed (with Gaussian density); the standard deviation isσ
=
V
(
f
)
/
n; MC is unbiased for all n (the expectation value is the real value of the integral); the estimator is consistent if V
(
f
)
<
∞ (the estimator converges to the true value of the integral for every large n); a sampling phase can be applied to compute the estimator if we do not know anything about the function f; it is just suitable for integration. The sampling phase can be expressed, in a stratified way, as
(6)
∫
a
b
‍
f
(
r
)
d
r
=
∫
a
r
1
‍
f
(
r
)
d
r
+
∫
r
1
r
2
‍
f
(
r
)
d
r
+
⋯
+
∫
r
n
b
‍
f
(
r
)
d
r
.MC estimations and MC event generators are necessary tools in most of HEP experiments being used at all their steps: experiments preparation, simulation running, and data analysis.An example of MC estimation is the Lorentz invariant phase space (LIPS) that describes the cross section for a typical HEP process withn particle in the final state.Consider(7)
σ
n
~
∫
‍
|
M
|
2
d
R
n
,
where M is the matrix describing the interaction between particles and d
R
n is the element of LIPS. We have the following estimation:
(8)
R
n
(
P
,
p
1
,
p
2
,
…
,
p
n
)
=
∫
‍
δ
(
4
)
(
P
-
∑
k
=
1
n
‍
p
k
)
∏
k
=
1
n
‍
(
δ
(
p
k
2
-
m
k
2
)
Θ
(
p
k
0
)
d
4
p
k
)
,
where P is total four-momentum of the n-particle system; p
kand m
k are four-momenta and mass of the final state particles; δ
(
4
)
(
P
-
∑
k
=
1
n
‍
p
k
) is the total energy momentum conservation; δ
(
p
k
2
-
m
k
2
) is the on-mass-shell condition for the final state system. Based on the integration formula
(9)
∫
‍
δ
(
p
k
2
-
m
k
2
)
Θ
(
p
k
0
)
d
4
p
k
=
d
3
p
k
2
p
k
0
,
obtain the iterative form for cross section:
(10)
R
n
(
P
,
p
1
,
p
2
,
…
,
p
n
)
=
∫
‍
R
n
-
1
(
P
-
p
n
,
p
1
,
p
2
,
…
,
p
n
-
1
)
d
3
p
n
2
p
n
0
,
which can be numerical integrated by using the recurrence relation. As result, we can construct a general MC algorithm for particle collision processes.Example 1.
Let us consider the interaction:e
+
e
-
→
μ
+
μ
- where Higgs boson contribution is numerically negligible. Figure 3 describes this interaction (Φ is the azimuthal angle, θ the polar angle, and p
1
,
p
2
,
q
1
,
q
2 are the four-momenta for particles).
The cross section is(11)
d
σ
=
α
2
4
s
[
W
1
(
s
)
(
1
+
cos
⁡
2
θ
)
+
W
2
(
s
)
cos
⁡
θ
]
d
Ω
,
where d
Ω
=
d
cos
⁡
θ
d
Φ, α
=
e
2
/
4
π (fine structure constant), s
=
(
p
1
0
+
p
2
0
)
2 is the center of mass energy squared, and W
1
(
s
) and W
2
(
s
) are constant functions. For pure processes we have W
1
(
s
)
=
1 and W
2
(
s
)
=
0, and the total cross section becomes
(12)
σ
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
d
2
σ
d
Φ
d
cos
⁡
θ
.Figure 3
Example of particle interaction:e
(
p
1
)
+
e
-
(
p
2
)
→
μ
+
(
q
1
)
μ
-
(
q
2
).We introduce the following notation:(13)
ρ
(
cos
⁡
θ
,
Φ
)
=
d
2
σ
d
Φ
d
cos
⁡
θ
,
and let us consider ρ
~
(
cos
⁡
θ
,
Φ
) an approximation of ρ
(
cos
⁡
θ
,
Φ
). Then σ
~
=
∬
d
Φ
d
cos
⁡
θ
ρ
~. Now, we can compute
(14)
σ
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
ρ
(
cos
⁡
θ
,
Φ
)
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
w
(
cos
⁡
θ
,
Φ
)
ρ
~
(
cos
⁡
θ
,
Φ
)
≈
〈
w
〉
ρ
~
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
ρ
~
(
cos
⁡
θ
,
Φ
)
=
σ
~
〈
w
〉
ρ
~
,
where w
(
cos
⁡
θ
,
Φ
)
=
ρ
(
cos
⁡
θ
,
Φ
)
/
ρ
~
(
cos
⁡
θ
,
Φ
) and 〈
w
〉
ρ
~ is the estimation of w based on ρ
~. Here, the MC estimator is
(15)
〈
w
〉
MC
=
1
n
∑
i
=
1
n
‍
w
i
,
and the standard deviation is
(16)
s
MC
=
(
1
n
(
n
-
1
)
∑
i
=
1
n
‍
(
w
i
-
〈
w
〉
MC
)
2
)
1
/
2
.The final numerical result based on MC estimator is(17)
σ
MC
=
σ
~
〈
w
〉
MC
±
σ
~
s
MC
.As we can show, the principle of a Monte Carlo estimator in physics is to simulate the cross section in interaction and radiation transport knowing the probability distributions (or an approximation) governing each interaction of elementary particles.Based on this result, the Monte Carlo algorithm used to generate events is as follows. It takes as inputρ
~
(
cos
⁡
θ
,
Φ
) and in a main loop considers the following steps: (
1
) generate (
cos
⁡
θ
,
Φ
) peer from ρ
~; (
2
) compute four-momenta p
1
,
p
2
,
q
1
,
q
2; (
3
) compute w
=
ρ
/
ρ
~. The loop can be stopped in the case of unweighted events, and we will stay in the loop for weighted events. As output, the algorithm returns four-momenta for particle for weighted events and four-momenta and an array of weights for unweighted events. The main issue is how to initialize the input of the algorithm. Based on d
σ formula (for W
1
(
s
)
=
1 and W
2
(
s
)
=
0), we can consider as input ρ
~
(
cos
⁡
θ
,
Φ
)
=
(
α
2
/
4
s
)
(
1
+
cos
⁡
2
θ
). Then σ
~
=
4
π
α
2
/
3
s.In HEP theoretical predictions used for particle collision processes modeling (as shown in presented example) should be provided in terms of Monte Carlo event generators, which directly simulate these processes and can provide unweighted (weight = 1) events. A good Monte Carlo algorithm should be used not only for numerical integration [7] (i.e., provide weighted events) but also for efficient generation of unweighted events, which is very important issue for HEP.


#### 2.2.2.  Markovian Monte-Carlo Chains

A classical Monte Carlo method estimates a functionF with F
^ by using a random variable. The main problem with this approach is that we cannot predict any value in advance for a random variable. In HEP simulation experiments the systems are described in states [8]. Let us consider a system with a finite set of possible states S
1
,
S
2
,
…, and S
t the state at the moment t. The conditional probability is defined as
(18)
P
(
S
t
=
S
j
∣
S
t
1
=
S
i
1
,
S
t
2
=
S
i
2
,
…
,
S
t
n
=
S
i
n
)
,
where the mappings (
t
1
,
i
1
)
,
…
,
(
t
n
,
i
n
) can be interpreted as the description of system evolution in time by specifying a specific state for each moment of time.The system is a Markov chain if the distribution ofS
t depends only on immediate predecessor S
t
-
1 and it is independent of all previous states as follows:
(19)
P
(
S
t
=
S
j
∣
S
t
-
1
=
S
i
t
-
1
,
…
,
S
t
2
=
S
i
2
,
S
t
1
=
S
i
1
)
=
P
(
S
t
=
S
j
∣
S
t
-
1
=
S
i
t
-
1
)
.To generate the time steps(
t
1
,
t
2
,
…
,
t
n
) we use the probability of a single forward Markovian step given by p
(
t
∣
t
n
) with the property ∫
t
n
∞
‍
p
(
t
∣
t
n
)
d
t
=
1 and we define p
(
t
)
=
p
(
t
∣
0
). The 1-dimensional Monte Carlo Markovian Algorithm used to generate the time steps is presented in Algorithm 3.Algorithm 3:1-Dimensional Monte Carlo Markovian Algorithm.
(1)   Generatet
1 according with p
(
t
1
)
=
p
(
t
1
∣
t
0
=
0
)
(2)if
t
1
<
t
max
⁡
then             ▹  Generate the initial state.
(3)P
N
≥
1
=
∫
0
t
max
⁡
‍
p
(
t
1
∣
t
0
)
d
t
1;        ▹  Compute the initial probability.
(4)Retain
t
1;
(5)end if
(6)if
t
1
>
t
max
⁡
then
▹  Discard all generated and computed data.
(7)N
=
0; P
0
=
∫
t
max
⁡
∞
‍
p
(
t
1
∣
t
0
)
d
t
1
=
e
-
t
max
⁡;
(8)Delete
t
1;
(9)EXIT.               ▹  The algorithm ends here.
(10)end if
(11)i
=
2;
(12)while (1) do            ▹  Infinite loop until a successful EXIT.
(13)   Generatet
i according with p
(
t
i
∣
t
i
-
1
)
(14)if
t
i
<
t
max
⁡
then
▹  Generate a new state and new probability.
(15)P
N
≥
i
=
∫
t
i
t
max
⁡
‍
p
(
t
i
∣
t
i
-
1
)
d
t
i;
(16)Retain
t
i;
(17)end if
(18)if
t
i
>
t
max
⁡
then
▹  Discard all generated and computed data.
(19)N
=
i
-
1; P
i
=
∫
t
max
⁡
∞
‍
p
(
t
i
∣
t
i
-
1
)
d
t
i;
(20)Retain
(
t
1
,
t
2
,
…
,
t
i
-
1
);Delete
t
i;
(21)EXIT.                  ▹  The algorithm ends here.
(22)end if
(23)i
=
i
+
1;
(24)end whileThe main result of Algorithm3 is that P
(
t
max
⁡
) follows a Poisson distribution:
(20)
P
N
=
∫
0
t
max
⁡
‍
p
(
t
1
∣
t
0
)
d
t
1
×
∫
t
1
t
max
⁡
‍
p
(
t
2
∣
t
1
)
d
t
2
×
⋯
×
∫
t
N
-
1
t
max
⁡
‍
p
(
t
N
∣
t
N
-
1
)
d
t
N
×
∫
t
max
⁡
∞
‍
p
(
t
N
+
1
∣
t
N
)
d
t
N
+
1
=
1
N
!
(
t
max
⁡
)
N
e
-
t
max
⁡
.We can consider the 1-dimensional Monte Carlo Markovian Algorithm as a method used to iteratively generate the systems’ states (codified as a Markov chain) in simulation experiments. According to the Ergodic Theorem for Markov chains, the chain defined has a unique stationary probability distribution [9, 10].Figures4 and 5 present the running of Algorithm 3. According to different values of parameter s used to generate the next step, the results are very different, for 1000 iterations. Figure 4 for s
=
1 shows a profile of the type of noise. For s
=
10,100,1000 profile looks like some of the information is filtered and lost. The best results are obtained for s
=
0.01 and s
=
0.1 and the generated values can be easily accepted for MC simulation in HEP experiments.Example of 1-dimensional Monte Carlo Markovian algorithm.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
Analysis of acceptance rate for 1-dimensional Monte Carlo Markovian algorithm for differents values.Figure5 shows the acceptance rate of values generated with parameter s used in the algorithm. And parameter values are correlated with Figure 4. Results in Figure 5 show that the acceptance rate decreases rapidly with increasing value of parameter s. The conclusion is that values must be kept small to obtain meaningful data. A correlation with the normal distribution is evident, showing that a small value for the mean square deviation provides useful results.


#### 2.2.3. Performance of Numerical Algorithms Used in MC Simulations

Numerical methods used to compute MC estimator use numerical quadratures to approximate the value of the integral for functionf on a specific domain by a linear compilation of function values and weights {
w
i
}
1
≤
i
≤
m as follows:
(21)
∫
a
b
‍
f
(
r
)
d
r
=
∑
i
=
1
m
‍
w
i
f
(
r
i
)
.We can consider a consistent MC estimatora a classical numerical quadrature with all w
i
=
1. Efficiency of integration methods for 1 dimension and for d dimensions is presented in Table 1. We can conclude that quadrature methods are difficult to apply in many dimensions for variate integration domains (regions) and the integral is not easy to be estimated.Table 1
Efficiency of integration methods for 1 dimension and ford dimensions.
Method
1 dimension
d dimensions
Monte Carlo
n
-
1
/
2
n
-
1
/
2
Trapezoidal rule
n
-
2
n
-
2
/
d
Simpson’s rule
n
-
4
n
-
4
/
d
m-points Gauss rule (m
<
n)
n
-
2
m
n
-
2
m
/
dAs practical example, in a typical high-energy particle collision there can be many final-state particles (even hundreds). If we haven final state particle, we face with d
=
3
n
-
4 dimensional phase space. As numerical example, for n
=
4 we have d
=
8 dimensions, which is very difficult approach for classical numerical quadratures.Full decomposition integration volume for one double number (10 Bytes) per volume unit isn
d
×
10 Bytes. For the example considered with d
=
8 and n
=
10 divisions for interval [
0,1
] we have, for one numerical integration,
(22)
n
d
×
10
Bytes
=
1
0
8
×
10
102
4
3
G
Bytes
≈
0.93
G
Bytes
.
Considering 1
0
6 events per second, one integration per event, the data produced in one hour will be ≈3197.4 P Bytes.The previous assumption is only for multidimensional arrays. But due to the factorization assumption,p
(
r
1
,
r
2
,
…
,
r
n
)
=
p
(
r
1
)
p
(
r
2
)
⋯
p
(
r
n
), we obtain for one integration
(23)
n
×
d
×
10
Bytes
=
800
Bytes
,
which means ≈2.62 T Bytes of data produce for one hour of simulations.


### 2.3. Unfolding Processes in Particle Physics and Kernel Estimation in HEP

In particle physics analysis we have two types of distributions: true distribution (considered in theoretical models) and measured distribution (considered in experimental models, which are affected by finite resolution and limited acceptance of existing detectors). A HEP interaction process starts with a true knows distribution and generate a measured distribution, corresponding to an experiment of a well-confirmed theory. An inverse process starts with a measured distribution and tries to identify the true distribution. These unfolding processes are used to identify new theories based on experiments [11].


#### 2.3.1. Unfolding Processes in Particle Physics

The theory of unfolding processes in particle physics is as follows [12]. For a physics variable t we have a true distribution f
(
t
) mapped in x and an n-vector of unknowns and a measured distribution g
(
s
) (for a measured variable s) mapped in an m-vector of measured data. A response matrix A
∈
R
m
×
n encodes a Kernel function K
(
s
,
t
) describing the physical measurement process [12–15]. The direct and inverse processes are described by the Fredholm integral equation [16] of the first kind, for a specific domain Ω,
(24)
∫
Ω
‍
K
(
s
,
t
)
f
(
t
)
d
t
=
g
(
s
)
.
In particle physics the Kernel function K
(
s
,
t
) is usually known from a Monte Carlo sample obtained from simulation. A numerical solution is obtained using the following linear equation: A
x
=
b. Vectors x and y are assumed to be 1-dimensional in theory, but they can be multidimensional in practice (considering multiple independent linear equations). In practice, also the statistical properties of the measurements are well known and often they follow the Poisson statistics [17]. To solve the linear systems we have different numerical methods.First method is based on linear transformationx
=
A
#
y. If m
=
n then A
#
=
A
-
1 and we can use direct Gaussian methods, iterative methods (Gauss-Siedel, Jacobi or SOR), or orthogonal methods (based on Householder transformation, Givens methods, or Gram-Schmidt algorithm). If m
>
n (the most frequent scenario) we will construct the matrix A
#
=
(
A
T
A
)
-
1
A
T (called pseudoinverse Penrose-Moore). In these cases the orthogonal methods offer very good and stable numerical solutions.Second method considers the singular value decomposition:(25)
A
=
U
Σ
V
T
=
∑
i
=
1
n
‍
σ
i
u
i
v
i
T
,
where U
∈
R
m
×
n and V
∈
R
n
×
n are matrices with orthonormal columns and the diagonal matrix Σ
=
diag
⁡
{
σ
1
,
…
,
σ
n
}
=
U
T
A
V. The solution is
(26)
x
=
A
#
y
=
V
Σ
-
1
(
U
T
y
)
=
∑
i
=
1
n
‍
1
σ
i
(
u
i
T
y
)
v
i
=
∑
i
=
1
n
‍
1
σ
i
c
i
v
i
,
where c
i
=
u
i
T
y, i
=
1
,
…
,
n, are called Fourier coefficients.


#### 2.3.2.  Random Matrix Theory

Analysis of particle spectrum (e.g., neutrino spectrum) faces with Random Matrix Theory (RMT), especially if we consider anarchic neutrino masses. The RMT means the study of the statistical properties of eigenvalues of very large matrices [18]. For an interaction matrix A (with size N), where A
i
j is an independent distributed random variable and A
H is the complex conjugate and transpose matrix, we define M
=
A
+
A
H, which describes a Gaussian Unitary Ensemble (GUE). The GUE properties are described by the probability distribution P
(
M
)
d
M: (
1
) it is invariant under unitary transformation, P
(
M
)
d
m
=
P
(
M
′
)
d
M
′, where M
′
=
U
H
M
U, U is a Hermitian matrix (U
H
U
=
I); (
2
) the elements of M matrix are statistically independent, P
(
M
)
=
∏
i
≤
j
‍
P
i
j
(
M
i
j
); and (
3
) the matrix M can be diagonalized as M
=
U
D
U
H, where U
=
diag
⁡
{
λ
1
,
…
,
λ
N
}, λ
i is the eigenvalue of M and λ
i
≤
λ
j if i
<
j
(27)
Propability
(
2
)

:
P
(
M
)
d
M
~
d
M
exp
⁡
{
-
N
2
T
r
(
M
H
M
)
}
;
Propability
(
3
)

:
P
(
M
)
d
M
~
d
U
∏
i
d
λ
i
∏
i
<
j
(
λ
i
-
λ
j
)
2
×
exp
⁡
{
-
N
2
∑
i
‍
(
λ
i
2
)
}
.The numerical methods used for eigenvalues computation are the QR method and Power methods (direct and indirect). The QR method is a numerical stable algorithm and Power method is an iterative one. The RMT can be used for many body systems, quantum chaos, disordered systems, quantum chromodynamics, and so forth.


#### 2.3.3.  Kernel Estimation in HEP

Kernel estimation is a very powerful solution and relevant method for HEP when it is necessary to combine data from heterogeneous sources like MC datasets obtained by simulation and from Standard Model expectation, obtained from real experiments [19]. For a set of data {
x
i
}
1
≤
i
≤
n with a constant bandwidth h (the difference between two consecutive data values), called the smoothing parameter, we have the estimation
(28)
f
^
(
x
)
=
1
n
h
∑
i
=
1
n
‍
K
(
x
-
x
i
h
)
,
where K is an estimator. For example, a Gauss estimator with mean μ and standard deviation σ is
(29)
K
(
x
)
=
1
σ
2
π
exp
⁡
{
-
(
x
-
μ
)
2
2
σ
2
}
,
and has the following properties: positive definite and infinitely differentiable (due to the exp function), and it can be defined for an infinite supports (n
→
∞). The kernel is a nonparametric method, which means that h is independent of dataset and for large amount of normally distributed data we can find a value for h that minimizes the integrated squared error of f
^
(
x
). This value for bandwidth is computed as
(30)
h
*
=
(
4
3
n
)
1
/
5
σ
.The main problem in Kernel Estimation is that the set of data{
x
i
}
1
≤
i
≤
n is not normally distributed and in real experiments the optimal bandwidth it is not known. An improvement of presented method considers adaptive Kernel Estimation proposed by Abramson [20], where h
i
=
h
/
f
(
x
i
) and σ are considered global qualities for dataset. The new form is
(31)
f
^
a
(
x
)
=
1
n
∑
i
=
1
n
‍
1
h
i
K
(
x
-
x
i
h
i
)
,
and the local bandwidth value that minimizes the integrated squared error of f
^
a
(
x
) is
(32)
h
i
*
=
(
4
3
n
)
1
/
5
σ
f
^
(
x
i
)
,
where f
^ is the normal estimator.Kernel estimation is used for event selection to confidence level evaluation, for example, in Markovian Monte Carlo chains or in selection of neural network output used in experiments for reconstructed Higgs mass. In general, the main usage of Kernel estimation in HEP is searching for new particle, by finding relevant data in a large dataset.A method based on Kernel estimation is the graphical representation of datasets using advanced shifted histogram algorithm (ASH). This is a numerical interpolation for large datasets with the main aim of creating a set ofn
b
i
n histograms H
=
{
H
i
}, with the same bin-width h. Algorithm 4 presents the steps of histograms generation starting with a specific interval [
a
,
b
], a number of points n in this interval, and a number of bins and a number of values used for kernel estimation, m. Figure 6 shows the results of kernel estimation if function f
=
-
(
1
/
2
)
x
2 on [
0,1
] and graphical representation with a different number of bins. The values on vertical axis are aggregated in step 17 of Algorithm 4 and increase with the number of bins.Algorithm 4:Advanced shifted histogram (1D algorithm).
(1)procedure
ASH
(
a
,
b
,
n
,
x
,
n
bin
,
m
,
f
m
)
(2)δ
=
(
b
-
a
)
/
n
bin; h
=
m
δ;
(3)for
k
=
1
…
n
bin
do
(4)v
k
=
0; y
k
=
0;
(5)end for
(6)for
i
=
1
…
n
do
(7)k
=
(
x
i
-
a
)
/
δ
+
1;
(8)if
k
∈
[
1
,
n
bin
]
then
(9)v
k
=
v
k
+
1;
(10)end if
(11)end for
(12)for
k
=
1
…
n
bin
do
(13)if
v
k
=
0
then
(14)k
=
k
+
1;
(15)end if
(16)for
i
=
max
⁡
{
1
,
k
-
m
+
1
}
…
min
⁡
{
n
bin
,
k
+
m
-
1
}
do
(17)y
i
=
y
i
+
v
k
f
m
(
i
-
k
);
(18)end for
(19)end for
(20)for
k
=
1
…
n
bin
do
(21)y
k
=
y
k
/
(
n
h
);
(22)t
k
=
a
+
(
k
-
0.5
)
δ;
(23)end for
(24)return
{
t
k
}
1
≤
k
≤
n
bin, {
y
k
}
1
≤
k
≤
n
bin.
(25)end procedureExample of advanced shifted histogram algorithm running for different bins: 10, 100, and 1000.
(a)
(b)
(c)


#### 2.3.4.  Performance of Numerical Algorithms Used in Particle Physics

All operations used in presented methods for particle physics (Unfolding Processes, Random Matrix Theory, and Kernel Estimation) can be reduced to scalar products, matrix-vector products, and matrix-matrix products. In [21] the design of new standard for the BLAS (Basic Linear Algebra Subroutines) in C language by extension of precision is described. This permits higher internal precision and mixed input/output types. The precision allows implementation of some algorithms that are simpler, more accurate, and sometimes faster than possible without these features. Regarding the precision of numerical computing, Dongarra and Langou established in [22] an upper bound for the residual check for A
x
=
y system, with A
∈
R
n
×
n a dense matrix. The residual check is defined as
(33)
r
∞
=
∥
A
x
-
y
∥
∞
n
ϵ
(
∥
A
∥
∞
∥
x
∥
∞
+
∥
y
∥
∞
)
<
16
,
where ϵ is the relative machine precision for the IEEE representation standard; ∥
y
∥
∞ is the infinite norm of a vector: ∥
y
∥
∞
=
max
⁡
1
≤
i
≤
n
{
|
y
i
|
}; and ∥
A
∥
∞ is the infinite norm of a matrix ∥
A
∥
∞
=
max
⁡
1
≤
i
≤
n
{
∑
j
=
1
n
‍
|
A
i
j
|
}.Figure7 presents the graphical representation of Dongarras result (using logarithmic scales) for simple and double precision. For simple precision, ϵ
s
=
2
-
24, for all n
≥
1.05
×
1
0
6 the residual check is always lower than imposed upper bound, similarly for double precision with ϵ
d
=
2
-
53, for all n
≥
5.63
×
1
0
14. If matrix size is greater than these values, it will not be possible to detect if the solution is correct or not. These results establish upper bounds for data volume in this model.Figure 7
Residual check analysis for solvingA
x
=
y system in HPL2.0 using simple and double precision representation.In a single-processor system, the complexity of algorithms depends only on the problem size,n. We can assume T
(
n
)
=
Θ
(
f
(
n
)
), where f
(
n
) is a fundamental function (f
(
n
)
∈
{
1
,
n
α
,
a
n
,
log
⁡
n
,
n
,
…
}). In parallel systems (multiprocessor systems, with p processors) we have the serial processing time T
*
(
n
)
=
T
1
(
n
) and parallel processing time T
p
(
n
). The performance of parallel algorithms can be analyzed using speed-up, efficiency, and isoefficiency metrics.(i)
Thespeed-up, S
(
p
), represents how a parallel algorithm is faster than a corresponding sequential algorithm. The speed-up is defined as S
(
p
)
=
T
1
(
n
)
/
T
p
(
n
). There are special bounds for speed-up [23]: S
(
p
)
≤
p
p
~
/
(
p
+
p
~
-
1
), where p
~
=
T
1
/
T
∞ is the average parallelism (the average number of busy processors given unbounded number of processors). Usually S
(
p
)
≤
p, but under special circumstances the speed-up can be S
(
p
)
>
p [24]. Another upper bound is established by the Amdahls law: S
(
p
)
=
(
s
+
(
(
1
-
s
)
/
p
)
)
1
/
2
≤
1
/
s where s is the fraction of a program that is sequential. The upper bound is considered for a 0 time of parallel fraction.
(ii)
Theefficiency is the average utilization of p processors: E
(
p
)
=
S
(
p
)
/
p.
(iii)
Theisoefficiency is the growth rate of workload W
p
(
n
)
=
p
T
p
(
n
) in terms of number of processors to keep efficiency fixed. If we consider W
1
(
n
)
-
E
W
p
(
n
)
=
0 for any fixed efficiency E we obtain p
=
p
(
n
). This means that we can establish a relation between needed number of processors and problem size. For example for the parallel sum of n numbers using p processors we have n
≈
E
(
n
+
p
log
⁡
p
), so n
=
Θ
(
p
log
⁡
p
).Numerical algorithms use for implementation a hypercube architecture. We analyze the performance of different numerical operations using the isoefficiency metric. For the hypercube architecture a simple model for intertask communication considersT
com
=
t
s
+
L
t
w where t
s is the latency (the time needed by a message to cross through the network), t
w is the time needed to send a word (1
/
t
w is called bandwidth), and L is the message length (expressed in number of words). The word size depends on processing architecture (usually it is two bytes). We define t
c as the processing time per word for a processor. We have the following results.(i)
External product  x
y
T. The isoefficiency is written as
(34)
t
c
n
≈
E
(
t
c
n
+
(
t
s
+
t
w
)
p
log
⁡
p
)
⟹
n
=
Θ
(
p
log
⁡
p
)
.
Parallel processing time is T
p
=
t
c
n
/
p
+
(
t
s
+
t
w
)
log
⁡
p. The optimality is computed using
(35)
d
T
p
d
p
=
0
⟹
-
t
c
n
p
2
+
t
s
+
t
w
p
=
0
⟹
p
≈
t
c
n
t
s
+
t
w
.
(ii)
Scalar product(internal product) x
T
y
=
∑
i
=
1
n
‍
x
i
y
i. The isoefficiency is written as
(36)
t
c
n
2
≈
E
(
t
c
n
2
+
t
s
2
p
log
⁡
p
+
t
w
2
n
p
log
⁡
p
)
⟹
n
=
Θ
(
p
(
log
⁡
p
)
2
)
.
(iii)
Matrix-vector product
y
=
A
x
,
y
i
=
∑
j
=
1
n
‍
A
i
j
x
j. The isoefficiency is written as
(37)
t
c
n
2
≈
E
(
t
c
n
2
+
t
s
p
log
⁡
p
+
t
w
n
p
log
⁡
p
)
⟹
n
=
Θ
(
p
(
log
⁡
p
)
2
)
.Table2 presented the amount of data that can be processed for a specific size. The cases that meet the upper bound n
≥
1.05
×
1
0
6 are marked with (*). To keep the efficiency high for a specific parallel architecture, HPC algorithms for particle physics introduce upper limits for the amount of data, which means that we have also an upper bound for Big Data volume in this case.Table 2
Isoefficiency for a hypercube architecture:n
=
Θ
(
p
log
⁡
p
) and n
=
Θ
(
p
(
log
⁡
p
)
2
). We marked with (*) the limitations imposed by Formula (33).
Scenario
Architecture size (p)
n
=
Θ
(
p
log
⁡
p
)
n
=
Θ
(
p
(
log
⁡
p
)
2
)
1
1
0
1
1.0
×
1
0
1
1.00
×
1
0
1
2
1
0
2
2.0
×
1
0
2
8.00
×
1
0
2
3
1
0
3
3.0
×
1
0
3
2.70
×
1
0
4
4
1
0
4
4.0
×
1
0
4
6.40
×
1
0
5
*
5
1
0
5
5.0
×
1
0
5
*
1.25
×
1
0
7
6
1
0
6
6.0
×
1
0
6
2.16
×
1
0
8
7
1
0
7
7.0
×
1
0
7
3.43
×
1
0
9
8
1
0
8
8.0
×
1
0
8
5.12
×
1
0
10
9
1
0
9
9.0
×
1
0
9
7.29
×
1
0
11The factors that determine the efficiency of parallel algorithms aretask balancing (work-load distribution between all used processors in a system → to be maximized); concurrency (the number/percentage of processors working simultaneously → to be maximized); and overhead (extra work for introduce by parallel processing that does not appear in serial processing → to be minimized).


## 2.1. General Approach of Event Generation

The most important aspect in simulation for HEP experiments is event generation. This process can be split into multiple steps, according to physical models. For example, structure of LHC (Large Hadron Collider) events:(
1
) hard process; (
2
) parton shower; (
3
) hadronization; (
4
) underlying event. According to official LHC website (http://home.web.cern.ch/about/computing): “approximately 600 million times per second, particles collide within the LHC
…
Experiments at CERN generate colossal amounts of data. The Data Centre stores it, and sends it around the world for analysis.” The analysis must produce valuable data and the simulation results must be correlated with physical experiments.Figure1 presents the general approach of event generation, detection, and reconstruction. The physical model is used to create simulation process that produces different type of events, clustered in vector of events (e.g., the fourth type of events in LHC experiments).Figure 1
General approach of event generation, detection, and reconstruction.In parallel, the real experiments are performed. The detectors identify the most relevant events and, based on reconstruction techniques, vector of events is created. The detectors can be real or simulated (software tools) and the reconstruction phase combines real events with events detected in simulation. At the end, the final result is compared with the simulation model (especially with generated vectors of events). The model can be corrected for further experiments. The goal is to obtain high accuracy and precision of measured and processed data.Software tools for event generation are based on random number generators. There are three types of random numbers: truly random numbers (from physical generators), pseudorandom numbers (from mathematical generators), and quasirandom numbers (special correlated sequences of numbers, used only for integration). For example, numerical integration using quasirandom numbers usually gives faster convergence than the standard integration methods based on quadratures. In event generation pseudorandom numbers are used most often.The most popular HEP application uses Poisson distribution combined with a basic normal distribution. The Poisson distribution can be formulated as(1)
P
[
X
=
k
]
=
μ
k
k
!
exp
⁡
{
-
μ
}
,
k
=
0,1
,
…
,
with E
(
k
)
=
V
(
k
)
=
μ (V is variance and E is expectation value). Having a uniform random number generator called RND() (Random based on Normal Distribution) we can use the following two algorithms for event generation techniques.The result of running Algorithms1 and 2 to generate around 1
0
6 random numbers is presented in Figure 2. In general, the second algorithm has better result for Poisson distribution. General recommendation for HEP experiments indicates the use of popular random number generators like TRNG (True Random Number Generators), RANMAR (Fast Uniform Random Number Generator used in CERN experiments), RANLUX (algorithm developed by Luscher used by Unix random number generators), and Mersenne Twister (the “industry standard”). Random number generators provided with compilers, operating system, and programming language libraries can have serious problem because they are based on system clock and suffer from lack of uniformity of distribution for large amounts of generated numbers and correlation of successive values.Algorithm 1:Random number generation for Poisson distribution using many random generated numbers with normal distribution (RND).
(1)procedureRa
ndom_Generator_Poisson(μ)
(2)n
u
m
b
e
r
←
-
1;
(3)a
c
c
u
m
u
l
a
t
o
r
←
1.0;
(4)q
←
exp
⁡
{
-
μ
};
(5)while
a
c
c
u
m
u
l
a
t
o
r
>
q
do
(6)r
n
d
_
n
u
m
b
e
r
←
R
N
D
(
);
(7)a
c
c
u
m
u
l
a
t
o
r
←
a
c
c
u
m
u
l
a
t
o
r
*
r
n
d
_
n
u
m
b
e
r;
(8)n
u
m
b
e
r
←
n
u
m
b
e
r
+
1;
(9)end while
(10)return number;
(11)end procedureAlgorithm 2:Random number generation for Poisson distribution using one random generated number with normal distribution.
(1)procedure Ra
ndom_Generator_Poisson_RND(
μ
,
r
)
(2)n
u
m
b
e
r
←
0;
(3)q
←
exp
⁡
{
-
μ
};
(4)a
c
c
u
m
u
l
a
t
o
r
←
q;
(5)p
←
q;
(6)while
r
>
a
c
c
u
m
u
l
a
t
o
r
do
(7)n
u
m
b
e
r
←
n
u
m
b
e
r
+
1;
(8)p
←
p
*
μ
/
n
u
m
b
e
r;
(9)a
c
c
u
m
u
l
a
t
o
r
←
a
c
c
u
m
u
l
a
t
o
r
+
p;
(10)end while
(11)return number;
(12)end procedureGeneral approach of event generation, detection, and reconstruction.
(a)
(b)The art of event generation is to use appropriate combinations of various random number generation methods in order to construct an efficient event generation algorithm being solution to a given problem in HEP.


## 2.2.  Monte Carlo Simulation and Markovian Monte Carlo Chains in HEP

In general, a Monte Carlo (MC) method is any simulation technique that uses random numbers to solve a well-defined problem,P. If F is a solution of the problem P (e.g., F
∈
R
n or F has a Boolean value), we define F
^, an estimation of F, as F
^
=
f
(
{
r
1
,
r
2
,
…
,
r
n
}
,
…
), where {
r
i
}
1
≤
i
≤
n is a random variable that can take more than one value and for which any value that will be taken cannot be predicted in advance. If ρ
(
r
) is the probability density function, ρ
(
r
)
d
r
=
P
[
r
<
r
′
<
r
+
d
r
], the cumulative distributed function is
(2)
C
(
r
)
=
∫
-
∞
r
‍
ρ
(
x
)
d
x
⟹
ρ
(
r
)
=
d
C
(
r
)
d
r
.C
(
r
) is a monotonically nondecreasing function with all values in [
0,1
]. The expectation value is
(3)
E
(
f
)
=
∫
‍
f
(
r
)
d
C
(
r
)
=
∫
‍
f
(
r
)
ρ
(
r
)
d
r
.
And the variance is
(4)
V
(
f
)
=
E
[
f
-
E
(
f
)
]
2
=
E
(
f
2
)
-
E
2
(
f
)
.


### 2.2.1. Monte Carlo Event Generation and Simulation

To define a MC estimator the “Law of Large Numbers (LLN)” is used. LLN can be described as follows: let one choosen numbers r
i randomly, with the probability density function uniform on a specific interval (
a
,
b
), each r
i being used to evaluate f
(
r
i
). For large n (consistent estimator),
(5)
1
n
∑
i
=
1
n
‍
f
(
r
i
)
⟶
E
(
f
)
=
1
b
-
a
∫
a
b
‍
f
(
r
)
d
r
.The properties of a MC estimator are being normally distributed (with Gaussian density); the standard deviation isσ
=
V
(
f
)
/
n; MC is unbiased for all n (the expectation value is the real value of the integral); the estimator is consistent if V
(
f
)
<
∞ (the estimator converges to the true value of the integral for every large n); a sampling phase can be applied to compute the estimator if we do not know anything about the function f; it is just suitable for integration. The sampling phase can be expressed, in a stratified way, as
(6)
∫
a
b
‍
f
(
r
)
d
r
=
∫
a
r
1
‍
f
(
r
)
d
r
+
∫
r
1
r
2
‍
f
(
r
)
d
r
+
⋯
+
∫
r
n
b
‍
f
(
r
)
d
r
.MC estimations and MC event generators are necessary tools in most of HEP experiments being used at all their steps: experiments preparation, simulation running, and data analysis.An example of MC estimation is the Lorentz invariant phase space (LIPS) that describes the cross section for a typical HEP process withn particle in the final state.Consider(7)
σ
n
~
∫
‍
|
M
|
2
d
R
n
,
where M is the matrix describing the interaction between particles and d
R
n is the element of LIPS. We have the following estimation:
(8)
R
n
(
P
,
p
1
,
p
2
,
…
,
p
n
)
=
∫
‍
δ
(
4
)
(
P
-
∑
k
=
1
n
‍
p
k
)
∏
k
=
1
n
‍
(
δ
(
p
k
2
-
m
k
2
)
Θ
(
p
k
0
)
d
4
p
k
)
,
where P is total four-momentum of the n-particle system; p
kand m
k are four-momenta and mass of the final state particles; δ
(
4
)
(
P
-
∑
k
=
1
n
‍
p
k
) is the total energy momentum conservation; δ
(
p
k
2
-
m
k
2
) is the on-mass-shell condition for the final state system. Based on the integration formula
(9)
∫
‍
δ
(
p
k
2
-
m
k
2
)
Θ
(
p
k
0
)
d
4
p
k
=
d
3
p
k
2
p
k
0
,
obtain the iterative form for cross section:
(10)
R
n
(
P
,
p
1
,
p
2
,
…
,
p
n
)
=
∫
‍
R
n
-
1
(
P
-
p
n
,
p
1
,
p
2
,
…
,
p
n
-
1
)
d
3
p
n
2
p
n
0
,
which can be numerical integrated by using the recurrence relation. As result, we can construct a general MC algorithm for particle collision processes.Example 1.
Let us consider the interaction:e
+
e
-
→
μ
+
μ
- where Higgs boson contribution is numerically negligible. Figure 3 describes this interaction (Φ is the azimuthal angle, θ the polar angle, and p
1
,
p
2
,
q
1
,
q
2 are the four-momenta for particles).
The cross section is(11)
d
σ
=
α
2
4
s
[
W
1
(
s
)
(
1
+
cos
⁡
2
θ
)
+
W
2
(
s
)
cos
⁡
θ
]
d
Ω
,
where d
Ω
=
d
cos
⁡
θ
d
Φ, α
=
e
2
/
4
π (fine structure constant), s
=
(
p
1
0
+
p
2
0
)
2 is the center of mass energy squared, and W
1
(
s
) and W
2
(
s
) are constant functions. For pure processes we have W
1
(
s
)
=
1 and W
2
(
s
)
=
0, and the total cross section becomes
(12)
σ
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
d
2
σ
d
Φ
d
cos
⁡
θ
.Figure 3
Example of particle interaction:e
(
p
1
)
+
e
-
(
p
2
)
→
μ
+
(
q
1
)
μ
-
(
q
2
).We introduce the following notation:(13)
ρ
(
cos
⁡
θ
,
Φ
)
=
d
2
σ
d
Φ
d
cos
⁡
θ
,
and let us consider ρ
~
(
cos
⁡
θ
,
Φ
) an approximation of ρ
(
cos
⁡
θ
,
Φ
). Then σ
~
=
∬
d
Φ
d
cos
⁡
θ
ρ
~. Now, we can compute
(14)
σ
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
ρ
(
cos
⁡
θ
,
Φ
)
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
w
(
cos
⁡
θ
,
Φ
)
ρ
~
(
cos
⁡
θ
,
Φ
)
≈
〈
w
〉
ρ
~
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
ρ
~
(
cos
⁡
θ
,
Φ
)
=
σ
~
〈
w
〉
ρ
~
,
where w
(
cos
⁡
θ
,
Φ
)
=
ρ
(
cos
⁡
θ
,
Φ
)
/
ρ
~
(
cos
⁡
θ
,
Φ
) and 〈
w
〉
ρ
~ is the estimation of w based on ρ
~. Here, the MC estimator is
(15)
〈
w
〉
MC
=
1
n
∑
i
=
1
n
‍
w
i
,
and the standard deviation is
(16)
s
MC
=
(
1
n
(
n
-
1
)
∑
i
=
1
n
‍
(
w
i
-
〈
w
〉
MC
)
2
)
1
/
2
.The final numerical result based on MC estimator is(17)
σ
MC
=
σ
~
〈
w
〉
MC
±
σ
~
s
MC
.As we can show, the principle of a Monte Carlo estimator in physics is to simulate the cross section in interaction and radiation transport knowing the probability distributions (or an approximation) governing each interaction of elementary particles.Based on this result, the Monte Carlo algorithm used to generate events is as follows. It takes as inputρ
~
(
cos
⁡
θ
,
Φ
) and in a main loop considers the following steps: (
1
) generate (
cos
⁡
θ
,
Φ
) peer from ρ
~; (
2
) compute four-momenta p
1
,
p
2
,
q
1
,
q
2; (
3
) compute w
=
ρ
/
ρ
~. The loop can be stopped in the case of unweighted events, and we will stay in the loop for weighted events. As output, the algorithm returns four-momenta for particle for weighted events and four-momenta and an array of weights for unweighted events. The main issue is how to initialize the input of the algorithm. Based on d
σ formula (for W
1
(
s
)
=
1 and W
2
(
s
)
=
0), we can consider as input ρ
~
(
cos
⁡
θ
,
Φ
)
=
(
α
2
/
4
s
)
(
1
+
cos
⁡
2
θ
). Then σ
~
=
4
π
α
2
/
3
s.In HEP theoretical predictions used for particle collision processes modeling (as shown in presented example) should be provided in terms of Monte Carlo event generators, which directly simulate these processes and can provide unweighted (weight = 1) events. A good Monte Carlo algorithm should be used not only for numerical integration [7] (i.e., provide weighted events) but also for efficient generation of unweighted events, which is very important issue for HEP.


### 2.2.2.  Markovian Monte-Carlo Chains

A classical Monte Carlo method estimates a functionF with F
^ by using a random variable. The main problem with this approach is that we cannot predict any value in advance for a random variable. In HEP simulation experiments the systems are described in states [8]. Let us consider a system with a finite set of possible states S
1
,
S
2
,
…, and S
t the state at the moment t. The conditional probability is defined as
(18)
P
(
S
t
=
S
j
∣
S
t
1
=
S
i
1
,
S
t
2
=
S
i
2
,
…
,
S
t
n
=
S
i
n
)
,
where the mappings (
t
1
,
i
1
)
,
…
,
(
t
n
,
i
n
) can be interpreted as the description of system evolution in time by specifying a specific state for each moment of time.The system is a Markov chain if the distribution ofS
t depends only on immediate predecessor S
t
-
1 and it is independent of all previous states as follows:
(19)
P
(
S
t
=
S
j
∣
S
t
-
1
=
S
i
t
-
1
,
…
,
S
t
2
=
S
i
2
,
S
t
1
=
S
i
1
)
=
P
(
S
t
=
S
j
∣
S
t
-
1
=
S
i
t
-
1
)
.To generate the time steps(
t
1
,
t
2
,
…
,
t
n
) we use the probability of a single forward Markovian step given by p
(
t
∣
t
n
) with the property ∫
t
n
∞
‍
p
(
t
∣
t
n
)
d
t
=
1 and we define p
(
t
)
=
p
(
t
∣
0
). The 1-dimensional Monte Carlo Markovian Algorithm used to generate the time steps is presented in Algorithm 3.Algorithm 3:1-Dimensional Monte Carlo Markovian Algorithm.
(1)   Generatet
1 according with p
(
t
1
)
=
p
(
t
1
∣
t
0
=
0
)
(2)if
t
1
<
t
max
⁡
then             ▹  Generate the initial state.
(3)P
N
≥
1
=
∫
0
t
max
⁡
‍
p
(
t
1
∣
t
0
)
d
t
1;        ▹  Compute the initial probability.
(4)Retain
t
1;
(5)end if
(6)if
t
1
>
t
max
⁡
then
▹  Discard all generated and computed data.
(7)N
=
0; P
0
=
∫
t
max
⁡
∞
‍
p
(
t
1
∣
t
0
)
d
t
1
=
e
-
t
max
⁡;
(8)Delete
t
1;
(9)EXIT.               ▹  The algorithm ends here.
(10)end if
(11)i
=
2;
(12)while (1) do            ▹  Infinite loop until a successful EXIT.
(13)   Generatet
i according with p
(
t
i
∣
t
i
-
1
)
(14)if
t
i
<
t
max
⁡
then
▹  Generate a new state and new probability.
(15)P
N
≥
i
=
∫
t
i
t
max
⁡
‍
p
(
t
i
∣
t
i
-
1
)
d
t
i;
(16)Retain
t
i;
(17)end if
(18)if
t
i
>
t
max
⁡
then
▹  Discard all generated and computed data.
(19)N
=
i
-
1; P
i
=
∫
t
max
⁡
∞
‍
p
(
t
i
∣
t
i
-
1
)
d
t
i;
(20)Retain
(
t
1
,
t
2
,
…
,
t
i
-
1
);Delete
t
i;
(21)EXIT.                  ▹  The algorithm ends here.
(22)end if
(23)i
=
i
+
1;
(24)end whileThe main result of Algorithm3 is that P
(
t
max
⁡
) follows a Poisson distribution:
(20)
P
N
=
∫
0
t
max
⁡
‍
p
(
t
1
∣
t
0
)
d
t
1
×
∫
t
1
t
max
⁡
‍
p
(
t
2
∣
t
1
)
d
t
2
×
⋯
×
∫
t
N
-
1
t
max
⁡
‍
p
(
t
N
∣
t
N
-
1
)
d
t
N
×
∫
t
max
⁡
∞
‍
p
(
t
N
+
1
∣
t
N
)
d
t
N
+
1
=
1
N
!
(
t
max
⁡
)
N
e
-
t
max
⁡
.We can consider the 1-dimensional Monte Carlo Markovian Algorithm as a method used to iteratively generate the systems’ states (codified as a Markov chain) in simulation experiments. According to the Ergodic Theorem for Markov chains, the chain defined has a unique stationary probability distribution [9, 10].Figures4 and 5 present the running of Algorithm 3. According to different values of parameter s used to generate the next step, the results are very different, for 1000 iterations. Figure 4 for s
=
1 shows a profile of the type of noise. For s
=
10,100,1000 profile looks like some of the information is filtered and lost. The best results are obtained for s
=
0.01 and s
=
0.1 and the generated values can be easily accepted for MC simulation in HEP experiments.Example of 1-dimensional Monte Carlo Markovian algorithm.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
Analysis of acceptance rate for 1-dimensional Monte Carlo Markovian algorithm for differents values.Figure5 shows the acceptance rate of values generated with parameter s used in the algorithm. And parameter values are correlated with Figure 4. Results in Figure 5 show that the acceptance rate decreases rapidly with increasing value of parameter s. The conclusion is that values must be kept small to obtain meaningful data. A correlation with the normal distribution is evident, showing that a small value for the mean square deviation provides useful results.


### 2.2.3. Performance of Numerical Algorithms Used in MC Simulations

Numerical methods used to compute MC estimator use numerical quadratures to approximate the value of the integral for functionf on a specific domain by a linear compilation of function values and weights {
w
i
}
1
≤
i
≤
m as follows:
(21)
∫
a
b
‍
f
(
r
)
d
r
=
∑
i
=
1
m
‍
w
i
f
(
r
i
)
.We can consider a consistent MC estimatora a classical numerical quadrature with all w
i
=
1. Efficiency of integration methods for 1 dimension and for d dimensions is presented in Table 1. We can conclude that quadrature methods are difficult to apply in many dimensions for variate integration domains (regions) and the integral is not easy to be estimated.Table 1
Efficiency of integration methods for 1 dimension and ford dimensions.
Method
1 dimension
d dimensions
Monte Carlo
n
-
1
/
2
n
-
1
/
2
Trapezoidal rule
n
-
2
n
-
2
/
d
Simpson’s rule
n
-
4
n
-
4
/
d
m-points Gauss rule (m
<
n)
n
-
2
m
n
-
2
m
/
dAs practical example, in a typical high-energy particle collision there can be many final-state particles (even hundreds). If we haven final state particle, we face with d
=
3
n
-
4 dimensional phase space. As numerical example, for n
=
4 we have d
=
8 dimensions, which is very difficult approach for classical numerical quadratures.Full decomposition integration volume for one double number (10 Bytes) per volume unit isn
d
×
10 Bytes. For the example considered with d
=
8 and n
=
10 divisions for interval [
0,1
] we have, for one numerical integration,
(22)
n
d
×
10
Bytes
=
1
0
8
×
10
102
4
3
G
Bytes
≈
0.93
G
Bytes
.
Considering 1
0
6 events per second, one integration per event, the data produced in one hour will be ≈3197.4 P Bytes.The previous assumption is only for multidimensional arrays. But due to the factorization assumption,p
(
r
1
,
r
2
,
…
,
r
n
)
=
p
(
r
1
)
p
(
r
2
)
⋯
p
(
r
n
), we obtain for one integration
(23)
n
×
d
×
10
Bytes
=
800
Bytes
,
which means ≈2.62 T Bytes of data produce for one hour of simulations.


## 2.2.1. Monte Carlo Event Generation and Simulation

To define a MC estimator the “Law of Large Numbers (LLN)” is used. LLN can be described as follows: let one choosen numbers r
i randomly, with the probability density function uniform on a specific interval (
a
,
b
), each r
i being used to evaluate f
(
r
i
). For large n (consistent estimator),
(5)
1
n
∑
i
=
1
n
‍
f
(
r
i
)
⟶
E
(
f
)
=
1
b
-
a
∫
a
b
‍
f
(
r
)
d
r
.The properties of a MC estimator are being normally distributed (with Gaussian density); the standard deviation isσ
=
V
(
f
)
/
n; MC is unbiased for all n (the expectation value is the real value of the integral); the estimator is consistent if V
(
f
)
<
∞ (the estimator converges to the true value of the integral for every large n); a sampling phase can be applied to compute the estimator if we do not know anything about the function f; it is just suitable for integration. The sampling phase can be expressed, in a stratified way, as
(6)
∫
a
b
‍
f
(
r
)
d
r
=
∫
a
r
1
‍
f
(
r
)
d
r
+
∫
r
1
r
2
‍
f
(
r
)
d
r
+
⋯
+
∫
r
n
b
‍
f
(
r
)
d
r
.MC estimations and MC event generators are necessary tools in most of HEP experiments being used at all their steps: experiments preparation, simulation running, and data analysis.An example of MC estimation is the Lorentz invariant phase space (LIPS) that describes the cross section for a typical HEP process withn particle in the final state.Consider(7)
σ
n
~
∫
‍
|
M
|
2
d
R
n
,
where M is the matrix describing the interaction between particles and d
R
n is the element of LIPS. We have the following estimation:
(8)
R
n
(
P
,
p
1
,
p
2
,
…
,
p
n
)
=
∫
‍
δ
(
4
)
(
P
-
∑
k
=
1
n
‍
p
k
)
∏
k
=
1
n
‍
(
δ
(
p
k
2
-
m
k
2
)
Θ
(
p
k
0
)
d
4
p
k
)
,
where P is total four-momentum of the n-particle system; p
kand m
k are four-momenta and mass of the final state particles; δ
(
4
)
(
P
-
∑
k
=
1
n
‍
p
k
) is the total energy momentum conservation; δ
(
p
k
2
-
m
k
2
) is the on-mass-shell condition for the final state system. Based on the integration formula
(9)
∫
‍
δ
(
p
k
2
-
m
k
2
)
Θ
(
p
k
0
)
d
4
p
k
=
d
3
p
k
2
p
k
0
,
obtain the iterative form for cross section:
(10)
R
n
(
P
,
p
1
,
p
2
,
…
,
p
n
)
=
∫
‍
R
n
-
1
(
P
-
p
n
,
p
1
,
p
2
,
…
,
p
n
-
1
)
d
3
p
n
2
p
n
0
,
which can be numerical integrated by using the recurrence relation. As result, we can construct a general MC algorithm for particle collision processes.Example 1.
Let us consider the interaction:e
+
e
-
→
μ
+
μ
- where Higgs boson contribution is numerically negligible. Figure 3 describes this interaction (Φ is the azimuthal angle, θ the polar angle, and p
1
,
p
2
,
q
1
,
q
2 are the four-momenta for particles).
The cross section is(11)
d
σ
=
α
2
4
s
[
W
1
(
s
)
(
1
+
cos
⁡
2
θ
)
+
W
2
(
s
)
cos
⁡
θ
]
d
Ω
,
where d
Ω
=
d
cos
⁡
θ
d
Φ, α
=
e
2
/
4
π (fine structure constant), s
=
(
p
1
0
+
p
2
0
)
2 is the center of mass energy squared, and W
1
(
s
) and W
2
(
s
) are constant functions. For pure processes we have W
1
(
s
)
=
1 and W
2
(
s
)
=
0, and the total cross section becomes
(12)
σ
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
d
2
σ
d
Φ
d
cos
⁡
θ
.Figure 3
Example of particle interaction:e
(
p
1
)
+
e
-
(
p
2
)
→
μ
+
(
q
1
)
μ
-
(
q
2
).We introduce the following notation:(13)
ρ
(
cos
⁡
θ
,
Φ
)
=
d
2
σ
d
Φ
d
cos
⁡
θ
,
and let us consider ρ
~
(
cos
⁡
θ
,
Φ
) an approximation of ρ
(
cos
⁡
θ
,
Φ
). Then σ
~
=
∬
d
Φ
d
cos
⁡
θ
ρ
~. Now, we can compute
(14)
σ
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
ρ
(
cos
⁡
θ
,
Φ
)
=
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
w
(
cos
⁡
θ
,
Φ
)
ρ
~
(
cos
⁡
θ
,
Φ
)
≈
〈
w
〉
ρ
~
∫
0
2
π
‍
d
Φ
∫
-
1
1
‍
d
cos
⁡
θ
ρ
~
(
cos
⁡
θ
,
Φ
)
=
σ
~
〈
w
〉
ρ
~
,
where w
(
cos
⁡
θ
,
Φ
)
=
ρ
(
cos
⁡
θ
,
Φ
)
/
ρ
~
(
cos
⁡
θ
,
Φ
) and 〈
w
〉
ρ
~ is the estimation of w based on ρ
~. Here, the MC estimator is
(15)
〈
w
〉
MC
=
1
n
∑
i
=
1
n
‍
w
i
,
and the standard deviation is
(16)
s
MC
=
(
1
n
(
n
-
1
)
∑
i
=
1
n
‍
(
w
i
-
〈
w
〉
MC
)
2
)
1
/
2
.The final numerical result based on MC estimator is(17)
σ
MC
=
σ
~
〈
w
〉
MC
±
σ
~
s
MC
.As we can show, the principle of a Monte Carlo estimator in physics is to simulate the cross section in interaction and radiation transport knowing the probability distributions (or an approximation) governing each interaction of elementary particles.Based on this result, the Monte Carlo algorithm used to generate events is as follows. It takes as inputρ
~
(
cos
⁡
θ
,
Φ
) and in a main loop considers the following steps: (
1
) generate (
cos
⁡
θ
,
Φ
) peer from ρ
~; (
2
) compute four-momenta p
1
,
p
2
,
q
1
,
q
2; (
3
) compute w
=
ρ
/
ρ
~. The loop can be stopped in the case of unweighted events, and we will stay in the loop for weighted events. As output, the algorithm returns four-momenta for particle for weighted events and four-momenta and an array of weights for unweighted events. The main issue is how to initialize the input of the algorithm. Based on d
σ formula (for W
1
(
s
)
=
1 and W
2
(
s
)
=
0), we can consider as input ρ
~
(
cos
⁡
θ
,
Φ
)
=
(
α
2
/
4
s
)
(
1
+
cos
⁡
2
θ
). Then σ
~
=
4
π
α
2
/
3
s.In HEP theoretical predictions used for particle collision processes modeling (as shown in presented example) should be provided in terms of Monte Carlo event generators, which directly simulate these processes and can provide unweighted (weight = 1) events. A good Monte Carlo algorithm should be used not only for numerical integration [7] (i.e., provide weighted events) but also for efficient generation of unweighted events, which is very important issue for HEP.


## 2.2.2.  Markovian Monte-Carlo Chains

A classical Monte Carlo method estimates a functionF with F
^ by using a random variable. The main problem with this approach is that we cannot predict any value in advance for a random variable. In HEP simulation experiments the systems are described in states [8]. Let us consider a system with a finite set of possible states S
1
,
S
2
,
…, and S
t the state at the moment t. The conditional probability is defined as
(18)
P
(
S
t
=
S
j
∣
S
t
1
=
S
i
1
,
S
t
2
=
S
i
2
,
…
,
S
t
n
=
S
i
n
)
,
where the mappings (
t
1
,
i
1
)
,
…
,
(
t
n
,
i
n
) can be interpreted as the description of system evolution in time by specifying a specific state for each moment of time.The system is a Markov chain if the distribution ofS
t depends only on immediate predecessor S
t
-
1 and it is independent of all previous states as follows:
(19)
P
(
S
t
=
S
j
∣
S
t
-
1
=
S
i
t
-
1
,
…
,
S
t
2
=
S
i
2
,
S
t
1
=
S
i
1
)
=
P
(
S
t
=
S
j
∣
S
t
-
1
=
S
i
t
-
1
)
.To generate the time steps(
t
1
,
t
2
,
…
,
t
n
) we use the probability of a single forward Markovian step given by p
(
t
∣
t
n
) with the property ∫
t
n
∞
‍
p
(
t
∣
t
n
)
d
t
=
1 and we define p
(
t
)
=
p
(
t
∣
0
). The 1-dimensional Monte Carlo Markovian Algorithm used to generate the time steps is presented in Algorithm 3.Algorithm 3:1-Dimensional Monte Carlo Markovian Algorithm.
(1)   Generatet
1 according with p
(
t
1
)
=
p
(
t
1
∣
t
0
=
0
)
(2)if
t
1
<
t
max
⁡
then             ▹  Generate the initial state.
(3)P
N
≥
1
=
∫
0
t
max
⁡
‍
p
(
t
1
∣
t
0
)
d
t
1;        ▹  Compute the initial probability.
(4)Retain
t
1;
(5)end if
(6)if
t
1
>
t
max
⁡
then
▹  Discard all generated and computed data.
(7)N
=
0; P
0
=
∫
t
max
⁡
∞
‍
p
(
t
1
∣
t
0
)
d
t
1
=
e
-
t
max
⁡;
(8)Delete
t
1;
(9)EXIT.               ▹  The algorithm ends here.
(10)end if
(11)i
=
2;
(12)while (1) do            ▹  Infinite loop until a successful EXIT.
(13)   Generatet
i according with p
(
t
i
∣
t
i
-
1
)
(14)if
t
i
<
t
max
⁡
then
▹  Generate a new state and new probability.
(15)P
N
≥
i
=
∫
t
i
t
max
⁡
‍
p
(
t
i
∣
t
i
-
1
)
d
t
i;
(16)Retain
t
i;
(17)end if
(18)if
t
i
>
t
max
⁡
then
▹  Discard all generated and computed data.
(19)N
=
i
-
1; P
i
=
∫
t
max
⁡
∞
‍
p
(
t
i
∣
t
i
-
1
)
d
t
i;
(20)Retain
(
t
1
,
t
2
,
…
,
t
i
-
1
);Delete
t
i;
(21)EXIT.                  ▹  The algorithm ends here.
(22)end if
(23)i
=
i
+
1;
(24)end whileThe main result of Algorithm3 is that P
(
t
max
⁡
) follows a Poisson distribution:
(20)
P
N
=
∫
0
t
max
⁡
‍
p
(
t
1
∣
t
0
)
d
t
1
×
∫
t
1
t
max
⁡
‍
p
(
t
2
∣
t
1
)
d
t
2
×
⋯
×
∫
t
N
-
1
t
max
⁡
‍
p
(
t
N
∣
t
N
-
1
)
d
t
N
×
∫
t
max
⁡
∞
‍
p
(
t
N
+
1
∣
t
N
)
d
t
N
+
1
=
1
N
!
(
t
max
⁡
)
N
e
-
t
max
⁡
.We can consider the 1-dimensional Monte Carlo Markovian Algorithm as a method used to iteratively generate the systems’ states (codified as a Markov chain) in simulation experiments. According to the Ergodic Theorem for Markov chains, the chain defined has a unique stationary probability distribution [9, 10].Figures4 and 5 present the running of Algorithm 3. According to different values of parameter s used to generate the next step, the results are very different, for 1000 iterations. Figure 4 for s
=
1 shows a profile of the type of noise. For s
=
10,100,1000 profile looks like some of the information is filtered and lost. The best results are obtained for s
=
0.01 and s
=
0.1 and the generated values can be easily accepted for MC simulation in HEP experiments.Example of 1-dimensional Monte Carlo Markovian algorithm.
(a)
(b)
(c)
(d)
(e)
(f)Figure 5
Analysis of acceptance rate for 1-dimensional Monte Carlo Markovian algorithm for differents values.Figure5 shows the acceptance rate of values generated with parameter s used in the algorithm. And parameter values are correlated with Figure 4. Results in Figure 5 show that the acceptance rate decreases rapidly with increasing value of parameter s. The conclusion is that values must be kept small to obtain meaningful data. A correlation with the normal distribution is evident, showing that a small value for the mean square deviation provides useful results.


## 2.2.3. Performance of Numerical Algorithms Used in MC Simulations

Numerical methods used to compute MC estimator use numerical quadratures to approximate the value of the integral for functionf on a specific domain by a linear compilation of function values and weights {
w
i
}
1
≤
i
≤
m as follows:
(21)
∫
a
b
‍
f
(
r
)
d
r
=
∑
i
=
1
m
‍
w
i
f
(
r
i
)
.We can consider a consistent MC estimatora a classical numerical quadrature with all w
i
=
1. Efficiency of integration methods for 1 dimension and for d dimensions is presented in Table 1. We can conclude that quadrature methods are difficult to apply in many dimensions for variate integration domains (regions) and the integral is not easy to be estimated.Table 1
Efficiency of integration methods for 1 dimension and ford dimensions.
Method
1 dimension
d dimensions
Monte Carlo
n
-
1
/
2
n
-
1
/
2
Trapezoidal rule
n
-
2
n
-
2
/
d
Simpson’s rule
n
-
4
n
-
4
/
d
m-points Gauss rule (m
<
n)
n
-
2
m
n
-
2
m
/
dAs practical example, in a typical high-energy particle collision there can be many final-state particles (even hundreds). If we haven final state particle, we face with d
=
3
n
-
4 dimensional phase space. As numerical example, for n
=
4 we have d
=
8 dimensions, which is very difficult approach for classical numerical quadratures.Full decomposition integration volume for one double number (10 Bytes) per volume unit isn
d
×
10 Bytes. For the example considered with d
=
8 and n
=
10 divisions for interval [
0,1
] we have, for one numerical integration,
(22)
n
d
×
10
Bytes
=
1
0
8
×
10
102
4
3
G
Bytes
≈
0.93
G
Bytes
.
Considering 1
0
6 events per second, one integration per event, the data produced in one hour will be ≈3197.4 P Bytes.The previous assumption is only for multidimensional arrays. But due to the factorization assumption,p
(
r
1
,
r
2
,
…
,
r
n
)
=
p
(
r
1
)
p
(
r
2
)
⋯
p
(
r
n
), we obtain for one integration
(23)
n
×
d
×
10
Bytes
=
800
Bytes
,
which means ≈2.62 T Bytes of data produce for one hour of simulations.


## 2.3. Unfolding Processes in Particle Physics and Kernel Estimation in HEP

In particle physics analysis we have two types of distributions: true distribution (considered in theoretical models) and measured distribution (considered in experimental models, which are affected by finite resolution and limited acceptance of existing detectors). A HEP interaction process starts with a true knows distribution and generate a measured distribution, corresponding to an experiment of a well-confirmed theory. An inverse process starts with a measured distribution and tries to identify the true distribution. These unfolding processes are used to identify new theories based on experiments [11].


### 2.3.1. Unfolding Processes in Particle Physics

The theory of unfolding processes in particle physics is as follows [12]. For a physics variable t we have a true distribution f
(
t
) mapped in x and an n-vector of unknowns and a measured distribution g
(
s
) (for a measured variable s) mapped in an m-vector of measured data. A response matrix A
∈
R
m
×
n encodes a Kernel function K
(
s
,
t
) describing the physical measurement process [12–15]. The direct and inverse processes are described by the Fredholm integral equation [16] of the first kind, for a specific domain Ω,
(24)
∫
Ω
‍
K
(
s
,
t
)
f
(
t
)
d
t
=
g
(
s
)
.
In particle physics the Kernel function K
(
s
,
t
) is usually known from a Monte Carlo sample obtained from simulation. A numerical solution is obtained using the following linear equation: A
x
=
b. Vectors x and y are assumed to be 1-dimensional in theory, but they can be multidimensional in practice (considering multiple independent linear equations). In practice, also the statistical properties of the measurements are well known and often they follow the Poisson statistics [17]. To solve the linear systems we have different numerical methods.First method is based on linear transformationx
=
A
#
y. If m
=
n then A
#
=
A
-
1 and we can use direct Gaussian methods, iterative methods (Gauss-Siedel, Jacobi or SOR), or orthogonal methods (based on Householder transformation, Givens methods, or Gram-Schmidt algorithm). If m
>
n (the most frequent scenario) we will construct the matrix A
#
=
(
A
T
A
)
-
1
A
T (called pseudoinverse Penrose-Moore). In these cases the orthogonal methods offer very good and stable numerical solutions.Second method considers the singular value decomposition:(25)
A
=
U
Σ
V
T
=
∑
i
=
1
n
‍
σ
i
u
i
v
i
T
,
where U
∈
R
m
×
n and V
∈
R
n
×
n are matrices with orthonormal columns and the diagonal matrix Σ
=
diag
⁡
{
σ
1
,
…
,
σ
n
}
=
U
T
A
V. The solution is
(26)
x
=
A
#
y
=
V
Σ
-
1
(
U
T
y
)
=
∑
i
=
1
n
‍
1
σ
i
(
u
i
T
y
)
v
i
=
∑
i
=
1
n
‍
1
σ
i
c
i
v
i
,
where c
i
=
u
i
T
y, i
=
1
,
…
,
n, are called Fourier coefficients.


### 2.3.2.  Random Matrix Theory

Analysis of particle spectrum (e.g., neutrino spectrum) faces with Random Matrix Theory (RMT), especially if we consider anarchic neutrino masses. The RMT means the study of the statistical properties of eigenvalues of very large matrices [18]. For an interaction matrix A (with size N), where A
i
j is an independent distributed random variable and A
H is the complex conjugate and transpose matrix, we define M
=
A
+
A
H, which describes a Gaussian Unitary Ensemble (GUE). The GUE properties are described by the probability distribution P
(
M
)
d
M: (
1
) it is invariant under unitary transformation, P
(
M
)
d
m
=
P
(
M
′
)
d
M
′, where M
′
=
U
H
M
U, U is a Hermitian matrix (U
H
U
=
I); (
2
) the elements of M matrix are statistically independent, P
(
M
)
=
∏
i
≤
j
‍
P
i
j
(
M
i
j
); and (
3
) the matrix M can be diagonalized as M
=
U
D
U
H, where U
=
diag
⁡
{
λ
1
,
…
,
λ
N
}, λ
i is the eigenvalue of M and λ
i
≤
λ
j if i
<
j
(27)
Propability
(
2
)

:
P
(
M
)
d
M
~
d
M
exp
⁡
{
-
N
2
T
r
(
M
H
M
)
}
;
Propability
(
3
)

:
P
(
M
)
d
M
~
d
U
∏
i
d
λ
i
∏
i
<
j
(
λ
i
-
λ
j
)
2
×
exp
⁡
{
-
N
2
∑
i
‍
(
λ
i
2
)
}
.The numerical methods used for eigenvalues computation are the QR method and Power methods (direct and indirect). The QR method is a numerical stable algorithm and Power method is an iterative one. The RMT can be used for many body systems, quantum chaos, disordered systems, quantum chromodynamics, and so forth.


### 2.3.3.  Kernel Estimation in HEP

Kernel estimation is a very powerful solution and relevant method for HEP when it is necessary to combine data from heterogeneous sources like MC datasets obtained by simulation and from Standard Model expectation, obtained from real experiments [19]. For a set of data {
x
i
}
1
≤
i
≤
n with a constant bandwidth h (the difference between two consecutive data values), called the smoothing parameter, we have the estimation
(28)
f
^
(
x
)
=
1
n
h
∑
i
=
1
n
‍
K
(
x
-
x
i
h
)
,
where K is an estimator. For example, a Gauss estimator with mean μ and standard deviation σ is
(29)
K
(
x
)
=
1
σ
2
π
exp
⁡
{
-
(
x
-
μ
)
2
2
σ
2
}
,
and has the following properties: positive definite and infinitely differentiable (due to the exp function), and it can be defined for an infinite supports (n
→
∞). The kernel is a nonparametric method, which means that h is independent of dataset and for large amount of normally distributed data we can find a value for h that minimizes the integrated squared error of f
^
(
x
). This value for bandwidth is computed as
(30)
h
*
=
(
4
3
n
)
1
/
5
σ
.The main problem in Kernel Estimation is that the set of data{
x
i
}
1
≤
i
≤
n is not normally distributed and in real experiments the optimal bandwidth it is not known. An improvement of presented method considers adaptive Kernel Estimation proposed by Abramson [20], where h
i
=
h
/
f
(
x
i
) and σ are considered global qualities for dataset. The new form is
(31)
f
^
a
(
x
)
=
1
n
∑
i
=
1
n
‍
1
h
i
K
(
x
-
x
i
h
i
)
,
and the local bandwidth value that minimizes the integrated squared error of f
^
a
(
x
) is
(32)
h
i
*
=
(
4
3
n
)
1
/
5
σ
f
^
(
x
i
)
,
where f
^ is the normal estimator.Kernel estimation is used for event selection to confidence level evaluation, for example, in Markovian Monte Carlo chains or in selection of neural network output used in experiments for reconstructed Higgs mass. In general, the main usage of Kernel estimation in HEP is searching for new particle, by finding relevant data in a large dataset.A method based on Kernel estimation is the graphical representation of datasets using advanced shifted histogram algorithm (ASH). This is a numerical interpolation for large datasets with the main aim of creating a set ofn
b
i
n histograms H
=
{
H
i
}, with the same bin-width h. Algorithm 4 presents the steps of histograms generation starting with a specific interval [
a
,
b
], a number of points n in this interval, and a number of bins and a number of values used for kernel estimation, m. Figure 6 shows the results of kernel estimation if function f
=
-
(
1
/
2
)
x
2 on [
0,1
] and graphical representation with a different number of bins. The values on vertical axis are aggregated in step 17 of Algorithm 4 and increase with the number of bins.Algorithm 4:Advanced shifted histogram (1D algorithm).
(1)procedure
ASH
(
a
,
b
,
n
,
x
,
n
bin
,
m
,
f
m
)
(2)δ
=
(
b
-
a
)
/
n
bin; h
=
m
δ;
(3)for
k
=
1
…
n
bin
do
(4)v
k
=
0; y
k
=
0;
(5)end for
(6)for
i
=
1
…
n
do
(7)k
=
(
x
i
-
a
)
/
δ
+
1;
(8)if
k
∈
[
1
,
n
bin
]
then
(9)v
k
=
v
k
+
1;
(10)end if
(11)end for
(12)for
k
=
1
…
n
bin
do
(13)if
v
k
=
0
then
(14)k
=
k
+
1;
(15)end if
(16)for
i
=
max
⁡
{
1
,
k
-
m
+
1
}
…
min
⁡
{
n
bin
,
k
+
m
-
1
}
do
(17)y
i
=
y
i
+
v
k
f
m
(
i
-
k
);
(18)end for
(19)end for
(20)for
k
=
1
…
n
bin
do
(21)y
k
=
y
k
/
(
n
h
);
(22)t
k
=
a
+
(
k
-
0.5
)
δ;
(23)end for
(24)return
{
t
k
}
1
≤
k
≤
n
bin, {
y
k
}
1
≤
k
≤
n
bin.
(25)end procedureExample of advanced shifted histogram algorithm running for different bins: 10, 100, and 1000.
(a)
(b)
(c)


### 2.3.4.  Performance of Numerical Algorithms Used in Particle Physics

All operations used in presented methods for particle physics (Unfolding Processes, Random Matrix Theory, and Kernel Estimation) can be reduced to scalar products, matrix-vector products, and matrix-matrix products. In [21] the design of new standard for the BLAS (Basic Linear Algebra Subroutines) in C language by extension of precision is described. This permits higher internal precision and mixed input/output types. The precision allows implementation of some algorithms that are simpler, more accurate, and sometimes faster than possible without these features. Regarding the precision of numerical computing, Dongarra and Langou established in [22] an upper bound for the residual check for A
x
=
y system, with A
∈
R
n
×
n a dense matrix. The residual check is defined as
(33)
r
∞
=
∥
A
x
-
y
∥
∞
n
ϵ
(
∥
A
∥
∞
∥
x
∥
∞
+
∥
y
∥
∞
)
<
16
,
where ϵ is the relative machine precision for the IEEE representation standard; ∥
y
∥
∞ is the infinite norm of a vector: ∥
y
∥
∞
=
max
⁡
1
≤
i
≤
n
{
|
y
i
|
}; and ∥
A
∥
∞ is the infinite norm of a matrix ∥
A
∥
∞
=
max
⁡
1
≤
i
≤
n
{
∑
j
=
1
n
‍
|
A
i
j
|
}.Figure7 presents the graphical representation of Dongarras result (using logarithmic scales) for simple and double precision. For simple precision, ϵ
s
=
2
-
24, for all n
≥
1.05
×
1
0
6 the residual check is always lower than imposed upper bound, similarly for double precision with ϵ
d
=
2
-
53, for all n
≥
5.63
×
1
0
14. If matrix size is greater than these values, it will not be possible to detect if the solution is correct or not. These results establish upper bounds for data volume in this model.Figure 7
Residual check analysis for solvingA
x
=
y system in HPL2.0 using simple and double precision representation.In a single-processor system, the complexity of algorithms depends only on the problem size,n. We can assume T
(
n
)
=
Θ
(
f
(
n
)
), where f
(
n
) is a fundamental function (f
(
n
)
∈
{
1
,
n
α
,
a
n
,
log
⁡
n
,
n
,
…
}). In parallel systems (multiprocessor systems, with p processors) we have the serial processing time T
*
(
n
)
=
T
1
(
n
) and parallel processing time T
p
(
n
). The performance of parallel algorithms can be analyzed using speed-up, efficiency, and isoefficiency metrics.(i)
Thespeed-up, S
(
p
), represents how a parallel algorithm is faster than a corresponding sequential algorithm. The speed-up is defined as S
(
p
)
=
T
1
(
n
)
/
T
p
(
n
). There are special bounds for speed-up [23]: S
(
p
)
≤
p
p
~
/
(
p
+
p
~
-
1
), where p
~
=
T
1
/
T
∞ is the average parallelism (the average number of busy processors given unbounded number of processors). Usually S
(
p
)
≤
p, but under special circumstances the speed-up can be S
(
p
)
>
p [24]. Another upper bound is established by the Amdahls law: S
(
p
)
=
(
s
+
(
(
1
-
s
)
/
p
)
)
1
/
2
≤
1
/
s where s is the fraction of a program that is sequential. The upper bound is considered for a 0 time of parallel fraction.
(ii)
Theefficiency is the average utilization of p processors: E
(
p
)
=
S
(
p
)
/
p.
(iii)
Theisoefficiency is the growth rate of workload W
p
(
n
)
=
p
T
p
(
n
) in terms of number of processors to keep efficiency fixed. If we consider W
1
(
n
)
-
E
W
p
(
n
)
=
0 for any fixed efficiency E we obtain p
=
p
(
n
). This means that we can establish a relation between needed number of processors and problem size. For example for the parallel sum of n numbers using p processors we have n
≈
E
(
n
+
p
log
⁡
p
), so n
=
Θ
(
p
log
⁡
p
).Numerical algorithms use for implementation a hypercube architecture. We analyze the performance of different numerical operations using the isoefficiency metric. For the hypercube architecture a simple model for intertask communication considersT
com
=
t
s
+
L
t
w where t
s is the latency (the time needed by a message to cross through the network), t
w is the time needed to send a word (1
/
t
w is called bandwidth), and L is the message length (expressed in number of words). The word size depends on processing architecture (usually it is two bytes). We define t
c as the processing time per word for a processor. We have the following results.(i)
External product  x
y
T. The isoefficiency is written as
(34)
t
c
n
≈
E
(
t
c
n
+
(
t
s
+
t
w
)
p
log
⁡
p
)
⟹
n
=
Θ
(
p
log
⁡
p
)
.
Parallel processing time is T
p
=
t
c
n
/
p
+
(
t
s
+
t
w
)
log
⁡
p. The optimality is computed using
(35)
d
T
p
d
p
=
0
⟹
-
t
c
n
p
2
+
t
s
+
t
w
p
=
0
⟹
p
≈
t
c
n
t
s
+
t
w
.
(ii)
Scalar product(internal product) x
T
y
=
∑
i
=
1
n
‍
x
i
y
i. The isoefficiency is written as
(36)
t
c
n
2
≈
E
(
t
c
n
2
+
t
s
2
p
log
⁡
p
+
t
w
2
n
p
log
⁡
p
)
⟹
n
=
Θ
(
p
(
log
⁡
p
)
2
)
.
(iii)
Matrix-vector product
y
=
A
x
,
y
i
=
∑
j
=
1
n
‍
A
i
j
x
j. The isoefficiency is written as
(37)
t
c
n
2
≈
E
(
t
c
n
2
+
t
s
p
log
⁡
p
+
t
w
n
p
log
⁡
p
)
⟹
n
=
Θ
(
p
(
log
⁡
p
)
2
)
.Table2 presented the amount of data that can be processed for a specific size. The cases that meet the upper bound n
≥
1.05
×
1
0
6 are marked with (*). To keep the efficiency high for a specific parallel architecture, HPC algorithms for particle physics introduce upper limits for the amount of data, which means that we have also an upper bound for Big Data volume in this case.Table 2
Isoefficiency for a hypercube architecture:n
=
Θ
(
p
log
⁡
p
) and n
=
Θ
(
p
(
log
⁡
p
)
2
). We marked with (*) the limitations imposed by Formula (33).
Scenario
Architecture size (p)
n
=
Θ
(
p
log
⁡
p
)
n
=
Θ
(
p
(
log
⁡
p
)
2
)
1
1
0
1
1.0
×
1
0
1
1.00
×
1
0
1
2
1
0
2
2.0
×
1
0
2
8.00
×
1
0
2
3
1
0
3
3.0
×
1
0
3
2.70
×
1
0
4
4
1
0
4
4.0
×
1
0
4
6.40
×
1
0
5
*
5
1
0
5
5.0
×
1
0
5
*
1.25
×
1
0
7
6
1
0
6
6.0
×
1
0
6
2.16
×
1
0
8
7
1
0
7
7.0
×
1
0
7
3.43
×
1
0
9
8
1
0
8
8.0
×
1
0
8
5.12
×
1
0
10
9
1
0
9
9.0
×
1
0
9
7.29
×
1
0
11The factors that determine the efficiency of parallel algorithms aretask balancing (work-load distribution between all used processors in a system → to be maximized); concurrency (the number/percentage of processors working simultaneously → to be maximized); and overhead (extra work for introduce by parallel processing that does not appear in serial processing → to be minimized).


## 2.3.1. Unfolding Processes in Particle Physics

The theory of unfolding processes in particle physics is as follows [12]. For a physics variable t we have a true distribution f
(
t
) mapped in x and an n-vector of unknowns and a measured distribution g
(
s
) (for a measured variable s) mapped in an m-vector of measured data. A response matrix A
∈
R
m
×
n encodes a Kernel function K
(
s
,
t
) describing the physical measurement process [12–15]. The direct and inverse processes are described by the Fredholm integral equation [16] of the first kind, for a specific domain Ω,
(24)
∫
Ω
‍
K
(
s
,
t
)
f
(
t
)
d
t
=
g
(
s
)
.
In particle physics the Kernel function K
(
s
,
t
) is usually known from a Monte Carlo sample obtained from simulation. A numerical solution is obtained using the following linear equation: A
x
=
b. Vectors x and y are assumed to be 1-dimensional in theory, but they can be multidimensional in practice (considering multiple independent linear equations). In practice, also the statistical properties of the measurements are well known and often they follow the Poisson statistics [17]. To solve the linear systems we have different numerical methods.First method is based on linear transformationx
=
A
#
y. If m
=
n then A
#
=
A
-
1 and we can use direct Gaussian methods, iterative methods (Gauss-Siedel, Jacobi or SOR), or orthogonal methods (based on Householder transformation, Givens methods, or Gram-Schmidt algorithm). If m
>
n (the most frequent scenario) we will construct the matrix A
#
=
(
A
T
A
)
-
1
A
T (called pseudoinverse Penrose-Moore). In these cases the orthogonal methods offer very good and stable numerical solutions.Second method considers the singular value decomposition:(25)
A
=
U
Σ
V
T
=
∑
i
=
1
n
‍
σ
i
u
i
v
i
T
,
where U
∈
R
m
×
n and V
∈
R
n
×
n are matrices with orthonormal columns and the diagonal matrix Σ
=
diag
⁡
{
σ
1
,
…
,
σ
n
}
=
U
T
A
V. The solution is
(26)
x
=
A
#
y
=
V
Σ
-
1
(
U
T
y
)
=
∑
i
=
1
n
‍
1
σ
i
(
u
i
T
y
)
v
i
=
∑
i
=
1
n
‍
1
σ
i
c
i
v
i
,
where c
i
=
u
i
T
y, i
=
1
,
…
,
n, are called Fourier coefficients.


## 2.3.2.  Random Matrix Theory

Analysis of particle spectrum (e.g., neutrino spectrum) faces with Random Matrix Theory (RMT), especially if we consider anarchic neutrino masses. The RMT means the study of the statistical properties of eigenvalues of very large matrices [18]. For an interaction matrix A (with size N), where A
i
j is an independent distributed random variable and A
H is the complex conjugate and transpose matrix, we define M
=
A
+
A
H, which describes a Gaussian Unitary Ensemble (GUE). The GUE properties are described by the probability distribution P
(
M
)
d
M: (
1
) it is invariant under unitary transformation, P
(
M
)
d
m
=
P
(
M
′
)
d
M
′, where M
′
=
U
H
M
U, U is a Hermitian matrix (U
H
U
=
I); (
2
) the elements of M matrix are statistically independent, P
(
M
)
=
∏
i
≤
j
‍
P
i
j
(
M
i
j
); and (
3
) the matrix M can be diagonalized as M
=
U
D
U
H, where U
=
diag
⁡
{
λ
1
,
…
,
λ
N
}, λ
i is the eigenvalue of M and λ
i
≤
λ
j if i
<
j
(27)
Propability
(
2
)

:
P
(
M
)
d
M
~
d
M
exp
⁡
{
-
N
2
T
r
(
M
H
M
)
}
;
Propability
(
3
)

:
P
(
M
)
d
M
~
d
U
∏
i
d
λ
i
∏
i
<
j
(
λ
i
-
λ
j
)
2
×
exp
⁡
{
-
N
2
∑
i
‍
(
λ
i
2
)
}
.The numerical methods used for eigenvalues computation are the QR method and Power methods (direct and indirect). The QR method is a numerical stable algorithm and Power method is an iterative one. The RMT can be used for many body systems, quantum chaos, disordered systems, quantum chromodynamics, and so forth.


## 2.3.3.  Kernel Estimation in HEP

Kernel estimation is a very powerful solution and relevant method for HEP when it is necessary to combine data from heterogeneous sources like MC datasets obtained by simulation and from Standard Model expectation, obtained from real experiments [19]. For a set of data {
x
i
}
1
≤
i
≤
n with a constant bandwidth h (the difference between two consecutive data values), called the smoothing parameter, we have the estimation
(28)
f
^
(
x
)
=
1
n
h
∑
i
=
1
n
‍
K
(
x
-
x
i
h
)
,
where K is an estimator. For example, a Gauss estimator with mean μ and standard deviation σ is
(29)
K
(
x
)
=
1
σ
2
π
exp
⁡
{
-
(
x
-
μ
)
2
2
σ
2
}
,
and has the following properties: positive definite and infinitely differentiable (due to the exp function), and it can be defined for an infinite supports (n
→
∞). The kernel is a nonparametric method, which means that h is independent of dataset and for large amount of normally distributed data we can find a value for h that minimizes the integrated squared error of f
^
(
x
). This value for bandwidth is computed as
(30)
h
*
=
(
4
3
n
)
1
/
5
σ
.The main problem in Kernel Estimation is that the set of data{
x
i
}
1
≤
i
≤
n is not normally distributed and in real experiments the optimal bandwidth it is not known. An improvement of presented method considers adaptive Kernel Estimation proposed by Abramson [20], where h
i
=
h
/
f
(
x
i
) and σ are considered global qualities for dataset. The new form is
(31)
f
^
a
(
x
)
=
1
n
∑
i
=
1
n
‍
1
h
i
K
(
x
-
x
i
h
i
)
,
and the local bandwidth value that minimizes the integrated squared error of f
^
a
(
x
) is
(32)
h
i
*
=
(
4
3
n
)
1
/
5
σ
f
^
(
x
i
)
,
where f
^ is the normal estimator.Kernel estimation is used for event selection to confidence level evaluation, for example, in Markovian Monte Carlo chains or in selection of neural network output used in experiments for reconstructed Higgs mass. In general, the main usage of Kernel estimation in HEP is searching for new particle, by finding relevant data in a large dataset.A method based on Kernel estimation is the graphical representation of datasets using advanced shifted histogram algorithm (ASH). This is a numerical interpolation for large datasets with the main aim of creating a set ofn
b
i
n histograms H
=
{
H
i
}, with the same bin-width h. Algorithm 4 presents the steps of histograms generation starting with a specific interval [
a
,
b
], a number of points n in this interval, and a number of bins and a number of values used for kernel estimation, m. Figure 6 shows the results of kernel estimation if function f
=
-
(
1
/
2
)
x
2 on [
0,1
] and graphical representation with a different number of bins. The values on vertical axis are aggregated in step 17 of Algorithm 4 and increase with the number of bins.Algorithm 4:Advanced shifted histogram (1D algorithm).
(1)procedure
ASH
(
a
,
b
,
n
,
x
,
n
bin
,
m
,
f
m
)
(2)δ
=
(
b
-
a
)
/
n
bin; h
=
m
δ;
(3)for
k
=
1
…
n
bin
do
(4)v
k
=
0; y
k
=
0;
(5)end for
(6)for
i
=
1
…
n
do
(7)k
=
(
x
i
-
a
)
/
δ
+
1;
(8)if
k
∈
[
1
,
n
bin
]
then
(9)v
k
=
v
k
+
1;
(10)end if
(11)end for
(12)for
k
=
1
…
n
bin
do
(13)if
v
k
=
0
then
(14)k
=
k
+
1;
(15)end if
(16)for
i
=
max
⁡
{
1
,
k
-
m
+
1
}
…
min
⁡
{
n
bin
,
k
+
m
-
1
}
do
(17)y
i
=
y
i
+
v
k
f
m
(
i
-
k
);
(18)end for
(19)end for
(20)for
k
=
1
…
n
bin
do
(21)y
k
=
y
k
/
(
n
h
);
(22)t
k
=
a
+
(
k
-
0.5
)
δ;
(23)end for
(24)return
{
t
k
}
1
≤
k
≤
n
bin, {
y
k
}
1
≤
k
≤
n
bin.
(25)end procedureExample of advanced shifted histogram algorithm running for different bins: 10, 100, and 1000.
(a)
(b)
(c)


## 2.3.4.  Performance of Numerical Algorithms Used in Particle Physics

All operations used in presented methods for particle physics (Unfolding Processes, Random Matrix Theory, and Kernel Estimation) can be reduced to scalar products, matrix-vector products, and matrix-matrix products. In [21] the design of new standard for the BLAS (Basic Linear Algebra Subroutines) in C language by extension of precision is described. This permits higher internal precision and mixed input/output types. The precision allows implementation of some algorithms that are simpler, more accurate, and sometimes faster than possible without these features. Regarding the precision of numerical computing, Dongarra and Langou established in [22] an upper bound for the residual check for A
x
=
y system, with A
∈
R
n
×
n a dense matrix. The residual check is defined as
(33)
r
∞
=
∥
A
x
-
y
∥
∞
n
ϵ
(
∥
A
∥
∞
∥
x
∥
∞
+
∥
y
∥
∞
)
<
16
,
where ϵ is the relative machine precision for the IEEE representation standard; ∥
y
∥
∞ is the infinite norm of a vector: ∥
y
∥
∞
=
max
⁡
1
≤
i
≤
n
{
|
y
i
|
}; and ∥
A
∥
∞ is the infinite norm of a matrix ∥
A
∥
∞
=
max
⁡
1
≤
i
≤
n
{
∑
j
=
1
n
‍
|
A
i
j
|
}.Figure7 presents the graphical representation of Dongarras result (using logarithmic scales) for simple and double precision. For simple precision, ϵ
s
=
2
-
24, for all n
≥
1.05
×
1
0
6 the residual check is always lower than imposed upper bound, similarly for double precision with ϵ
d
=
2
-
53, for all n
≥
5.63
×
1
0
14. If matrix size is greater than these values, it will not be possible to detect if the solution is correct or not. These results establish upper bounds for data volume in this model.Figure 7
Residual check analysis for solvingA
x
=
y system in HPL2.0 using simple and double precision representation.In a single-processor system, the complexity of algorithms depends only on the problem size,n. We can assume T
(
n
)
=
Θ
(
f
(
n
)
), where f
(
n
) is a fundamental function (f
(
n
)
∈
{
1
,
n
α
,
a
n
,
log
⁡
n
,
n
,
…
}). In parallel systems (multiprocessor systems, with p processors) we have the serial processing time T
*
(
n
)
=
T
1
(
n
) and parallel processing time T
p
(
n
). The performance of parallel algorithms can be analyzed using speed-up, efficiency, and isoefficiency metrics.(i)
Thespeed-up, S
(
p
), represents how a parallel algorithm is faster than a corresponding sequential algorithm. The speed-up is defined as S
(
p
)
=
T
1
(
n
)
/
T
p
(
n
). There are special bounds for speed-up [23]: S
(
p
)
≤
p
p
~
/
(
p
+
p
~
-
1
), where p
~
=
T
1
/
T
∞ is the average parallelism (the average number of busy processors given unbounded number of processors). Usually S
(
p
)
≤
p, but under special circumstances the speed-up can be S
(
p
)
>
p [24]. Another upper bound is established by the Amdahls law: S
(
p
)
=
(
s
+
(
(
1
-
s
)
/
p
)
)
1
/
2
≤
1
/
s where s is the fraction of a program that is sequential. The upper bound is considered for a 0 time of parallel fraction.
(ii)
Theefficiency is the average utilization of p processors: E
(
p
)
=
S
(
p
)
/
p.
(iii)
Theisoefficiency is the growth rate of workload W
p
(
n
)
=
p
T
p
(
n
) in terms of number of processors to keep efficiency fixed. If we consider W
1
(
n
)
-
E
W
p
(
n
)
=
0 for any fixed efficiency E we obtain p
=
p
(
n
). This means that we can establish a relation between needed number of processors and problem size. For example for the parallel sum of n numbers using p processors we have n
≈
E
(
n
+
p
log
⁡
p
), so n
=
Θ
(
p
log
⁡
p
).Numerical algorithms use for implementation a hypercube architecture. We analyze the performance of different numerical operations using the isoefficiency metric. For the hypercube architecture a simple model for intertask communication considersT
com
=
t
s
+
L
t
w where t
s is the latency (the time needed by a message to cross through the network), t
w is the time needed to send a word (1
/
t
w is called bandwidth), and L is the message length (expressed in number of words). The word size depends on processing architecture (usually it is two bytes). We define t
c as the processing time per word for a processor. We have the following results.(i)
External product  x
y
T. The isoefficiency is written as
(34)
t
c
n
≈
E
(
t
c
n
+
(
t
s
+
t
w
)
p
log
⁡
p
)
⟹
n
=
Θ
(
p
log
⁡
p
)
.
Parallel processing time is T
p
=
t
c
n
/
p
+
(
t
s
+
t
w
)
log
⁡
p. The optimality is computed using
(35)
d
T
p
d
p
=
0
⟹
-
t
c
n
p
2
+
t
s
+
t
w
p
=
0
⟹
p
≈
t
c
n
t
s
+
t
w
.
(ii)
Scalar product(internal product) x
T
y
=
∑
i
=
1
n
‍
x
i
y
i. The isoefficiency is written as
(36)
t
c
n
2
≈
E
(
t
c
n
2
+
t
s
2
p
log
⁡
p
+
t
w
2
n
p
log
⁡
p
)
⟹
n
=
Θ
(
p
(
log
⁡
p
)
2
)
.
(iii)
Matrix-vector product
y
=
A
x
,
y
i
=
∑
j
=
1
n
‍
A
i
j
x
j. The isoefficiency is written as
(37)
t
c
n
2
≈
E
(
t
c
n
2
+
t
s
p
log
⁡
p
+
t
w
n
p
log
⁡
p
)
⟹
n
=
Θ
(
p
(
log
⁡
p
)
2
)
.Table2 presented the amount of data that can be processed for a specific size. The cases that meet the upper bound n
≥
1.05
×
1
0
6 are marked with (*). To keep the efficiency high for a specific parallel architecture, HPC algorithms for particle physics introduce upper limits for the amount of data, which means that we have also an upper bound for Big Data volume in this case.Table 2
Isoefficiency for a hypercube architecture:n
=
Θ
(
p
log
⁡
p
) and n
=
Θ
(
p
(
log
⁡
p
)
2
). We marked with (*) the limitations imposed by Formula (33).
Scenario
Architecture size (p)
n
=
Θ
(
p
log
⁡
p
)
n
=
Θ
(
p
(
log
⁡
p
)
2
)
1
1
0
1
1.0
×
1
0
1
1.00
×
1
0
1
2
1
0
2
2.0
×
1
0
2
8.00
×
1
0
2
3
1
0
3
3.0
×
1
0
3
2.70
×
1
0
4
4
1
0
4
4.0
×
1
0
4
6.40
×
1
0
5
*
5
1
0
5
5.0
×
1
0
5
*
1.25
×
1
0
7
6
1
0
6
6.0
×
1
0
6
2.16
×
1
0
8
7
1
0
7
7.0
×
1
0
7
3.43
×
1
0
9
8
1
0
8
8.0
×
1
0
8
5.12
×
1
0
10
9
1
0
9
9.0
×
1
0
9
7.29
×
1
0
11The factors that determine the efficiency of parallel algorithms aretask balancing (work-load distribution between all used processors in a system → to be maximized); concurrency (the number/percentage of processors working simultaneously → to be maximized); and overhead (extra work for introduce by parallel processing that does not appear in serial processing → to be minimized).


## 3. New Challenges for Big Data Science

There are a lot of applications that generate Big Data, like social networking profiles, social influence, SaaS & Cloud Apps, public web information, MapReduce scientific experiments and simulations (especially HEP simulations), data warehouse, monitoring technologies, and e-government services. Data grow rapidly, since applications produce continuously increasing volumes of both unstructured and structured data. The impact on the approach to data processing, transfer, and storage is the need to reevaluate the way and solutions to better answer the users’ needs [25]. In this context, scheduling models and algorithms for data processing have an important role becoming a new challenge for Big Data Science.HEP applications consider both experimental data (that are application with TB of valuable data) and simulation data (with data generated using MC based on theoretical models). The processing phase is represented by modeling and reconstruction in order to find properties of observed particles (see Figure8). Then, the data are analyzed a reduced to a simple statistical distribution. The comparison of results obtained will validate how realistic is a simulation experiment and validate it for use in other new models.Figure 8
Processing flows for HEP experiments.Since we face a large variety of solutions for specific applications and platforms, a thorough and systematic analysis of existing solutions for scheduling models, methods, and algorithms used in Big Data processing and storage environments is needed. The challenges for scheduling impose specific requirements in distributed systems: the claims of the resource consumers, the restrictions imposed by resource owners, the need to continuously adapt to changes of resources’ availability, and so forth. We will pay special attention to Cloud Systems and HPC clusters (datacenters) as reliable solutions for Big Data [26]. Based on these requirements, a number of challenging issues are maximization of system throughput, sites’ autonomy, scalability, fault-tolerance, and quality of services.When discussing Big Data we have in mind the 5 Vs: Volume, Velocity, Variety, Variability, and Value. There is a clear need of many organizations, companies, and researchers to deal with Big Data volumes efficiently. Examples include web analytics applications, scientific applications, and social networks. For these examples, a popular data processing engine for Big Data is Hadoop MapReduce [27]. The main problem is that data arrives too fast for optimal storage and indexing [28]. There are other several processing platforms for Big Data: Mesos [29], YARN (Hortonworks, Hadoop YARN: A next-generation framework for Hadoop data processing, 2013 (http://hortonworks.com/hadoop/yarn/)), Corona (Corona, Under the Hood: Scheduling MapReduce jobs more efficiently with Corona, 2012 (Facebook)), and so forth. A review of various parallel and distributed programming paradigms, analyzing how they fit into the Big Data era is presented in [30]. The challenges that are described for Big Data Science on the modern and future Scientific Data Infrastructure are presented in [31]. The paper introduces the Scientific Data Life-cycle Management (SDLM) model that includes all the major stages and reflects specifics in data management in modern e-Science. The paper proposes the SDI generic architecture model that provides a basis for building interoperable data or project centric SDI using modern technologies and best practices. This analysis highlights in the same time performance and limitations of existing solutions in the context of Big Data. Hadoop can handle many types of data from disparate systems: structured, unstructured, logs, pictures, audio files, communications records, emails, and so forth. Hadoop relies on an internal redundant data structure with cost advantages and is deployed on industry standard servers rather than on expensive specialized data storage systems [32]. The main challenges for scheduling in Hadoop are to improve existing algorithms for Big Data processing: capacity scheduling, fair scheduling, delay scheduling, longest approximate time to end (LATE) speculative execution, deadline constraint scheduler, and resource aware scheduling.Data transfer scheduling in Grids, Cloud, P2P, and so forth represents a new challenge that is the subject to Big Data. In many cases, depending on applications architecture, data must be transported to the place where tasks will be executed [33]. Consequently, scheduling schemes should consider not only the task execution time, but also the data transfer time for finding a more convenient mapping of tasks [34]. Only a handful of current research efforts consider the simultaneous optimization of computation and data transfer scheduling. The big-data I/O scheduler [35] offers a solution for applications that compete for I/O resources in a shared MapReduce-type Big Data system [36]. The paper [37] reviews Big Data challenges from a data management respective and addresses Big Data diversity, Big Data reduction, Big Data integration and cleaning, Big Data indexing and query, and finally Big Data analysis and mining. On the opposite side, business analytics, occupying the intersection of the worlds of management science, computer science, and statistical science, is a potent force for innovation in both the private and public sectors. The conclusion is that the data is too heterogeneous to fit into a rigid schema [38].Another challenge is the scheduling policies used to determine the relative ordering of requests. Large distributed systems with different administrative domains will most likely have different resource utilization policies. For example, a policy can take into consideration the deadlines and budgets, and also the dynamic behavior [39]. HEP experiments are usually performed in private Clouds, considering dynamic scheduling with soft deadlines, which is an open issue.The optimization techniques for the scheduling process represent an important aspect because the scheduling is a main building block for making datacenters more available to user communities, being energy-aware [40] and supporting multicriteria optimization [41]. An example of optimization is multiobjective and multiconstrained scheduling of many tasks in Hadoop [42] or optimizing short jobs [43]. The cost effectiveness, scalability, and streamlined architectures of Hadoop represent solutions for Big Data processing. Considering the use of Hadoop in public/private Clouds; a challenge is to answer the following questions: what type of data/tasks should move to public cloud, in order to achieve a cost-aware cloud scheduler? And is public Cloud a solution for HEP simulation experiments?The activities for Big Data processing vary widely in a number of issues, for example, support for heterogeneous resources, objective function(s), scalability, coscheduling, and assumptions about system characteristics. The current research directions are focused on accelerating data processing, especially for Big Data analytics (frequently used in HEP experiments), complex task dependencies for data workflows, and new scheduling algorithms for real-time scenarios.


## 4. Conclusions

This paper presented general aspects about methods used in HEP: Monte Carlo methods and simulations of HEP processes, Markovian Monte Carlo, unfolding methods in particle physics, kernel estimation in HEP, Random Matrix Theory used in analysis of particles spectrum. For each method the proper numerical method had been identified and analyzed. All of identified methods produce data-intensive applications, which introduce new challenges and requirements for Big Data systems architecture, especially for processing paradigms and storage capabilities. This paper puts together several concepts: HEP, HPC, numerical methods, and simulations. HEP experiments are modeled using numerical methods and simulations: numerical integration, eigenvalues computation, solving linear equation systems, multiplying vectors and matrices, interpolation. HPC environments offer powerful tools for data processing and analysis. Big Data was introduced as a concept for a real problem: we live in a data-intensive world, produce huge amount of information, we face with upper bound introduced by theoretical models.



---

*Source: 507690-2014-02-26.xml*