# A novel method for the line-of-response and time-of-flight reconstruction in
TOF-PET detectors based on a library of synchronized model signals

P. Moskal  N. Zoń  T. Bednarski  P. Białas  E. Czerwiński  A. Gajos  D.
Kamińska  Ł. Kapłon  A. Kochanowski  G. Korcyl  J. Kowal  P. Kowalski  T.
Kozik  W. Krzemień  E. Kubicz  Sz. Niedźwiecki  M. Pałka  L. Raczyński  Z.
Rudy  O. Rundel  P. Salabura  N.G. Sharma  M. Silarski  A. Słomski  J. Smyrski
A. Strzelecki  A. Wieczorek  W. Wiślicki  M. Zieliński  Faculty of Physics,
Astronomy and Applied Computer Science, Jagiellonian University, 30-348
Cracow, Poland  Institute of Metallurgy and Materials Science of Polish
Academy of Sciences,
30-059 Cracow, Poland.  Faculty of Chemistry, Jagiellonian University, 30-060
Cracow, Poland  Świerk Computing Centre, National Centre for Nuclear Research,
05-400 Otwock-Świerk, Poland

###### Abstract

A novel method of hit time and hit position reconstruction in scintillator
detectors is described. The method is based on comparison of detector signals
with results stored in a library of synchronized model signals registered for
a set of well-defined positions of scintillation points. The hit position is
reconstructed as the one corresponding to the signal from the library which is
most similar to the measurement signal. The time of the interaction is
determined as a relative time between the measured signal and the most similar
one in the library. A degree of similarity of measured and model signals is
defined as the distance between points representing the measurement- and
model-signal in the multi-dimensional measurement space. Novelty of the method
lies also in the proposed way of synchronization of model signals enabling
direct determination of the difference between time-of-flights (TOF) of
annihilation quanta from the annihilation point to the detectors. The
introduced method was validated using experimental data obtained by means of
the double strip prototype of the J-PET detector and 22Na sodium isotope as a
source of annihilation gamma quanta. The detector was built out from plastic
scintillator strips with dimensions of 5 mm x 19 mm x 300 mm, optically
connected at both sides to photomultipliers, from which signals were sampled
by means of the Serial Data Analyzer. Using the introduced method, the spatial
and TOF resolution of about 1.3 cm ($\sigma$) and 125 ps ($\sigma$) were
established, respectively.

###### keywords:

Scintillator Detectors , Positron Emission Tomography , Time-of-Flight , J-PET

††journal: Journal of LaTeX Templates

##  1 Introduction

Positron Emission Tomography (PET) is at present one of the most
technologically advanced diagnostics methods that allows for non-invasive
imaging of physiological processes occurring in a patient’s body. In the PET
tomography the information about the distribution of annihilation points, and
hence about the density distribution of the administered radiopharmaceuticals
inside the patient’s body, is carried out by pairs of gamma quanta which are
registered in detectors surrounding the patient. All commercial PET devices
use inorganic scintillator materials as radiation detectors \- usually these
are the LBS (BGO) (GE Healthcare), LSO (Siemens) or LYSO (Philips) crystals
[1, 2, 3, 4]. Determination of the interaction point of gamma quanta in PET
detectors is based on the measurement of charge of signals generated by
photomultipliers or avalanche photodiodes (APD) connected optically to
inorganic crystal blocks cut into array of smaller elements. The spatial
resolution achievable with this method is equal approximately to the
dimensions of the small elements of the crystal block. Determination of
interaction points for both annihilation quanta enables reconstruction of the
line-of-response (LOR). In turn, the measurement of the difference between the
arrival times of gamma quanta to the detectors, referred to as time-of-flight
(TOF) difference, allows to calculate position of the annihilation point along
the LOR. The TOF resolution of about of 400 ps achievable with LSO crystals
[5], allows for a substantial improvement of a signal to noise ratio in
reconstruction of PET images [1, 3, 6].

Although detectors used in Positron Emission Tomography are presently at the
highly advanced stage of development there is still a large room for
improvement, and there is ongoing research especially aiming at (i) refinement
of time resolution by search and adaptation of new inorganic crystals [7, 8,
9, 10, 11], (ii) reduction of parallax errors due to the unknown depth of
interaction (DOI) e.g. by application of new geometrical configurations of
crystals and APD and photomultipliers [12, 13, 14, 15, 16], (iii) finding
cost-effective solutions which would allow for construction of large detectors
enabling single-bed whole-body PET imaging as e.g. straw tubes drift chambers
[17, 18] or large area Resistive Plate Chambers (RPC) [19, 20], and (iv)
adaptation of PET detectors for their simultaneous usage together with MRI and
CT modalities [12, 21, 22, 23, 24, 25, 26].

Recently a new concept of large acceptance Jagiellonian PET (J-PET) system
(see Fig. 1.) based on strips of polymer scintillators arranged in a large
acceptance detectors was proposed [27, 28, 29, 30, 31]. The J-PET detector
allows to solve the challenges discussed above in an utterly new way. It
offers improvement of TOF resolution due to the usage of fast plastic
scintillators, it enables a fusion with MRI and CT modalities in a way
allowing for simultaneous morphological and functional imaging [32, 33], it
permits to determine the depth of interaction [34], and constitutes a
promising solution for single-bed whole-body PET imaging. At present it is
however in its early stage of development and requires elaborations of new
hit-position [35] and TOF reconstruction methods which would allow to make use
of the potential it offers. This article is devoted to the presentation of a
reconstruction method that allows to exploit the advantages of the J-PET
detector but it may also be applied to other types of scintillator detectors.

In scintillator detectors, amplitude and shape of signals change strongly with
distance of the hit position to the converter, leading to a deterioration of
the spatial and time resolution. The proposed method of position and time
reconstruction turns this disadvantage into an advantage, and makes use of the
signal shape variation in hit position reconstruction. The method is based on
determination of the degree of similarity between measured signals and
standard signals stored in the data base and on a novel concept of signals
synchronization.

In the following, for the sake of completeness, the J-PET concept is briefly
described. Next, in order to facilitate a clear explanation of the
reconstruction method we introduce a way of representing signals and describe
an example of the creation of the library of model signals. Further on we
describe the invented method of signals’ synchronization, which is crucial for
the reconstruction of LOR and TOF. Finally, the experimental results are
presented in the last section of this article.

##  2 The J-PET detector system

The J-PET test chamber is built out of strips of organic scintillator, forming
a cylinder. One of the possible arrangements of strips is visualized
schematically in Fig. 1. Light signals from each strip are converted to
electrical signals by two photomultipliers placed at opposite ends of the
strip. The position and time of reaction of gamma quanta in the detector
material can be determined based on the time of arrival of light signals to
the ends of the scintillator strips. In article [30] we argued that
disadvantages of polymer scintillators due to the low detection efficiency and
negligible probability for photoelectric effect can be compensated by a large
acceptance, significantly improved time resolution and possibility of usage of
several independent detection layers. Especially promising is the possibility
of extension of the diagnostic chamber in the J-PET detector which does not
entail an increase in the number of photomultipliers and electronics channels
when increasing the axial field of view (AFOV). This feature, in contrast to
crystal-based PET scanners, allows for building single-bed, whole-body PET
scanners without significant increase of costs with respect to scanners with
short AFOV.

Figure 1: Left: Schematic view of the two detection modules of the PET
detector referred to as J-PET [28, 29]. A single detection module consists of
a scintillator strip read out by two photomultipliers labeled with letters PM.
In the first approximation the hit distance from the center of the
scintillator ($\Delta l$) is determined based on time difference measured at
both ends of the scintillator strip, and the position ($\Delta x$) along the
line-of-response is determined from time difference measured between two
modules. Right: An example of the two layers version of the J-PET detection
chamber.

The shape (distribution of number of photons as a function of time) and the
amplitude of the light signal reaching the photomultiplier changes as a
function of the distance between photomultiplier and the place where the light
signal was created. Variations of shapes and amplitudes of light signals
become stronger with the increasing size of the scintillator, and they
constitute a limitation in an achievable time resolution with presently used
electronics readout systems utilizing single threshold constant-fraction or
constant-level discriminators. In the case of J-PET modality with long polymer
scintillator strips, this time resolution determines also the uncertainty of
reconstruction of ionization point. Moreover, distribution of amplitude of
light signals induced by the gamma quanta is continuous due to the fact that
in practice for annihilation quanta only the Compton scattering plays a role
in polymer scintillators and the probability for the photoelectric process is
negligible. As a consequence the amplitude of signals used for the J-PET image
reconstruction varies even if they originate from the same interaction point.
Therefore, a new hit positions reconstruction method is required.

##  3 Signal representation

In the current TOF-PET detectors the reconstruction of line-of-response and of
TOF values is based on the charge and time distributions measured for each
annihilation event without referring to the external sets of model signals. In
this article we present a novel method for the reconstruction of the
interaction point in PET detectors and for the reconstruction of time
differences between the arrival of the annihilation quanta to PET detectors.

The description of the proposed method is based on the example of the J-PET
detector (Fig. 1), consisting of scintillator strips connected optically at
two ends to photomultipliers. Fig. 2B shows schematically signals (voltage as
a function of time) from a single detection module for a few irradiation
positions indicated in Fig. 2A. Figs. 2A and 2B illustrate qualitative changes
of the signal shape and amplitude as a function of the hit position of the
gamma quantum along the scintillator strip.

Figure 2: (A) Pictorial illustration of model signals generation for the J-PET
detector. (B) Example of generated signals and (C) synchronized signals.
Synchronization procedure is explained in section 5. Figure (B) shows
situation when signal from the left photomultiplier was used to define the
trigger. Letters L and R in figure (A) denote left and right photomultiplier,
respectively. The time of the trigger is denoted by $t_{trig}$, and beginnings
of signals for left and right side are indicated by $t_{L}$ and $t_{R}$,
respectively. In figure (C) index “s” is added to indicate times after
synchronization. More details are given in the text.

Signals from photomultipliers are processed by the read-out electronics
enabling determination of their charges and times at which they pass through
given reference voltages [36].

The method described below may in general be used in PET modalities in which
signals are sampled in the voltage domain by means of multi-threshold
constant-level discriminators, or in the domain of the fractions of amplitude
by means of constant-fraction discriminators. Preferably both kinds of
discriminators should be applied for sampling since they deliver complementary
information. Constant level discriminators are used to determine a moment of
time in which the detector signal crosses a defined reference voltage, whereas
the constant fraction discriminators allow to determine the time when signal
crosses voltage level equal to a certain fraction of the signal’s amplitude.

The shape and amplitude of signals corresponding to the registration of gamma
quantum in the scintillator detector changes from event to event and depends
on many factors such as e.g. statistical character of light emission, energy
absorbed in the scintillator, the location of the interaction point in the
scintillator strip and many others. Therefore, in order to reach high
precision of the time and position measurement, signals in the J-PET detector
are split and read by means of the multi-threshold constant level and multi-
threshold constant fraction discriminators.

The registration of each gamma quantum may be represented as a point in a
measurement space $\Omega_{m}$ with the number of dimensions equal to the
number of measured parameters such as times or charges of signals induced in a
single detection module. Each time or charge measurement increases the
$\Omega_{m}$ space by one dimension. Further on we denote number of
measurements done on signals induced by a single gamma quantum by
$N_{m}=2(N_{f}+N_{cr}+N_{cf}+1)$, where $N_{f}$ denotes the number of
thresholds at the constant fraction discriminators, $N_{cr}$ is the number of
time measurements with a constant-level discriminator at the leading edge of
the signal, $N_{cf}$ stands for the number of measurements with a constant
level discriminator at the trailing edge, and $(+1)$ corresponds to the charge
measurement. Factor of „2” before the parenthesis reflects the fact that each
scintillator is read out by two photomultipliers. Hence, a result of the
registration of a single gamma quantum corresponds to a point $P$ in an
$N_{m}$ dimensional measurement space $\Omega_{m}$. First $N_{m}/2$
coordinates of the point correspond to the measurement at one side of the
strip and the next $N_{m}/2$ coordinates to the measurements at the other
side.

Various coordinates of the point and their mutual relations are sensitive in a
different way to the changes of amplitude, time or the shape of the signal.
Therefore, based on the measured signals, it is possible to disentangle
information of the time, position and energy deposited in the scintillator.
For example, values of $P(i)$, for $i=1,...,N_{f}$, corresponding to the
measurement of the time by the constant-fraction discriminators, are sensitive
to variation of signal’s shape, but are not sensitive to the changes of the
signals amplitude provided that the shape of the signal and its time of origin
are not changing. On the other hand values of $P(i)$ corresponding to the
measurement of time with the constant-level discriminators depend on the
signal amplitude even if the shape and time of the origin of the signal are
not changing. Moreover, time differences measured for a given reference
voltage at different sides of the scintillator strip strongly depend on the
place of the gamma quantum interaction. In general, coordinate $P(i)$
representing the time measurement may be expressed as:
$P(i)=t_{measurement}(i)+t_{delay}(i)+t_{trig}$, where $t_{trig}$ acquires the
same value for all coordinates $i$, $t_{measurement}(i)$ denotes the time at
which signal crosses a reference voltage at discriminator corresponding to the
ith dimension in the space $\Omega_{m}$, and $t_{delay}(i)$ stands for the
constant which is subject to calibration, and which denotes the time elapsed
to the moment of the measurement from the moment at which electronic signal
would be created if it was created at the edge of the scintillator without
delays due to the photomultiplier, cables and read-out electronics. The time
offsets $t_{delay}(i)$ may be determined for each detector module with respect
to the reference detector utilizing “beta plus” radioactive isotopes rotating
inside a scanner (see e.g. Fig. 3) or by other methods [37]. Therefore, we
assume that the $t_{delay}(i)$ constants are known for each detection module
and for simplicity, and without loss of generality, we will skip them in the
further considerations.

##  4 Generation of the library of model signals

The reconstruction method described in this article requires generation of a
data base of synchronized model signals for various interaction points. The
library of model signals is generated by scanning the scintillator strip with
a collimated beam of annihilation quanta with profile smaller than the spatial
resolution required for the hit position reconstruction. For example a beam
with the profile width of FWHM equal to 1 mm can be used. Scanning may be
performed using a source of annihilation gamma quanta placed inside a
collimator which may rotate around the axis of the detector and which can be
moved simultaneously along the strips, giving possibility to irradiate each
place of the detector’s inner surface as it is depicted in Fig. 3. Movements
of the collimator must be synchronized with the data acquisition system in
order to assign a place of irradiation to each measured signal. The
information about the position of irradiation is added to each signal in the
library. For each position a high statistics of signals is collected.

Figure 3: Schematic cross sections of the J-PET detector with a collimated
beam of annihilation quanta rotating around the axis and moving along the
detector, as indicated by the arrows. RD denotes a reference detector
permitting to select signals corresponding to the annihilation quanta and
improving the collimation.

As a next step, signals in the generated data base are synchronized in such a
way that for all of them a time corresponding to the reaction of the gamma
quantum amounts to zero, as it is demonstrated in Fig. 2C. The synchronized
library of model signals constitutes a set of baseline points in the space
$\Omega_{m}$, and it is used in the reconstruction procedure for determining
sets of LOR and TOF values from the PET measurements. The synchronization
procedure is described in the next section.

##  5 Method of the model signals synchronization

A raw data base determined as described in the previous section constitutes a
set of points $P$ in a space $\Omega_{m}$. Coordinates of points $P$
correspond to moments of time in which signals pass the discriminator
thresholds with respect to the time of the triggering signal. Few exemplary
signals from the raw data base are shown in Fig. 2B. The signals may be
synchronized with respect to the signal of the reference detector, which
allows for the precise determination of the calibration offsets ($t_{delay}$),
since in this case an average value from the large statistics sample of
signals can be used. However, the reference detector introduces smearing of
time when the single signals from the library are used for the reconstruction
of time of the interaction. Therefore, we introduce a synchronization
technique independent of the performance of the reference detector. The main
idea of synchronization of signals in the data base lies in shifting them in
time such that the moment when gamma quantum hits the detector would be the
same for all events in the library. The absolute value of this time is not
relevant, therefore for simplicity the signals in the library will be
synchronized to the hit-time equal to zero. With the appropriate choice of the
calibration constants, the time at which gamma quantum undergoes scattering in
the detector with respect to the time of the trigger signal is determined by
the mean value: $t_{hit}=(t_{L}+t_{R})/2$, where $t_{L}$, and $t_{R}$ denote
the beginning of the signal measured at the left and right side of the strip,
respectively. The beginning of the signal, represented by a given point P may
be determined based on the coordinates of this point by fitting to them a
function which describes the shape of the signal at its leading edge. A
function describing the shape of signals for a given hit-position may be
determined experimentally e.g. by averaging signals in a high statistical
sample collected at this position. Synchronization is realized by transforming
point $P$ into point $P_{s}$ using the following prescription:
$(P_{s}(i)=P(i)+t_{synch},i=1,...,N_{m})$, where a value of $t_{synch}$, is
chosen such that after the transformation $t_{Ls}+t_{Rs}=0$. Thus,
$t_{synch}=-(t_{L}+t_{R})/2$, and it needs to be determined separately for
each point from the data base. An example of synchronized signals is shown in
Fig. 2C. Note that synchronization of signals in the library to the value of
$t_{hit}=0$ implies that always one of the times from the pair
$(t_{Ls},t_{Rs})$ is negative and the other positive. The above described
synchronization procedure enables to determine not only LOR but also TOF for
each registered event.

##  6 The reconstruction method of LOR and TOF

The reconstruction of the time and position of the interaction of gamma quanta
in the detector is based on comparison of measurement signals for a given
event with synchronized model signals stored in the library. The hit position
is reconstructed as this which corresponds to the signal from the library
which is most similar to the measurement signal, and hit time of the
interaction is reconstructed as a relative time between the measured signal
and the most similar one in the library. A degree of similarity is defined as
a distance between points representing the measurement- and model-signal in
the multi-dimensional measurement space $\Omega_{m}$. The distance is
determined taking into account measurement uncertainties of charge and times
at various reference voltages and correlations between these measurements.

In order to determine the time and position of the gamma quantum interaction
in a given scintillator strip, the algorithm searches through the set of
points in the library of synchronized model signals to find a point $P_{s0}$
which is closest to the point $P$ representing the measurement signal.
Measurements, and as a consequence coordinates of points in the measurement
space $\Omega_{m}$, are burdened with uncertainties which may be correlated
with each other. These uncertainties are described by the covariance matrix
which should be determined for each detection module separately. An inverse
covariance matrix constitutes a metric determining distance in the measurement
space $\Omega_{m}$. Such defined distance, which takes into account
measurement uncertainties and their correlations is in the literature referred
to as Mahalanobis distance [38]. In general a measure of the distance between
points, and thus the measure of the degree of similarity between signals
represented by these points, may be defined in many manners, as for example:
(i) probability that two compared signals are the same (applicable in the case
of the maximum likelihood method), (ii) Chi-square ($\chi^{2}$) value used in
the case of the minimum square method, or (iii) Hausdorff distance used as a
degree of resemblance between two signals [39]. In order to compare a
measurement signal (from the diagnosis of the patient) represented by point
$P$ with the model signal from the synchronized library represented by point
$P_{s}$ one has to perform minimization of the distance between points $P$ and
$P_{s}$ varying the relative time ($t_{rel}$) between the synchronized basis
and the signal $P$ from the diagnosis of the patient. This first step may be
understood as superimposing of signals $P$ and $P_{s}$ on each other. Thus,
the degree of similarity, e.g. a Mahalanobis distance between points $P$ and
$P_{s}$ is expressed as a function of $t_{rel}$:
$Mahanalobis(P+T_{rel},P_{s})$, with $T_{rel}=(t_{rel},...,t_{rel})$, where
$t_{rel}$ is a fit parameter. For each point $P_{s}$ from the synchronized
library a minimum value of $min[Mahalanobis(P+T_{rel},P_{s})]$ is determined
with respect to $t_{rel}$, and next as a point $P_{sfit}$, being closest to
the point $P$, such point $P_{s}$ is chosen, for which a value of
$min[Mahanalobis(P+T_{rel},P_{s})]$ is the smallest. Finally, a point of
interaction of the gamma quantum is determined as a place at which a beam of
annihilation quanta was directed at the moment when a signal represented in
the library by $P_{sfit}$ was registered, and as the time of the interaction
the value of $t_{rel}$ is chosen for which $Mahanalobis(P+T_{rel},P_{sfit})$
is minimal. Such choice of the value of time of interaction constitutes one of
the crucial ideas of the described reconstruction method. It ensures that the
difference between times of interactions reconstructed in different detectors
for the quanta from the same annihilation process, correspond to the true
difference (TOF) between times of arrival of these quanta to the detectors.
This feature is proven below by the reasoning illustrated in Fig. 4.

Figure 4: Pictorial illustration of the method for TOF reconstruction based on
the comparison of the measured signal (solid line) with the synchronized model
signals (dotted lines). For the detailed description see text. $tA_{L}$,
$tA_{R}$ denote the beginning of signals measured at the left and right side
of the detector A, and $tA$ stands for the hit-time calculated as a mean value
of times $tA_{L}$ and $tA_{R}$. Analogously, beginnings of signals from
detector B and signals from the library are shown. In the right-lower corner
of the figure, an event corresponding to the measurement signal is
illustrated. AL, AR, BL, BR denote photomultipliers on the left and right
sides of the A and B detectors, respectively. A cross indicates a center of
the line-of-response, and the dot on the line denotes the point of
annihilation which is by $\Delta x$ away from the center of the LOR.

In order to focus the attention of the reader, without loss of generality, we
assume that the gamma quanta were registered in detectors A and B (see right-
lower corner of Fig. 4). Then, $tA$, a time of the reaction of gamma quantum
in the detector A with respect to the time of the trigger, may be determined
as:

$tA=(tA_{L}+tA_{R})/2$,

and analogously:

$tB=(tB_{L}+tB_{R})/2$,

where $tA_{L}$ denotes the time of the beginning of the signal generated in
the left photomultiplier of detector A measured with respect to the time of
the trigger, and $tA_{R}$, $tB_{L}$ and $tB_{R}$ denote correspondingly
beginning of signals in the right side of detector A and beginning of signals
in left and right sides of detector B. Solid line in Fig. 4 represents pulses
registered in left and right side of detectors A and B for an exemplary event
where annihilation process occurred by $\Delta x$ away from the center of LOR.
The dotted lines indicate signals from the synchronized library of model
signals which were determined by means of reconstruction procedure as most
similar to the measurement signal indicated by solid line. The reconstruction
procedure described above returns $t_{rel}A$ and $t_{rel}B$ as times at which
gamma quanta hit detector A and B, respectively. Thus, Fig. 4 clearly shows
that:

$t_{rel}A=tA+t_{trig}$ and $t_{rel}B=tB+t_{trig}$,

and hence

$t_{rel}B-t_{rel}A=tB-tA=TOF$.

It is important to note that the above result is independent of the time of
the trigger. The result of the above reasoning proves that the synchronization
and reconstruction methods presented in this article allows for the direct
determination of LOR and TOF once the most similar signal to the measurement
signal was found in the library of synchronized model signals.

##  7 Double-strip J-PET prototype

The J-PET detector system shown in Fig. 1 is axially symmetric and its
performance may be tested using a double strip prototype which allows for
simultaneous registration of two annihilation quanta and reconstruction of
both LOR and TOF. Therefore, the functioning of the J-PET detector and
validation of the reconstruction method proposed in this article was verified
using the double strip prototype outlined in Fig. 5. The prototype is built
out from BC-420 scintillator strips [40] with dimensions of 5 mm x 19 mm x 300
mm wrapped with the 3M Vikuiti specular reflector foil [41]. The strips are
read out at both sides by Hamamatsu R4998 and R5320 photomultipliers [42]. Two
different kinds of available photomulitpliers R4998 and R5320 were used.
However, they differ only in quantum efficiency for the registration of
photons in the ultra-violet region not relevant for the emission spectra of
BC-420 scintillator.

Figure 5: A schematic view of the double-strip J-PET prototype. Detailed
description is given in the text.

The source of ${}^{22}Na$ with its active part in the form of cylinder with
diameter of 3 mm and thickness of 1 mm was located within a lead collimator
with a 1.5 mm wide and 20 cm long slit providing a well collimated beam of
annihilation quanta with the spatial profile of about 1.5 mm (FWHM). A
dedicated mechanical system and a step-motor allowed to move the collimator
along ”z” axis with the precision of a fraction of millimeter so as to permit
irradiation of chosen point within the detector. Signals from photomultipliers
were sampled with 100 ps intervals by means of the Serial Data Analyzer
(Lecroy SDA6000A). A library of model signals was created by moving the
collimator in steps of 3 mm and collecting 5000 events for each irradiated
position. The information about the position of irradiation is added to each
signal in the library. Exemplary signals measured for three different
positions are shown in Fig. 6. Left and right panels of Fig. 6 present signals
measured closer to the left and right photomultiplier, respectively. As
expected, signals measured by the photomultiplier nearer to the interaction
site are larger and arrive earlier than signals from the other
photumultiplier, whereas signals collected in the center of the strip are
characterized by the same shape and the same onset time (the signals in the
center may not ideally overlap due to the possible differences in the
photomultiplier transit times and different delays of the SDA channels).

Figure 6: Exemplary signals measured closer to the left photomultiplier at
positions z = 51 mm (left), 150 mm (center) and closer to the right
photomultiplier at position z = 249 mm (right).

Plastic scintillators such as BC-420 consist of carbon and hydrogen and due to
the low atomic number of these elements the probability for the photoelectric
effect for the 511 keV annihilation quanta is negligible. In practice
interactions of annihilation quanta in plastic scintillators occur only via
Compton scattering [43, 44], and the spectrum of energy deposition and hence
distribution of charge of registered signals is continuous and ranges from
zero to 0.341 MeV (2/3 of electron mass). The example charge spectrum of
signals registered by irradiating the middle of the scintillator is shown as a
black solid line in the left panel of Fig. 7. In order to avoid large
fluctuations in shape of signals consisting of small numbers of photoelectrons
for the further analysis we have selected only these events for which energy
depositions were larger than 0.2 MeV in both scintillator strips. In order to
find relation between the measured charge and deposited energy the Klein-
Nishina formula [45] convoluted with the detector resolution was fitted to the
experimental data with energy calibration constant and normalisation as free
parameters [44]. An example of result of such fit is shown as dashed red
histogram in Fig. 7. It is worth to stress that in a reconstruction of the
tomographic image such filtering of signals will be performed and only signals
with energy deposited larger than 0.2 MeV will be considered in order to
suppress events originating from scattering of the annihilation quanta in the
patient’s body [30].

Figure 7: Left panel: (Solid black histogram) Charge distribution of signals
from one of the photomultipliers determined by irradiating the center of the
scintillator strip with collimated beam of annihilation gamma quanta. (Dashed
red histogram) Theoretical distribution of energy of scattered electrons [45]
convoluted with the resolution of the detector [44] and fitted to the
experimental data with normalization and energy calibration constants as the
free parameters. Details of the fitting procedure are described in reference
[44]. The charge of the signal is expressed in the number of photoelectrons
estimated using a method described in reference [46]. The lower range of the
experimental spectrum is cut by the threshold set at SDA. Right panel: An
example of the distribution of the noise of the measured signals. The figure
shows distribution of voltage of a single signal for times lower than -2.5 ns,
i.e. before the onset of real pulse (examine Fig. 6). A superimposed line
indicates a result of the fit of the Gaussian distribution.

###  7.1 LOR and TOF reconstruction

According to the description included in section 5 we have synchronized model
signals in the library such that the time of the interaction of gamma quantum
corresponding to each signal is equal to zero. To this end each model signal
was shifted in time by the value of $t_{synch}=-(t_{L}+t_{R})/2$. For the sake
of simplicity, we have determined $t_{L}$ and $t_{R}$ as times at which the
signal cross the threshold voltage of 80 mV. A value of 80 mV was chosen to
optimise smearing of time due to the noise and due to the time walk effect. In
order to decrease the influence of the time walk effect the threshold should
be as low as possible but on the other hand it should be sufficiently high to
decrease the influence of the electronic noise which typically amounts to
about 10 mV to 20 mV(sigma) depending on the applied voltage. An example of
the electronic noise spectrum is presented in Fig. 7. It is also important to
stress that before the determination of $t_{L}$ and $t_{R}$ each signal was
corrected for the pedestal which in the example shown in Fig. 7 amounts to 9.9
mV.

For the test of the reconstruction method introduced in section 6 we have
chosen events measured when the collimated beam was irradiating strips at the
following positions: z = 51 mm, z = 99 mm, z = 150 mm, z = 201 mm and z = 249
cm. The position and TOF resultions of the J-PET prototype are determined from
distributions of the differences between true and reconstructed values of
position and TOF, respectively. As $true$ positions we denote real positions
of irradiation.

For the purpose of this demonstrative analysis each signal is represented as
an array with coordinates P(i), where i = 1,…,42. According to the description
from section 3 first 21 coordinates describe results of measurement with left
signal and next 21 corresponds to the right signal. In particular:

P(1),…,P(10) corresponds to times at which a left signal is crossing ith
threshold voltage $V_{i}$ = 60 mV + (i-1) $\times$ 50 mV,

P(11),…,P(20) corresponds to times at which a left signal is crossing ith
fraction of its amplitue $f_{i}$ = 0.1 + (i-11) $\times$ 0.05,

P(21) corresponds to the signal charge. And analogously P(22),…,P42 are
defined for the right signal.

As introduced in section 6, $P$ and $P_{s}$ denote signals from the tested
subset (P) and from the library of synchronized model signals ($P_{s}$),
respectively. In order to reconstruct a place of gamma interaction
corresponding to a given signal P, this signal is compared with all signals in
the library. Next, position assigned to the most similar model signal is taken
as the reconstructed position. As a measure of similarity a $\chi^{2}$ like
variable is used which is defined as follows:

$\begin{array}[]{c}\chi^{2}(P,P_{s},t_{rel})=\\\
\sum_{i=1}^{20}(P(i)-P_{s}(i)-t_{rel})^{2}/\sigma^{2}(t)\\\
+~{}\sum_{i=22}^{41}(P(i)-P_{s}(i)-t_{rel})^{2}/\sigma^{2}(t)\\\
+~{}(P(21)-P_{s}(21))^{2}/\sigma^{2}(Q)\\\
+~{}(P(42)-P_{s}(42))^{2}/\sigma^{2}(Q),\end{array}$ (1)

where $\sigma(t)$ varies between $\sim 13$ ps and $\sim 40$ ps depending on
the thereshold, and it was determined by the measurement of distributions of
time differences of the same signal split into two different SDA channels. As
regards the charge, the studies described in reference [44] revealed that the
uncertainty of the measurement of a signal’s charge is dominated by the
statistical fluctuation of the number of photoelectrons. Therefore, we express
charge in units of photolectrons $N_{phe}$ and estimate its uncertainty as
$\sqrt{N_{phe}}$. Thus, the $\sigma^{2}(Q)$ denoting the variance of the
difference of the measured charges is equal to the sum of the number of
photoelectrons from the compared signals. Fig. 8 presents an example plot of
minimum values of $\chi^{2}$ determined during the reconstruction process for
one of the P signals measured at position z = 150 mm. Each point at this plot
corresponds to a minimum value of $\chi^{2}$ resulting from the comparison of
signal $P$ with a model signal $P_{s}$. A minimum value of $\chi^{2}$ is found
with respect to $t_{rel}$.

Figure 8:  $\chi^{2}_{min}$ as a function of position assigned to the model
signals. Superimpsed solid line indicates average value of $\chi^{2}_{min}$
determined for each position separately.

Figure 9:  (Left) Distribution of differences between the true and
reconstructed position for signals measured at z = 150 mm. (Right) Position
resolution as a function of the place of gamma quantum interaction.

It should be noted that for each position many signals in the library are
stored, and hence many points are visible in Fig. 8 at each position. Finally,
as a most similar signal to the processed $P$ signal, such model signal
$P_{sfit}$ is chosen for which $\chi^{2}_{min}$ is the smallest. Position
assigned to $P_{sfit}$ is taken as reconstructed position corresponding to
signal $P$. Left panel of Fig. 9 shows distribution of differences between the
true and reconstructed position for signals measured at z = 150 mm, where by
fitting a Gaussian function the resolution of $\sigma$ $\approx$ 13 mm was
established. Right panel indicates that this resolution does not change with
the position.

Figure 10:  (Left) Distribution of differences between the true and
reconstructed TOF for signals measured at z = 150 mm. The non-zero value of
the mean is mainly due to the difference in the cable lengths and it is matter
of calibration. (Right) TOF resolution as a function of the place of gamma
quantum interaction.

As the last step of the analysis a TOF value is reconstructed as $t_{rel}A$ \-
$t_{rel}B$ according to the procedure described in section 6 and illustrated
in Fig. 4, where $t_{rel}A$ and $t_{rel}B$ stand for the $t_{rel}$ values for
which the global minima of $\chi^{2}$ were found for first and second detector
strip, respectively. As regards the true TOF value, it is equal to zero since
the source was positioned in the middle between the detection modules.
However, due to the differences in delays caused by different electronics
channels and cables, the reconstructed TOF may differ from zero but it should
be the same for each event. The result of the TOF distribution reconstructed
for signals measured at z = 150 mm is shown in the left panel of Fig. 10, and
the right panel shows the TOF resolution as a function of the position. The
determined resolution is equal to about 125 ps($\sigma$) over the full length
of 30 cm long detector. It is important to stress that the TOF resolution
includes contribution from the spread of the time of interaction due to the
unknown depth of interaction of the gamma quantum within 19 mm thick
scintillators (see Fig. 5) and due to the size of the source (3 mm in
diameter). These two effects together cause a spread of about 27 ps($\sigma$).

##  8 Summary

A method enabling reconstruction of hit time and hit position of gamma quanta
in scintillator detectors was described and validated based on the
experimental data collected with the double-module prototype of the J-PET
detector. The method is based on a comparison of measured signal probed in the
voltage or time domains with synchronized model signals from the library. The
hit time and hit position are reconstructed as these which correspond to the
signal from the library which is most similar to the measured signal. A
measure of similarity is defined as the distance between points representing
the measured- and model-signal in the multi-dimensional measurement space. In
order to compare a measured signal (from the diagnosis of the patient) with
the model signal from the synchronized library, minimization of the distance
between points representing these signals is performed as a function of
relative time between them. The relative time resultant from the minimization
is taken as the time at which gamma quantum interacted in the detector. Such
choice of the value of time of interaction constitutes one of the crucial
ideas of the described reconstruction method, and ensures that the difference
between times of interactions reconstructed in different detectors for the
quanta from the same annihilation process, correspond to the true difference
(TOF) between times of arrival of these quanta to the detectors. The novelty
of the method lies also in synchronization of the model signals in a way
enabling determination of TOF and in a manner of determining the time of the
interaction of gamma quantum in detectors. In the article an exemplary
procedure for generating a library of model signals was also presented which
is based on scanning the scintillator strip with a collimated beam of
annihilation quanta with profile smaller than the spatial resolution required
for the hit-position reconstruction.

The introduced method was validated by means of the experimental data
collected by the double strip prototype of the J-PET detector built out from
plastic scintillator strips with dimensions of 5 mm x 19 mm x 300 mm read out
at both sides by photomultipliers. The strips were irradiated by the
annihilation quanta from the ${}^{22}Na$ source placed in the middle of a lead
collimator. A library of model signals was created using a dedicated electro-
mechanical system permitting to move the collimator along the scintillators in
a way synchronized with the data acquisition system. Signals from
photomultipliers were sampled with 100 ps intervals by means of the Serial
Data Analyzer. Applying the method introduced in this article a spatial
resolution of about 1.3 cm ($\sigma$) for the hit-position reconstruction and
TOF resolution of about 125 ps ($\sigma$) were established.

The obtained result for the TOF resolution for the detector of 30 cm length is
better by about a factor of two with respect to the current TOF-PET tomographs
characterized by typical field of views of about 16 cm and TOF resolution of
about 230 ps ($\sigma$) [47].

The result presented in this article can still be improved in the future by
more elaborated method of the determination of the onset of the signals used
for the synchronization of the library of model signals (e.g. by utilising
more than one value of time determined at different threshold levels) and by
application of the measure of similarities which would account for the
possible correlations between the times measured at different thresholds.

##  9 Acknowledgements

We acknowledge technical and administrative support by T. Gucwa-Ryś, A.
Heczko, M. Kajetanowicz, G. Konopka-Cupiał, W. Migdał, K. Wołek and the
financial support by the Polish National Center for Development and Research
through grant No. INNOTECH-K1/IN1/64/159174/NCBR/12, the Foundation for Polish
Science through MPD programme, the EU, MSHE Grant No. POIG.02.03.00-161
00-013/09, and Doctus - the Lesser Poland PhD Scholarship Fund.

## References

  * [1] M. Conti, Phys. Med. 25 (2009) 1-11.
  * [2] J. L. Humm, A. Rosenfeld, A. Del Guerra, Eur. J. Nucl. Med. Mol. Imaging 30 (2003) 1574.
  * [3] J. S. Karp et al., J. Nucl. Med. 49 (2008) 462.
  * [4] D. W. Townsend, Ann. Acad. Med. Singapore 22 (2004) 133.
  * [5] W. W. Moses, S.E. Derenzo, IEEE Trans. Nucl. Sci. 46 (1999) 474.
  * [6] W. W. Moses, IEEE Trans. Nucl. Sci. 50 (2003) 1325.
  * [7] M. Conti, L. Eriksson, H. Rothfuss, C. L. Melcher, IEEE Trans. Nucl. Sci. 56 (2009) 926.
  * [8] M. Moszynski, T. Szczesniak, Acta Phys. Pol. B. Proc. Supp. 4 (2011) 59.
  * [9] R. D. Schaart et al., Phys. Med. Biol. 55 (2010) N179.
  * [10] R. D. Schaart et al., Phys. Med. Biol. 54 (2009) 3501.
  * [11] A. Kuhn et al., IEEE Trans. Nucl. Sci. 53 (2006) 1090.
  * [12] D. B. Crosetto, Int. J. Med. Tech. 1 (2003) 1.
  * [13] R. S. Miyaoka et al., IEEE Trans. Nucl. Sci. 45 (1998) 1069.
  * [14] W. W. Moses, S. E. Derenzo, IEEE Trans. Nucl. Sci. 41 (1994) 1441.
  * [15] A. Saoudi et al., IEEE Trans. Nucl. Sci. 46 (1999) 462.
  * [16] J. J. Vaquero et al., Conf. Rec. IEEE Nucl. Sci. Symp. Med. Imaging Conf. (1998) M6-29.
  * [17] J. L. Lacy, Ch. S. Martin, L. P. Armendarez, Nucl. Instrum. and Meth. A 471 (2001) 88.
  * [18] N. N. Shehad et al., Nuclear Science Symposium Conference Record IEEE 4 (2005) 2895.
  * [19] G. Belli et al., J. Phys. Conf. Ser. 41 (2006) 555.
  * [20] A. Blanco et al., Nucl. Instrum. and Meth. A 602 (2009) 780.
  * [21] K. M. Gilbert, W. B. Handler, T. J. Scholl, J. W. Odegaard, B. A. Chronik, Phys. Med. Biol. 51 (2006) 2825.
  * [22] M. S. Judenhofer et al., Nat. Med. 14 (2008) 459.
  * [23] P. K. Marsden et al., Brit. J. Radiology 75 (2002) S53.
  * [24] B. Pichler, H. F. Wehrl, A. Kolb, M. S. Judenhofer, Semin. Nucl. Med. 38 (2008) 199.
  * [25] H. Quick et al., MAGNETOM Flash 1/2011 (2011) 88.
  * [26] D. W. Townsend, Phys. Med. Biol. 53 (2008) R1–R39.
  * [27] P. Moskal, Patent (2014) No. EP2454612B1, WO2011008119.
  * [28] P. Moskal et al., Bio-Algorithms and Med-Systems 7 (2011) 73. [arXiv:1305.5187 [physics.med-ph]].
  * [29] P. Moskal et al., Nuclear Medicine Review 15 (2012) C68. [arXiv:1305.5562 [physics.ins-det]].
  * [30] P. Moskal et al., Nuclear Medicine Review 15 (2012) C81. [arXiv:1305.5559 [physics.ins-det]].
  * [31] P. Moskal et al., Radiotheraphy and Oncology 110 (2014) S69.
  * [32] P. Moskal, Patent Application No. PCT/EP2014/068373 (2014).
  * [33] P. Moskal, Patent Application No. PCT/EP2014/068363 (2014).
  * [34] J. Smyrski et al., Bio-Algorithms and Med-Systems 10 (2014) 59. [arXiv:1311.6129 [physics.ins-det]].
  * [35] L. Raczyński et al., Nucl. Instrum. and Meth. A764 (2014) 186-192. [arXiv:1407.8293 [physics.ins-det]].
  * [36] M. Pałka et al., Bio-Algorithms and Med-Systems 10 (2014) 41. [arXiv:1311.6127 [physics.ins-det]].
  * [37] M. Silarski et al., Bio-Algorithms and Med-Systems 10 (2014) 19. [arXiv:1311.6152 [physics.ins-det]].
  * [38] P. C. Mahalanobis, Proceedings of the National Institute of Sciences of India 2 (1936) 49–55
  * [39] D. P. Huttenlocher et al., IEEE Trans. Pat. Ana. Mach. Int. 15 (1993) 850-863
  * [40] Saint Gobain Crystals, http://www.crystals.saint-gobain.com
  * [41] 3M Optical Systems, www.3M.com/Vikuiti
  * [42] Hamamatsu, http://www.hamamatsu.com
  * [43] K. Szymański et al., Bio-Algorithms and Med-Systems 10 (2014) 71-77. [arXiv:1312.0250 [physics.ins-det]].
  * [44] P. Moskal et al., Nucl. Instrum. and Meth. A764 (2014) 317-321. [arXiv:1407.7395 [physics.ins-det]].
  * [45] O. Klein, Y. Nishima, Z. Phys. 52 (1929) 853
  * [46] T. Bednarski et al., Bio-Algorithms and Med-Systems 10 (2014) 13. [arXiv:1312.2744 [physics.ins-det]].
  * [47] V. Bettinardi et al., Medical Physics 38 (2011) 5394.


∎

emailEAe-mail<EMAIL_ADDRESS>emailVBe-mail<EMAIL_ADDRESS>emailSFe-mail<EMAIL_ADDRESS>emailKGe-mail:
<EMAIL_ADDRESS>emailHMe-mail<EMAIL_ADDRESS>emailDSe-
mail<EMAIL_ADDRESS>emailHSe-mail<EMAIL_ADDRESS>emailPSe-mail:
<EMAIL_ADDRESS>emailKTe-mail<EMAIL_ADDRESS>
11institutetext:  Department of Physics, Brookhaven National Laboratory,
Upton, New York 11973 22institutetext:  National Centre for Nuclear Research
(NCBJ), Pasteura 7, 02-093 Warsaw, Poland 33institutetext:  University of
Calabria & INFN-Cosenza, Italy 44institutetext:  University of Glasgow,
Glasgow G12 8QQ, United Kingdom 55institutetext:  IRFU, CEA, Université Paris-
Saclay, F-91191 Gif-sur-Yvette, France 66institutetext:  PRISMA+ Cluster of
Excellence, Institut für Physik, Johannes Gutenberg-Universität, D-55099
Mainz, Germany

# EpIC: novel Monte Carlo generator for exclusive processes

E. C. AschenaueremailEA,addressBNL  V. BatozskayaemailVB,addressNCBJ  S.
FazioemailSF,addressCAL  K. GatesemailKG,addressGlasgow  H.
MoutardeemailHM,addressCEA  D. SokhanemailDS,addressCEA,addressGlasgow  H.
SpiesbergeremailHS,addressMainz  P. SznajderemailPS,addressNCBJ  K.
TezginemailKT,addressBNL

(Received: date / Accepted: date)

###### Abstract

We present the EpIC Monte Carlo event generator for exclusive processes
sensitive to generalised parton distributions. EpIC utilises the PARTONS
framework, which provides a flexible software architecture and a variety of
modelling options for the partonic description of the nucleon. The generator
offers a comprehensive set of features, including multi-channel capabilities
and radiative corrections. It may be used both in analyses of experimental
data, as well as in impact studies, especially for future electron-ion
colliders.

††journal: Eur. Phys. J. C

##  1 Introduction

The objective of this study is to develop a Monte Carlo event generator for
exclusive processes involving hadrons remaining coherent during interactions
with high-energy leptons. The factorisation theorems developed in the
framework of quantum chromodynamics (QCD) allow one to describe such processes
in terms of non-perturbative generalised parton distribution functions (GPDs)
Muller:1994ses ; Ji:1996ek ; Ji:1996nm ; Radyushkin:1996ru ; Radyushkin:1997ki
convoluted with perturbatively calculable hard scattering subprocesses.

GPDs are universal, process-independent functions that parametrise the off-
forward nucleon matrix elements of quark and gluon bilinear operators with
light-like separations. In case there is no momentum transfer to the nucleon,
i.e. in the forward limit, certain GPDs become equivalent to PDFs.
Additionally, the first Mellin moments of GPDs are related to elastic form
factors. In this regard, GPDs may be viewed as a unified concept of elastic
form factors studied via elastic scattering processes and one-dimensional
parton distribution functions studied via (semi-) inclusive scattering
processes. Another key aspect of GPDs is their relation to nucleon tomography.
The Fourier transform of GPDs are related to the impact parameter space
distributions when there is no collinear, but finite transverse momentum
transfer to the nucleon Burkardt:2000za ; Burkardt:2002hr ; Burkardt:2004bv .
This relation enables nucleon tomography, which involves the correlation
between impact space distribution functions, parton polarization, and the
longitudinal momentum fractions carried by partons. GPDs also exhibit a unique
relationship with energy-momentum tensor (EMT) form factors, which encode
fundamental properties of the nucleon, such as mass and spin decomposition of
the nucleon into its constituent parts (including orbital angular momentum)
Ji:1996ek ; Ji:1996nm ; Ji:1997gm  as well as its internal structure based on
the so-called “mechanical” forces Polyakov:2002yz ; Polyakov:2018zvc ;
Lorce:2018egm . For more information on GPDs, we refer to the available
reviews on the subject, like Ref. Diehl:2003ny .

The extraction of GPDs from data on exclusive processes is not an easy task,
mainly due to the difficulty of the inverse problem one must solve. Namely,
typically many types of GPDs contribute to a given process, and it is
necessary to deconvolute them all from the process amplitude. In order to
accomplish this, measurements of many different types of processes are needed
in a wide range of kinematic domains. Therefore, measurements of exclusive
processes have been performed in several facilities, such as DESY, JLAB, and
CERN. In addition to that, GPDs are pillars of scientific programmes for a new
generation of machines. This includes both electron-ion colliders, such as EIC
Accardi:2012qut ; AbdulKhalek:2021gbh , EIcC Anderle:2021wcy  and LHeC
LHeCStudyGroup:2012zhm , and fixed target experiments, such as AMBER at CERN
Adams:2676885  and JLAB12 Arrington:2021alx .

As part of this global effort, we have developed a Monte Carlo (MC) generator
called EpIC, whose logo appears in Fig. 1. The EpIC generator features a novel
architecture based on the modularity programming paradigm. In this way, the
code structure is kept as simple as possible and the addition of new
developments, such as new channels or algorithms to generate random numbers,
is made as easy as possible. The architecture of EpIC, including the
nomenclature for naming its elements, is based on the PARTONS framework
Berthou:2015oaw . PARTONS is also used to evaluate the Born cross-section for
a given process, which is used to generate MC events after the inclusion of
radiative corrections (RCs).

Figure 1: Logo of the project.

EpIC can generate events for a number of exclusive processes. The following
are available at the time of this writing: deeply virtual Compton scattering
(DVCS), time-like Compton scattering (TCS), and deeply virtual meson
production (DVMP) of $\pi^{0}$ and $\pi^{+}$ mesons. In this article, we
choose the case of DVCS whenever an exemplary process is needed. Adapting to
other processes is straightforward and typically only requires replacing the
acronym, e.g. the equivalence of DVCSGeneratorService for TCS is
TCSGeneratorService. EpIC, as in the case of PARTONS, follows the so-called
Trento convention Bacchetta:2004jz  for the definition of e.g. scattering
angles.

The purpose of this article is to provide a reference to be used for a better
understanding of EpIC’s structure and its functioning. We will mostly focus on
the explanation of the architecture, starting with a brief description of the
PARTONS framework in Sect. 2. Topics covered in that section are crucial to
understand the technical aspects of EpIC. In Sect. 3, we present the
architecture of EpIC. In particular, we list all types of modules and show the
relation between them. The user interface is introduced in Sect. 4, while the
demonstration of the generator’s performance is shown in Sect. 5. The summary
is given in Sect. 6. In the appendix of the article (A), we review the basics
of radiative corrections and explain how to compute them in the collinear
approximation. The purpose of this appendix is to help users familiarise
themselves with the parameters of RC used in EpIC.

EpIC is written in the C++ programming language. The code can be accessed at
the GitHub platform www:epic_github  and is distributed under the GPL 3.0
licence. The project’s webpage, Ref. www:epic , includes information such as a
technical documentation of C++ classes, a guide on how to compile the code,
and example code for users.

##  2 PARTONS framework

PARTONS (PARtonic Tomography Of Nucleon Software) Berthou:2015oaw  is a
software framework for studying the 3D structure of hadrons. It provides
essential tools for QCD phenomenology and allows one to compute observables
from the models of non-perturbative objects. For instance, it includes various
models of GPDs and methods to perform their pQCD evolution. It also includes
methods to evaluate amplitudes and cross-sections for a variety of exclusive
processes, like DVCS and TCS, described at various orders of pQCD precision.
The framework can be used in analyses to predict observables from existing
models, see e.g. Refs. Grocholski:2019pqj ; Anderle:2021wcy ;
AbdulKhalek:2021gbh , but also to constrain new models from available
experimental data Moutarde:2018kwr ; Moutarde:2019tqa . Because of its
versatility, the framework can also act as a laboratory for studying new
concepts for phenomenology, like new ways of modelling, see e.g. Refs.
Bertone:2021yyz ; Chavez:2021llq ; Dutrieux:2021wll .

What distinguishes PARTONS from other software projects used in the field of
particle physics is its architecture utilising the modular programming
paradigm. The central role in this architecture, which has been also adopted
in EpIC, is played by _the modules_. A module is a single encapsulated
development of a given type, for instance a single GPD model. The only purpose
of modules is to provide data for a specific input, like values of GPDs for a
given GPD kinematics. The construction of modules utilises both the
inheritance and polymorphism mechanisms of C++. Therefore, the development of
new modules is remarkably straightforward, as developers need only to provide
their code related to the developments they are interested in and adapt it to
the predefined classes in the framework. We demonstrate this with the example
of GPDGK11 module, which is used to implement the Goloskokov-Kroll (GK) GPD
model Goloskokov:2005sd ; Goloskokov:2007nt ; Goloskokov:2009ia , from the
PARTONS library. The following is an excerpt from the header file: namespace
PARTONS

class GPDGK11: public GPDModule

public:

//Unique ID to automatically register the class in the registry. static const
unsigned int classId;

//Constructor. //Class name is passed to GPDModule constructor. GPDGK11(const
std::string& className);

//Used for automation. virtual void configure(const ElemUtils::Parameters
parameters);

//Function containing the evaluation of GPD $H$. virtual PartonDistribution
computeH();

… ; where only the most relevant members of the class are shown, some of which
will be discussed further in this section. The purpose of the GPDGK11 module
is to evaluate the values of GPDs. This task is predefined in GPDModule, being
in terms of inheritance the parent class of GPDGK11. The predefinition means
here the declaration of a virtual function, so that each derived classes
(GPDGK11, GPDGK16, …) have a function of the same signature among its members,
but the actual implementation of this function in those classes is in general
different. The evaluation of GPDs takes place, for instance, in the following
function: using namespace PARTONS;

PartonDistribution GPDGK11::computeH()

//Protected variables that store input GPD kinematics to be used by
developers. m_x; m_xi; m_t; m_MuF2, m_MuR2;

…

//Container to be returned. PartonDistribution partonDistribution;

//Store values of $H$, $H^{(+)}$ and $H^{(-)}$, respectively, for up quarks.
//Results for other flavours can be added analogously.
partonDistribution.addQuarkDistribution( QuarkDistribution( QuarkFlavor::UP,
…, …, …));

//Store value of $H$ for gluons. partonDistribution.setGluonDistribution(…);

return partonDistribution; The function GPDGK11::computeH() returns a
container that stores the numerical values of GPDs of type $H$, including
singlet, $H^{(+)}$, and non-singlet, $H^{(-)}$, combinations, as defined in
Ref. Diehl:2003ny . The input GPD kinematics is accessible in the body of the
function via the m_x, m_xi, m_t, m_MuF2 and m_MuR2 variables, which are
defined in the GPDModule class. Those variables correspond to $x$ (average
longitudinal momentum fraction carried by the active parton), $\xi$ (skewness
parameter), $t$ (four-momentum transfer to the hadron target),
$\mu_{\mathrm{F}}^{2}$ (factorisation scale squared) and
$\mu_{\mathrm{R}}^{2}$ (renormalisation scale squared) GPDs depend on. They
are set when one evaluates GPDs via PartonDistribution compute(const
GPDKinematic& kinematic, GPDType::Type gpdType), which is also defined in
GPDModule. As one may see, in our example quark flavours and GPD types are
distinguished by QuarkFlavor and GPDType enums, respectively, while inputs and
outputs are encapsulated in the QuarkDistribution, PartonDistribution and
GPDKinematic containers.

A single instance of each module is loaded into memory at the beginning of the
program execution, namely, when static constant variables are initialised. The
addresses of such instances are stored by _the registry_ , which later acts as
a phone book, allowing one to access the instances by either their IDs or
defined names. We demonstrate this mechanism with the following example: using
namespace PARTONS;

const unsigned int GPDGK11::classId = BaseObjectRegistry::getInstance()-¿
registerBaseObject(new GPDGK11(”GPDGK11”)); When the program starts,
GPDGK11::classId is initialised with a unique ID assigned by the registry
(BaseObjectRegistry class being a singleton). During this process, the
registry takes and stores the address of each new instance. All modules
inherit from the same basic class called BaseObject, allowing the registry to
store addresses of the same type of objects. The constructor of BaseObject
requires a unique string of characters that is used as a human-readable name
of the module.

The registry stores addresses of single instances of pre-configured modules.
Users may get copies of those instances by using _the factory_ (also being a
singleton), for instance: using namespace PARTONS;

GPDModule* pGPDModel = nullptr;

//Get by ID. pGPDModel = Partons::getInstance()-¿ getModuleObjectFactory()-¿
newGPDModule(GPDGK11::classId);

//Get by name pGPDModel = Partons::getInstance()-¿ getModuleObjectFactory()-¿
newGPDModule(”GPDGK11”); The registry and factory are the key architectural
ingredients, which allow users to add an unlimited number of new modules in
the PARTONS architecture, without having to modify existing modules.

There are many types of modules, but one may _configure_ all of them, i.e.
pass additional information to the modules, in a generic way, avoiding, for
instance, the static casting. This is possible thanks to the void
configure(const ElemUtils::Parameters& parameters) virtual functions, which
developers may use to allow transferring additional parameters to modules. For
instance: using namespace PARTONS;

void configure(const ElemUtils::Parameters& parameters)

//Check if parameter is available. if (parameters.isAvailable(”myKey”))

//Set. int myVariable = parameters.getLastAvailable().toInt(); can be used to
set myVariable variable to $1$ in the following way: //Container for multiple
parameters. ElemUtils::Parameters parameters;

//Set. parameters.add(ElemUtils::Parameter(”myKey”, 1));

//Pass. pMyModule-¿configure(parameters); where ”myKey” is the unique name of
the parameter. The way of supplying modules with additional parameters
presented here plays a crucial role in constructing user interface (UI).

To avoid “manually” repeating some tasks, _the services_ are used. In this
way, users can perform complex tasks in a straightforward and robust manner by
hiding the complexity of low-level functions. The services link such tasks as
parsing input files, running many computations, employing multi-threading
computing, filling the database, printing the output to the screen etc.
Without services, all those tasks would have to be done explicitly by users,
requiring a profound knowledge of the functioning of the project. This way of
working would also be time consuming and vulnerable to mistakes.

##  3 EpIC’s architecture

The generic architecture of the generator, including the flow of data
transferred between the modules, is shown in Fig. 2. Each module takes input
data, processes it according to predefined tasks, and returns the
corresponding output. The actual implementation of some types of modules
depends on the simulated process. For instance, RCModule is the parent class
for DVCSRCModule and TCSRCModule, which are used for DVCS and TCS processes,
respectively. These two types of modules process data according to the same
task predefined in RCModule, but they use different inputs. It is either DVCS
or TCS kinematics. Because of this difference, RCModule has been developed as
a C++ template, allowing DVCSRCModule and TCSRCModule to handle input
containers of various types. With those templates, the addition of new
processes to the generator is as easy as possible. DVCSRCModule and
TCSRCModule are used as base classes for actual modules, where specific
physics developments are encoded. In Fig. 3 we present the entire tree of
inheritance for modules of RCModule type.

Below, we provide a brief description of each type of module and its purpose.
Additionally, we indicate whether the actual implementation of modules of a
given type depends on the simulated process. This information is useful in
recognising which parts of the generator need to be extended when a new
process is added to the generator.

Figure 2: Sketch showing the generic architecture of the project, i.e. without
specifying the corresponding physics process for some modules. Details are
provided in the text. Figure 3: Inheritance tree of modules used for the
implementation of radiative corrections. Each box represents a single C++
class. Two first classes are the part of PARTONS library, the rest are the
part of EpIC project. Classes acting as modules to be used be users are
indicated by boldface fonts.

###  3.1 RCModule

_Purpose:_ Simulation of radiative corrections.

_Does it depend on the generated process?_ Yes.

_Description:_ In EpIC events are generated according to the probability
distribution given by:

$\mathrm{d}\sigma(X^{\prime},Z)=\mathrm{d}\sigma_{0}(X)~{}R(Z)\,.$ (1)

Here, $d\sigma_{0}(X)$ is the Born cross-section (given by PARTONS), while
$R(Z)$ is the radiator function. The set of variables that the Born cross-
section depends on is denoted by $X$, for instance
$X=\\{x_{\mathrm{Bj}},t,Q^{2},\phi,\phi_{S},E_{l}\\}$ for DVCS, where
$x_{\mathrm{Bj}}$ and $Q^{2}$ are the usual DIS variables, $\phi$ is the angle
between the lepton scattering plane and the production plane, $\phi_{S}$ is
the angle between the lepton scattering plane and the target spin component
perpendicular to the direction of the virtual photon, and $E_{l}$ is the
lepton energy in the fixed target frame Berthou:2015oaw . Because of the
radiative corrections, the values of those variables may be different from
those for a generated event, $X^{\prime}$. The radiator function, on the other
hand, depends on the set of variables $Z$. The number of variables making this
set depends on the approximation one uses. For instance, if one considers the
collinear approximation and only the initial and final state radiations of
single photons from the incoming and scattered leptons, two variables make the
set $Z$, namely $Z=\\{z_{1},z_{3}\\}$, see A for more details.

Taking the above in mind, modules of RCModule type perform three actions: _i_)
they define a number of variables that the radiative corrections depend on,
i.e. they define the set $Z$, _ii_) they implement a function used for the
evaluation of $R(Z)$, but also for the evaluation of “true” kinematics
entering the Born cross-section, $X=f(X^{\prime},Z)$, _iii_) they implement a
function that stores the four-momenta of radiated photons in the final event
records.

###  3.2 ProcessModule

_Purpose:_ Evaluation of Born cross-section.

_Does it depend on the generated process?_ Yes.

_Description:_ The Born cross-section is evaluated in EpIC using the modules
of ProcessModule type, which are defined in the PARTONS library. In the most
general case, where the Born cross-section is directly calculated from GPDs,
ProcessModule requires modules of other types for proper functioning:
ConvolCoeffFunctionModule for the evaluation of the process amplitudes,
GPDModule for the evaluation of GPDs at a reference scale, GPDEvolutionModule
for the evolution of GPDs, XiConverterModule for the evaluation of GPD
skewness variable, $\xi$, from the process kinematics, and ScalesModule for
the evaluation of renormalisation and factorisation scales from the process
kinematics. However, this general case can hardly be used in the Monte Carlo
generation because of the nested integrations evaluating the Born cross-
section. Instead, one may start the evaluation, for instance, from the level
of amplitudes parametrised in lookup tables. In this case, GPDModule and
GPDEvolutionModule are not used as there is no need for a convolution of GPDs
with hard scattering coefficient functions. More details about the
ProcessModule modules can be found in Ref. Berthou:2015oaw .

###  3.3 EventGeneratorModule

_Purpose:_ Generation of kinematic configurations according to the probability
distribution given by the product of the Born cross-section and the radiator
function.

_Does it depend on the generated process?_ No.

_Description:_ Modules of this type play a central role in the generation of
events. Three actions are implemented in those modules: _i_) Initialisation,
where the algorithm implemented in the module probes the probability
distribution, for instance to “memorise” it using an interpolation method, or
to capture its key features, like the location of places in which the
probability distribution tends to change rapidly. For some complicated
probability distribution functions, the initialisation process can be time-
consuming, but the results of the process can be stored in a file for use in
subsequent jobs. _ii_) Generation of kinematic configurations according to the
probability distribution, based on the data collected during the
initialisation stage. _iii_) Evaluation of the total cross-section, which is
needed to normalise the event distributions to a given integrated luminosity.

###  3.4 KinematicModule

_Purpose:_ Evaluation of four-momenta for a given kinematic configuration.

_Does it depend on the generated process?_ Yes.

_Description:_ Modules of this type are responsible for the evaluation of
four-momenta for a given kinematic configuration. In addition, they implement
a predefined function that determines if a given kinematic configuration is
physical, i.e. if it does not break any kinematic limit. This function is
utilised by GeneratorModule to probe only valid kinematics.

###  3.5 WriterModule

_Purpose:_ Composing and saving event records.

_Does it depend on the generated process?_ No.

_Description:_ The purpose of modules of this type is to create and store
event records in output files. Different formats may be implemented by
different modules. Additional information, like the integrated cross-section,
can also be stored in output files, e.g. in their headers.

##  4 User interface

Both compilation and linking of the project make use of the CMake tool
www:cmake . In the current version ($1.0.0$) EpIC requires the following
external libraries: PARTONS Berthou:2015oaw  (providing elements of the
architecture and used for the evaluation of Born cross-sections), ROOT
Brun:1997pa  (used in one of the EventGeneratorModule modules), HepMC3
Buckley:2019xhk  (used in one of the WriterModule modules), GSL galassi2009gnu
(used for the generation of random numbers). For more detailed and always up-
to-date information we refer to the online documentation www:epic .

The executable of the project, epic, must be invoked with two arguments, like
this: ./epic –seed=SEED –scenario=SCENARIO_PATH where SEED is the random seed
(unsigned integer) to be used in the initialisation of modules that deal with
random numbers, and SCENARIO_PATH is the relative or absolute path to _the
scenario_ containing all options used in the generation.

The EpIC scenarios are written in XML markup language which, thanks to the
usage of tags based on standard words, is an easy-to-read format for humans
while being straightforward for machines to process. EpIC assumes input data
to be provided in units of GeV (and its powers, whenever applicable), and
radians for angles. The general structure of a single scenario is the
following: ¡!– XML header –¿ ¡?xml version=”1.0” encoding=”UTF-8”
standalone=”yes” ?¿

¡!– Definition of scenario –¿ ¡!– For bookkeeping it includes date and
description –¿ ¡scenario date=”2022-01-01” description=”My first scenario”¿

¡!– Selection of service and its method–¿ ¡task service=”DVCSGeneratorService”
method=”generate”¿

¡!– General configuration –¿ ¡general_configuration¿ …
¡/general_configuration¿

¡!– Selection of kinematic ranges –¿ ¡kinematic_range¿ … ¡/kinematic_range¿

¡!– Indication of experimental conditions –¿ ¡experimental_conditions¿ …
¡/experimental_conditions¿

¡!– Configuration of ProcessModule –¿ ¡computation_configuration¿ …
¡/computation_configuration¿

¡!– Configuration of GeneratorModule –¿ ¡generator_configuration¿ …
¡/generator_configuration¿

¡!– Configuration of KinematicModule –¿ ¡kinematic_configuration¿ …
¡/kinematic_configuration¿

¡!– Configuration of RCModule –¿ ¡rc_configuration¿ … ¡/rc_configuration¿

¡!– Configuration of WriterModule –¿ ¡writer_configuration¿ …
¡/writer_configuration¿

¡/task¿ ¡/scenario¿ Here, we only expose blocks containing specific
information, like ¡general_configuration¿ … ¡/general_configuration¿. Each
block is described in one of the subsequent subsections. In our demonstration
we only show an exemplary scenario for the DVCS case. More examples, including
those for other processes, can be found online www:epic .

###  4.1 general_configuration

¡general_configuration¿

¡!– Number of events to be generated –¿ ¡param name=”number_of_events”
value=”10” /¿

¡!– Subprocess. Possible values: ”DVCS” (pure DVCS), ”BH” (pure Bethe-
Heitler), ”DVCS—BH” (pure DVCS and BH), ”DVCS—BH—INT” (pure DVCS, BH and the
interference between the two) –¿ ¡param name=”subprocess_type” value=”DVCS” /¿

¡/general_configuration¿

###  4.2 kinematic_range

¡kinematic_range¿

¡!– Range of $x_{\text{B}}$ (Bjorken variable)–¿ ¡param name=”range_xB”
value=”0.—1.” /¿

¡!– Range of $t$ –¿ ¡param name=”range_t” value=”-1.—0.” /¿

¡!– Range of $Q^{2}$ –¿ ¡param name=”range_Q2” value=”1.—10.” /¿

¡!– Range of $\phi$ –¿ ¡param name=”range_phi” value=”0.—2*pi” /¿

¡!– Range of $\phi_{S}$ –¿ ¡param name=”range_phiS” value=”0.—2*pi” /¿

¡!– Range of $y$ –¿ ¡param name=”range_y” value=”0.01—0.95” /¿

¡/kinematic_range¿

###  4.3 experimental_conditions

¡experimental_conditions¿

¡!– Energy of the lepton beam –¿ ¡param name=”lepton_energy” value=”5.” /¿

¡!– Type of the lepton beam, here electron –¿ ¡param name=”lepton_type”
value=”e-” /¿

¡!– Polarisation of the lepton beam –¿ ¡param name=”lepton_helicity” value=”1”
/¿

¡!– Energy of the hadron beam –¿ ¡!– For a target in rest frame use:
value=”fixed_target” –¿ ¡param name=”hadron_energy” value=”10.” /¿

¡!– Type of the hadron beam, here proton –¿ ¡param name=”hadron_type”
value=”p” /¿

¡!– Polarisation of the hadron beam –¿ ¡!– Possible values: ”0—0—0” for
unpolarised target, ”0—0—$\pm 1$” for long. polarised target, ”$\pm 1$—0—0” or
”0—$\pm 1$—0” for trans. polarised target –¿ ¡param name=”hadron_polarisation”
value=”0.—0.—0.” /¿

¡/experimental_conditions¿

###  4.4 computation_configuration

¡computation_configuration¿

¡!– Selection of module –¿ ¡!– DVCSProcessBMJ12 module encodes the BH, DVCS,
and interference terms of the Born cross-section according to Refs.
Belitsky:2001ns  and Belitsky:2012ch  –¿ ¡module type=”DVCSProcessModule”
name=”DVCSProcessBMJ12”¿

¡!– Selection of module used for evaluation of factorisation and
renormalisation scales from DVCS kinematics –¿ ¡!– DVCSScalesQ2Multiplier
module identifies both scales in terms of the $Q^{2}$ variable –¿ ¡module
type=”DVCSScalesModule” name=”DVCSScalesQ2Multiplier”¿ ¡/module¿

¡!– Selection of module used for evaluation of the GPD skewness variable,
$\xi$, from DVCS kinematics –¿ ¡!– DVCSXiConverterXBToXi module use the
conversion $\xi=x_{\mathrm{Bj}}/(2-x_{\mathrm{Bj}})$ –¿ ¡module
type=”DVCSXiConverterModule” name=”DVCSXiConverterXBToXi”¿ ¡/module¿

¡!– Selection of module used for the evaluation DVCS Compton form factors
(CFFs) –¿ ¡!– DVCSCFFCMILOU3DTables module does not evaluate CFFs from a GPD
model during EpIC run-time. Instead, it uses look-up tables of CFFs evaluated
ahead of time for a particular GPD module. A number of these tables are
included with EpIC.–¿ ¡module type=”DVCSConvolCoeffFunctionModule”
name=”DVCSCFFCMILOU3DTables”¿

¡!– pQCD order of evaluation –¿ ¡param name=”qcd_order_type” value=”LO” /¿

¡!– Path to a look-up table, here the one based on the GK GPD model is
specified –¿ ¡param name=”cff_set_file” value=”PATH/epic/data/
DVCSCFFCMILOU3DTables/tables_GK.root” /¿ ¡/module¿

¡/module¿ ¡/computation_configuration¿

###  4.5 generator_configuration

¡generator_configuration¿

¡!– Selection of module –¿ ¡!– EventGeneratorFOAM module uses mini-FOAM
library Jadach:2005ex  that is issued with ROOT Brun:1997pa  –¿ ¡module
type=”EventGeneratorModule” name=”EventGeneratorFOAM”¿

¡!– Parameters of the FOAM algorithm –¿

¡!– Maximum number of cells –¿ ¡param name=”nCells” value=”10000” /¿

¡!–Number of MC events when exploring a cell –¿ ¡param name=”nSamples”
value=”2000” /¿

¡!– Number of bins in edge histogram for a cell –¿ ¡param name=”nBins”
value=”2000” /¿

¡!– After the initialisation process, the state of FOAM will be saved in
PATH/state.root. To use this file in another EpIC run, so as to skip the
initialisation in that run, use the read_state_file_path option. –¿ ¡param
name=”save_state_file_path” value=”PATH/state.root” /¿ ¡/module¿
¡/generator_configuration¿

###  4.6 kinematic_configuration

¡kinematic_configuration¿

¡!– Selection of module –¿ ¡!– DVCSKinematicDefault module provides default
evaluation of four-momenta from DVCS kinematics –¿ ¡module
type=”DVCSKinematicModule” name=”DVCSKinematicDefault”¿ ¡/module¿
¡/kinematic_configuration¿

###  4.7 rc_configuration

¡rc_configuration¿

¡!– Selection of module –¿ ¡!– DVCSRCCollinear module provides evaluation of
initial and final state radiative corrections from lepton lines using the
collinear approximation, see Sect. A for more details –¿ ¡module
type=”DVCSRCModule” name=”DVCSRCCollinear”¿ ¡/module¿ ¡/rc_configuration¿

###  4.8 writer_configuration

¡writer_configuration¿

¡!– Selection of module –¿ ¡!– WriterHepMC3 module is used to save event
records in HepMC3 format Buckley:2019xhk  –¿ ¡module type=”WriterModule”
name=”WriterHepMC3”¿

¡!– Path to output file containing event records –¿ ¡param
name=”output_file_path” value=”test.txt” /¿

¡!– Format of output file containing event records: ”ascii” for text format,
or ”root” for binary ROOT format–¿ ¡param name=”HepMC3_writer_type”
value=”ascii” /¿ ¡/module¿ ¡/writer_configuration¿

##  5 Demonstration of EpIC

Here, we demonstrate EpIC in a detailed manner by using specific components of
the PARTONS framework. The demonstration is based on one million MC events
generated for the DVCS sub-process, i.e. exclusive leptoproduction of a single
photon without the Bethe-Heitler contribution, at $10\,\mathrm{GeV}$ positive
helicity electron and $100\,\mathrm{GeV}$ unpolarized proton energies. The
following modules are used during the generation process:
DVCSCFFCMILOU3DTables for the parameterisation of CFFs obtained from the GK
GPD model Goloskokov:2005sd ; Goloskokov:2007nt ; Goloskokov:2009ia  and LO
coefficients functions, DVCSProcessBMJ12 for the evaluation of DVCS cross-
section based on the set of expressions published in Ref. Belitsky:2012ch ,
DVCSRCNull for the simulation without radiative corrections,
EventGeneratorFOAM for the generation of DVCS kinematics based on cross-
section with the FOAM algorithm Jadach:2005ex , DVCSKinematicDefault for the
default evaluation of four-momenta from DVCS kinematics, and WriterHepMC3 to
save four-momenta in HepMC3 format.

The following cuts are used in the generation process: $0.0001\leq
x_{\text{B}}\leq 0.6$, $0.01\leq y\leq 0.95$ (here, $y$ is the inelasticity
parameter), $1\leq Q^{2}\leq 100$ $\text{GeV}^{2}$, $0\leq|t|\leq 1$
$\text{GeV}^{2}$, $0\leq\phi\leq 2\pi$, and $0\leq\phi_{\text{S}}\leq 2\pi$.
On the basis of these cuts, the distribution of events is shown in Fig. 4. The
following FOAM library parameters are used in producing those events
Jadach:2005ex : $\texttt{nCells}=3000$, $\texttt{nSamples}=600$,
$\texttt{nBins}=600$. For the FOAM library, the initialisation lasted about
$40$ minutes, whereas the generation time after initialisation took around
$0.0052$ seconds per event at BNL computer farms. The initialisation time can
be shortened by changing the parameters of the FOAM library, which may,
however, lead to a longer time needed by EpIC to generate a single event. We
remind that the result of the initialisation can be stored in a file, allowing
to avoid this step in subsequent MC runs.

To ensure our results are consistent, we also plot the theory expectation
values, as solid lines, on top of each histogram. The number of events
associated with the theory curves is calculated as follows: at each bin, the
cross-section for the given bin width is calculated, and then the result is
divided by the total cross-section of the phase space. One can then multiply
this ratio by the total number of events to find the number of events in that
specific bin. To obtain those curves in Fig. 4, first we assign the number of
events to the middle point of each bin, then interpolate all these results. As
a result, we observe that EpIC generates events in accordance with the theory
values.

Figure 4: Distribution of one million DVCS events produced by the
DVCSProcessBMJ12 module utilising GPDGK16 GPDs within the PARTONS framework.
Events are binned with respect to
$\text{log}_{10}(x_{\text{B}}),y,Q^{2},|t|,\phi$ and $\phi_{\text{S}}$. The
solid curves on top of the histograms depict the expected number of events for
the given model. The curves are normalised according to the integrated
luminosity estimated by EpIC.

##  6 Summary

In this work, we described in detail the overall structure and essential
features of the novel MC generator EpIC. This generator is based on the
PARTONS framework in the model selection. It is designed to generate events
for exclusive processes, such as DVCS, TCS, DVMP and more broadly to all
exclusive processes implemented in the PARTONS framework. As these processes
have the potential to reveal the 3D structure of the nucleon through GPDs,
they have been studied extensively at facilities such as DESY, JLab, and CERN,
and they will remain crucial for future EIC programmes.

The EpIC MC generator stands out with its flexible architecture that utilises
a modular programming paradigm. Consequently, the code is easier to
understand, as each element of the architecture is designed to perform a
specific task. The distribution of events generated by EpIC is highly precise
since it utilises an accurate representation of the cross-section of the
underlying process in the FOAM library. Besides generating accurate events,
EpIC also allows users to generate events with radiative corrections as well
as to run multi-channel analyses. On top of this, EpIC is released through
GitHub as an open-source project under a GPL licence. In particular it means
that users can trace code modifications and new contributors can easily join
the development team. To the best of our knowledge, a tool such as EpIC is
unique in the GPD community by its very nature. The comprehensive set of
features and the variety of models available in EpIC make it ideal for the
systematic study of the 3D structure of the nucleon.

The generator was designed in a way that facilitates easy expansion. The
encapsulation of elements is one of the fundamental characteristics of the
generator when it comes to maintaining the project for a long time and with a
changing team of developers. This allows us to easily adapt the architecture
to the latest developments. Finally, to make this project as useful to the
particle physics community as possible, we promote open-source standards and
offer easy access to its development. This perspective is well suited to the
experimental timeline of the 3D hadron structure community.

###### Acknowledgements.

The authors thank Alexander Jentsch, Krešimir Kumerički, and Kornelija Passek-
Kumerički for fruitful discussions. This work was supported by the Grant No.
2019/35/D/ST2/00272 of the National Science Center, Poland; by the European
Union’s Horizon 2020 research and innovation programme under grant agreement
No 824093; in the framework of the GLUODYNAMICS project funded by the ”P2IO
LabEx (ANR-10-LABX-0038)” in the framework ”Investissements d’Avenir”
(ANR-11-IDEX-0003-01) managed by the Agence Nationale de la Recherche (ANR),
France; by the Ile-de-France region via the Blaise Pascal Chair of
International Excellence. KT was supported by U.S. Department of Energy under
Contract No. de-sc0012704.

##  Appendix A Radiative corrections

Radiative corrections can have a significant impact on the interpretation of
experimental data, as they introduce uncertainties in the reconstruction of
“true” kinematics of the measured process. On the other hand, neglecting
radiative corrections in simulations typically leads to unrealistic
distributions of generated events, making e.g. impact studies to be much less
robust. To address both of these problems it is highly desirable to have
radiative corrections implemented in generators of MC events.

In the literature, it has been discussed that the leptonic initial and final
state radiative corrections, i.e. when photons are emitted by either incoming
or outgoing leptons, respectively, can be computed with a good accuracy in the
collinear approximation Mo:1968cg ; Kripfganz:1990vm . In this approximation
the transverse component of the momenta of emitted photons are neglected.
Therefore, the radiated photons are restricted to move along the direction of
the source leptons.

The collinear approximation allows one to introduce a single parameter,
$z_{1}$ or $z_{3}$, for the initial or final state radiation, respectively.
The parameter describes the energy of the emitted photon:

$z_{1}=\frac{E_{e}-E_{\gamma}}{E_{e}}\,,\qquad
z_{3}=\frac{E_{e^{\prime}}}{E_{e^{\prime}}+E_{\gamma^{\prime}}}\,,$ (2)

where $E_{e}$ ($E_{e^{\prime}}$) and $E_{\gamma}$ ($E_{\gamma^{\prime}}$) are
the energies of the incoming (outgoing) lepton and the initial (final) state
emitted photon, respectively. The DIS cross section can be then expressed by
Kripfganz:1990vm :

$\displaystyle\frac{\mathrm{d}^{2}\sigma}{\mathrm{d}x\,\mathrm{d}y}=\int_{z_{1}^{\text{min}}}^{1}\frac{\mathrm{d}z_{1}}{z_{1}}D(z_{1})\int_{z_{3}^{\text{min}}}^{1}\frac{\mathrm{d}z_{3}}{z_{3}^{2}}\overline{D}(z_{3})\,\frac{y}{\hat{y}}\,\frac{\mathrm{d}^{2}\hat{\sigma}_{\text{Born}}}{\mathrm{d}\hat{x}\,\mathrm{d}\hat{y}}\,,$
(3)

where $\mathrm{d}^{2}\hat{\sigma}_{\text{Born}}$ is the differential Born
cross-section evaluated for “true” variables, which are related to the
observed ones by Kripfganz:1990vm :

$\displaystyle\hat{x}=\frac{z_{1}xy}{z_{1}z_{3}+y-1}\,,~{}\hat{y}=\frac{z_{1}z_{3}+y-1}{z_{1}z_{3}}\,,$
(4)

and where due to kinematic limits one has:

$z_{1}^{\text{min}}=\frac{1-y}{1-xy},\qquad z_{3}^{\text{min}}=1-y(1-x)\,.$
(5)

In Eq. (3) the radiator functions, $D(z_{1})$ and $\overline{D}(z_{3})$,
describe the weight of the rescaled cross-section in the presence of the
initial and final state radiations, respectively. At the leading order
$D(z_{1})$ and $\overline{D}(z_{3})$ are equal to each other and given by:

$\displaystyle D(z)$
$\displaystyle=\Bigg{[}\delta(1-z)\Bigg{[}1+\frac{\alpha}{2\pi}L\Bigg{(}2\ln\,\epsilon+\frac{3}{2}\Bigg{)}\Bigg{]}$
$\displaystyle+\theta(1-\epsilon-z)\frac{\alpha}{2\pi}L\frac{1+z^{2}}{1-z}\Bigg{]}\,,$
(6)

with $L=\ln(Q^{2}/m_{l}^{2})$. Here, $\alpha$ is the fine-structure constant,
and $m_{l}$ is the lepton mass. The parameter $\epsilon$ acts as a cutoff
avoiding the generation of a very soft photons characterised by energies
$E_{\gamma}\leq\epsilon E_{e}$ and $E_{\gamma^{\prime}}\leq\epsilon
E_{e^{\prime}}$. It is crucial to choose the cutoff carefully so that
$\epsilon E_{e}$ and $\epsilon E_{e^{\prime}}$ are much smaller than the
energy resolution of the experiment. Due to the cutoff Eq. (6) consists of two
parts: the term proportional to the Dirac delta function collectively
describes the contribution of very soft photons and virtual corrections, while
the term proportional to the step function describes the effects of real
radiation of photons with sizeable energies.

The generalisation of leading order radiative corrections for the DVCS case
can be achieved as follows:

$\displaystyle\frac{\mathrm{d}^{5}\sigma}{\mathrm{d}x\,\mathrm{d}Q^{2}\,\mathrm{d}t\,\mathrm{d}\phi\,\mathrm{d}\phi_{S}}$
$\displaystyle=\int_{z_{1}^{\text{min}}}^{1}\mathrm{d}z_{1}z_{1}D(z_{1})\int_{z_{3}^{\text{min}}}^{1}\frac{\mathrm{d}z_{3}}{z_{3}^{2}}\overline{D}(z_{3})$
$\displaystyle\times\frac{y}{\hat{y}}\,\frac{\mathrm{d}^{5}\hat{\sigma}_{\text{Born}}}{\mathrm{d}\hat{x}\,\mathrm{d}\hat{Q}^{2}\,\mathrm{d}t\,\mathrm{d}\,\phi\,\mathrm{d}\phi_{S}}\,.$
(7)

## References

  * (1) D. Müller, D. Robaschik, B. Geyer, F.M. Dittes, J. Hořejši, Fortsch. Phys. 42, 101 (1994).  DOI 10.1002/prop.2190420202
  * (2) X.D. Ji, Phys. Rev. Lett. 78, 610 (1997).  DOI 10.1103/PhysRevLett.78.610
  * (3) X.D. Ji, Phys. Rev. D 55, 7114 (1997).  DOI 10.1103/PhysRevD.55.7114
  * (4) A.V. Radyushkin, Phys. Lett. B 385, 333 (1996).  DOI 10.1016/0370-2693(96)00844-1
  * (5) A.V. Radyushkin, Phys. Rev. D 56, 5524 (1997).  DOI 10.1103/PhysRevD.56.5524
  * (6) M. Burkardt, Phys. Rev. D 62, 071503 (2000).  DOI 10.1103/PhysRevD.62.071503.  [Erratum: Phys.Rev.D 66, 119903 (2002)]
  * (7) M. Burkardt, Int. J. Mod. Phys. A 18, 173 (2003).  DOI 10.1142/S0217751X03012370
  * (8) M. Burkardt, Phys. Lett. B 595, 245 (2004).  DOI 10.1016/j.physletb.2004.05.070
  * (9) X.D. Ji, W. Melnitchouk, X. Song, Phys. Rev. D 56, 5511 (1997).  DOI 10.1103/PhysRevD.56.5511
  * (10) M.V. Polyakov, Phys. Lett. B 555, 57 (2003).  DOI 10.1016/S0370-2693(03)00036-4
  * (11) M.V. Polyakov, P. Schweitzer, Int. J. Mod. Phys. A 33(26), 1830025 (2018).  DOI 10.1142/S0217751X18300259
  * (12) C. Lorcé, H. Moutarde, A.P. Trawiński, Eur. Phys. J. C 79(1), 89 (2019).  DOI 10.1140/epjc/s10052-019-6572-3
  * (13) M. Diehl, Phys. Rept. 388, 41 (2003).  DOI 10.1016/j.physrep.2003.08.002
  * (14) A. Accardi, et al., Eur. Phys. J. A 52(9), 268 (2016).  DOI 10.1140/epja/i2016-16268-9
  * (15) R. Abdul Khalek, et al., Science Requirements and Detector Concepts for the Electron-Ion Collider: EIC Yellow Report (2021).  arXiv:physics.ins-det/2103.05419
  * (16) D.P. Anderle, et al., Front. Phys. (Beijing) 16(6), 64701 (2021).  DOI 10.1007/s11467-021-1062-0
  * (17) J.L. Abelleira Fernandez, et al., J. Phys. G 39, 075001 (2012).  DOI 10.1088/0954-3899/39/7/075001
  * (18) B. Adams, et al., COMPASS++/AMBER: Proposal for Measurements at the M2 beam line of the CERN SPS Phase-1: 2022-2024.  Tech. rep., CERN, Geneva (2019).  CERN-SPSC-2019-022, SPSC-P-360
  * (19) J. Arrington, et al., Physics with CEBAF at 12 GeV and Future Opportunities (2021).  arXiv:nucl-ex/2112.00060
  * (20) B. Berthou, et al., Eur. Phys. J. C 78(6), 478 (2018).  DOI 10.1140/epjc/s10052-018-5948-0
  * (21) A. Bacchetta, U. D’Alesio, M. Diehl, C.A. Miller, Phys. Rev. D 70, 117504 (2004).  DOI 10.1103/PhysRevD.70.117504
  * (22) EpIC project on GitHub.  https://github.com/pawelsznajder/epic
  * (23) EpIC project webpage.  https://pawelsznajder.github.io/epic
  * (24) O. Grocholski, H. Moutarde, B. Pire, P. Sznajder, J. Wagner, Eur. Phys. J. C 80(2), 171 (2020).  DOI 10.1140/epjc/s10052-020-7700-9
  * (25) H. Moutarde, P. Sznajder, J. Wagner, Eur. Phys. J. C 78(11), 890 (2018).  DOI 10.1140/epjc/s10052-018-6359-y
  * (26) H. Moutarde, P. Sznajder, J. Wagner, Eur. Phys. J. C 79(7), 614 (2019).  DOI 10.1140/epjc/s10052-019-7117-5
  * (27) V. Bertone, H. Dutrieux, C. Mezrag, H. Moutarde, P. Sznajder, Phys. Rev. D 103(11), 114019 (2021).  DOI 10.1103/PhysRevD.103.114019
  * (28) J.M.M. Chavez, V. Bertone, F.D.S. Borrero, M. Defurne, C. Mezrag, H. Moutarde, J. Rodríguez-Quintero, J. Segovia, Pion GPDs: A path toward phenomenology (2021).  arXiv:hep-ph/2110.06052
  * (29) H. Dutrieux, O. Grocholski, H. Moutarde, P. Sznajder, Eur. Phys. J. C 82(3), 252 (2022).  DOI 10.1140/epjc/s10052-022-10211-5
  * (30) S.V. Goloskokov, P. Kroll, Eur. Phys. J. C 42, 281 (2005).  DOI 10.1140/epjc/s2005-02298-5
  * (31) S.V. Goloskokov, P. Kroll, Eur. Phys. J. C 53, 367 (2008).  DOI 10.1140/epjc/s10052-007-0466-5
  * (32) S.V. Goloskokov, P. Kroll, Eur. Phys. J. C 65, 137 (2010).  DOI 10.1140/epjc/s10052-009-1178-9
  * (33) CMake project webpage.  https://cmake.org
  * (34) R. Brun, F. Rademakers, Nucl. Instrum. Meth. A 389, 81 (1997).  DOI 10.1016/S0168-9002(97)00048-X
  * (35) A. Buckley, P. Ilten, D. Konstantinov, L. Lönnblad, J. Monk, W. Pokorski, T. Przedzinski, A. Verbytskyi, Comput. Phys. Commun. 260, 107310 (2021).  DOI 10.1016/j.cpc.2020.107310
  * (36) M. Galassi, et al., _GNU Scientific Library: Reference Manual_ (Network Theory, 2009).  ISBN 0954612078
  * (37) A.V. Belitsky, D. Mueller, A. Kirchner, Nucl. Phys. B 629, 323 (2002).  DOI 10.1016/S0550-3213(02)00144-X
  * (38) A.V. Belitsky, D. Müller, Y. Ji, Nucl. Phys. B 878, 214 (2014).  DOI 10.1016/j.nuclphysb.2013.11.014
  * (39) S. Jadach, P. Sawicki, Comput. Phys. Commun. 177, 441 (2007).  DOI 10.1016/j.cpc.2007.02.112
  * (40) L.W. Mo, Y.S. Tsai, Rev. Mod. Phys. 41, 205 (1969).  DOI 10.1103/RevModPhys.41.205
  * (41) J. Kripfganz, H.J. Mohring, H. Spiesberger, Z. Phys. C 49, 501 (1991).  DOI 10.1007/BF01549704