.. _sec_techno: 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 :ref:`sec_ps`. Next, in Section :ref:`sec_fksfor`, we discuss how the :term:`FKS` scheme :cite:`Frixione:1995ms, Frederix:2009yq, Ulrich:2019fks, Engel:2019nfw, Ulrich:2020phd`. 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 :ref:`sec_ps`. Next, in Section :ref:`sec_fksfor`, we discuss how the :term:`FKS` scheme :cite:`Frixione:1995ms, Frederix:2009yq, Ulrich:2019fks, Engel:2019nfw, Ulrich:2020phd` (cf. Appendix :ref:`sec_fks` 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 :ref:`sec_rng`. Finally, we give an account of how the statefiles work and how they are used to store distributions in Section :ref:`sec_vegasff`. .. _sec_particle_string: Particle string framework ------------------------- .. note:: This section needs to be completed, link to `issue `_ .. _sec_ps: Phase-space generation ---------------------- We use the ``vegas`` algorithm for numerical integration :cite:`Lepage:1980jk`. As ``vegas`` only works on the hypercube, we need a routine that maps :math:`[0,1]^{3n-4}` to the momenta of an :math:`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 :numref:`lst_psn`. .. literalinclude:: listings/psn.f90 :language: fortran :caption: 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 :math:`E_i \ge m_i`. :name: lst_psn As soon as we start using :term:`FKS`, we cannot use this simplistic approach any longer. The :math:`c`-distributions of :term:`FKS` require the photon energies :math:`\xi_i` to be variables of the integration. We can fix this by first generating the photon explicitly as .. math:: :label: eq_kdef k_1 = p_{n+1} = \frac{\sqrt s}2 \xi_1 (1,\sqrt{1-y_1 ^2}\vec e_\perp,y_1 ) where :math:`\vec e_\perp` is a :math:`(d-2)` dimensional unit vector and the ranges of :math:`y_1` (the cosine of the angle) and :math:`\xi_1` (the scaled energy) are :math:`-1\le y_1 \le 1` and :math:`0\le\xi_1 \le\xi_{\text{max}}`, respectively. The upper bound :math:`\xi_{\text{max}}` depends on the masses of the outgoing particles. Following :cite:`Frederix:2009yq` we find .. math:: \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 :term:`PCS`\ s 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 :math:`\mu\to\nu\bar\nu e\gamma`. The :term:`PCS` comes from .. math:: \mathcal{M}_{n+1}^{(\ell )} \propto \frac{1}{q\cdot k} = \frac1{\xi^2}\frac1{1-y\beta} where :math:`y` is the angle between photon (:math:`k`) and electron (:math:`q`). For large velocities :math:`\beta` (or equivalently small masses), this becomes almost singular as :math:`y\to1`. If now :math:`y` is a variable of the integration this can be mediated. An example implementation is shown in :numref:`lst_psmu`. .. literalinclude:: listings/psmu.f90 :language: fortran :caption: Example implementation of a so-called :term:`FKS` phase-space where the fifth particle is an :term:`FKS` photon that may becomes soft. Not shown are checks whether :math:`E_i \ge m_i` . :name: lst_psmu 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 :math:`\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 :math:`\mu`-:math:`e`-scattering .. math:: 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 :eq:`eq_kdef`. The outgoing electron :math:`p_3` is more complicated. Writing the :math:`p_4`-phase-space four- instead of three-dimensional .. math:: \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 :math:`\delta` function for :math:`p_4` and proceed for the generation :math:`p_3` and :math:`p_5` almost as for the muon decay above. Doing this we obtain for the final :math:`\delta` function .. math:: :label: eq_psdel \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 :math:`E_3`, we need to take care to avoid extraneous solutions of this radical equation :cite:`Gkioulekas:2018`. We have now obtained our phase-space parametrisation, albeit with one caveat: for anti-collinear photons, i.e. :math:`-1 s_{35} \big) + \theta\big( s_{15} < s_{35} \big) with :math:`s_{ij} = 2p_i\cdot p_j` as usual. The integrand of the first :math:`\theta` function has a final-state :term:`PCS` and hence we use the parametrisation obtained by solving :eq:`eq_psdel`. The second :math:`\theta` function, on the other hand, has an initial-state :term:`PCS` which can be treated by just directly parametrising the photon in the centre-of-mass frame as per :eq:`eq_kdef`. This automatically makes :math:`s_{15}\propto(1-\beta_{\text{in}}y_1)` a variable of the integration. For the double-real corrections of :math:`\mu`-:math:`e` scattering, we proceed along the same lines except now the argument of the :math:`\delta` function is more complicated. .. _sec_fksfor: Implementation of FKS schemes ----------------------------- Now that we have a phase-space routine that has :math:`\xi_i` as variables of the integration, we can start implementing the relevant :math:`c`-distributions :eq:`eq_xidist` .. math:: :nowrap: :label: eq_event_counterevent \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 :term:`event` and the second as the :term:`counter-event`. Note that, due to the presence of :math:`\delta(\xi_1)` in the counter-event (that is implemented through the eikonal factor :math:`\mathcal{E}`) the momenta generated by the phase-space :math:`\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 :term:`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 :math:`\xi` and sometimes also on a parameter that encodes the :term:`PCS` such as :math:`y={\tt y2}` in :eq:`eq_kdef` and :numref:`lst_psmu`. Events that have values of :math:`\xi` smaller than this :term:`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 :term:`soft cut` has to be chosen small enough to not influence the result. An example implementation can be found in :numref:`lst_fks`. .. literalinclude:: listings/fks.f90 :language: fortran :caption: An example implementation of the :term:`FKS` scheme in Fortran. Not shown are various checks performed, the binning as well as initialisation blocks. :name: lst_fks .. _sec_calling: 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 :f:subr:`init_piece` and then called throughout :f:mod:`integrands` and :f:mod:`phase_space`. The pointers for the phase-space generator and integrand are just assigned using the ``=>`` operator, i.e. .. code:: fortran ps => psx2 ; fxn => sigma_0 The relevant abstract interface for the integrand ``fxn`` is .. code:: fortran 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 .. code:: fortran 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:: Arguments for ``set_func`` :widths: auto +--------------+--------------+-------------------------------+ | 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``. .. _sec_integrandpara: Optional parameters for integrands ---------------------------------- The integration is configured during the :f:subr:`initpiece` routine. Additionally to identifying what is to be integrated (cf. Section :ref:`sec_calling`), one also configures other parameters such as the dimensionality or the masses involved. .. list-table:: Frozen Delights! :header-rows: 1 * - 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 :math:`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 :math:`\xi_c` for the first subtraction - for real corrections * - ``xicut2`` - ``real`` - the value of :math:`\xi_c` for the second subtraction - for double-real corrections * - ``xieik1`` - ``real`` - the value of :math:`\xi_c` for the first eikonal - for virtual or real-virtual corrections * - ``xieik2`` - ``real`` - the value of :math:`\xi_c` for the second eikonal - for double-virtual corrections * - ``polarised`` - ``integer`` - the number of polarised particles - no, defaults to 0 * - ``symmfac`` - ``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 :term:`NTS` switching point - only for NTS matrix elemnts :math:`\xi_c` parameters ~~~~~~~~~~~~~~~~~~~~~~~~ For the :math:`\xi_c` parameters, the user enters a value between zero (exclusive) and one (inclusive). However, the :term:`FKS` procedure requires the bounds of :eq:`eq_xic` 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 :math:`\xi` (:math:`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 :eq:`eq_event_counterevent` .. math:: :nowrap: \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. .. _sec_rng: Random number generation ------------------------ A Monte Carlo integrator relies on a (pseudo) *random number generator* (:term:`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 :term:`RNG` needs to be simple enough to be called many billion times without being a significant source of runtime. :term:`RNG`\ s 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 :term:`RNG`, also known as ``minstd``\ :cite:`Park:1988RNG`. As a linear congruential generator, its :math:`(k+1)`\ th output :math:`x_{k+1}` can be found as .. math:: 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 :math:`m` is a large, preferably prime, number and :math:`2` ``1B1C`` ``0014`` ``integer`` :math:`n_q` number of histograms ``integer`` :math:`n_b` number of bins ``integer`` :math:`n_s` len. histogram name ``1B30`` :math:`10n_q+8` ``real(n_q)`` ``minv`` lower bounds ``real(n_q)`` ``maxv`` upper bounds :math:`n_sn_q+8` ``character(n_s,n_q)`` ``names`` names of :math:`S` :math:`10n_q(n_b+2)+8` ``real(n_q,n_b+2)`` ``quantsum`` accu. histograms ``real(n_q,n_b+2)`` ``quantsumsq`` accu. histograms squared ``-0144`` ``0010`` ``real`` ``time`` current runtime in seconds ``-0134`` ``0134`` ``character(300)`` ``msg`` any message ``-0000`` ``EOF`` ========= ====================== ====================== ============== ==================================================== .. _sec_docker: Basics of containerisation -------------------------- McMule is Docker-compatible. Production runs should be performed with Docker :cite:`Merkel:2014`, or its user-space complement *udocker* :cite:`Gomes:2017hct`, 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 :term:`container` needs to be generated. Essentially, this involves uncompressing all layers into a directory an ``chroot``\ ing into said directory. It is important to note, that :term:`containers` are ephemeral, i.e. changes made to the :term:`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 :term:`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 .. code:: docker 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 .. code:: docker 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 .. code:: bash 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, :term:`containers` are usually created and run in one command .. code:: bash $ docker run --rm $imagename $cmd The flag ``–rm`` makes sure the :term:`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 .. code:: bash $ udocker create $imagename # this prints the container id $ udocker run $containerid $cmd # work in container $ udocker rm $containerid or in one step .. code:: bash $ udocker run --rm $imagename $cmd Running :term:`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 :math:`\xi=0`. .. [2] Note that, because of the simple recursion the :term:`RNG` will not repeat any number until the full period is complete .. [3] If :math:`p` is prime, for any integer :math:`a`, :math:`a^p-a` is a multiple of :math:`p`. .. [4] An infamous example is ``randu`` that used :math:`a=2^{16}+3` and :math:`m=2^{31}` that in three dimension produces only 15 planes instead of the maximum 2344. .. [5] To be precise, the actual dimensions are :math:`(n_b+2)\times n_q` to accommodate under- and overflow bins