Master thesis
Fast Magnetic Field Query Algorithm for
the ALICE O2 Project
Shuto Yamasaki
M164026
Hiroshima University, Graduate School of Science,
Department of Physical Science, Quark Physics Laboratory
Supervisor: Prof. Toru Sugitate
Chief Examiner: Prof. Toru Sugitate
Vice-Examiner: Associate Prof. Tsunefumi Mizuno
February 28, 2018
Abstract
ALICE experiment studies strongly interacting matters in the LHC
energy by analyzing unique physics in heavy ion collisions. The large
amount of data (1TB/s) that will be collected to measure the rare physics
events in ALICE Upgrade project makes record every collision nearly impossible. To address this issue, O2 （Online-Oﬄine）project was planned
to build a new event collecting and processing infrastructure.
Though it is expected to reduce the amount of data to be transferred
and recorded by selectively recording only the physical events excluding
the background by reconstruction of particle tracks. The conventional
data collection program has a issue of speed to calculate the magnetic
ﬁeld applied to the particles inside the magnet.
In this thesis, we proposed a faster implementation for the spacesolenoid electromagnet applied, and for the space dipole magnet applied, respectively in ALICE detector. As a result, the solenoid region
achieved 8.2 times speedup and the dipole region achieved 4.0 times
speedup.
Conventional implementation has used corrected measured ﬁelds as
its sources. It has divided the detector space in order to balance speed
and precision. It has used approximation algorithm which can approximate any function in arbitrary precision.
In the proposed implementation, the solenoid region is approximatingby at most 3rd order polynomial and by using simpliﬁed space division
scheme to reduced search time.

1

In the dipole region, the approximation polynomial forms are restoredfrom the array of coeﬃcients given to the algorithm of the conventional
implementation, and it is modiﬁed and pre-compiled to minimize thecomputation cost, thereby speeding up without decreasing much precision. Segment search in dipole region was done conventionally withbinary search. We replaced it with time complexity O(1) algorithm at thecost of memory to gain speed.
We compared error and speed between the conventional and the proposed implementations. Comparison was conducted to conﬁrm practicality. Proposed implementation will be used in both data collection andsimulation for the ALICE experiment in the future.

Contents
1 Introduction
1.1 Quark Gluon Plasma (QGP) . . . . . . . . . . . . . . . . . . . . .
1.2 Physics program with the upgraded ALICE detector . . . . . .
1.3 Technical background . . . . . . . . . . . . . . . . . . . . . . . .
1.4 Implementation diﬃculties . . . . . . . . . . . . . . . . . . . . .
1.5 Purpose . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5
6
8
11
15
16

2 Magnetic ﬁeld calculation inside ALICE detector
16
2.1 Source measurement . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3 Methods
18
3.1 Solenoid Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1.1 Approximation with linear regression model . . . . . . . 18
3.1.2 Space division toward improved precision . . . . . . . . 19
3.2 Dipole Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2.1 Loop and Branch elimination . . . . . . . . . . . . . . . . 21
3.2.2 Extra multiplication reduction and Fused Multiply Add . 22
3.2.3 Constant time complexity O(1) segment search . . . . . 23
3.2.4 Quick segment search: Preparation . . . . . . . . . . . . 27
4 Results
29
4.1 Result of Solenoid Region . . . . . . . . . . . . . . . . . . . . . . 29
4.2 Result of Dipole Region . . . . . . . . . . . . . . . . . . . . . . . 30
5 Discussions
30
5.1 Trade-oﬀ between polynomials calculation speed and precision 30
5.2 Comparison of regression or approximation methods other thanmultiple linear regression . . . . . . . . . . . . . . . . . . . . . . 32
5.3 Memory usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
6 Conclusion

34

2

References

35

A Deﬁnition of the Chebyshev polynomials Tn

36

B Indices of precision evaluation

37

C Programs

37

List of Figures
1.1 Schematic generation mechanism of the QGP. Quarks are conﬁned by strong force in room tempreature. Quarks are deconﬁned in high tempreature and high density environment andundergo phase transition to QGP phase (ﬁg.1.2 top.) . . . . . .
1.2 QCD phase diagram. . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 The ALICE detector after upgrade (Run3) . . . . . . . . . . . . .
1.4 Data ﬂow and processing pipeline of the O2 system. ITS/TPC
track ﬁndings are performed online (synchronously) at the Event
Processing Nodes (EPNs). Track ﬁnding is also performed in
the Muon CHember (MCH) and Muon Forward Tracker (MFT)
detector. [12] . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.5 Possible hardware implementation of the logical data ﬂow.
The EPNs compress 500GB/s input data into 90GB/s.[12] . .
1.6 Particle tracks, projected on a plane perpendicular to beamaxis (Z), of pp collisions measured with TPC detector. Redlines are physics tracks and black lines are background tracks.
Nearly half tracks are unnecessary because these are comefrom gas molecules in the detector and electrons in the beam
pipe. In the ﬁeld of L3 solenoid magnet, some particles with
about 50M eV /c of momentum are moving circularly. . . . . . .
2.1 Schematic ﬂow of the existing B-ﬁeld calculation. Given a position (x, y, z) in the detector volume, searches corresponding
subvolume (segment) ﬁrst. Then calculate B-ﬁeld with the
Chebyshev polynomials using the coeﬃcients from the segment struct. The last used segment is cached for speed because it is likely near points are queried continuously in track
ﬁnding. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.1 Partitioning of the solenoid volume into chambers. (i) Oneof the chambers. (ii) The volume is divided by the detector
boundary radius and quadrants. (iii) The volume is divided by
the sign of the Z axis. . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Bx of conventional implementation near Z = 0 (−25 ≤ Z ≤ 30)
on the X-Y plane, plots of 5cm step. The ﬁeld is ﬂipped near
Z = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

7
8
9

12
13

14

17

20

21

3.3 By of conventional implementation near Z = 0 (−25 ≤ Z ≤30) on the X-Y plane, plots of 5cm step. The ﬁeld is ﬂipped
near Z = 0. Additionally, ϕ -dependency is remarkable at Z =−5, 0, 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.4 Code snippet of a generated Bz function. (shortened to 6 digits
for visibility) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.5 Comparison of assemblies w/ and w/o Fused Multiply Add(FMA) intrinsics. Left: w/o FMA, right: w/ FMA. Instructions
which starts with ʼfmaʼ computes a = ax + b in 5 cycles of latency. On the other hand, mulss and addss takes 5, 3 cyclesof latency respectively, the former is 1.6 times faster in simple
calculation. Used compiler is Clang-900.0.39.2 on Intel Core
i7 Haswell Refresh. . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.6 Split dipole region −1760 ≤ Z ≤ −558.3 into 89 segments by
maximum and minimum Z values of all the chambers. Calla segment with minimum Z as zid = 0, and maximum Z as
zid = 88. A crossing line is placed to show possibility of linear
relationship between Z positions and zids. The gap between
the line and the segments implies Interpolation search works
less eﬃciently. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.7 Schematic ﬁgure of Quick Segment Search algorithm. The
points are just an example and nothing to do with real data.
Prepare evenly spaced slices (orange arithmetic sequence)
which contain only one or zero Z value of a chamber border(blue point) in advance. The id of one of evenly spaced slices
is immediately calculated from given z (e.g. 6). Every slice
knows a z position of the blue point and a zid of the segment
which is smaller than the blue point. The zid of larger side
is always (the zid of smaller side + 1). Then by comparing z
and z of the blue point, a zid for queried z is obtained (e.g. 3).
Thus, zid is used to lookup the other value (data of X,Y axis oran index of a polynomial). . . . . . . . . . . . . . . . . . . . . . . 27
3.8 Parametrized part of dipole magnet applied region. Coloredcubes are the chambers of dipole ﬁeld parametrization from
AliMagF class. This parametrization contains 1482 chambers. 28
4.1 Positional dependency of mean precision per magnetic ﬁeld
component in solenoid region. . . . . . . . . . . . . . . . . . . . 30
5.1 Trade-oﬀ between polynomials calculation speed and precision 31

4

List of Tables
1

2

3
4
5
6

7

Comparison of the physics reach, minimum accessible pT andrelative statistical uncertainty, for selected observables between the approved (current) scenario and the proposed upgrade scenario. [3] . . . . . . . . . . . . . . . . . . . . . . . . . .
Estimated signal per event (S/ev), signal-to-background ratio(S/B) and number of background candidates per event (Bʼ/ev)√
for central Pb-Pb collisions at sNN = 5.5 TeV[12]. These
values are derived in the Conceptual Design Report and in
the Technical Design Report for the Inner Tracking System
upgrade documents[5], where the geometrical conditions arealso described. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Data reduction plan for the Pb-Pb collision 50kHz[12]. . . . . .
Consumed CPU time for single query of one random point. . .
Importance rank obtained with Recursive Feature Elimination
technique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Fitting results of regression methods available in scikit-learn(Python) library sorted on coeﬃcient of determination. ʼ-ʼ indicates NODATA. Percentages in the table is ratio against target
precision 10−3 ×Bz = 5e−03; less than 100% means the ﬁtting
successfully reached target precision. . . . . . . . . . . . . . .
Memory usage comparison between the conventional and the
proposed implementation. Both 2 kilo Gauss and 5 kilo Gausscases are shown. The values are for reference. See text for
detail. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

12
13
29
31

33

34

List of Algorithms
1
2
3
4

1

Prepopulated Quick Segment Search (search part). See text
for details. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Prepopulated Quick Segment Search (construction part). See
text for details. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Mathematica reimplementation of the AliMagF class. It is used
to extract polynomials from existing parametrizations. . . . . .
Mathematica C-source generator of the B-ﬁeld parametrization in HornerForm. . . . . . . . . . . . . . . . . . . . . . . . . . .

27
28
38
39

Introduction

ALICE experiments are studying matters that interact with strong force bysearching for physics peculiar to heavy ion collisions at the energy of LHC
accelerator. In the prediction of lattice Quantum ChromoDynamics (lQCD),

5

quarks are released from conﬁnement due to strong interaction in a hightemperature and high energy density environment as in the early universe,
and a new state of matter Quark Gluon Plasma (QGP) is generated. Studies on the speciﬁc viscosity and transport coeﬃcient of QGP reproduced byenergy heavy ion collision and the mechanism of formation of quark of the
QGP derived from QGP are progressing. If heavy quark interacts on medium,its mass and ﬂavor are hardly changed, so QGP in the ALICE upgrade plan,by raising the frequency of heavy ion collision to 50 kHz by 2021, it is con− +sidered that D0 → K − π + ，Λ+c → pK π , and other rare physical events. In
the collection of scarce physical events scheduled to be measured in ALICE
experiments, it is diﬃcult to record everything because data of 1 TB/s ormore is generated. Therefore, the O2 (Online-Oﬄine) plan was drafted forthe purpose of building a new event collection and processing infrastructure.
It is expected to reduce the amount of data to be transferred and recorded by
selectively recording only the physical events excluding the background byreconstruction of particle tracks, but with the conventional data collection
program, To 5% was spent in the calculation of the magnetic ﬁeld applied
to the particles inside the magnet, which was regarded as a problem. If
the calculation speed can not catch up with the occurrence of data, there isa possibility that it is necessary to reduce the collected data and respond.
Therefore, speeding up reconﬁguration has become important. The current
implementation of the magnetic ﬁeld calculation above has already been undertaken by a general-purpose optimization approach, and further speeding
up method is not trivial. In this paper, we propose a method that can speed
up while maintaining the precision of the current implementation to a certain extent (or perfectly), without relying on this special speed hardware,
And compare the speed with the current implementation.

1.1

Quark Gluon Plasma (QGP)

A nucleus is a ball of protons and neutrons, generally referred to as nucleons, which is consisted of three quarks. These quarks are bound by gluons
so that they cannot be taken out from a hadron at ordinary environment (ﬁg.
1.1). This is called conﬁnement. Quantum Chromodynamics (QCD) calculations suggests strongly interacting particles (hadrons) should cause phasetransition at extreme environment that is high temperature (150-200MeV)
and high density (> 1GeV/fm3 ). In the new state of matter, quarks and gluons move freely without conﬁnement. We call the hot and dense medium
as Quark-Gluon Plasma (QGP). Some predictions says critical temperature
for phase transition is around 200 MeV (∼ 2 × 1012 K) and the QGP existed
at t = 106 to 105 seconds from the Big Bang. Thus, study of the characteristics of the QGP lead to our knowledge of the early universe. Quark-Gluon
Plasma is also important for understanding mass generation mechanism of
hadrons. The mass of the quark constituting the nucleon is estimated tobe ∼ 20 × 10−30 kg, and there is a big diﬀerence between triple the quark

6

Figure 1.1: Schematic generation mechanism of the QGP. Quarks are conﬁned by strong force in room tempreature. Quarks are deconﬁned in hightempreature and high density environment and undergo phase transition to
QGP phase (ﬁg.1.2 top.)mass and proton mass 1700 × 10−30 kg. Most of our mass is thought to begenerated by ”spontaneous break of chiral symmetry” caused by the strong
interaction of QCD, and QGP is a mass gain mechanism as a ﬁeld wherechiral symmetry is also restored It will bring important insights to.
High energy heavy ion collisions are used to generate Quark Gluon Plasma
in laboratory and to investigate its properties. In this way, heavy ion like Pb
are collided in high speed so that a hot and dense matter is created there.
Such kind of QGP generation experiments has undergone since 1980s: anexperiment using the BEVALAC accelerator of Lawrence Berkeley Laboratory in the United States, the AGS accelerator of Brookhaven National Laboratory (BNL) of the United States, and the SPS accelerator of the European Common Nuclear Research Organization (CERN). In the past, BNLAGS gold beams with 10 GeV per nucleus, CERN-SPS with 200 GeV sulfur
beams per nucleus and 160 GeV per nucleon beam ﬁxed target experiments
have been conducted. The worldʼs ﬁrst hadron collider BNL-RHIC accelerator has started since 2000, and the European CERN-LHC accelerator has
started to operate since 2009. Veriﬁcation of the QGP with 10 to 100 times
the collision energy was done until now. The experimental goals are below.
• Discovery of the QCD phase diagram• Elucidation of properties of hot and dense QCD multi-body system• Precise veriﬁcation of QCD in non-perturbative regions• Elucidation of mass acquisition mechanism of quarks• Understanding the conﬁnement mechanism in hadron

7Temperature (MeV)

250
200
150

l
Primordia -6
1(< 0 s)
Universe

Quark-Gluon Plasma

Hadronic gas

Ph

p

100

as

e
Tr
an
sit
ion
"Solid" state

K

50

Neutron Sta

rs"Ordinary" state

Net baryonic density
0

5

1

8

(normalised, d/d0)

Figure 1.2: QCD phase diagram.
• Understanding the QCD phase transition happened in the early universe after the Big Bang
Several evidence of the QGP was found as a result of the previous experiments (Run 1 and Run 2). Representative results include discovery of hightransverse momentum (high-pT ) hadron yield suppression (aka. jet quenching), azimuthal anisotropy and ﬂow, thermal photon observation.
Currently it goes beyond the phase of discovery of QGP and researcheson physical properties of QGP are in progress. In the future, the following
topics are to be understood.
• Temperature dependency of QGP quantities (Speciﬁc viscosity, transportation coeﬃcients, and so on)• QGP response (behavior and propagation of energy lost by jet)
• Initial condition of collisions, early thermalization mechanism• Hadron mass acquisition mechanism• QCD phase structure
For above advanced study, more accurate measurements and new measurements that were impossible precisely are indispensable.

1.2

Physics program with the upgraded ALICE detector

LHC will be upgraded by 2022 and develop into high brightness LHC with increased collision frequency (HL - LHC). To further study the physical properties of QGP in this second generation LHC, we focus on rare probes of8

Figure 1.3: The ALICE detector after upgrade (Run3)
heavy ion physics in experiments after Long Shutdown 2 (LS 2) (Run 3, Run
4). The research subject is the coupling and the hadronization process of
QGP and the medium, and it consists of heavy ﬂavor particles, quartzone
state (probe of temperature of QGP), real photon, virtual photon, jet including heavy quark and jet with high lateral momentum , Correlation between
jet and other probes, measurement of lepton pair, etc. For example, in orderto enable these high precision measurements, the LHC-ALICE experiment isplanning to construct a large-scale new measuring instrument. Comparisonof measurement precision before and after upgrade [3] is shown in the table1.
The upgraded ALICE detector will target for an integrated luminosity of
10nb−1 at full magnetic ﬁeld (B = 0.5 T) in the ALICE solenoid and 3nb−1
with reduced ﬁeld (B = 0.2 T) for Pb-Pb collisions. The requirement for an
integrated luminosity of 13nb−1 is motivated mainly by the following performance ﬁgures as described in the ALICE upgrade document[12].
• Heavy ﬂavour measurements‒ The nuclear modiﬁcation factor RAA and elliptic ﬂow v2 of strangecharmed mesons (Ds ) down to a transverse momentum pT of atleast 2 GeV/c with a statistical precision better than 10% for bothobservables. This will allow a precise comparison of strange andnon-strange charm meson dynamics;‒ Λc baryon RAA and v2 down to 2 GeV/c and 3 GeV/c, respectively,
with a precision of about 20% and baryon/meson ratio for charm(Λc /D) down to 2 GeV/c with the same precision. This will allow
to address the charm quark hadronization mechanisms at low andintermediate momentum;‒ RAA and v2 of beauty-decay particles via non-prompt D0 , nonprompt J/ψ and beauty-decay leptons (the two latter at both cen9

Observable

D meson RAA
D meson from B decays RAA
D meson elliptic flow
D meson from B elliptic flow
Charm baryon-to-meson ratio
Ds meson RAA
J/ψ RAA (forward rapidity)
J/ψ RAA (mid rapidity)
J/ψ elliptic flow
ψ(2S) yield
Temperature (intermediate mass)
Elliptic flow
Low-mass spectral function
Hyper(anti)nuclei4Λ H yield
Hyper(anti)nucleiΛΛ4 H yield

ApprovedpAmin
statistical
T(GeV/c)uncertainty
Heavy Flavour1
10 % at pAmin
T3
30 % at pAmin
T1
50 % at pAmin
T
not accessible
not accessible
4
15 % at pAmin
T
Charmonia
0
1 % at 1 GeV/c
05 % at 1 GeV/c
015 % at 2 GeV/c
030 %Dielectrons
not accessible
not accessible
not accessible
Heavy Nuclear States
35 %not accessible

pUmin
T(GeV/c)Upgradestatisticaluncertainty

0
2
0
2
2
1

0.3 % at pAmin
T
1 % at pAmin
T
2.5 % at pAmin
T
20 % at pUmin
T
15 % at pUmin
T
1 % at pAmin
T

0
0
0
0

0.3 % at 1 GeV/c
0.5 % at 1 GeV/c
5 % at 2 GeV/c
10 %

0.3

10 %
10 %
20 %
3.5 %
20 %

Table 1: Comparison of the physics reach, minimum accessible pT and relative statistical uncertainty, for selected observables between the approved(current) scenario and the proposed upgrade scenario. [3]10tral and forward rapidity) down to 1‒2 GeV/c with precisions froma few percent to 10%. This will allow a detailed assessment of theb quark transport properties in the medium;‒ B meson fully-reconstructed decays (B + → D0 π + ) down to 3
GeV/c with a precision of about 10%. This will provide an important direct measurement of beauty production;‒ Λb baryon production for pT > 7GeV/c. This will be a unique measurement in heavy ion collisions and should allow for the determination of the nuclear modiﬁcation factor of the beauty baryon,
which is sensitive to the b quark hadronization mechanism.
• Charmonium measurements‒ RAA of J/ψdown to pT = 0 with statistical precision better than
1%, at both central and forward rapidity;‒ RAA of ψ(2S) down to pT = 0 with precision of about 10%, at both
central and forward rapidity;‒ v2 of J/ψ down to pT = 0 with precision of about 0.05 (absolute
uncertainty), at both central and forward rapidity;
these measurements will allow a detailed investigation of the mechanisms of dissociation and regeneration for charmonium states inthe deconﬁned medium.
• Low-mass dileptons‒ the additional sample of 3nb−1 with a reduced magnetic ﬁeld value(0.2 T) in the central barrel is essential for low-mass dielectron
analysis to obtain the projected precision of about 10% on theslope of the high-invariant-mass region and of about 10% on the
dielectron elliptic ﬂow. This measurement will make it possible toassess the time-evolution of the thermal radiation emitted by the
hot medium.
Most of these analysis (with the exception of the exclusive reconstruction of
beauty hadron decays) are characterized by a very small signal-to-backgroundratio. It implies signal candidate events in Pb-Pb collisions are indistinguishable from background events. This means that it is not possible to use somededicated triggers to select collisions online for oﬄine analysis. Instead, aminimum bias trigger must be used to “record all collisions of Pb-Pb”. A
minimum bias trigger is the most basic trigger used to record all the events.
Table 2 summarize such branches.

1.3

Technical backgroundTo tackle the data recording challenge, the new computing infrastructures
has planned: Online-Oﬄine (O2 ) Project. Fig. 1.5 shows O2 hardware and11Analysis0
D → K − π++
Ds → K + K − π +− +Λ+c → pK π+− +Λc → pK π (pT > 2GeV/c)

S/ev [5]7.6 · 10−3
2.3 · 10−3
6.5 · 10−4
3.7 · 10−4

S/B
10−2
< 2 · 10−3
< 10−42 · 10−4

B ′ /ev
3.0
> 4.6
> 26
7.4

Table 2: Estimated signal per event (S/ev), signal-to-background ratio (S/B)
and number of background candidates per event (Bʼ/ev) for central Pb-Pb√
collisions at sNN = 5.5 TeV[12]. These values are derived in the Conceptual
Design Report and in the Technical Design Report for the Inner Tracking
System upgrade documents[5], where the geometrical conditions are also
described.

Adjusted accounting for current luminosity

procedures
Adjusted withmultiplicity

multiplicity
reﬁttng

Vertexing

Local processing

DCS data

calorimeters

matching
PID calibrations

clusterization
ﬁnding…
trigg
association)

CTF

Step 0

Step 1

AOD

Step 2

Step 3

Step 4

Figure 1.4: Data ﬂow and processing pipeline of the O2 system. ITS/TPC
track ﬁndings are performed online (synchronously) at the Event Processing
Nodes (EPNs). Track ﬁnding is also performed in the Muon CHember (MCH)and Muon Forward Tracker (MFT) detector. [12]12Detectors

250 FLPs

1.1 TB/s
8300 Readout
Links

Switching
Network

1500 EPNs

Storage 34 Storage
Storage
Network Servers

500 GB/s

90 GB/s

Input: 250 ports
Output : 1500 ports

1500 x 60MB/s

Figure 1.5: Possible hardware implementation of the logical data ﬂow. The
EPNs compress 500GB/s input data into 90GB/s.[12]Detector
TPC
ITS
TRD
TOF
MCH
MFT
EMCal
PHOS
(Total)Raw rate (GB/s)
1,000
40
20
2.5
2.2
10
4
2
1080.7

Compressed (GB/s)
50
26
3
2
0.7
5
1
0.5
88.2

Reduction (1/x)
20
1.5
6.6
1.25
2.9
2
4
4
12.3

Table 3: Data reduction plan for the Pb-Pb collision 50kHz[12].

13Figure 1.6: Particle tracks, projected on a plane perpendicular to beam axis
(Z), of pp collisions measured with TPC detector. Red lines are physics
tracks and black lines are background tracks. Nearly half tracks are unnecessary because these are come from gas molecules in the detector andelectrons in the beam pipe. In the ﬁeld of L3 solenoid magnet, some particles with about 50M eV /c of momentum are moving circularly.

14

network communication between them. The data ﬂowing from the detector to the ﬁrst stage processor (First Level Processor) group reaches 1.1
terabytes per second, which is diﬃcult to record as it is. Therefore, in O2 ,background information irrelevant to physics is immediately discarded, andother data is compressed and recorded. The value read from the detector
which is the original data is basically discarded. As a result, data amountreduction as shown in the table 3 is scheduled. The original data is organized and compressed from 1.1 TB/s to 90 GB/s. An element which is very
important in this data compression process is calculation of particle trajectory. Figure 1.6 lists the particle trajectory plot of the TPC detector locatedin the solenoid region. Half of that is the background, giving the opportunity
to reduce the data volume by half. As a function of the detector software
used to obtain the momentum from the curvature of the charged particle forthe calculation of the trajectory there is an inquiry of the static magnetic
ﬁeld vector applied inside the detector. Although this function is necessary,it took a little less than 5% to acquire the magnetic ﬁeld out of the CPU time
involved in the reconstruction of the track, which was regarded as a problem. In this paper, the calculation of this magnetic ﬁeld has been speeded
up.

1.4

Implementation diﬃculties

There are some restrictions when implementing acceleration of magnetic
ﬁeld calculation.
First, since this calculation is performed very frequently, it is required tocomplete at a speed of 0.1 μs order. Therefore, it is not possible to use
inference by a neural network which depends on many tensor products.
Next, the precision of the calculated magnetic ﬁeld must maintain highprecsion of ±10−3 Bz (Bz = 2, 5 kG). Its domain covers a wide range of 11 +12 m in length and 6 - 10 m in diameter.
This magnetic ﬁeld is a non-uniform and nonlinear magnetic ﬁeld basedon measured values, and it is diﬃcult to handle with a single equation.
Attempts to obtain precision by directly applying several statistical machine learning methods (this is an analytical approach diﬀerent from socalled AI) failed.
Therefore, one can think of an approach to lower the amount of calculation by changing the algorithm of existing implementation code. However,as described later, existing implementations have already been considerably
optimized, including consideration for caching, so the improvement is nottrivial. For example, the binary search method used in existing implementations will be the fastest category as a list search method for real-valued
divisions where the distribution is biased.
If coping with software is diﬃcult, one will consider using dedicatedhardware such as GPU, FPGA, ASIC. In O2 , we plan to install FPGA in FLP
and GPU in EPN. Since track ﬁnding is performed by EPN from ﬁg.1.4, there15

is a possibility that GPU can be used only when using O2 . However, magnetic ﬁeld calculation is also required for Monte Carlo simulation outside
the ALICE computer room and CERN computer center. Therefore, the existence of special hardware can not be assumed, and it must operate on acommercially available ordinal x86̲64 machines.
As with the lattice QCD simulation, it may be possible to arrange a large
amount of lattice points that have previously calculated the magnetic ﬁeld in
the detector volume. However, the amount of available memory is also limited. The magnetic ﬁeld calculation is performed in a distributed computing
environment and runs in many instances. To cope with other physical calculations, only a few MB of memory can be used for magnetic ﬁeld calculation.
According to rough estimation, at least 70 GB per instance will be required
to obtain suﬃcient precision with the above method.

1.5

PurposeThe purpose of this study is to solve the above problem and to facilitate themeasurement of new physical channels by supporting data collection in the
coming Run 3. And to document explicitly one of the high-speed implementation which is faster than binary search which is considered the fastest incommon sense (this is special case, though).

2

Magnetic ﬁeld calculation inside ALICE detector

2.1

Source measurement

Analysis of the magnetic ﬁeld applied by ALICEʼs L3 solenoid electromagnet was made by Ruben Shahoyan[10] and the analysis of the magnetic ﬁeld
applied by the dipole electromagnet was done by Ruben and Morsch[11]. In
the conventional implementation, internally the magnetic ﬁeld is a set of parameterization (3 dimensional to 3 dimensional mapping) based on Chebyshev polynomial in the cylindrical coordinate system (see appendix for its
deﬁnition), deﬁned for each divided chember of the magnetic ﬁeld.
The degree of the polynomial is automatically adjusted so as to guarantee
the required precision (the maximum value of the diﬀerence from the measured ﬁeld). The precision varies depending on the location, but it is about
2 × 10−4 ∼ 5 × 10−3 (solenoid region, in the case of applying 5 k gauss).
Where the data was missing, supplementation by a spline curve is performed. There are two sets of data: 12 kA (low, 2 k Gauss) and 30 kA (high,
5 k Gauss). These dataset are the most accurate data set available.
When referring to conventional (existing) implementation in this paper,
it refers to AliMagF class[2] implemented by Ruben and these data sets.

16

(x, y, z)
Check cached segment
Mismatch
Match
Binary Search (Z)Segment SearchLinear Search (X, Y)
Cache the segment
Chebyshev polynomials

Field Calculation(Bx, By, Bz)
Figure 2.1: Schematic ﬂow of the existing B-ﬁeld calculation. Given a position (x, y, z) in the detector volume, searches corresponding subvolume(segment) ﬁrst. Then calculate B-ﬁeld with the Chebyshev polynomials
using the coeﬃcients from the segment struct. The last used segment iscached for speed because it is likely near points are queried continuously intrack ﬁnding.

2.2

Overview

When a magnetic ﬁeld inquiry of a certain point inside the detector occurs,
the ﬂow until the magnetic ﬁeld is calculated by the existing implementation is shown in the ﬁgure 2.1. The given point (x, y, z) is checked to see if
it is contained in the ﬁrst cached chamber. If it matches, skip search and
go to calculation of Chebyshev polynomial. If it does not match, it ﬁnds thesegment to which z belongs by binary searching the array sorted by the maximum and minimum values (boundary values) of the Z axis of the chamber.
Subsequently, an array obtained by sorted the boundary values of the X axis
associated with that division is searched by linear search. Likewise, the Yaxis is linearly searched to ﬁnd chambers.
Binary search is considered to be the fastest as a search method for arbitrary length sorted real number list (Z axis boundary values) with unequal
distribution. The found chamber is cached for the next time. Thereafter, a
Chebyshev polynomial is calculated. Chebyshev polynomials are orthogonal
systems, and any function can be expressed when the number of terms is
inﬁnite. Unlike the Fourier transform it is very fast because division is not17

used in it.
Thus, existing implementations are already optimized for speed.

3

Methods

Bottleneck of conventional method are chamber search and calculation ofcomplete system (Chebyshev polynomial) describing the magnetic ﬁeld of
the chamber. Both took the same degree of time. In the conventional implementation, the solenoid region and the dipole region were calculated in thesame way, but in the proposed implementation, the optimization was performed in diﬀerent ways due to the diﬀerence in the gradient of the magnetic
ﬁeld.

3.1

Solenoid Region

To decrease the ﬁeld query overhead, High Level Trigger (HLT) has usedless precise and faster parametrization from before. In this method, only Bz
which is the main component of the solenoid region is computed as:bz = c[0] + c[1] ∗ z + c[2] ∗ r + c[3] ∗ z ∗ z + c[4] ∗ z ∗ r + c[5] ∗ r ∗ r;
where c is an array of coeﬃcients, z is distance in the beam line directionin the ALICE local coordinate system with the Interaction Point as origin, andr is the distance from the beam line. As seen above, φ-dependency of Bz andtransverse components (Bx , By ) are ignored. This is not a issue because
current HLT reconstruction does not require much precision, but the result
cannot be used in ﬁnal ﬁtting. Therefore, we extended the parametrizationof HLT, consider the φ-dependency and transverse components and createdan implementation that can be used for the ﬁnal ﬁtting. The target precision
is set to 10−3 ×Bz for the value returned by the conventional implementation.
3.1.1

Approximation with linear regression model

In the solenoid region, we regard the conventional implementation as a blackbox and create expressions that reproduce its behavior (aka. supervised
learning). Approximate conventional implementation with low order real coeﬃcient polynomial by regarding it as pure function B(x, y, z) . The new
implementation of the solenoid region is practically a ﬁtting to the conventional implementation. The existing implementation adopted the cylindricalcoordinate system (r, ϕ, z) for the solenoid region, so this ﬁtting makes it
possible to omit the coordinate system transformation which includes time
consuming ϕ calculation. The problem of ﬁnding a function that reproduces
data is called a regression problem. The linear regression model is the most18

basic model. First, the objective variable Bx is represented by a linear sumof coordinates (explanatory variables).
Bx = w0 + w1 x1 + w2 x2 + w3 x3 + ... + wm xm = wT x
Where wi ∈ R is weight, w is a vector (w0 , w1 , ..., wm ) and x is a vector(x1 , x2 , ..., xm ). In normal linear regression, x = (x, y, z), but this time theexplanatory variable contains not only the given coordinates x, y, z but also
the product of coordinates (x2 , xy, z 2 x , etc). This method is called polynomial regression. The number of explanatory variables must be chosen suchthat suﬃcient precision is obtained and there is no speed degradation. As
a result of ﬁtting and measuring the precision, it was found that at mosta third order real coeﬃcient polynomial is suﬃcient for the gradient of thesolenoid region. Therefore, the expression of Bx is given as follows.
Bx (x, y, z) = a0 + a1 x + a2 y + a3 z + a4 xx + a5 xy + a6 xz + a7 yy + a8 yz + a9 zz+ a10 xxx + a11 xxy + a12 xxz + a13 xyy + a14 xyz + a15 xzz + a16 yyy + a17 yyz + a18 yzz + a19 zzz
Where w = (a0 , a1 , ..., a19 ) . To reduce multiplications and additions,transformed as follows. However, do not change x.
Bx (x, y, z) = a0 + x(a1 + x(a4 + xa10 + ya11 + za12 ) + y(a5 + za14 ))+ y(a2 + y(a7 + xa13 + ya16 + za17 ) + z(a8 ))+ z(a3 + z(a9 + xa15 + ya18 + za19 ) + x(a6 ))Next, randomly collect sample points from existing implementations N
points. Let the pair of explanatory variable vector and Bx corresponding
to the i-th point be (xi , Bx,i ). At this time, care must be taken not to bias
the distribution of the sample points. Simply taking sample points in a grid
shape has scale dependence. Therefore, it is better to sample them at random. Also, when random numbers are obtained according to the geometry(cylindrical shape) of the detector, it is necessary to check whether points
are evenly distributed over the detector volume.
Finally, ﬁnd the weight w using the method of least squares. For calculation, we used sklearn.linear̲model.LinearRegression() of scikit-learn[8].
The approximation precision of this approach is described in the results section.
3.1.2

Space division toward improved precision

When tried to approximated the entire solenoid region by one equation, theprecision is not suﬃcient. So in order to improve it, whole space was divided
into chambers as shown in ﬁg. 3.1. The considerations on the division are
as follows.

19

Figure 3.1: Partitioning of the solenoid volume into chambers. (i) One of thechambers. (ii) The volume is divided by the detector boundary radius andquadrants. (iii) The volume is divided by the sign of the Z axis.
• For R axis
‒ Divided along the detectorsʼ material boundary taking advantageof the fact that the detectors in the barrel are centered around thepoint at R = 0. The gradient is greatly diﬀerent inside and outsideof detector boundaries.
‒ More speciﬁcally, the volume was divided into ﬁve regions: ITS
detector (R ≤ 80 cm), TPC detector (80 < R ≤ 250), TOF detector(250 < R ≤ 400), a gap which is just outside of the TOF detector(400 < R ≤ 423), the volume where there are calorimeters such as
PHOS, EMCal, DCal detector (423 < R ≤ 500).
√‒ r = x2 + y 2 is used for branching depending on R, but by squaring the above constants beforehand, extra calculation of squareroot can be reduced by comparing it with x2 + y 2 .
• For ϕ axis
‒ In some cases ϕ dependency is large like ﬁg. 3.3, and approximation precision improves by dividing in ϕ axis direction.
‒ Simpliﬁed conditional branching by splitting just at 90 degrees.
• For Z axis
‒ Bx , By has symmetry around Z = 0 as in ﬁg. 3.2,3.3, therefore, aswith the ϕ axis, splitting with Z = 0 improves precision.

20

Figure 3.2: Bx of conventional implementation near Z = 0 (−25 ≤ Z ≤ 30)
on the X-Y plane, plots of 5cm step. The ﬁeld is ﬂipped near Z = 0.
‒ −550 ≤ z (cm) ≤ 550 to divide 22 at 50 cm intervals to ensure
suﬃcient precision in the wide region.
‒ Since it is divided at regular intervals, an index of the array can be
speciﬁed with one division, and a penalty (pipeline stall in CPU)due to failure of conditional branch does not occur.
The entire solenoid region is divided into 5 × 4 × 22 = 440 chambers. Since a
simple division method is used, it is possible to identify the chamber whichcontains the point to be queried quickly. 100,000 points were sampled andﬁtted for each chamber.

3.2

Dipole Region

Unlike the solenoid region, the gradient of the dipole region is remarkably
large, and if ﬁtted in the same way, the number of divisions becomes verylarge. Therefore, we did not dare ﬁt in the dipole region, but speed up by
extracting polynomial from existing parametrization and modifying it.
3.2.1

Loop and Branch elimination

The Chebyshev polynomial calculation was implemented with a code having
a triple loop and a branch within the loop. Since branch prediction is mainly
21

Figure 3.3: By of conventional implementation near Z = 0 (−25 ≤ Z ≤ 30) on
the X-Y plane, plots of 5cm step. The ﬁeld is ﬂipped near Z = 0. Additionally,
ϕ -dependency is remarkable at Z = −5, 0, 5.
performed based on the previous results that passed through the branch, if
the for loop satisﬁes the termination condition, the branch prediction fails
and the penalty is received. Therefore, it is considered that speed can beincreased by removing the loop. Other branches should be removed as muchas possible.
Parametrizations of the magnetic ﬁeld using the conventional implementation is stored in the ROOT ﬁle format. From this ﬁle one can write parametrization in plain text format. In order to make it easier to handle, we createdsoftware that parses this text and converts it into portable JSON text format. Furthermore, using Mathematica, parametrizations using the Chebyshev polynomials in the dipole region are extracted from the converted JSONﬁle, and obtained the polynomials consisting of just x, y, z with no loop or
branch by simpliﬁcation for each chamber and components Bx , By , Bz .
3.2.2

Extra multiplication reduction and Fused Multiply Add

Then converted the obtained polynomials into C code and made it compilable
beforehand (Fig.3.4). For the magnetic ﬁeld query, two function is requried.
One is to calculate all of (Bx , By , Bz ) and another calculate only Bz . The
calculation of (Bx , By , Bz ) and Bz was separately deﬁned so as not to include extra jump instructions in the generated code. Figure 3.4 shows one
Bz (x, y, z) function of a chamber describing a magnetic ﬁeld of 5k gauss.
22

f l o a t dip5k1bz ( const f l o a t p [ 3 ] ) {const f l o a t x = p [ 0 ] , y = p [ 1 ] , z = p [ 2 ] ;r e t u r n 5484.90 f +−0.00734 f * y+x *( −0.01580 f +0.00008 f * y+−0.00004 f * z+x *(0.00008 f +1.67927e−7f * x+−1.06067e−6f * y +3.12794e−7f * z ))+z *(40.0776 f+z *(0.10987 f+z *(0.00013 f +6.10621e−8f * z ) ) ) ;};Figure 3.4: Code snippet of a generated Bz function. (shortened to 6 digits
for visibility)
This function takes the Cartesian coordinates of the point to query the magnetic ﬁeld as an argument and returns Bz at that point.
A polynomial in the Figure3.4 was made from the expanded form by theHornerʼs method. Hornerʼs method is like following.
1 + 2x + 3x2 + 4x3 → 1 + x(2 + x(3 + 4x))
As shown above, the variable x is extracted so as to reduce the numberof multiplications. Hornerʼs method has been proved to have the smallest
number of multiplications when used for polynomials[7]. Also, the numberof additions is also minimum[6]. In addition, recent CPUs have dedicatedinstructions for calculating ax+b in Horner format. In the Intel Haswell architecture and after, this Fused Multiply Add (FMA) instructions was extended
to handle ﬂoating point numbers only from integers. O2 computers havenot been purchased yet, so there is no problem that old CPUs mix in online
processing. For non-compliant CPUs used in places other than O2 , they arecompiled to use normal AVX instructions or similar instructions (using the-march compiler option). Figure 3.5 shows the diﬀerence in assembly depending on the presence or absence of FMA. Instructions starting with fma
could calculate a = ax + b with 5 cycles of latency. For mulss and addss
has 5 and3 cycles of latency respectively, it takes 8 cycles in total, 1.6 times
faster with fma.
Since the proposed implementation of the dipole region requires compiling in advance, it can be said that it is an implementation specialized for
calculation speed, discarding versatility.
3.2.3

Constant time complexity O(1) segment search

The computational complexity of the binary search of ﬁgure 2.1 is O(log n)(n is the number of Z segments to be searched + 1 = 90), but about 7 array accesses (90 < 27 ) and three cache miss opportunities (if the number
is single precision ﬂoating point number and the cache line is 64 bytes).
Penalties due to cache miss have been improved year by year, but still there
is a big diﬀerence in the time taken when L1 cache hit and main memory is23

_dip5k1bz:
pushq %rbp
movq %rsp, %rbp
vmovss (%rdi), %xmm0
vmovss 4(%rdi), %xmm1
vmovss 8(%rdi), %xmm2
vmovss 2538974(%rip), %xmm3
vfmadd213ss
2538969(%rip), %xmm1, %xmm3
vmovss 2538965(%rip), %xmm4
vfmadd213ss
2538960(%rip), %xmm1, %xmm4
vfmadd231ss
2538955(%rip), %xmm2, %xmm4
vmovss 2538951(%rip), %xmm5
vfmadd213ss
2538946(%rip), %xmm0, %xmm5
vfmadd231ss
2538941(%rip), %xmm1, %xmm5
vfmadd231ss
2538936(%rip), %xmm2, %xmm5
vfmadd213ss%xmm4, %xmm0, %xmm5
vfmadd213ss%xmm3, %xmm0, %xmm5
vmovss 2538922(%rip), %xmm0
vfmadd213ss
2538917(%rip), %xmm2, %xmm0
vfmadd213ss
2538912(%rip), %xmm2, %xmm0
vfmadd213ss
2538907(%rip), %xmm2, %xmm0
vfmadd213ss%xmm5, %xmm2, %xmm0
popq %rbp
retq
nopw %cs:(%rax,%rax)_dip5k1bz:
pushq %rbp
movq %rsp, %rbp
movss (%rdi), %xmm0
movss 4(%rdi), %xmm3
movss 8(%rdi), %xmm1
movss 3297658(%rip), %xmm2
mulss %xmm3, %xmm2
addss 3297650(%rip), %xmm2
movss 3297646(%rip), %xmm4
mulss %xmm3, %xmm4
addss 3297638(%rip), %xmm4
movss 3297634(%rip), %xmm5
mulss %xmm1, %xmm5
addss %xmm4, %xmm5
movss 3297622(%rip), %xmm4
mulss %xmm0, %xmm4
addss 3297614(%rip), %xmm4
mulss 3297610(%rip), %xmm3
addss %xmm4, %xmm3
movss 3297602(%rip), %xmm4
mulss %xmm1, %xmm4
addss %xmm3, %xmm4
mulss %xmm0, %xmm4
addss %xmm5, %xmm4
mulss %xmm0, %xmm4
addss %xmm2, %xmm4
movss 3297574(%rip), %xmm0
mulss %xmm1, %xmm0
addss 3297566(%rip), %xmm0
mulss %xmm1, %xmm0
addss 3297558(%rip), %xmm0
mulss %xmm1, %xmm0
addss 3297550(%rip), %xmm0
mulss %xmm1, %xmm0
addss %xmm4, %xmm0
popq %rbp
retq

Figure 3.5: Comparison of assemblies w/ and w/o Fused Multiply Add (FMA)
intrinsics. Left: w/o FMA, right: w/ FMA. Instructions which starts withʼfmaʼ computes a = ax + b in 5 cycles of latency. On the other hand, mulss
and addss takes 5, 3 cycles of latency respectively, the former is 1.6 times
faster in simple calculation. Used compiler is Clang-900.0.39.2 on Intel Core
i7 Haswell Refresh.
24referenced, improving cache misses may aﬀect performance[4]. It is considered that the Intel Haswell architecture is about 40 times faster in this case.
Therefore, it is necessary not only to reduce the computation complexity, butalso to make an access such that main memory access is reduced or cache
prefetch succeeds. For linear search, there are few data to be searched and
if it is in the cache, the order of the computation time may be better thanthe binary search despite O(n).
Interpolation search[9] is a search algorithm that may be faster than binary search if the data to be searched does not ﬁt in the cache line. In
Interpolation search, a linear search is performed by determining the indexs to start searching by linear expression, expecting that the values are
distributed linearly in the sorted real number array. This is similar to how
to ﬁnd a speciﬁc person from the phone book. Interpolation search completes in O (log log n) steps on average on uniformly distributed arrays on
average[9]. However, if the data is skewed, it is substantially the sameas the linear search and it becomes the worst calculation complexity O (n).
Fig.3.6 is a plot of the variation in values in order to investigate whether
Interpolation search can be used for searching Z section. Dipole region−1760 ≤ Z ≤ −558.3 is divided into 89 segments without duplication using the Z coordinate maximum value minimum value of the chambers. Letʼscall zid = 0 for the segment with the smallest z and zid = 88 for the segment with the biggest z. The straight line in the ﬁgure tries to associate theZ position with zid by a linear equation, and the gap between this straightline and the segment indicates the possibility that the interpolation searchdoes not operate eﬃciently.
The proposed implementation is a new algorithm that corrects the latter part of this interpolation search so that the search always ends at one
branch. Since there is no loop, both the average computational complexity
and the worst computational complexity are O(1). Such ideas have been
known at least within ALICE for a long time but have not been clariﬁed.
We call this method ”Prepopulated Quick Segment Search”. The fast segment search consists of a slow build part that creates a table from the data
to be searched and a fast search part that uses a table. A schematic diagram of the search part is shown in ﬁgures 3.7 and pseudocode in thealgorithm1. In the ﬁgure, since it is a one-dimensional search method, the Z
axis is taken as an example. The black dot z is the z coordinate of the point
where the magnetic ﬁeld is desired to be acquired. Black line segments represent divisions. The distance of the minimum value zmin and the maximum
value zmax width = zmax −zmin , the number of divisions number of divisions,
The diﬀerence between the smallest element of the sequence and the origin is characterized by of f set = zmin . Indices starting from 0 are given toslices created by this arithmetic progression. The number of divisions is determined beforehand so that the maximum value of only one Z axis value(boundary value of segments) of the end of the chamber is included. You
can follow the procedure later. The ﬁrst step in the search is to identify
this index from z (the ﬁrst line of pseudocode). In the example shown in25

- 1800- 1600- 1400- 1200- 1000

- 800

- 600

Figure 3.6: Split dipole region −1760 ≤ Z ≤ −558.3 into 89 segments by
maximum and minimum Z values of all the chambers. Call a segment withminimum Z as zid = 0, and maximum Z as zid = 88. A crossing line is placed
to show possibility of linear relationship between Z positions and zids. The
gap between the line and the segments implies Interpolation search works
less eﬃciently.

26

Z (cm)

0
０

１

zid=０

2

3

4

5

１

２

z

6

7
３

8
４

9
5

10

11

6

Figure 3.7: Schematic ﬁgure of Quick Segment Search algorithm. The points
are just an example and nothing to do with real data. Prepare evenly spaced
slices (orange arithmetic sequence) which contain only one or zero Z valueof a chamber border (blue point) in advance. The id of one of evenly spaced
slices is immediately calculated from given z (e.g. 6). Every slice knows az position of the blue point and a zid of the segment which is smaller thanthe blue point. The zid of larger side is always (the zid of smaller side + 1).
Then by comparing z and z of the blue point, a zid for queried z is obtained(e.g. 3). Thus, zid is used to lookup the other value (data of X,Y axis or an
index of a polynomial).
Algorithm 1 Prepopulated Quick Segment Search (search part). See text for
details.
index ← ⌊(z − of f set) ∗ (number of divisions/width)⌋.
if index > number of divisions then search f ailed.
slice ← slices[index].
id ← f ormer zid of slice +(if z < segment end position of slice then 0 else 1) .
return id.
the ﬁgure, index = 6. 0 ≤ index < number of divisions, continue, else thesearch is unsuccessful because it is out of range. By the work in the construction part, an array slices is prepared in which slice boundary values
segment end position and the set of zid (f ormer zid) of the division smaller
than the boundary value are stored for all the slices. If there is no boundary
value in the slice, segment end position = +∞ is stored. On the other hand,since the divisions are sorted, the zid of the division larger than the boundary value is always zid + 1. On the fourth line of the pseudo code, we use
this property to obtain the zid of the segment containing z. In the example
shown in the ﬁgure, zid = 3. You can obtain a search table for other axes
by subtracting zid obtained in this way with another sequence. As with the
Z axis, search in the X axis and Y axis direction can be performed, and ﬁnally the number of the chambers can be obtained, and an approximation
function pointer is found.
3.2.4

Quick segment search: Preparation

Describe the construction part of the high-speed segment search algorithm.
Pseudocode is shown in the algorithm2.
27

7

Z

Figure 3.8: Parametrized part of dipole magnet applied region. Coloredcubes are the chambers of dipole ﬁeld parametrization from AliMagF class.
This parametrization contains 1482 chambers.

Algorithm 2 Prepopulated Quick Segment Search (construction part). See
text for details.
n←1while(true){
maxSegments ← 0lastzid ← 0
f or(i ← 0; i < n; + + i){
slice start ← i ∗ (width/n) + of f set
slice end ← (i + 1) ∗ (width/n) + of f set
nSegs ← 0
slices[i] ← (lastzid, +inf )f or(j ← 0; j < segends.length; + + j)
if (slice start < segends[j] < slice end)
nSegs + +slices[i] ← (j, segends[j])lastzid ← j}}maxSegments ← max(maxSegments, nSegs)}
if (maxSegments == 1)return(n, slices)
n++
}

28

µs/callSolenoid region
Dipole region

Conventional
0.43（exact）
0.61（exact）

Proposed
0.053（approx.）0.15（exact.）

Conventional/Proposed
8.2
4.0

Table 4: Consumed CPU time for single query of one random point.
The basic idea is to increase n from 1 in order until the maximum value ofthe number of segment boundaries contained in the slice becomes 1 when
divided into n.
Since this preliminary preparation takes a relatively long time, it is suitable for inquiries about data whose data to be queried is rarely changed.
In this way, the implementation of the dipole region discards versatility andspecializes in the search speed.

4

Results

Improvement of speed by proposed implementation is shown in the table 4.
A speedup of 8.2 times in the solenoid region and 4.0 times in the dipole
region was observed. For the compiler, Clang-900.0.39.2 was used. The
environment is Intel Core i7 Haswell Refresh.

4.1

Result of Solenoid Region

The magnetic ﬁeld approximation precision of the solenoid region is shownin Fig.4.1. Focusing on the center column (Z dependency) of this ﬁgure, theerror exceeds 0.005 kG of the target precision in the region of|Z| > 240 cm.
In this region, the gradient is large and it is diﬃcult to approximate witha three-dimensional polynomial. Due to this inﬂuence, some errors (φ, Zdependence) exceeding target precision are also seen. Fortunately, this Z
range region is outside the acceptance region, so it has no direct eﬀect on
physics. There is a possibility of small inﬂuence on background particles, but
it is not very important that there is high precision when considering them.
Therefore, the proposed implementation seems to have suﬃcient precision
for use in real time reconstruction and simulation. Therefore, we decided toallow errors on |Z| > 250 cm. One can fall back to existing implementations
whenever more accurate results are needed, and you can also examine theeﬀect of errors by comparing the results using existing implementations andthose using suggested implementations.
By approximation based on simple linear regression and simple splitting,precision that is not problematic in practical use could be secured.

29

0

0

− 0.0005− 0.0005− 0.001− 0.001− 0.0015− 0.0015−2

ΔBy, kGauss

0.0015

−1

0

1

2

3

ϕ

ΔBy

0.001
0.00050− 0.0005− 0.001− 0.001
−2

−1

0

1

2

3

0.0015ΔBz

0.001
0.0005− 0.0015

ϕ

ΔBz, kGauss

ΔBz, kGauss−3− 0.0005− 0.001− 0.001
−3

−2

−1

0

1

2

3

ϕ

− 0.0015

ALICE performance

100

200

0

Z

100

200

300

400

500
R, cm

200

300

400

500
R, cm

200

300

400

500
R, cm

By

0.004
0.002
0
− 0.002− 0.004

−200

−100

0

100

200

0

Z

ΔBz

0.001

− 0.0005

ALI−PERF−135528

0ΔBy

0.00050− 0.0015−100

0.001500.002− 0.004
−200

0.001

− 0.0005

Bx

0.004

− 0.0020.00050− 0.0015

ALICE performance00.0015

ΔBy, kGauss

−3

ΔBx

0.001
0.0005ΔBx , kGauss

0.0015

ΔBy, kGauss

ALICE performance

ΔBz, kGauss

ΔBx

0.001
0.0005ΔBx , kGauss

ΔBx , kGauss

0.0015

100

Bz

0.004
0.002
0
− 0.002− 0.004−200

ALI−PERF−135524

−100

0

100

200

Z

0

ALI−PERF−135520

100

Figure 4.1: Positional dependency of mean precision per magnetic ﬁeld component in solenoid region.

4.2

Result of Dipole Region

Since the calculation of the dipole region is implemented so as to be compatible with the conventional implementation including the behavior at chamberinterface, there should be no error.
The proposed implementation using loop and multiplication elimination,
precompilation and fast search algorithm was not as fast as the solenoid region using approximation, but since it was able to secure the same precision
as in the conventional implementation, we used the conventional implementation in all use cases It can safely be replaced.

5

Discussions

5.1

Trade-oﬀ between polynomials calculation speed andprecision

Consider the trade-oﬀ between precision and speed by reducing the numberof polynomial terms. As the general tendency, the average approximation
precision is higher (lower) as the number of terms of the approximate polynomial is larger (smaller), and the time taken for calculation is increased(decreased). Therefore, Recursive Feature Elimination was used to extract
terms contributing to high precision. The ranking of importance of termsby RFE is shown up to the 4th order term tab 5. RFE is performed in thefollowing procedure.

30

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

z
y
x
zz
yy
xx
xy
xz
yz
xyz
xzz
yzz
yyy
xyy
xxy
xxx
xxz

p
pp
ppp
pppp
z
r
zp
zpp
zppp
rp
rpp
rppp
zz
zr
rr
rrp
zzp

18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

yyz
zzz
xyzz
xxyy
xxxz
xxxx
yyyy
xzzz
yzzz
yyyz
xyyy
xyyz
xxyz
xxxy
yyzz
zzzz
xxzz

rrpp
zrp
zrpp
zzr
rrr
zzpp
rrrp
zzz
zrrp
zzzp
zrr
zzrp
rrrr
zzzr
zzrr
zrrr
zzzz

Table 5: Importance rank obtained with Recursive Feature Elimination technique.

1D

2D

3D

4D

Req. prec(TPC 2kG, 0<z<250cm, 90°<φ<180°)

Figure 5.1: Trade-oﬀ between polynomials calculation speed and precision31

1. Create a model starting from all features (terms) and delete featureswith the lowest importance (absolute value of coeﬃcient)
2. Create the model again and repeat deletion of the feature quantity witha small importance until there is no item
In the table, r and p represent the radius and the azimuth angle ϕ in cylindrical coordinates.
The relationship between the number of terms used and the maximum
deviation / average deviation sorted in this order is shown in 5.1.
The departure from the conventional implementation decreases in the
Cartesian coordinate system due to the increase in order. With the precisionof cylindrical coordinate display, it is inferior to Cartesian coordinate display
up to item 23.
Considering the above-mentioned trade-oﬀ between precision and speed,
we made technical decisions to adopt polynomials of at most 3 in the Cartesian coordinate system.

5.2

Comparison of regression or approximation methods
other than multiple linear regression

What other approaches to eﬃciently satisfy the required precision? The
table 6 is the result of using Pythonʼs scikit-learn package which contains
many regression algorithms. The four regression algorithms reached the target precision with respect to the maximum deviation amount. The most systematic was the linear regression and the Bayesian ridge regression. Since
the ﬁtting speed was higher in the linear regression, the Bayesian ridge regression algorithm was not selected. We tried mini (DFP method), conjugate
gradient method, Levenberg-Marquardt method, BFGS method, SVM, decision tree, original algorithm etc. as a method not limited to the scikit-learn
regression algorithm, but suﬃcient precision was obtained I did not use itor because I did not ﬁnish the calculation.

5.3

Memory usage

You should not allocate more memory to calculate the magnetic ﬁeld. The
table 7 is a comparison of memory consumption of the conventional implementation and the proposed implementation. Root [0] new AliMagF () //Finish class loading Similarly, 2492 usage measurement was made using
GetProcInfo () of ROOT, but it should be regarded as a reference value because the blur is large for each measurement. Conventional implementationloads two parameterizations besides solenoid and dipole, so it is reasonable
to compare with half of the values in the table.
In addition, although the dipole region of the proposed implementation
is implemented by dynamic loading of shared libraries, a reference count ismade with global variables, and when all instances are destroyed, unloading32

33

Max deviation
3.13e-03 (62.6%)
3.13e-03 (62.6%)
3.00e-03 (60.0%)
3.34e-03 (66.7%)
1.97e-02 (395%)
3.61e-02 (721%)
1.07e-01 (2140%)
3.03e+00 (60630%)
3.91e+00 (78169%)
8.10e+00 (1.6e5%)
1.20e+02 (2.4e6%)
3.34e+23 (6.7e27%)

Mean deviation
2.57e-04 (5.1%)
2.57e-04 (5.1%)
2.85e-04 (5.7%)
2.69e-04 (5.4%)
4.59e-03 (91.7%)
1.30e-02 (260%)
2.95e-02 (591%)
3.56e-01 (7119%)7.12e-01 (14237%)
1.99e+00 (39804%)
2.02e+01 (4.4e5%)
6.40e+22 (1.2e27%)RMS
3.84e-04 (7.7%)
3.84e-04 (7.7%)
4.11e-04 (8.2%)
4.35e-04 (8.7%)
5.57e-03 (111%)
1.51e-02 (301.8%)
3.55e-02 (709.5%)-

Table 6: Fitting results of regression methods available in scikit-learn (Python) library sorted on coeﬃcient of determination. ʼ-ʼ indicates NODATA. Percentages in the table is ratio against target precision 10−3 × Bz = 5e − 03; less than
100% means the ﬁtting successfully reached target precision.

LinearRegression()
BayesianRidge()
Lasso(alpha = 0.1)
TheilSenRegressor()
ElasticNetCV(cv=5)
OrthogonalMatchingPursuit()
LassoLars(alpha=.1)
Lars()
HuberRegressor()
RidgeCV(alphas=[0.1,1,10])PassiveAggressiveRegressor()
SGDRegressor()coeﬀ. of det. R2
0.99988
0.99988
0.99987
0.999850.97536
0.819100-171.50-795.85-4603-726315-5.7e48

param\impl (KB)
5kGauss
2kGauss

AliMagF（conv.）
38,724
33,104

AliMagFast（prop.）
4,724
6,208

Table 7: Memory usage comparison between the conventional and the proposed implementation. Both 2 kilo Gauss and 5 kilo Gauss cases are shown.
The values are for reference. See text for detail.
is performed and the memory is freed It has become. Therefore, the memory consumption of the suggested implementation should be larger than thesize of the shared library to be dynamically loaded, but the diﬀerence ofthe value returned by GetProcInfo () only increased or decreased by around1000 KB. Therefore, the size of the shared library is listed as a reference
value in the right column of the table 7. Considering the above, looking at
the table, there is a possibility that the memory consumption of the proposed
implementation is less than half that of the conventional implementation.

6

Conclusion

In the heavy ion collision physics, the discovery of the quark gluon plasma
is now an important subject of study for its physical properties. For rarephysical events scheduled to be measured by the ALICE experimental advancement plan to reach unexplored physics, the use of online triggers is
not eﬃcient for the low signal / background ratio. Therefore, it is necessary
to collect all data with lead-lead collision. However, since data of 1 TB/s ormore is generated from the detector, it is diﬃcult to record everything. For
data compression, it is required to calculate the momentum of the particles
from the magnetic ﬁeld inside the detectors, reconstruct their track, andeliminate the background tracks. In the conventional data acquisition program, ~ 5% of the reconstruction time was spent in calculating the magnetic
ﬁeld applied to the particles inside the magnet.
For each space applied by solenoid electromagnet and dipole electromagnet, we proposed a faster implementation than before. In the solenoid
region, approximation by a polynomial of at most 3rd order is performed, andby simplifying the division method of the space, the time spent on the search
can be greatly shortened. The solenoid region is divided into 5×4×22 = 440,
and the precision around Bz ×10−3 is ensured, especially in the center where
precision is required.
In the dipole region, approximation polynomials that record magnetic
ﬁelds recorded using conventional implementation are extracted, modiﬁed
and precompiled so as to minimize the number of multiplications and additions, thereby speeding up without decreasing precision. The search of
the divided space in the dipole region was conventionally binary search andlinear search combined, but by replacing the algorithm with the time com34

plexity O(1) at a price of memory, the speed was increased. In this algorithm,
in order to associate a one-dimensional position with a segment, the interval
between the maximum value and the minimum value is divided into slices.
The number of divisions is chosen so that the number of the boundaries ineach division is at most 1 and the number is minimum. Such divided slices
hold the boundary position of the segment and the id of the section smaller
than the boundary. As a result, the division is deﬁnitely found by dividing
once, lookup of the array, and comparing the position to be inquired with theboundary value. A high-speed segment search is performed for each of theZ axis, the X axis, and the Y axis, so that a polynomial corresponding to thequeried point can be obtained.
For both solenoid and dipole regions, practicality was conﬁrmed by comparing the measurement of the deviation from the conventional implementation and the speed. The proposed implementation will be used both fordata collection and physical simulation of future ALICE experiments.

References[1] P. L. Chebyshev. Théorie des mécanismes connus sous le nom de parallélogrammes. Mémoires des Savants étrangers présentés à lʼAcadémiede Saint-Pétersbourg, vol. 7:539‒586, 1854.

[2] R.Shahoyan (ALICE collaboration). https://github.com/alisw/AliRoot/blob/master/STEER/STEERBa[3] B Abelev et al (The ALICE Collaboration). Upgrade of the ALICE experiment: Letter of intent. Journal of Physics G: Nuclear and Particle
Physics, 41(8):087001, August 2014.
[4] Intel Corporation. Intel® 64 and IA-32 Architectures Optimization Reference Manual, December 2017.
[5] Musa, L (CERN). Conceptual Design Report for the Upgrade of the ALICE ITS. Technical Report CERN-LHCC-2012-005. LHCC-G-159, CERN,
Geneva, Mar 2012.
[6] A. Ostrowski. On two problems in abstract algebra connected withhornerʼs rule. In Studies in Mathematics and Mechanics presented
to Richard von Mises, pages 40‒48. Academic Press, 1954.
[7] V Ya Pan. ”METHODS OF COMPUTING VALUES OF POLYNOMIALS”.
Russ Math Surv, 21(1):105‒136, feb 1966.
[8] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel,
M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas,
A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.
Scikit-learn: Machine learning in Python. Journal of Machine Learning
Research, 12:2825‒2830, 2011.
35

[9] Yehoshua Perl, Alon Itai, and Haim Avni. Interpolation search̶a log
logN search. Communications of the ACM, 21(7):550‒553, jul 1978.
[10] R.Shahoyan. Summary of the l3 magnet ﬁeld analysis. ALICE Internal
Note ALICE-INT-2007-012 v.1, CERN, 2007.
[11] A.Morsch R.Shahoyan. Summary of the alice dipole magnet ﬁeld analysis. ALICE Internal Note ALICE-INT-2008-019 version 1.0, CERN, 2008.
[12] The ALICE Collaboration. Technical Design Report for the Upgrade of
the Online-Oﬄine Computing System. Technical Report CERN-LHCC2015-006, ALICE-TDR-019, CERN, 2015.

A Deﬁnition of the Chebyshev polynomials Tn
The Chebyshev polynomials[1] of the ﬁrst kind are deﬁned by the recurrence
relation:T0 (x) = 1
T1 (x) = x
Tn+1 (x) = 2xTn (x) − Tn−1 (x)Here is a example of the Chebyshev polynomials of the ﬁrst kind.
T0 (x) = 1
T1 (x) = x
T2 (x) = 2x2 − 1
T3 (x) = 4x3 − 3x
T4 (x) = 8x4 − 8x2 + 1
T5 (x) = 16x5 − 20x3 + 5x
T6 (x) = 32x6 − 48x4 + 18x2 − 1
T7 (x) = 64x7 − 112x5 + 56x3 − 7x
T8 (x) = 128x8 − 256x6 + 160x4 − 32x2 + 1
T9 (x) = 256x9 − 576x7 + 432x5 − 120x3 + 9x
T10 (x) = 512x10 − 1280x8 + 1120x6 − 400x4 + 50x2 − 1
T11 (x) = 1024x11 − 2816x9 + 2816x7 − 1232x5 + 220x3 − 11x
The Chebyshev polynomials Tn are polynomials with the largest possible
leading coeﬃcient, but subject to the condition that their absolute value on
the interval [−1,1] is bounded by 1. The interpolation built with the Chebyshev polynomial minimizes the problem of Rungeʼs phenomenon and provides an approximation that is close to the polynomial of best approximation
to a continuous function under n → ∞.
36

B

Indices of precision evaluation

Several indices of precision evaluation are considered where a ﬁeld calculated by the conventional implementation at given point pj is Ai,j (i = x, y, z)
and approximate ﬁeld Bi,j . The deﬁnitions are:• Mean deviation

1 ∑|Bi,j − Ai,j |N j=1
N⟨∆Bi ⟩ =• Root Mean Square (RMS)

v
u
N
u1 ∑
2
RM S(∆Bi ) = t(Bi,j − Ai,j )N j=1• Maximum deviationmax (∆B) =maxj∈Samples|Bi,j − Ai,j |Compared above and target precision Bz × 10−3 (=0.002 or 0.005kG)

C Programs37

Algorithm 3 Mathematica reimplementation of the AliMagF class. It is used
to extract polynomials from existing parametrizations.
(* Cheb: A piece of scalar parametrization *)Cheb[p_ , {r_ , ϕ_ , z_ }] :=Length[ p ]Length p jLength p j,k



j=1

k=1

l=1

- p[[j, k, l]] ChebyshevT[l - 1, z] ChebyshevT[k - 1, ϕ ]ChebyshevT[j - 1, r];(* Inbox: Point (r,φ,z) ∈ b or not *)InBox[b_ ] := And @@ Thread[Between[{r, ϕ, z}, MapThread[List, b ]]];(* Cheb3D: A pair of 3 params (Bx,By,Bz) and its interpolation region *)Cheb3D[block_] :={Cheb[block [["interpolationOutputs", # , "chebyshevPolynomialCoeffs"]],RescalingTransform[block [["interpolationRegion"]] // Transpose, {{-1, 1}, {-1, 1}, {-1, 1}}][{r, ϕ, z}]], InBox[block [["interpolationRegion"]]]} & /@ {1, 2, 3};(* AliMagWrapCheb: Create a 3D-field parametrization *)AliMagWrapCheb[cheb_] := Piecewise /@ (Cheb3D /@ cheb // Transpose);(* AliMagF: Naive reimplementation of AliMagF#Field() *)AliMagF[measurement_ , param_] :=If[param == "dipoleParams", # /. {r → x, ϕ → y} &,TransformedField["Cylindrical" → "Cartesian", # , {r, ϕ, Null} → {x, y, z}] &] /@AliMagWrapCheb[Import[measurement , "RawJSON"][[2, param, All, 2]]];(* Usage *)sol = AliMagF["Sol30_Dip6_Hole.json", "solenoidParams"];sol /. {x → 0, y → 0, z → 0}
dip = AliMagF["Sol30_Dip6_Hole.json", "dipoleParams"];
dip /. {x → 0, y → 0, z → -1000}

38

Algorithm 4 Mathematica C-source generator of the B-ﬁeld parametrizationin HornerForm.
(* Cheb: A piece of scalar parametrization *)Cheb[p_ , {r_ , ϕ_ , z_ }] :=Length[ p ]Length p j

Length p j,k







j=1

k=1

l=1

- p[[j, k, l]] ChebyshevT[l - 1, z] ChebyshevT[k - 1, ϕ ]ChebyshevT[j - 1, r];(* Cheb3D: A pair of 3 params (Bx,By,Bz) and its interpolation region *)Cheb3D[block_] :={Cheb[block [["interpolationOutputs", # , "chebyshevPolynomialCoeffs"]],RescalingTransform[(block [["interpolationRegion"]] // Transpose), {{-1, 1}, {-1, 1}, {-1, 1}}][{x, y, z}]], block [["interpolationRegion"]]} & /@ {1, 2, 3};(* ExpandPower: Replace Power[y,2] to (y*y) *)ExpandPower[e_ ] := StringReplace[ToString[CForm[e //. Power[y_, n_] ⧴ StringJoin["(", Riffle[Table[SymbolName[y], n], "*"], ")"]]], {" " → "", "\"" → ""}]; (* remove unneeded characters *)(* CHornerForm: Encode a polynomial to C-like efficient form *)CHornerForm[e_ ] := ExpandPower[HornerForm[e]];fieldTemplate = StringTemplate["void ``(const float p[3], float b[3]) { const float x = p[0], y = p[1], z = p[2];b[0] = ``; b[1] = ``; b[2] = ``bz(p); };"];
bzTemplate = StringTemplate["float ``bz(const float p[3]) { const float x = p[0], y = p[1], z = p[2]; return ``; };"];(* ChebCodeGen: Write C functions to be used with Field() and GetBz() *)ChebCodeGen[hornered_, prefix_ ] := StringRiffle[MapIndexed[
Module[{funcname = prefix <> ToString[First[#2 ] - 1]},bzTemplate[funcname, #1 [[2, 3]]] <> "\n"<> fieldTemplate[funcname, #1 [[2, 1]], #1 [[2, 2]], funcname]] &,hornered ], "\n\n"];(* PrepareFastDipoleData: Take a jsonified AliMagWrapCheb data and write dipole parametrizations *)PrepareFastDipoleData[jsonFilePath_, exportFilePath_, prefix_ ] :=Module[{data = Import[jsonFilePath , "RawJSON"]},
Module[{dipole = Cheb3D /@ data[[2, "dipoleParams", All, 2]]},Module[{hornered = Map[{(*region*) # [[1, 2]], (*Bx,By,Bz*)CHornerForm /@ # [[All, 1]]} &, dipole]},Export[exportFilePath , ChebCodeGen[hornered, prefix], "String"]]]];
AliRootGit = "alice/AliRoot";
AliPhysicsGit = "alice/AliPhysics";
SourceFile[name_] := FileNameJoin[{$HomeDirectory, AliPhysicsGit, "PWGPP/FieldParam", name}];DestFile[name_] := FileNameJoin[{$HomeDirectory, AliRootGit, "data/maps", name}];
PrepareFastDipoleData[SourceFile["Sol12_Dip6_Hole.json"], DestFile["dip2k.c"], "dip2k"]
PrepareFastDipoleData[SourceFile["Sol30_Dip6_Hole.json"], DestFile["dip5k.c"], "dip5k"]

39