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 \(e^-e^-\to e^-e^-\) could be implemented.

1. A new 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 Creating a new process group.

   In our case, \(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 LO and NLO: \(\mathcal{M}_{n}^{(0)}\) and \(\mathcal{M}_{n+1}^{(0)}\). This is relatively straightforward and – crucially – unambiguous as both are finite in \(d=4\). We will come back to an example calculation in Section Example calculations in Mathematica.

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 \(\mathcal{M}_{n}^{(0)}\).

   After each matrix element, the PID needs to be denoted in a comment. Further, all required masses as well as the centre-of-mass energy, called scms to avoid collisions with the function \({\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 \(\mathcal{M}_{n}^{(0)}\). For example, \(\mathcal{M}_{n}^{(0)}\) is implemented there as shown in Listing 13.

Listing 13 An example implementation of \(\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 \(8e^4 = 128\pi^2\alpha^2\) is included at the end.

FUNCTION EE2EE(p1, p2, p3, p4)
  !! e-(p1) e-(p2) -> e-(p3) e-(p4)
  !! for massive (and massless) electrons
implicit none
real(kind=prec), intent(in) :: p1(4), p2(4), p3(p4), p4(4), ee2ee
real(kind=prec) :: den1, den2, t, scms, m2
t = sq(p1-p3) ; scms = sq(p1+p2) ; m2 = sq(p1)
den1 = sq(p1-p3) ; den2 = sq(p1-p4)

ee2ee=(8**m2**2 - 8*m2*scms + 2*s**2 + 2*scms*t + t**2)/den1**2
ee2ee=ee2ee+2*(12*m2**2 - 8*m2*scms + scms**2) / den1 / den2
ee2ee=ee2ee+(24*m2**2 + scms**2 + t**2 - 8*m2*(s + t))/den2**2

ee2ee = ee2ee * 128*pi**2*alpha**2
END FUNCTION

4. Further, we need an interface file that also contains the soft limits. In our case this is called ee/ee.f95.

   The abstract interface partInterface (cf. Section Technical routines) can take care of the generation of all soft limits for a given particle string, as shown in Listing 14. [1] See also Section Calling procedures and function pointers for more details.

Listing 14 An example implementation of the soft limits for Møller scattering in the particle framework.

FUNCTION EE2EE_part(p1, p2, p3, p4)
  !! e-(p1) e-(p2) -> e-(p3) e-(p4)
  !! for massive (and massless) electrons
implicit none
real(kind=prec) :: p1(4), p2(4), p3(p4), p4(4)
type(particles) :: ee2ee_part

ee2ee_part = parts((/part(p1, 1, 1), part(p2, 1, 1), part(p3, 1, -1), part(p4, 1, -1)/))
END FUNCTION

5. Because \(\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 \(\mathcal{M}_{n}^{(1)}\), renormalised in the 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 \(c_{-1}\), \(c_0\), and \(c_1\)

   \[\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 \(c_{-1}\) and \(c_0\) this is equivalent to the conventions employed by Package-X [20] up to a factor \(1/16\pi^2\). While not strictly necessary, it is generally advisable to also include \(c_{-1}\) in the Fortran code.

   For NLO calculations, \(c_1\) does not enter. However, we wish to include Møller scattering up to 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

   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 \(c_0\) in ee2eel and, if present \(c_{-1}\) and \(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 Listing 15.

Listing 15 Test routine for \(ee\to ee\) matrix elements and integrands. The reference values for the integration are yet to be determined.

SUBROUTINE TESTEEMATEL
implicit none
real (kind=prec) :: x(2),y(5)
real (kind=prec) :: single, finite, lin
real (kind=prec) :: weight
integer ido

call blockstart("ee matrix elements")
scms = 40000.
musq = me
x = (/0.75,0.5/)
call ps_x2(x,scms,p1,me,p2,me,p3,me,p4,me,weight)
call check("ee2ee" ,ee2ee (p1,p2,p3,p4), 2.273983244890001e4, threshold=2e-8)
call check("ee2eel",ee2eel(p1,p2,p3,p4), 6.964297070440638e7, threshold=2e-8)

scms = 40000.
y = (/0.3,0.6,0.8,0.4,0.9/)
call ps_x3_fks(y,scms,p1,me,p2,me,p3,me,p4,me,p5,weight)
call check("ee2eeg",ee2eeg(p1,p2,p3,p4,p5),7.864297444955537e2, threshold=2e-8)

call blockend(3)
END SUBROUTINE

SUBROUTINE TESTMEEVEGAS
xinormcut1 = 0.2
xinormcut2 = 0.3

call blockstart("Moller VEGAS test")

call test_INT('ee2ee0', sigma_0, 2,10,10, NaN)
call test_INT('ee2eeF', sigma_0, 2,10,10, NaN)
call test_INT('ee2eeR', sigma_1, 5,10,10, NaN)
call blockend(3)
END SUBROUTINE

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

   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 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 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 which_piece in the initpiece().

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

      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 which_piece ee2ee0 that is calculated using the matrix element ee2ee. The phase space is generated with psx4() and integrated using sigma_0() (no subtraction). The process involves 4 particles and, since it is a \(2\to2\) process, the integration is two-dimensional. The masses of the involved particles are all Me.

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

      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 \(\xi_c\) is properly normalised. The user will enter a value from 0 to 1 which needs to be matched to \(\xi_c\) as defined in (12)

    • for real corrections we need to use a subtracting integrand, i.e. sigma_1() for single-real and sigma_2() for double-real corrections.

      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 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 \(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 Optional parameters for integrands.

    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 random seed, this is expected to be possible very precisely. Unfortunately, COLLIER [6] might produce slightly different results on different machines. Hence, integrands involving complicated loop functions are only required to agree up to \(\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 Phase-space generation.

12. Per default the stringent soft cut, that may be required to stabilise the numerical integration (cf. Section Implementation of FKS schemes), is set to zero. Study what the smallest value is that still permits integration.

13. Perform very precise \(\xi_{c}\) independence studies. Tips on how to do this can be found in Section Study of \xi_{c} dependence.

At this stage, the 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 NNLO precision be required, the following steps should be taken

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

15. Prepare the \(n\)-particle contribution \(\sigma_n^{(2)}\). In a pinch, massified results can be used also for \(\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 NLO, this might no longer work at NNLO due to the added complexity. Adapt and partition further if necessary, adding more test integrations in the process.

17. Perform yet more detailed \(\xi_{c}\) and soft cut analyses.

In the following we comment on a few aspects of this procedure such as the \(\xi_{c}\) study (Section Study of \xi_{c} dependence), the calculation of matrix elements (Section Example calculations in Mathematica), and a brief style guide for McMule code (Section Coding style and best practice).

Creating a new process group

Adding Møller scattering to McMule, the example discussed above, requires the addition of a new process group ee. For this we create a new folder in McMule called ee containing a makefile (Listing 16), a mat_el file (ee/ee_mat_el.f95, Listing 17) and a module file (ee/ee.f95, Listing 18). Finally, the name of the process group needs to be added to the GROUPS and WGROUPS variables of the makefile.

Listing 16 The bare makefile for the new process group ee. Large matrix elements that are stored in extra files such as ee/ee2eeg.f95 or ee/ee_ee2eel.f95 need to be added to the list of AUXFILES

group=ee
AUXFILES=ee_ee2eel.f95 ee_ee2eeg.f95

MAIN=$(group)_mat_el.f95 $(group).f95


include ../makefile.conf


all: $(group).a $(group).mod .obj/tree.sha

$(OBJ): ../.obj/functions.mod

.obj/$(group)_mat_el.o .obj/$(group)_mat_el.mod: \
    $(group)_mat_el.f95 $(MOD)
.obj/$(group).o .obj/$(group).mod: \
    $(group).f95 .obj/$(group)_mat_el.mod $(MOD)

$(group).mod: .obj/$(group).mod
    /+cp+/ $< $@

$(group).a:$(OBJ)
    @/+echo AR+/ $@
    @$(AR) $@ $^

clean:
	rm -f .obj/*.o .obj/*.mod .obj/*.gcda .obj/*.gcno *.mod *.a

Listing 17 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 \(\mathcal{M}_n^{(1)f}\) that is called from integrands.

                 !!!!!!!!!!!!!!!!!!!!!!!!!
                       MODULE EE_MAT_EL
                 !!!!!!!!!!!!!!!!!!!!!!!!!

  use functions
  use ee_ee2eel, only: ee2eel
  use ee_ee2eeg, only: ee2eeg
  implicit none

  contains

  FUNCTION EE2EE(p1,p2,q1,q2)
    !! e-(p1) e-(p2) -> e-(q1) e-(q2)
    !! for massive electrons
  ...
  END FUNCTION EE2EE

  FUNCTION EE2EEF(p1,p2,q1,q2)
    !! e-(p1) e-(p2) -> e-(q1) e-(q2)
    !! massive electrons
  real(kind=prec) :: p1(4),p2(4),q1(4),q2(4)
  real(kind=prec) :: ee2eef, mat0, Epart

  Epart = sqrt(sq(p1+p2))
  mat0 = ee2ee(p1,p2,q1,q2)

  ee2eef = ee2eel(p1,p2,q1,q2) + alpha / (2 * pi) * mat0 * (&
      - Ieik(xieik1,Epart,p1,p2) + Ieik(xieik1,Epart,p1,q1) &
      + Ieik(xieik1,Epart,p1,q2) + Ieik(xieik1,Epart,p2,q1) &
      + Ieik(xieik1,Epart,p2,q2) - Ieik(xieik1,Epart,q1,q2))

  END FUNCTION EE2EEF
                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
                       END MODULE EE_MAT_EL
                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!

Listing 18 The module file ee/ee.f95 which imports all matrix elements of ee_mat_el and defines the soft limits.

                 !!!!!!!!!!!!!!!!!!!!!!!!!
                       MODULE EE
                 !!!!!!!!!!!!!!!!!!!!!!!!!

  use functions
  use phase_space, only: ksoft, ksoftA, ksoftB
  use ee_mat_el
  implicit none
  contains

  FUNCTION EE2EE_part(p1, p2, p3, p4)
    !! e-(p1) e-(p2) --> e-(p3) e-(p4)
    !! both massive and massless electrons
  real (kind=prec) :: p1(4),p2(4),p3(4),p4(4)
  type(particles) :: ee2ee_part
  ee2ee_part = parts((/part(p1, 1, 1), part(p2, 1, 1), part(p3, 1, -1), part(p4, 1, -1)/))
  END FUNCTION EE2EE_part

                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!
                       END MODULE EE
                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!

Study of \(\xi_{c}\) dependence

When performing calculations with McMule, we need to check that the dependence of the unphysical \(\xi_{c}\) parameter introduced in the FKS scheme (cf. Appendix The FKS^2 scheme) actually drops out at NLO and 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 \(\xi_{c}\)) for production runs.

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

\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 \(\xi_{c}\) independence of \(\sigma^{(1)}\) of course requires

(7)\[ a_{0,1}+ a_{1,1} = 0\]

At NNLO we have

\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

\[a_{0,i} + a_{1,i} + a_{2,i} = 0\]

for \(i=1,2\). However, the IR structure allows for an even stronger statement for the \(a_{j,2}\) terms

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

Of course we cannot directly calculate any of the \(a_{1,i}\) or \(a_{2,i}\) because we use numerical integration to obtain the \(\sigma^{(\ell)}_{n+j}\). Still, knowing the coefficients can be extremely helpful when debugging the code or to just quantify how well the \(\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, \(\sigma^{(2)}_{n+2}(\xi_{c})\) is much worse than the one for \(\sigma^{(2)}_{n+1}(\xi_{c})\), a problem in the double-real corrections can be expected.

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 NLO corrections or the real-virtual contribution in 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 [25] 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 [20]. 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 \(\mu\to\nu\bar\nu e\gamma\) can be found in Listing 19. 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 \(d\) dimensions. RunQGraf has an option to keep this from happening.

Listing 19 An example on how to calculate the renormalised one-loop matrix element for \(\mu\to\nu\bar\nu e\) in fdf.

<<qgraf.wl
onshell = {
  p.p -> M^2, q.q -> m^2, p.q -> s/2
};

A0 = (4GF/Sqrt[2]) "diag1"/.RunQGraf[{"mum"},{"nu","elm"},0] //. {
  line[_, x_] -> x, p1->p, q1->p-q, q2->q, _δZ | δm -> 0
};
A1 = pref /. RunQGraf[{"mum"},{"nu","elm"},1] //. {
  line[_, x_] -> x, p1->p, q1->p-q, q2->q, _δZ | δm -> 0
};

M0=Block[{Dim=4},Simplify[Contract[
  1/2 Z2[m] Z2[M] FermionSpinSum[
    A0 /. DiracPL -> (Dirac1 - Z5 γ5)/2,
    A0 /. DiracPL -> (Dirac1 + Z5 γ5)/2
  ]
]] /. onshell]/.{
  Z2[M_] -> 1 + (α/(4π)) (-3/(2ε)-5/2 + 3/2 Log[M^2/Mu^2]),
  Z5     -> 1 - (α/(4π))
};
M1=Block[{Dim=4},Simplify[Contract[
  1/2 FermionSpinSum[
    A1/.γ.k1 -> γ. 4[k1]+I γ5 μ,
    A0
  ]
] /. onshell /. {
  μ^n_ /; EvenQ[n] -> μ2^(n/2), μ -> 0
}/.{
  4[k1]. 4[k1] -> k1.k1 + μ2, 4[k1] -> k1
}]];

M1bare = Simplify[KallenExpand[LoopRefine[LoopRelease[
  Pro2LoopIntegrate[
    Coefficient[M1, μ2, 0]/(16 π^2)
  ]
  + μIntegrate[
    Coefficient[M1, μ2, 1]/(64 π^3),
    1
  ],
  onshell
]]] /. e -> Sqrt[4 πα]];

There is a subtlety here that only arise for complicated matrix elements. Because the function Package-X uses for box integrals, ScalarD0IR6(), is so complicated, no native Fortran implementation exists in McMule. Instead, we are defaulting to COLLIER [6] 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 DiscB() and 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 DiscB() and ScalarC0IR6(), it makes sense to do so. Otherwise, it is often easier to directly link to COLLIER.

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

  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 menu files 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 menu files 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 \(e\mu\to e\mu\) scattering, the function 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 em2em_mm_part can be used for all soft limits from the muon line only. However, the function 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