.. _sec_implement:

Implementing new processes in McMule
====================================

In this section we will discuss how new processes can be added to McMule.
Not all of the points below might be applicable to any particular process.
Further, all points are merely guidelines that could be deviated from if necessary as long as proper precautions are taken.

As an example, we will discuss how Møller scattering :math:`e^-e^-\to e^-e^-` could be implemented.

1. A new :term:`process group` may need to be created if the process does not fit any of the presently implemented groups.
   This requires a new folder with a makefile as well as modifications to the main makefile as discussed in Section :ref:`sec_newpg`.

   In our case, :math:`ee\to ee` does not fit any of the groups, so we create a new group that we shall call ``ee``.

2. Calculate the tree-level matrix elements needed at :term:`LO` and :term:`NLO`: :math:`\mathcal{M}_{n}^{(0)}` and :math:`\mathcal{M}_{n+1}^{(0)}`.
   This is relatively straightforward and – crucially – unambiguous as both are finite in :math:`d=4`.
   We will come back to an example calculation in Section :ref:`sec_matel`.

3. A generic matrix element file is needed to store 'simple' matrix elements as well as importing more complicated matrix elements.
   Usually, this file should not contain matrix elements that are longer than a few dozen or so lines.
   In most cases, this applies to :math:`\mathcal{M}_{n}^{(0)}`.

   After each matrix element, the :term:`PID` needs to be denoted in a comment.
   Further, all required masses as well as the centre-of-mass energy, called :f:var:`scms` to avoid collisions with the function :math:`{\tt s(pi,pj)}=2{\tt pi}\cdot{\tt pj}`, need to be calculated in the matrix element to be as localised as possible.

   In the case of Møller scattering, a file ``ee/ee_mat_el.f95`` will contain :math:`\mathcal{M}_{n}^{(0)}`.
   For example, :math:`\mathcal{M}_{n}^{(0)}` is implemented there as shown in :numref:`lst_mollerlo`.

.. literalinclude:: listings/mollerlo.f90
   :language: fortran
   :caption: An example implementation of :math:`\mathcal{M}_n^{(0)}` for Møller scattering. Note that the electron mass and the centre-of-mass energy are calculated locally. A global factor of :math:`8e^4 = 128\pi^2\alpha^2` is included at the end.
   :name: lst_mollerlo

4.  Further, we need an interface file that implements the abstract interface :f:type:`partInterface` (cf. Section :ref:`sec_particle_string`).
    In our case this is called ``ee/ee.f95``, as shown in :numref:`lst_mollerpart`.

    In most cases, the :f:type:`particle` string framework can automatically define all soft limits and :term:`NTS` matrix elements.
    In rare cases, this needs to be done manually in ``ee/ee.f95``. [1]_

.. literalinclude:: listings/mollerpart.f90
   :language: fortran
   :caption: An example implementation of the soft limits for Møller scattering in the particle framework.
   :name: lst_mollerpart

5. Because :math:`\mathcal{M}_{n+1}^{(0)}` is border-line large, we will assume that it will be stored in an extra file, ``ee/ee2eeg.f95``.
   The required functions are to be imported in ``ee/ee_mat_el.f95``.

6. Calculate the one-loop virtual matrix element :math:`\mathcal{M}_{n}^{(1)}`, renormalised in the :term:`OS` scheme.
   Of course, this could be done in any regularisation scheme.
   However, results in McMule shall be in the FDH (or equivalently the FDF) scheme.
   Divergent matrix elements in McMule are implemented as :math:`c_{-1}`, :math:`c_0`, and :math:`c_1`

   .. math::

          \mathcal{M}_{n}^{(1)} = \frac{(4\pi)^\epsilon}{\Gamma(1-\epsilon)}\Bigg(
              \frac{c_{-1}}\epsilon + c_0 + c_1\epsilon+\mathcal{O}(\epsilon^2)
          \Bigg)\,.

   For :math:`c_{-1}` and :math:`c_0` this is equivalent to the conventions employed by Package-X :cite:`Patel:2015tea` up to a factor :math:`1/16\pi^2`.
   While not strictly necessary, it is generally advisable to also include :math:`c_{-1}` in the Fortran code.

   For :term:`NLO` calculations, :math:`c_1` does not enter.
   However, we wish to include Møller scattering up to :term:`NNLO` and hence will need it sooner rather than later anyway.

   In our case, we will create a file ``ee/ee_ee2eel.f95``, which defines a function

   .. code:: fortran

        FUNCTION EE2EEl(p1, p2, p3, p4, sing, lin)
          !! e-(p1) e-(p2) -> e-(p3) e-(p4)
          !! for massive electrons
        implicit none
        real(kind=prec), intent(in) :: p1(4), p2(4), p3(p4), p4(4)
        real(kind=prec) :: ee2eel
        real(kind=prec), intent(out), optional :: sing, lin
        ...
        END FUNCTION

   The function shall return :math:`c_0` in ``ee2eel`` and, if ``present`` :math:`c_{-1}` and :math:`c_1` in ``sing`` and ``lin``.

7. At this stage, a new subroutine in the program ``test`` with reference values for all three matrix elements should be written to test the Fortran implementation.
   This is done by generating a few points using an appropriate phase-space routine and comparing to as many digits as possible using the routine ``check``.

   In our case, we would construct a subroutine ``TESTEEMATEL`` as shown in :numref:`lst_mollertest`.

.. literalinclude:: listings/mollertest.f90
   :language: fortran
   :caption: Test routine for :math:`ee\to ee` matrix elements and integrands. The reference values for the integration are yet to be determined.
   :name: lst_mollertest

8. In addition, McMule provides built-in routines for testing the convergence of real-emission matrix elements to the corresponding soft limits, for ever smaller photon energies.

   In our case, we would construct a subroutine

   .. code:: fortran

	SUBROUTINE TESTEESOFTN1
	implicit none
	real(kind=prec) y0(5)

	call blockstart("e-e \xi->0")
	call initflavour("muone")
	xinormcut1 = 0.3
	y0 = (/0.01,0.6,0.8,0.999,0.01/)

	call test_softlimit(y0, ["ee2eeR"])
	END SUBROUTINE

   where :f:type:`test_soft_limit` compares the real matrix element (``ee2eeR``) with its soft limit implemented in ``ee/ee.f95``.
   The comparison starts at the energy defined by the phase-space point ``y0`` and proceeds with ever smaller photon energies.
   A flavour (``muone``) as well as ``xinormcut1`` are required in order to complete the phase-space generation.

9. Define a default observable in ``user`` for this process.
   This observable must be defined for any :f:var:`which_piece` that might have been defined and test all relevant features of the implementation such as polarisation if applicable.

10. Add the matrix elements to the integrands defined in ``integrands.f95``.
    This is done by adding a new ``case`` corresponding to the new :f:var:`which_piece` in the :f:subr:`initpiece`.

     - for a IR-finite, non-radiative piece (i.e. :term:`LO` but also :term:`VP`), add

        .. code:: fortran

            case('eb2eb0')
              call set_func('00000000', eb2eb)
              ps => psx2 ; fxn => sigma_0
              nparticle = 4 ; ndim = 2
              masses(1:4) = (/ Me, Me, Me, Me /)

       which adds a :f:var:`which_piece` ``ee2ee0`` that is calculated using the matrix element ``ee2ee``.
       The phase space is generated with :f:subr:`psx4` and integrated using :f:func:`sigma_0` (no subtraction).
       The process involves 4 particles and, since it is a :math:`2\to2` process, the integration is two-dimensional.
       The masses of the involved particles are all :f:var:`Me`.

     - for pieces with an IR cancellation between real and virtual corrections, i.e. calculations involving photon loops, we need to specify :f:var:`xieik1` (at one-loop) and/or :f:var:`xieik2` (at two-loop)

        .. code:: fortran

            case('eb2ebF')
              call set_func('00000000', eb2ebf)
              ps => psx2 ; fxn => sigma_0
              nparticle = 4 ; ndim = 2
              masses(1:4) = (/ Me, Me, Me, Me /)
              xieik1 = xinormcut*(1.-(2*me)**2/scms)

       One needs to take care that :math:`\xi_c` is properly normalised.
       The user will enter a value from 0 to 1 which needs to be matched to :math:`\xi_c` as defined in :eq:`eq_xic`

     - for real corrections we need to use a subtracting integrand, i.e. :f:func:`sigma_1` for single-real and :f:func:`sigma_2` for double-real corrections.

        .. code:: fortran

            case('eb2ebR')
              call set_func('00000000', eb2ebg)
              call set_func('00000001', eb2eb)
              call set_func('11111111', eb2eb_part)
              ps => psx3_fks ; fxn => sigma_1
              nparticle = 5 ; ndim = 5
              masses(1:5) = (/ Me, Me, Me, Me, 0._prec /)
              xicut1 = xinormcut*(1.-(2*me)**2/scms)


       Additionally to the real matrix element ``eb2ebg``, we also specified the reduced matrix element ``eb2eb`` and the :f:type:`particle` string function ``eb2eb_part``.
       Note further changes to the number of particles and phase space dimension to accommodate the extra photon.

    Additionally to these required parameters, there are number of optional parameters such as ``symmfac`` (which is set to 2 for :math:`e^-e^-\to e^-e^-` because of the two indistinguishable final state particles), ``polarised`` (whether to consider the process polarised), and ``softCut`` and ``collCut``.
    For a full list of parameters, see Section :ref:`sec_integrandpara`.

    Once ``integrands`` are implemented, a second test routine should be written that runs short integrations against a reference value.
    Because ``test_INT`` uses a fixed :term:`random seed`, this is expected to be possible very precisely.
    Unfortunately, COLLIER :cite:`Denner:2016kdg` might produce slightly different results on different machines.
    Hence, integrands involving complicated loop functions are only required to agree up to :math:`\mathcal{O}(10^{-8})`.

11. After some short test runs, it should be clear whether new phase-space routines are required.
    Add those, if need be, to ``phase_space`` as described in Section :ref:`sec_ps`.

12. Per default the stringent :term:`soft cut`, that may be required to stabilise the numerical integration (cf.  Section :ref:`sec_fksfor`), is set to zero.
    Study what the smallest value is that still permits integration.

13. Perform very precise :math:`\xi_{c}` independence studies.
    Tips on how to do this can be found in Section :ref:`sec_xicut`.

At this stage, the :term:`NLO` calculation is complete and may, after proper integration into McMule and adherence to coding style has been confirmed, be added to the list of McMule processes in a new release.
Should :term:`NNLO` precision be required, the following steps should be taken

14. Calculate the real-virtual and double-real matrix elements :math:`\mathcal{M}_{n+1}^{(1)}` and :math:`\mathcal{M}_{n+2}^{(0)}` and add them to the test routines as well as integrands.

15. Prepare the :math:`n`-particle contribution :math:`\sigma_n^{(2)}`.
    In a pinch, massified results can be used also for :math:`\hat{\mathcal{E}}(\xi_{c})\mathcal{M}_{n}^{(1)}` though of course one should default to the fully massive results.

16. Study whether the pre-defined phase-space routines are sufficient.
    Even if it was possible to use an old phase-space at :term:`NLO`, this might no longer work at :term:`NNLO` due to the added complexity.
    Adapt and partition further if necessary, adding more test integrations in the process.

17. Perform yet more detailed :math:`\xi_{c}` and :term:`soft cut` analyses.

In the following we comment on a few aspects of this procedure such as the :math:`\xi_{c}` study (Section :ref:`sec_xicut`), the calculation of matrix elements (Section :ref:`sec_matel`), and a brief style guide for McMule code (Section :ref:`sec_style`).

.. _sec_newpg:

Creating a new process group
----------------------------

Adding Møller scattering to McMule, the example discussed above, requires the addition of a new :term:`process group` ``ee``.
For this we create a new folder in McMule called ``ee`` containing a meson.build (:numref:`lst_eemeson`), a ``mat_el`` file (``ee/ee_mat_el.f95``, :numref:`lst_eematel`) and a module file (``ee/ee.f95``, :numref:`lst_eemod`).
Finally, the name of the :term:`process group` needs to be added to the ``GROUPS`` and ``WGROUPS`` variables of the makefile.

.. literalinclude:: listings/eemeson.build
   :language: meson
   :caption: The bare meson.build for the new :term:`process group` ``ee``. Large matrix elements are stored in extra files such as ``ee/ee2eeg.f95`` or ``ee/ee_ee2eel.f95``.
   :name: lst_eemeson

.. literalinclude:: listings/eematel.f90
   :language: fortran
   :caption: The file ``ee/ee_mat_el.f95`` imports the complicated matrix elements ``ee2eel`` and ``ee2eegl``, defines the simple matrix element ``ee2ee`` as per 16, and provides an interface for the :math:`\mathcal{M}_n^{(1)f}` that is called from `integrands`.
   :name: lst_eematel

.. literalinclude:: listings/eemod.f90
   :language: fortran
   :caption: The module file ``ee/ee.f95`` which imports all matrix elements of ``ee_mat_el`` and defines the soft limits.
   :name: lst_eemod

.. _sec_xicut:

Study of :math:`\xi_{c}` dependence
-----------------------------------

When performing calculations with McMule, we need to check that the dependence of the unphysical :math:`\xi_{c}` parameter introduced in the :term:`FKS` scheme (cf. Appendix :ref:`sec_fks`) actually drops out at :term:`NLO` and :term:`NNLO`.
In principle it is sufficient to do this once during the development phase.
However, we consider it good practice to also do this (albeit with a reduced range of :math:`\xi_{c}`) for production runs.

Because the :math:`\xi_{c}` dependence is induced through terms as :math:`\xi_{c}^{-2\epsilon}/\epsilon`, we know the functional dependence of :math:`\sigma^{(\ell)}_{n+j}`.
For example, at :term:`NLO` we have

.. math::
   :nowrap:

    \begin{align}
       \sigma^{(1)}_{n  }(\xi_{c}) &= a_{0,0} + a_{0,1}\log(\xi_{c})\\
       \sigma^{(1)}_{n+1}(\xi_{c}) &= a_{1,0} + a_{1,1}\log(\xi_{c})
    \end{align}

where :math:`\xi_{c}` independence of :math:`\sigma^{(1)}` of course requires

.. math::
   :label: eq_xinlo

    a_{0,1}+ a_{1,1} = 0

At :term:`NNLO` we have

.. math::
   :nowrap:

   \begin{align}
   \sigma^{(2)}_{n  }(\xi_{c}) &= a_{0,0} + a_{0,1}\log(\xi_{c}) + a_{0,2}\log(\xi_{c})^2\,,\\
   \sigma^{(2)}_{n+1}(\xi_{c}) &= a_{1,0} + a_{1,1}\log(\xi_{c}) + a_{1,2}\log(\xi_{c})^2\,,\\
   \sigma^{(2)}_{n+2}(\xi_{c}) &= a_{2,0} + a_{2,1}\log(\xi_{c}) + a_{2,2}\log(\xi_{c})^2\,.
   \end{align}

We require

   .. math::

      a_{0,i} + a_{1,i} + a_{2,i} = 0

for :math:`i=1,2`.
However, the :term:`IR` structure allows for an even stronger statement for the :math:`a_{j,2}` terms

   .. math::

      a_{0,2} = a_{2,2} = -\frac{a_{1,2}}2\,.

Of course we cannot directly calculate any of the :math:`a_{1,i}` or :math:`a_{2,i}` because we use numerical integration to obtain the :math:`\sigma^{(\ell)}_{n+j}`.
Still, knowing the coefficients can be extremely helpful when debugging the code or to just quantify how well the :math:`\xi_{c}` dependence vanishes.
Hence, we use a fitting routine to fit the Monte Carlo results *after* any phase-space partitioning has been undone.
Sometimes non of this is sufficient to pin-point the source of a problem to any one integrand.
However, if the goodness of, for example, :math:`\sigma^{(2)}_{n+2}(\xi_{c})` is much worse than the one for :math:`\sigma^{(2)}_{n+1}(\xi_{c})`, a problem in the double-real corrections can be expected.

.. _sec_matel:

Example calculations in Mathematica
-----------------------------------

A thorough understanding of one-loop matrix elements is crucial for any higher-order calculation.
In McMule, one-loop matrix elements either enter as the virtual contribution to :term:`NLO` corrections or the real-virtual contribution in :term:`NNLO` calculations.
In any case, a fast numerical routine is required that computes the matrix element.

We perform all one-loop calculations in FDF as this is arguably the simplest scheme available.
For theoretical background, we refer to :cite:`Ulrich:2020phd` and references therein.

We use Qgraf for the diagram generation.
Using the in-house Mathematica package ``qgraf`` we convert Qgraf’s output for manipulation with Package-X :cite:`Patel:2015tea`.
This package is available on request through the McMule collaboration

   https://gitlab.com/mule-tools/qgraf

An example calculation for the one-loop calculation of :math:`\mu\to\nu\bar\nu e\gamma` can be found in :numref:`lst_pkgx`.
Of course this example can be made more efficient by, for example, feeding the minimal amount of algebra to the loop integration routine.

When using ``qgraf`` for fdf some attention needs to be paid when considering diagrams with closed fermion loops.
By default, ``qgraf.wl`` evaluates these traces in :math:`d` dimensions.
``RunQGraf`` has an option to keep this from happening.

.. literalinclude:: listings/px.m
   :language: wl
   :caption: An example on how to calculate the renormalised one-loop matrix element for :math:`\mu\to\nu\bar\nu e` in fdf.
   :name: lst_pkgx

There is a subtlety here that only arise for complicated matrix elements.
Because the function Package-X uses for box integrals, :f:func:`ScalarD0IR6`, is so complicated, no native Fortran implementation exists in McMule.
Instead, we are defaulting to COLLIER :cite:`Denner:2016kdg` and should directly evaluate the finite part of the ``PVD`` function above.
The same holds true for the more complicated triangle functions.
In fact, only the simple :f:func:`DiscB` and :f:func:`ScalarC0IR6` are natively implemented without need for external libraries.
For any other functions, a judgement call is necessary of whether one should ``LoopRefine`` the finite part in the first place.
In general, if an integral can be written through logarithms and dilogs of simple arguments (resulting in real answers) or :f:func:`DiscB` and :f:func:`ScalarC0IR6`, it makes sense to do so.
Otherwise, it is often easier to directly link to COLLIER.

.. _sec_style:

Coding style and best practice
------------------------------

A large-scale code base like McMule cannot live without some basic agreements regarding coding style and operational best practice.
These range from a (recommended but not enforced) style guide over the management of the git repository to how to best run McMule in development scenarios.
All aspects have been discussed within the McMule collaboration.

Fortran code in McMule is (mostly) written in accordance with the following style guide.
If new code is added, compliance would be appreciated but deviation is allowed if necessary.
If in doubt, contact any member of the McMule collaboration.

-  Indentation width is two spaces.
   In Vim this could be implemented by adding the following to ``.vimrc``

   .. code:: vim

      autocmd FileType fortran set tabstop=8 softtabstop=0 expandtab shiftwidth=2 smarttab

-  Function and subroutine names are in all-upper case.

-  A function body is not indented beyond its definition.

-  When specifying floating point literals specify the precision when possible, i.e. ``1._prec``.

-  Integrands should have ``ndim`` specified.

-  Internal functions should be used where available.

-  Masses and other kinematic parameters must be calculated in the matrix elements as local variables; using the global parameters ``Mm`` and ``Me`` is strictly forbidden.

-  These rules also hold for matrix elements.

For python code, i.e. pymule as well as the analysis code, PEP8 compliance is strongly encouraged with the exception of ``E231`` (Missing whitespace after ``,``, ``;``, and ``:``), ``E731`` (Do not assign a lambda expression, use a ``def``) as well, in justified cases, i.e. if required by the visual layout, ``E272`` (Multiple spaces before keyword), and ``E131`` (Continuation line unaligned for hanging indent).

McMule uses a git repositories for version management.
Development usually happens on feature branches that are merged into the ``devel`` branch semi-frequently by the McMule collaboration after sufficient vetting was performed.
Finally, once a project has been finished, the ``devel`` branch gets merged into the ``release`` branch that is to be used by McMule's users.

In general, developers are encouraged to not commit wrong or unvetted code though this can obviously not be completely avoided in practice.
To avoid uncontrollable growth of the git repository, large files movements are strongly discouraged.
This also means that matrix elements should not be completely overhauled barring unanimous agreement.
Instead, developers are encouraged to add a new matrix element file and link to that instead.

Even when running McMule for development purposes the usage of :term:`menu files<menu file>` is strongly encouraged because the code will do its utmost to automatically document the run by storing the git version as well as any modification thereof.
This allows for easy and unique reconstruction of what was running.
For production runs this is not optional; these must be conducted with :term:`menu files<menu file>` after which the run folder must be stored with an analysis script and all data on the AFS as well as the user file library to ensure data retention.

.. [1]
    Further coding may be required if the user needs to isolate different gauge-invariant contributions to the process.
    For example, for :math:`e\mu\to e\mu` scattering, the function :f:type:`em2em_ee_part` = ``parts((/part(p1, 1, 1), part(p2, 1, 1), part(p3, 1, -1), part(p4, 1, -1)/), "e")`` can be used to generate all soft limits due to emissions from the electron line only.
    Similarly, the function :f:type:`em2em_mm_part` can be used for all soft limits from the muon line only.
    However, the function :f:type:`em2em_em_part` = ``parts((/part(p1, 1, 1), part(p2, 1, 1), part(p3, 1, -1), part(p4, 1, -1)/), "x")`` generates only a subset of all soft contributions due to emissions from both lepton lines.
    Thus, for this mixed case, the required soft limits are hard-coded in ``mue/mue.f95``