Technical aspects of McMule

In this section, we will review the very technical details of the implementation. This is meant for those readers, who wish to truly understand the nuts and bolts holding the code together. We begin by discussing the phase-space generation and potential pitfalls in Section Phase-space generation. Next, in Section Implementation of FKS schemes, we discuss how the FKS scheme [8, 23, 25, 28, 29]. This is meant for those readers, who wish to truly understand the nuts and bolts holding the code together. We begin by discussing the phase-space generation and potential pitfalls in Section Phase-space generation. Next, in Section Implementation of FKS schemes, we discuss how the FKS scheme [8, 23, 25, 28, 29] (cf. Appendix The FKS^2 scheme for a review) is implemented in Fortran code. This is followed by a brief review of the random number generator used in McMule in Section Random number generation. Finally, we give an account of how the statefiles work and how they are used to store distributions in Section Differential distributions and intermediary state files.

Phase-space generation

We use the vegas algorithm for numerical integration [15]. As vegas only works on the hypercube, we need a routine that maps \([0,1]^{3n-4}\) to the momenta of an \(n\)-particle final state, including the corresponding Jacobian. The simplest way to do this uses iterative two-particle phase-spaces and boosting the generated momenta all back into the frame under consideration. An example of how this is done is shown in Listing 10.

Listing 10 Example implementation of iterative phase-space. Not shown are the checks to make sure that all particles have at least enough energy for their mass, i.e. that \(E_i \ge m_i\).

! use a random number to decide how much energy should 
! go into the first particle
minv3 = ra(1)*energy

! use two random numbers to generate the momenta of 
! particles 1 and the remainder in the CMS frame
call pair_dec(ra(2:3),energy,q2,m2,qq3,minv3)

! adjust the Jacobian
weight = minv3*energy/pi
weight = weight*0.125*sq_lambda(energy**2,m2,minv3)/energy**2/pi

! use a random number to decide how much energy should 
! go into the second particle
minv4 = ra(4)*energy
! use two random numbers to generate the momenta of 
! particles 2 and the remainder in their rest frame
call pair_dec(ra(5:6),minv3,q3,m3,qq4,minv4)

! adjust the Jacobian
weight = weight*minv4*energy/pi
weight = weight*0.125*sq_lambda(minv3**2,m3,minv4)/minv3**2/pi

! repeat this process until all particles are generated

! boost all generated particles back into the CMS frame

q4 = boost_back(qq4, q4)
q5 = boost_back(qq4, q5)

q3 = boost_back(qq3, q3)
q4 = boost_back(qq3, q4)
q5 = boost_back(qq3, q5)

As soon as we start using FKS, we cannot use this simplistic approach any longer. The \(c\)-distributions of FKS require the photon energies \(\xi_i\) to be variables of the integration. We can fix this by first generating the photon explicitly as

(4)\[k_1 = p_{n+1} = \frac{\sqrt s}2 \xi_1 (1,\sqrt{1-y_1 ^2}\vec e_\perp,y_1 )\]

where \(\vec e_\perp\) is a \((d-2)\) dimensional unit vector and the ranges of \(y_1\) (the cosine of the angle) and \(\xi_1\) (the scaled energy) are \(-1\le y_1 \le 1\) and \(0\le\xi_1 \le\xi_{\text{max}}\), respectively. The upper bound \(\xi_{\text{max}}\) depends on the masses of the outgoing particles. Following [28] we find

\[\xi_{\text{max}} = 1-\frac{\Big(\sum_i m_i\Big)^2}{s}\]

Finally, the remaining particles are generated iteratively again. This can always be done and is guaranteed to work.

For processes with one or more PCSs this approach is suboptimal. The numerical integration can be improved by orders of magnitude by aligning the pseudo-singular contribution to one of the variables of the integration, as this allows vegas to optimise the integration procedure accordingly. As an example, consider once again \(\mu\to\nu\bar\nu e\gamma\). The PCS comes from

\[\mathcal{M}_{n+1}^{(\ell )} \propto \frac{1}{q\cdot k} = \frac1{\xi^2}\frac1{1-y\beta}\]

where \(y\) is the angle between photon (\(k\)) and electron (\(q\)). For large velocities \(\beta\) (or equivalently small masses), this becomes almost singular as \(y\to1\). If now \(y\) is a variable of the integration this can be mediated. An example implementation is shown in Listing 11.

Listing 11 Example implementation of a so-called FKS phase-space where the fifth particle is an FKS photon that may becomes soft. Not shown are checks whether \(E_i \ge m_i\) .

xi5 = ra(1)
y2 = 2*ra(2) - 1.

! generate electron q2 and photon q5 s.t. that the 
! photon goes into z diractions

eme = energy*ra(3)
pme = sqrt(eme**2-m2**2)
q2 = (/ 0., pme*sqrt(1. - y2**2), pme*y2, eme /)
q5 = (/ 0., 0.                  , 1.    , 1.  /)
q5 = 0.5*energy*xi5*q5

! generate euler angles and rotate all momenta

euler_mat = get_euler_mat(ra(4:6))

q2 = matmul(euler_mat,q2)
q5 = matmul(euler_mat,q5)

qq34 = q1-q2-q5
minv34 = sqrt(sq(qq34))

! The event weight, note that a factor xi5**2 has been ommited
weight = energy**3*pme/(4.*(2.*pi)**4)

        ! generate remaining neutrino momenta

call pair_dec(ra(7:8),minv34,q3,m3,q4,m4,enough_energy)
weight = weight*0.125*sq_lambda(minv34**2,m3,m4)/minv34**2/pi

q3 = boost_back(qq34, q3)
q4 = boost_back(qq34, q4)

The approach outlined above is very easy to do in the case of the muon decay as the neutrinos can absorb any timelike four-momentum. This is because the \(\delta\) function of the phase-space was solved through the neutrino’s pair_dec. However, for scattering processes where all final state leptons could be measured, this fails. Writing a routine for \(\mu\)-\(e\)-scattering

\[e(p_1)+\mu(p_2) \to e(p_3)+\mu(p_4) + \gamma(p_5)\]

that optimises on the incoming electron is rather trivial because its direction stays fixed s.t. the photon just needs to be generated according to (4). The outgoing electron \(p_3\) is more complicated. Writing the \(p_4\)-phase-space four- instead of three-dimensional

\[\mathrm{d}\Phi_5 = \delta^{(4)}(p_1+p_2-p_3-p_4-p_5) \delta(p_4^2-M^2) \Theta(E_4) \frac{\mathrm{d}^4\vec p_4}{(2\pi)^4} \frac{\mathrm{d}^3\vec p_3}{(2\pi)^32E_3} \frac{\mathrm{d}^3\vec p_5}{(2\pi)^32E_5}\]

we can solve the four-dimensional \(\delta\) function for \(p_4\) and proceed for the generation \(p_3\) and \(p_5\) almost as for the muon decay above. Doing this we obtain for the final \(\delta\) function

(5)\[\delta(p_4^2-M^2) = \delta\bigg( m^2-M^2+s(1-\xi) +E_3\sqrt{s}\Big[\xi-2-y\xi\beta_3(E_3)\Big] \bigg)\,.\]

When solving this for \(E_3\), we need to take care to avoid extraneous solutions of this radical equation [11]. We have now obtained our phase-space parametrisation, albeit with one caveat: for anti-collinear photons, i.e. \(-1<y<0\) with energies

\[\xi_1 = 1-\frac{m}{\sqrt{s}}+\frac{M^2}{\sqrt{s}(m-\sqrt{s}} < \xi < \xi_{\text{max}} = 1-\frac{(m+M)^2}{s}\]

there are still two solutions. One of these corresponds to very low-energy electron that are almost produced at rest. This is rather fortunate as most experiments will have an electron detection threshold higher that this. Otherwise, phase-spaces optimised this way also define a which_piece for this corner region.

There is one last subtlety when it comes to these type of phase-space optimisations. Optimising the phase-space for emission from one leg often has adverse effects on terms with dominant emission from another leg. In other words, the numerical integration works best if there is only one PCS on which the phase-space is tuned. As most processes have more than one PCS we need to resort to something that was already discussed in the original FKS paper [29]. Scattering processes that involve multiple massless particles have overlapping singular regions. The FKS scheme now mandates that the phase-space is partitioned in such a way as to isolate at most one singularity per region with each region having its own phase-space parametrisation. Similarly we have to split the phase-space to contain at most one PCS as well as the soft singularity. In McMule \(\mu\)-\(e\) scattering for instance is split as follows [1]

\[1 = \theta\big( s_{15} > s_{35} \big) + \theta\big( s_{15} < s_{35} \big)\]

with \(s_{ij} = 2p_i\cdot p_j\) as usual. The integrand of the first \(\theta\) function has a final-state PCS and hence we use the parametrisation obtained by solving (5). The second \(\theta\) function, on the other hand, has an initial-state PCS which can be treated by just directly parametrising the photon in the centre-of-mass frame as per (4). This automatically makes \(s_{15}\propto(1-\beta_{\text{in}}y_1)\) a variable of the integration.

For the double-real corrections of \(\mu\)-\(e\) scattering, we proceed along the same lines except now the argument of the \(\delta\) function is more complicated.

Implementation of FKS schemes

Now that we have a phase-space routine that has \(\xi_i\) as variables of the integration, we can start implementing the relevant \(c\)-distributions (11)

\begin{align} \mathrm{d} \sigma_h^{(1)}(\xi_c) &= \mathrm{d}\Upsilon_1 \mathrm{d}\Phi_{n,1} \Bigg(\frac1{\xi_1}\Bigg)_c \mathrm{d}\xi_1 \Big(\xi_1^2\mathcal{M}_{n+1}^{(0)}\Big)\\ &= \mathrm{d}\xi_1 \Bigg( \mathrm{d}\Upsilon_1 \mathrm{d}\Phi_{n,1} \Big(\xi_1^2\mathcal{M}_{n+1}^{(0)}\Big) - \mathrm{d}\Upsilon_1 \mathrm{d}\Phi_{n,1} \Big(\mathcal{E}\mathcal{M}_{n}^{(0)}\Big) \theta(\xi_c-\xi_1) \Bigg) \end{align}

We refer to the first term as the event and the second as the counter-event.

Note that, due to the presence of \(\delta(\xi_1)\) in the counter-event (that is implemented through the eikonal factor \(\mathcal{E}\)) the momenta generated by the phase-space \(\mathrm{d}\Upsilon_1\mathrm{d}\Phi_{n,1}\) are different. Thus, it is possible that the momenta of the event pass the cuts or on-shell conditions, while those of the counter event fail, or vice versa. This subtlety is extremely important to properly implement the FKS scheme and many problems fundamentally trace back to this.

Finally, we should note that, in order to increase numerical stability, we introduce cuts on \(\xi\) and sometimes also on a parameter that encodes the PCS such as \(y={\tt y2}\) in (4) and Listing 11. Events that have values of \(\xi\) smaller than this soft cut are discarded immediately and no subtraction is considered. The dependence on this slicing parameter is not expected to drop out completely and hence, the soft cut has to be chosen small enough to not influence the result.

An example implementation can be found in Listing 12.

Listing 12 An example implementation of the FKS scheme in Fortran. Not shown are various checks performed, the binning as well as initialisation blocks.

FUNCTION SIGMA_1(x, wgt, ndim)

! The first random number x(1) is xi.
arr = x

! Generate momenta for the event using the function pointer ps
call gen_mom_fks(ps, x, masses(1:nparticle), vecs, weight)

! Whether unphysical or not, take the value of xi
xifix = xiout

! Check if the event is physical ...
if(weight > zero ) then
  ! and whether is passes the cuts
  var = quant(vecs(:,1), vecs(:,2), vecs(:,3), vecs(:,4), ...)
  cuts = any(pass_cut)
  if(cuts) then
    ! Calculate the xi**2 * M_{n+1}^0 using the pointer matel
    mat = matel(vecs(:,1), vecs(:,2), vecs(:,3), vecs(:,4), ...)
    mat = xifix*weight*mat
    sigma_1 = mat
  end if
end if

! Check whether soft subtraction is required
if(xifix < xicut1) then
  ! Implement the delta function and regenerate events
  arr(1) = 0._prec
  call gen_mom_fks(ps, arr, masses(1:nparticle), vecs, weight)
  ! Check whether to include the counter event
  if(weight > zero) then
    var = quant(vecs(:,1), vecs(:,2), vecs(:,3), vecs(:,4), ...)
    cuts = any(pass_cut)
    if(cuts) then
      mat = matel_s(vecs(:,1), vecs(:,2), vecs(:,3), vecs(:,4), ...)
      mat = weight*mat/xifix
      sigma_1 = sigma_1 - mat
    endif
  endif
endif
END FUNCTION SIGMA_1

Calling procedures and function pointers

McMule uses function pointers to keep track of which functions to call for the integrand, phase-space routine, and matrix element(s). These pointers are assigned during init_piece() and then called throughout integrands and phase_space. The pointers for the phase-space generator and integrand are just assigned using the => operator, i.e.

ps => psx2 ; fxn => sigma_0

The relevant abstract interface for the integrand fxn is

abstract interface
  function integrand(x,wgt,ndim)
    import prec
    integer :: ndim
    real(kind=prec) :: x(ndim),wgt
    real(kind=prec) :: integrand
  end function integrand
end interface

Doing the same for the matrix elements is not possible as they do not have a consistent interface. Instead, we are using a C function set_func that is implemented in a separate file to assign the functions, ignoring the interface

call set_func('00000000', pm2enngav)
call set_func('00000001', pm2ennav)
call set_func('11111111', m2enn_part)

The first argument corresponds to the type of functions that is being set.

Table 2 Arguments for set_func

Bitmask	Name	Description
00000000	matel0	hard matrix element
00000001	matel1	reduced matrix element
00000010	matel2	doubly reduced matrix element
11111111	partfunc	particle string function
10000001	matel_s	single soft limit
10000010	matel_hs	hard-soft limit
10000100	matel_sh	soft-hard limit
10000110	matel_ss	double soft limit

If the soft limits are not assigned, they are auto-generated using the partfunc.

Optional parameters for integrands

The integration is configured during the initpiece() routine. Additionally to identifying what is to be integrated (cf. Section Calling procedures and function pointers), one also configures other parameters such as the dimensionality or the masses involved.

Table 3 Frozen Delights!

Variable	Type	Description	Required
nparticle	integer	the number of total particles (initial & final)	yes
ndim	integer	the dimensionality of the phase space. usually this is \(3n_f-4\), for calculations with extra integrations, these are included	yes
masses	real(:)	the masses of all particles	yes
xicut1	real	the value of \(\xi_c\) for the first subtraction	for real corrections
xicut2	real	the value of \(\xi_c\) for the second subtraction	for double-real corrections
xieik1	real	the value of \(\xi_c\) for the first eikonal	for virtual or real-virtual corrections
xieik2	real	the value of \(\xi_c\) for the second eikonal	for double-virtual corrections
polarised	integer	the number of polarised particles	no, defaults to 0
symfac	real	the symmetry factor for indistinguishable final states	no, defaults to 1
softcut	real	the soft cut parameter	no, but recommended, defaults to 0
collcut	real	the collinear cut parameter	no, but recommended, defaults to 0
ntsSwitch	real	the NTS switching point	only for NTS matrix elemnts

\(\xi_c\) parameters

For the \(\xi_c\) parameters, the user enters a value between zero (exclusive) and one (inclusive). However, the FKS procedure requires the bounds of (12) and the parameters hence need to be rescaled accordingly. In principle the user may enter two different values (xinormcut = xinormcut1 and xinormcut2) though this is rarely called for.

Soft and collinear cut parameter

To improve numerical stability, we set events that have a value of \(\xi\) (\(y\)) lower than softcut (collcut) to zero.

Warning

This introduces a systematic error that needs to be studied. For small values, the improvement in stability is generally worth a small error that is anyway drowned out by the statistical error

This means that we are changing the integration (6)

\begin{align} \mathrm{d} \sigma_h^{(1)}(\xi_c) \to \mathrm{d}&\xi_1 \times \theta(\xi-{\tt softcut}) \\&\times \Bigg( \mathrm{d}\Upsilon_1 \mathrm{d}\Phi_{n,1} \Big(\xi_1^2\mathcal{M}_{n+1}^{(0)}\Big) - \mathrm{d}\Upsilon_1 \mathrm{d}\Phi_{n,1} \Big(\mathcal{E}\mathcal{M}_{n}^{(0)}\Big) \theta(\xi_c-\xi_1) \Bigg) \end{align}

and similarly with collcut. We have found that values of softcut = 1e-10 and collcut = 1.e-11 give reliable results.

Random number generation

A Monte Carlo integrator relies on a (pseudo) random number generator (RNG or PRNG) to work. The pseudo-random numbers need to be of high enough quality, i.e. have no discernible pattern and a long period, to consider each point of the integration independent but the RNG needs to be simple enough to be called many billion times without being a significant source of runtime. RNGs used in Monte Carlo applications are generally poor in quality and often predictable s.t. they could not be used for cryptographic applications.

A commonly used trade-off between unpredictability and simplicity, both in speed and implementation, is the Park-Miller RNG, also known as minstd[19]. As a linear congruential generator, its \((k+1)\)th output \(x_{k+1}\) can be found as

\[z_{k+1} = a\cdot z_k \ \text{mod}\ m = a^{k+1} z_1 \ \text{mod}\ m \qquad\text{and}\qquad x_k = z_k / m \in (0,1)\]

where \(m\) is a large, preferably prime, number and \(2<a<m-1\) an integer. The initial value \(z_1\) is called the random seed and is chosen integer between 1 and \(m-1\). It can easily be seen that any such RNG has a fixed period [2] \(p<m\) s.t. \(z_{k+p} = z_k\) because any \(z_{k+1}\) only depends on \(z_k\) and there are finitely many possible \(z_k\). We call the RNG attached to \((m,a)\) to be of full period if \(p=m-1\), i.e. all integers between 1 and \(m-1\) appear in the sequence \(z_k\).

Assuming \(z_1=1\) then the existence of \(p\) s.t. \(z_{p+1}=1\) is guaranteed by Fermat’s Theorem [3]. This means that the RNG is of full period iff \(a\) is a primitive root modulo \(m\), i.e.

\[\forall g\ \text{co-prime to $m$}\quad \exists k\in\mathbb{Z} \quad\text{s.t.}\quad a^k\equiv g\ (\text{mod}\ m)\]

Park and Miller suggest to use the Mersenne prime \(m=2^{31}-1\), noting that there are 534,600,000 primitive roots of which 7 is the smallest. Because \(7^b\ \text{mod}\ m\) is also a primitive root as long as \(b\) is co-prime to \((m\!-\!1)\), [19] settled on \(b=5\), i.e. \(a=16807\) as a good choice for the multiplier that, per construction, has full period and passes certain tests of randomness.

The points generated by any such RNG will fall into \(\sqrt[n]{n!\cdot m}\) hyperplanes if scattered in an \(n\) dimensional space [16]. However, for bad choices of the multiplier \(a\) the number of planes can be a lot smaller [4].

Presently, the period length of \(p=m-1=2^{31}-2\) is believed to be sufficient though detailed studies quantifying this would be welcome.

Differential distributions and intermediary state files

Distributions are always calculated as histograms by binning each event according to its value for the observable \(S\). This is done by having an \((n_b\times n_q)\)-dimensional array [5] quant() where \(n_q\) is the number of histograms to be calculated (nr_q) and \(n_b\) is the number of bins used (nr_bins). The weight of each event \(\mathrm{d}\Phi\times\mathcal{M} \times w\) is added to the correct entry in bit_it where \(w={\tt wgt}\) is the event weight assigned by vegas.

After each iteration of vegas we add quant() (\({\tt quant}^2\)) to an accumulator of the same dimensions called quantsum (quantsumsq). After \(i\) iterations, we can calculate the value and error as

\[\frac{\mathrm{d}\sigma}{\mathrm{d}S} \approx \frac{\tt quantsum}{\Delta\times i} \qquad\text{and}\qquad \delta\bigg(\frac{\mathrm{d}\sigma}{\mathrm{d}S}\bigg)\approx \frac1\Delta \sqrt{ \frac{{\tt quantsumsq}-{\tt quantsum}^2/i}{i(i-1)} }\]

where \(\Delta\) is the bin-size.

Related to this discussion is the concept of intermediary state files. Their purpose is to record the complete state of the integrator after every iteration in order to recover should the program crash – or more likely be interrupted by a batch system. McMule uses a custom file format .vegas for this purpose which uses Fortran’s record-based (instead of stream- or byte-based) format. This means that each entry starts with 32bit unsigned integer, i.e. 4 byte, indicating the record’s size and ends with the same 32bit integer. As this is automatically done for each record, it minimises the amount of metadata that have to be written.

The current version (v3) must begin with the magic header and version self-identification shown in Table 4. The latter includes file version information and the first five characters the source tree’s SHA1 hash, obtained using make hash.

The header is followed by records describing the state of the integrator as shown in Table 5. Additionally to information required to continue integration such as the current value and grid information, this file also has 300 bytes for a message. This is usually set by the routine to store information on the fate of the integration such as whether it was so-far uninterrupted or whether there is reason to believe it to be inconsistent.

The latter point is particularly important. While McMule cannot read intermediary files from a different version of the file format, it will continue any integration for which it can read the state file. This also includes cases where the source tree has been changed. In this case McMule prints a warning but continues the integration deriving potentially inconsistent results.

Table 4 The magic header and version information used by \(v_3\). \(v_1\) indicates the current version number and \(v_2\) whether long integers are used (L) or not (N). \(s_1\)-\(s_5\) indicate the first five characters of the SHA1 hash produced by the source code at compile time (make hash).

offset	00	01	02	03	04	05	06	07	08	09	0A	0B	0C	0D	0E	0F
hex	09	00	00	00	20	4D	63	4D	75	6C	65	20	20	09	00	00
ASCII	\t				‘ ‘	M	c	M	u	l	e	‘ ‘	‘ ‘	\t
offset	10	11	12	13	14	15	16	17	18	19	1A	1B	1C	1D	1E	1F
hex	00	0A	00	00	00	76	xx	xx	20	20	20	20	20	20	20	0A
ASCII		\n				v	\(v_1\)	\(v_2\)	‘ ‘	‘ ‘	‘ ‘	‘ ‘	‘ ‘	‘ ‘	‘ ‘	\n
offset	20	21	22	23	24	25	26	27	28	29	2A	2B	2C	2D	2E	2F
hex	00	00	00	05	00	00	00	xx	xx	xx	xx	xx	05	00	00	00
ASCII								\(s_1\)	\(s_2\)	\(s_3\)	\(s_4\)	\(s_5\)

Table 5 The body of a .vegas file storing all important information. Each horizontal line indicates as dressed record. In the offset and length columns, all integers are in hexadecimal notation. Negative numbers count from the end of file (EOF).

Off	Len	Type	Var.	Comment
0030	000C	integer	it	the current iteration
003C	000C	integer	ndo	subdiv. on an axis
0048	0010	real	si	\(\sigma/(\delta\sigma)^2\)
0058	0010	real	swgt	\(1/(\delta\sigma)^2\)
0068	0010	real	schi	\((1-{\tt it})\chi + \sigma^2/(\delta\sigma)^2\)
0078	1A98	real(50,17)	xi	the integration grid
1B10	000C	integer	randy	the current random number seed
1B1C	0014	integer integer integer	\(n_q\) \(n_b\) \(n_s\)	number of histograms number of bins len. histogram name
1B30	\(10n_q+8\) \(n_sn_q+8\) \(10n_q(n_b+2)+8\)	real(n_q) real(n_q) character(n_s,n_q) real(n_q,n_b+2) real(n_q,n_b+2)	minv maxv names quantsum quantsumsq	lower bounds upper bounds names of \(S\) accu. histograms accu. histograms squared
-0144	0010	real	time	current runtime in seconds
-0134	0134	character(300)	msg	any message
-0000	EOF

Basics of containerisation

McMule is Docker-compatible. Production runs should be performed with Docker [17], or its user-space complement udocker [12], to facilitate reproducibility and data retention. On Linux, Docker uses chroot to simulate an operating system with McMule installed. In our case, the underlying system is Alpine Linux, a Linux distribution that is approximately 5MB in size.

Terminology

To understand Docker, we need to introduce some terms

• An image is a representation of the system’s ’hard disk’. One host system can have multiple images. In (u)Docker, the images can be listed with docker image ls (udocker images).

• Images can be have names, called tags, otherwise Docker assigns a name as the SHA256 hash.

• Because keeping multiple full file systems is rather wasteful, images are split into layers that can be shared among images. In uDocker, these are tar files containing the changes made to the file system.

• To execute an image, a container needs to be generated. Essentially, this involves uncompressing all layers into a directory an chrooting into said directory.

It is important to note, that containers are ephemeral, i.e. changes made to the container are not stored unless explicitly requested. This is usually not required anyway.

For external interfacing, folders of the host system are mounted into the container.

Building images

Docker images are built using Dockerfiles, a set of instruction on how to create the image from external information and a base image. To speed up building of the image, McMule uses a custom base image called mcmule-pre that is constructed as follows

FROM alpine:3.11

LABEL maintainer="yannick.ulrich@psi.ch"
LABEL version="1.0"
LABEL description="The base image for the full McMule suite"

# Install a bunch of things
RUN apk add py3-numpy py3-scipy ipython py3-pip git tar gfortran gcc make curl musl-dev
RUN echo "http://dl-8.alpinelinux.org/alpine/edge/community" >> /etc/apk/repositories && \
    apk add py3-matplotlib && \
    sed -i '$ d' /etc/apk/repositories

On top of this, McMule is build

FROM yulrich/mcmule-pre:1.0.0

LABEL maintainer="yannick.ulrich@psi.ch"
LABEL version="1.0"
LABEL description="The full McMule suite"
RUN pip3 install git+gitlab.com/mule-tools/pymule.git
COPY . /monte-carlo
WORKDIR /monte-carlo
RUN ./configure
RUN make

To build this image, run

mcmule$ docker  build -t $mytagname . # Using Docker
mcmule$ udocker build -t=$mytagname . # Using udocker

The CI system uses udocker to perform builds after each push. Note that using udocker for building requires a patched version of the code that is available from the McMule collaboration.

Creating containers and running

In Docker, containers are usually created and run in one command

$ docker run --rm $imagename $cmd

The flag –rm makes sure the container is deleted after it is completed. If the command is a shell (usually ash), the flag -i also needs to be provided.

For udocker, creation and running can be done in two steps

$ udocker create $imagename
# this prints the container id
$ udocker run $containerid $cmd
# work in container
$ udocker rm $containerid

or in one step

$ udocker run --rm $imagename $cmd

Running containers can be listed with udocker ps and docker ps. For further details, the reader is pointed to the manuals of Docker and udocker.

[1]

When implementing this, care must be taken to ensure that the split is also well defined if the photon is soft, i.e. if \(\xi=0\).

[2]

Note that, because of the simple recursion the RNG will not repeat any number until the full period is complete

[3]

If \(p\) is prime, for any integer \(a\), \(a^p-a\) is a multiple of \(p\).

[4]

An infamous example is randu that used \(a=2^{16}+3\) and \(m=2^{31}\) that in three dimension produces only 15 planes instead of the maximum 2344.

[5]

To be precise, the actual dimensions are \((n_b+2)\times n_q\) to accommodate under- and overflow bins