Fortran reference guide

McMule’s Fortran code has hundreds of functions and subroutine and we will not document all of them here. However, we will list the user-facing function that are intended to help construct user files.

User-modifiable parameters

The following parameters may be modified by the user though it might become necessary to completely recompile McMule ones done.

type  real (kind=prec) [fixed]

The real number type used in McMule. This cannot be changed at runtime by the user but should be used for all interactions with the code. It usually refers to double precision

real(kind=prec) musq

The renormalisation scale \(\mu^2\). This variable needs to be set by the user, otherwise McMule will fail.

integer nel

Set to 1 if electron VP loops are to be included, set to 0 otherwise. More options may be added later

integer nmu

Set to 1 if muon VP loops are to be included, set to 0 otherwise. More options may be added later

integer ntau

Set to 1 if tau VP loops are to be included, set to 0 otherwise. More options may be added later

integer nhad

Set to 1 if HVP loops are to be included, set to 0 otherwise. More options may be added later

real (kind=prec) pol1(4)

The polarisation of the first polarised particle

real (kind=prec) pol2(4)

The polarisation of the second polarised particle

function  real(kind=prec) sachs_gel(q2)

The electric Sachs form factor of the proton. In the dipole approximation this is

\[G_e(Q^2) = \frac1{(1+Q^2/\Lambda)^2}\]

Parameters:

q2 [real(kind=prec)] :: the value of \(Q^2\)

function  real(kind=prec) sachs_gmag(q2)

The magnetic Sachs form factor of the proton. In the dipole approximation this is

\[G_m(Q^2) = \frac{\kappa}{(1+Q^2/\Lambda)^2}\]

Parameters:

q2 [real(kind=prec)] :: the value of \(Q^2\)

subroutine  init_flavour(flavour)

The definitions of the flavour. Users may edit this to add new experiments etc.

real(kind=prec) GF

The Fermi constant. For predominantly historic reasons, this is set to 1._prec.

real(kind=prec) alpha

The fine-structure constant in the OS scheme. For predominantly historic reasons, this is set to 1._prec.

real(kind=prec) sw2

The weak mixing angle \(\sin(\theta_W)^2\). This can be changed by the user at runtime to modify the EW scheme that is used.

real(kind=prec) Mel

The numerical value of the electron mass in MeV, irregardless of the flavour

real(kind=prec) Mmu

The numerical value of the muon mass in MeV, irregardless of the flavour

real(kind=prec) Mtau

The numerical value of the tau mass in MeV, irregardless of the flavour

real(kind=prec) Mproton

The numerical value of the proton mass in MeV, irregardless of the flavour

real(kind=prec) MZ

The numerical value of the Z boson mass in MeV

real(kind=prec) Mm

The actual value of the m particle, usually the muon mass but if flavour is eg. tau-e the tau mass

real(kind=prec) Me

The actual value of the e particle, usually the electron mass but if flavour is eg. tau-mu the muon mass

real(kind=prec) Mt

The actual value of the t particle, usually the tau mass

real(kind=prec) scms

The numerical value of the centre-of-mass energy

real(kind=prec) lambda

The dipole coefficient in the Sachs form factors of the proton in MeV:sup:2

real(kind=prec) kappa

The magnetic moment of the proton in Sachs form factors

Technical parameters

The following parameters should not be modified by the user unless especially advised to do so

character which_piece(25)

The piece being integrated, cf. Section Available processes and which_piece

character flavour(15)

The flavour configuration being used

real(kind=prec) softcut

The value of \(\xi\) below which the integrand is set to zero without subtraction

real(kind=prec) colcut

The value of \(\cos\theta\) below which the integrand is set to zero

real(kind=prec) sSwitch

The value of \(\xi\) below which the matrix element is approximated at LP. This is only available for some matrix elements

real(kind=prec) ntsSwitch

The value of \(\xi\) below which the matrix element is approximated at NLP. This is only available for some matrix elements

User-facing functions

The following function are available for the user to construct observables. Momenta are of the form (/ px, py, pz, E /).

Scalar quantities

function  real(kind=prec) s(p1, p2)

The scalar product \(2p_1\cdot p_2\)

Parameters:

• p1 (4) [real(kind=prec)] :: the first momentum \(p_1\)

• p2 (4) [real(kind=prec)] :: the second momentum \(p_2\)

function  real(kind=prec) sq(p)

The Lorentz square \(p^2\)

Parameters:

p (4) [real(kind=prec)] :: the momentum \(p\)

function  real(kind=prec) asymtensor(p1, p2, p3, p4)

The total asymmetric tensor \(\varepsilon_{\mu\nu\rho\sigma} p_1^\mu p_2^\nu p_3^\rho p_4^\sigma\)

Parameters:

• p1 (4) [real(kind=prec)] :: the momentum \(p_1\)

• p2 (4) [real(kind=prec)] :: the momentum \(p_2\)

• p3 (4) [real(kind=prec)] :: the momentum \(p_3\)

• p4 (4) [real(kind=prec)] :: the momentum \(p_4\)

function  real(kind=prec) eta(p)

The pseudorapidity w.r.t. the \(z\) axis

\[\eta = \frac12\log\frac{|\vec p| + p_z}{|\vec p| - p_z}\]

Parameters:

p (4) [real(kind=prec)] :: the momentum \(p\)

function  real(kind=prec) rap(p)

The rapidity w.r.t. the \(z\) axis

\[y = \frac12\log\frac{E + p_z}{E - p_z}\]

Parameters:

p (4) [real(kind=prec)] :: the momentum \(p\)

function  real(kind=prec) pt(p)

The transverse momentum w.r.t. the \(z\) axis

\[p_T = \sqrt{p_x^2+p_z^2}\]

Parameters:

p (4) [real(kind=prec)] :: the momentum \(p\)

function  real(kind=prec) absvec(p)

The length of the three-vector part \(|\vec p\)

Parameters:

p (4) [real(kind=prec)] :: the momentum \(p\)

function  real(kind=prec) phi(p)

The azimuthal angle of \(p\), \(-\pi<\phi<\pi\)

Parameters:

p (4) [real(kind=prec)] :: the momentum \(p\)

Note

This may return \(100\pi\) if the calculation fails.

function  real(kind=prec) rij(p1, p2)

The jet distance \(R_{12}\) between the two momenta \(p_1\) and \(p_2\), normalised by \(D_{\text{res}} = 0.7\)

\[R_{12} = \frac{\Delta y_{12}^2 + \Delta\phi_{12}^2}{D_\text{res}^2}\]

Parameters:

• p1 (4) [real(kind=prec)] :: the first momentum \(p_1\)

• p2 (4) [real(kind=prec)] :: the second momentum \(p_2\)

function  real(kind=prec) cos_th(p1, p2)

The cosine of the angle between the two momenta \(p_1\) and \(p_2\)

\[\cos\theta_{12} = \frac{\vec p_1\cdot \vec p_2}{|\vec p_1|\ |\vec p_2|}\]

Parameters:

• p1 (4) [real(kind=prec)] :: the first momentum \(p_1\)

• p2 (4) [real(kind=prec)] :: the second momentum \(p_2\)

Note

This will return 0 if the computation fails

Transformations

function  boost_back(rec, mo)

boosts the momentum mo from the frame where rec is at rest to the frame where rec is specified, i.e.

boost_back(rec, (/ 0., 0., 0., sqrt(sq(rec)) /)) = rec

This function can be viewed as the inversion of boost_rf().

Parameters:

• rec (4) [real(kind=prec)] :: the system to boost into

• mo (4) [real(kind=prec)] :: the momentum to boost

Return:

boost_back (4) [real(kind=prec)] :: the boosted momentum

function  boost_rf(rec, mo)

boosts mo to (non-unique) rest frame of rec, i.e.

boost_rf(rec, rec) = (/ 0., 0., 0., sqrt(sq(rec)) /)

This function can be viewed as the inversion of boost_back().

Parameters:

• rec (4) [real(kind=prec)] :: the system to boost into

• mo (4) [real(kind=prec)] :: the momentum to boost

Return:

boost_back (4) [real(kind=prec)] :: the boosted momentum

function  euler_mat(a, b, c)

gives the Euler rotation matrix formed by rotation by \(\alpha\) around the current \(z\) axis, then by \(\beta\) around the current \(y\) axis, and the by \(\gamma\) around the current \(z\) axis.

\[\begin{split}\begin{pmatrix} c_{\alpha}c_{\beta }c_{\gamma}-s_{\alpha}s_{\gamma}&-c_{\alpha}c_{\beta } s_{\gamma}-c_{\gamma}s_{\alpha}&c_{\alpha}s_{\beta }&0\\ c_{\beta }c_{\gamma}s_{\alpha}+c_{\alpha}s_{\gamma}& c_{\alpha}c_{\gamma}-c_{\beta } s_{\alpha}s_{\gamma}&s_{\alpha}s_{\beta }&0\\ -c_{\gamma}s_{\beta } & s_{\beta }s_{\gamma} &c_{\beta } &0\\ 0 & 0 & 0 &1 \end{pmatrix}\end{split}\]

Parameters:

• a [real(kind=prec)] :: the angle \(\alpha\)

• b [real(kind=prec)] :: the angle \(\beta\)

• c [real(kind=prec)] :: the angle \(\gamma\)

Return:

euler_mat (4,4) [real(kind=prec)] :: the \(4\times4\) Euler matrix

The user file

Mandatory functions

The user must implement the following functions in the user file

nr_q [integer,parameter= n]

The number of distributions the user intends to calculate

nr_bins [integer,parameter= n]

The number of bins in the distributions the user intends to calculate

min_val (nr_q) [real(kind=prec)]

The lower bounds of the distributions

max_val (nr_q) [real(kind=prec)]

The upper bounds of the distributions

userdim [integer]

The number of integrations the user wishes to carry out to account eg. for beam effects

pass_cut (nr_q) [logical]

This controls whether the event is acceptable. If at least one entry of this array is .true. the event will be calculated and added to the cross section. If individual elements are .false., this event will not be added to the corresponding histogram.

Note

Even though it is possible to calculate multiple closely related cuts simultaneously, this can harm the speed of convergence as the VEGAS algorithm optimises for the cross section and not for the distributions.

userweight [real(kind=prec)]

The weight the user wishes to attach to a given event

names (nr_q) [character(len=namesLen)]

The names of the distributions the user wishes the calculate

filenamesuffix [character(len=filenamesuffixLen)]

The observable-specific suffix to the vegas file

subroutine  fix_mu()

The user needs to choose the renormalisation scale \(\mu^2\) by writing to the variable musq. This can be done on a per-event basis.

A common example would be

SUBROUTINE FIX_MU
musq = Mm**2
END SUBROUTINE FIX_MU

subroutine  inituser()

This is called without arguments once as soon as McMule starts and has read all other configuration, meaning that it can access which_piece and flavour. It may be used to read any further information (like cut configuration etc). The user does not have to print hashes – this is already taken care of – but is very much invited to include information of what it is they are doing.

If the user is using the cut channel of the menu, they may need to set the filenamesuffix variable which is appended to the name of the VEGAS file.

Example for reading a cut:

SUBROUTINE INITUSER
integer cut
read*,cut
write(filenamesuffix,'(I2)') cut
END SUBROUTINE INITUSER

with a global variable cut

function  quant(q1, q2, q3, q4, q5, q6, q7)

The measurement function the user wishes to calculate. This needs to at least set pass_cut but also returns the values of the observables that are to be computed. It usually also calls fix_mu() to fix the renormalisation scale though this can be done elsewhere. If the user wishes to consider polarised scattering, pol1 and pol2 need to be set.

A minimal example that accepts every event and does not calculate a distribution would be

FUNCTION QUANT(q1,q2,q3,q4,q5,q6,q7)
real (kind=prec), intent(in) :: q1(4),q2(4),q3(4),q4(4),q5(4),q6(4),q7(4)
real (kind=prec) :: quant(nr_q)
!! ==== keep the line below in any case ==== !!
call fix_mu
pol1 = 0.
pass_cut = .true.
END FUNCTION QUANT

Parameters:

qi (4) [real(kind=prec)] :: the momenta

Return:

quant (nr_q) [real(kind=prec)] :: the observables that are to be histogrammed

subroutine  userevent(x, ndim)

The user may use this routine in combination with userweight to integrate over further parameters, i.e. to calculate

\[\sigma \sim \int_0^1 {\rm d} x_1 \int_0^1 {\rm d} x_2 \cdots \int_0^1 {\rm d} x_m \times \int {\rm d}\Phi\,\, \left\vert\mathcal{M}_n\right\vert^2\,\, f(x_1,x_2,\cdots,x_n;p_1,\cdots,p_n)\]

with a generalised measurement function \(f\). A minimal example that does not include extra intgration is

SUBROUTINE USEREVENT(X, NDIM)
integer :: ndim
real(kind=prec) :: x(ndim)
userweight = 1.
END SUBROUTINE USEREVENT

Parameters:

• x (ndim) [real(kind=prec)] :: the values of the integration

• ndim [integer] :: the dimension of x, should equal userdim.

Tweaking parameters

In rare cases it may be necessary to tweak some parameters.

integer namesLen

The maximally allowed length of the histogram names.

integer filenamesuffixLen

The maximally allowed length of the observable name as specified in filenamesuffix.

integer bin_kind

The binning mechanism being used, 0 for \({\rm d}\sigma/{\rm d}Q\) and 1 for \(Q {\rm d}\sigma/{\rm d}Q\).

Warning

Note that the latter is not properly tested and should only be used with great care

Technical routines

The following types, variables, and routines are unlikely to be needed by the typical user and are instead aimed at McMule’s developers.

The particle framework

integer maxparticles

The maximal number of particles allowed

type  mlm

Type fields:

• % momentum (4) [real(kind=prec)] :: the momentum

type  particle

Type fields:

• % momentum (4) [real(kind=prec)] :: the momentum

• % effcharge [integer] :: the effective charge, corresponding to the +charge for incoming and -charge for outgoing particles.

• % charge [integer] :: the actual charge

• % incoming [logical] :: .true. for incoming particles

• % lepcharge [integer] :: the lepton family (1 for electrons, 2 for muons, 3 for taus), defaults to zero

type  particles

Type fields:

• % vec (maxparticles) [type(particle)] :: the constituent partciles

• % n [integer] :: the number of particles actually used

• % combo [character(len=1)] :: the flavour combination used, allowed values are * (any combination), x (only mixed), e (only electronic), m (only muonic), t (only tauonic)

function  make_mlm(qq)

Construct a mlm, i.e. a massless momentum

Parameters:

qq (4) [real(kind=prec),in] :: the momentum

function  part(qq, charge, inc[, lepcharge])

Construct a particle.

Parameters:

• qq (4) [real(kind=prec),in] :: the momentum

• charge [integer,in] :: the charge of the particle

• inc [integer,in] :: +1 for incoming, -1 for outgoing

Options:

lepcharge [integer,1] :: the lepton family number

function  parts(ps[, combo])

Construct particles from a list of particles

Parameters:

ps (*) [type(particle),in] :: a list of particle

Options:

combo [character(len=1)] :: the flavour combination used, allowed values are * (any combination), x (only mixed), e (only electronic), m (only muonic), t (only tauonic)

function  eik()

An interface to construct the eikonal factor. eik can be called with

• (kg, pp), using the type particles. The optional flavour combination combo restricts the emission to the desired set of fermion lines. If combo is set to x, all contributions but the self-eikonal are included.

• ({q1,k1}, kg, {q2,k2}), with an explicit call to the momenta of the {massive, massless} emitter, before (1) and after (2) the emission.

Parameters:

• pp [type(particles),in] :: the fermions involved in the photon emission

• qi (4) [real(kind=prec),in] :: the momenta of the massive emitter

• ki [type(mlm),in] :: the momenta of the massless emitter

• kg [type(mlm),in] :: the momentum of the photon

Return:

eik :: the eikonal factor

function  ieik()

An interface to construct the integrated eikonal factor [28]. ieik can be called with

• (xicut, epcmf, pp[, pole]), using the type particles. The optional flavour combination combo restricts the emission to the desired set of fermion lines. If combo is set to x, all contributions but the self-eikonal are included.

• (xicut, epcmf, q1, q2[, pole]), with an explicit call to the momenta of the massive emitter, before (1) and after (2) the emission.

Parameters:

• xicut [real(kind=prec),in] :: \(\xi_c\) (cf. Section Running at NLO and beyond)

• epcmf [real(kind=prec),in] :: square root of scms

• pp [type(particles),in] :: the fermions involved in the photon emission

• qi (4) [real(kind=prec),in] :: the momenta of the massive emitter

Options:

pole [real(kind=prec),out] :: the singular part of the integrated eikonal, as a coefficient of \(1/\epsilon\)

Return:

ieik :: the finite part of the integrated eikonal factor

function  ntssoft(pp, kk, pole)

The (universal) soft contribution to the LBK theorem at 1 loop [9], i.e. the NTS soft function. The optional flavour combination for the particles pp restricts the emission to the desired set of fermion lines. [1]

Parameters:

• pp [type(particles),in] :: the fermions involved in the photon emission

• kk (4) [real(kind=prec),in] :: the momentum of the photon

Options:

pole [real(kind=prec),out] :: the singular part of the NTS soft function, as a coefficient of \(1/\epsilon\)

Return:

ieik :: the finite part of the NTS soft function

Matrix element interface

function  partInterface(q1, q2, q3, q4, q5, q6, q7)

an abstract interface to construct particles for a given process.

Parameters:

qi (4) [real(kind=prec)] :: the momenta

Return:

partInterface [particles] :: the constructed particle string

Package-X function

Note

This section needs to be completed, link to issue

function  DiscB()

function  DiscB_cplx()

function  ScalarC0IR6()

function  ScalarC0IR6_cplx()

function  ScalarC0()

function  ScalarC0_cplx()

function  ScalarD0IR16()

function  ScalarD0IR16_cplx()

VP functions

Note

This section needs to be completed, link to issue

Phase spaces

McMule has implemented a number of phase routines that map from the hypercube to the physical momenta. Here is a list of currently used ones

subroutine  PSD3(ra, q1, m1, q2, m2, q3, m3, weight)

Generic phase space routine for \(1\to2\) decays

Parameters:

• ra (2) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD4(ra, q1, m1, q2, m2, q3, m3, q4, m4, weight)

Generic phase space routine for \(1\to3\) decays

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD4_FKS(ra, q1, m1, q2, m2, q3, m3, q4, weight)

FKS phase space routine for \(1\to3\) decays, requires \(m_4=0\). Tuned for \(\sphericalangle(p_2,q_4)\) and \(E_4\)

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD5(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, weight)

Generic phase space routine for \(1\to4\) decays

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD5_25(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, weight)

Phase space routine for \(1\to4\) decays, tuned for \(\sphericalangle(p_2,q_5)\) and \(E_5\), collinear limit is ra(2) -> 0

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD5_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

FKS phase space routine for \(1\to4\) decays, requires \(m_5=0\). Tuned for \(\sphericalangle(p_2,q_5)\) and \(E_5\)

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD6(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, m6, weight)

Generic phase space routine for \(1\to5\) decays

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD6_23_24_34(ra, q1, m1, q2, m2, q5, m5, q6, m6, q3, m3, q4, m4, weight)

Phase space routine for \(1\to5\) decays with FKS-ish tuning. This is designed for the decay \(\mu^+\to e^+\nu\bar\nu e^+e^-\). q2 should be the unique particle (electron) and q3 and q4 are the identical particles (postirons):

                   _/  q
                   /|   2
---<---*~~~~~~~~~~*
       |          |
       ^          ^  q
       | q        |   3
          4

The ‘spectator’ neutrinos are q5 and q6. Start by generating p2 and p3 at an angle * = arccos(y2):

              ^ p2
             |||
             |||
p3         __|||
 --__     /  |||
     --__/  *|||
         --__|||
             |||

Generate p4 at an angle * = arccos(y3) and rotating by an angle phi w.r.t. to p3:

             |||      / p4
          _--|||--_  /
         -_  |||  _-/
p3         --___-- /
 --__        |||__/
     --__    |||*/
         --__|||/
             |||

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD6_23_24_34_E56(ra, q1, m1, q2, m2, q5, m5, q6, m6, q3, m3, q4, m4, weight)

Phase space routine for \(1\to5\) decays with FKS-ish tuning, similar to PSD6_23_24_34() but with special tuning on the \(E_5+E_6\) tail.

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD6_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, weight)

FKS phase space routine for \(1\to5\) decays, requires \(m_6=0\). Tuned for \(\sphericalangle(p_2,q_6)\) and \(E_6\)

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD6_25_26_m50_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, weight)

FKS phase space routine for \(1\to5\) decays, requires \(m_5=m_6=0\). Tuned for \(\sphericalangle(p_2,q_{5,6})\) and \(E_{5,6}\)

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD6_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

Double-FKS phase space routine for \(1\to5\) decays, requires \(m_5=m_6=0\). Tuned for \(\sphericalangle(p_2,q_{5,6})\) and \(E_{5,6}\)

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD7(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, m6, q7, m7, weight)

Generic phase space routine for \(1\to6\) decays

Parameters:

• ra (14) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD7_27_37_47_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, m6, q7, weight)

FKS phase space routine for \(1\to6\) decays, tuned for \(\sphericalangle(p_{2,3,4},q_7)\) and \(E_7\)

Parameters:

• ra (14) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD7_27_37_47_E56_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, m6, q7, weight)

FKS phase space routine for \(1\to6\) decays, tuned for \(\sphericalangle(p_{2,3,4},q_7)\) and \(E_7\) and tuned for the \(E_5+E_6\) tail

Parameters:

• ra (14) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX2(ra, q1, m1, q2, m2, q3, m3, q4, m4, weight)

Generic phase space routine for \(2\to2\) cross sections

Parameters:

• ra (2) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

FKS phase space routine for \(2\to3\) cross sections, requires \(m_5=0\). Tuned for ISR and not FSR

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_35_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight[, sol])

FKS phase space routine for \(2\to3\) cross sections, requires \(m_5=0\). Tuned for \(\sphericalangle(q_3,q_5)\) and \(E_5\)

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

Options:

sol [integer,in] :: which solution to pick

subroutine  PSX4(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, m6, weight)

Generic phase space routine for \(2\to4\) cross sections

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

Double-FKS phase space routine for \(2\to4\) cross sections, requires \(m_5=m_6=0\). Tuned for ISR and not :term`FSR`

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_35_36_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight[, sol])

Double-FKS phase space routine for \(2\to4\) cross sections, requires \(m_5=m_6=0\). Tuned for \(\sphericalangle(q_3,q_{5,6})\) and \(E_{5,6}\)

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

Options:

sol [integer,in] :: which solution to pick

subroutine  PSD6_P_25_26_m50_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, weight)

FKS phase space routine for \(1\to5\) decays, requires \(m_5=m_6=0\). Tuned for \(\sphericalangle(p_2,q_{5,6})\) and \(E_{5,6}\) Partioning of PSD6_25_26_m50_FKS() with \(s_{26}<s_{25}\).

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD6_26_2x5(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, m6, weight)

Phase space routine for \(1\to5\) decays with FKS-ish tuning. Modification of PSD6_23_24_34() with \(2\leftrightarrow 5\). This is designed for the decay \(\mu^+\to e^+\nu\bar\nu e^+e^-\).

Parameters:

• ra (11) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSD7_27_37_47_2x5_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, m5, q6, m6, q7, weight)

FKS phase space routine for \(1\to6\) decays, tuned for \(\sphericalangle(p_{2,3,4},q_7)\) and \(E_7\) Modification of PSD7_27_37_47_FKS() with \(2\leftrightarrow 5\).

Parameters:

• ra (14) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_P_15_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

FKS phase space routine for \(2\to3\) cross sections, requires \(m_5=0\). Tuned for ISR and not FSR. Partioning of PSX3_FKS() with \(s_{15}<s_{35}\).

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_P13_35_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

FKS phase space routine for \(2\to3\) cross sections, requires \(m_5=0\). Tuned for \(\sphericalangle(q_3,q_5)\) and \(E_5\) Partioning of PSX3_35_FKS() with \(s_{15}>s_{35}\).

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_coP13_35_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

The corner piece to PSX3_P13_35_FKS()

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_P_15_25_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

FKS phase space routine for \(2\to3\) cross sections, requires \(m_5=0\). Tuned for ISR and not FSR. Partioning of PSX3_FKS() with \({\rm min}\big(s_{15},s_{25}\big)<{\rm min}\big(s_{35},s_{45}\big)\).

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_P_35_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

FKS phase space routine for \(2\to3\) cross sections, requires \(m_5=0\). Tuned for \(\sphericalangle(q_3,q_5)\) and \(E_5\) Partioning of PSX3_35_FKS() with \(s_{35}<{\rm min}\big(s_{15},s_{25},s_{45}\big)\).

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_P_45_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

FKS phase space routine for \(2\to3\) cross sections, requires \(m_5=0\). Tuned for \(\sphericalangle(q_3,q_5)\) and \(E_5\) Partioning of PSX3_35_FKS() with \(s_{45}<{\rm min}\big(s_{15},s_{25},s_{35}\big)\) and \(3\leftrightarrow4\)

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_coP_35_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

The corner piece to PSX3_P_35_FKS()

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX3_coP_45_FKS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, weight)

The corner piece to PSX3_P_45_FKS()

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_P_15_16_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

Double-FKS phase space routine for \(2\to4\) cross sections, requires \(m_5=m_6=0\). Tuned for ISR and not :term`FSR` Partioning of PSX4_FKSS() with \({\rm min}\big(s_{15}, s_{16}\big)<{\rm min}\big(s_{35},s_{36}\big)\).

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_P_35_36_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

Double-FKS phase space routine for \(2\to4\) cross sections, requires \(m_5=m_6=0\). Tuned for \(\sphericalangle(q_3,q_{5,6})\) and \(E_{5,6}\) Partioning of PSX4_35_36_FKSS() with \({\rm min}\big(s_{15}, s_{36}\big)<{\rm min}\big(s_{15},s_{25},s_{45},s_{16},s_{26},s_{46}\big)\).

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_coP_35_36_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

The corner piece to PSX4_P_35_36_FKSS()

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_P13_35_36_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

Double-FKS phase space routine for \(2\to4\) cross sections, requires \(m_5=m_6=0\). Tuned for \(\sphericalangle(q_3,q_{5,6})\) and \(E_{5,6}\) Partioning of PSX4_35_36_FKSS() with \({\rm min}\big(s_{15}, s_{16}\big)>{\rm min}\big(s_{35},s_{36}\big)\).

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_coP13_35_36_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

The corner piece to PSX4_P13_35_36_FKSS()

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_P_15_16_25_26_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

Double-FKS phase space routine for \(2\to4\) cross sections, requires \(m_5=m_6=0\). Tuned for ISR and not :term`FSR` Partioning of PSX4_FKSS() with \({\rm min}\big(s_{15}, s_{16},s_{25}, s_{26}\big)<{\rm min}\big(s_{35},s_{36},s_{54},s_{46}\big)\).

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_P_45_46_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

Double-FKS phase space routine for \(2\to4\) cross sections, requires \(m_5=m_6=0\). Tuned for \(\sphericalangle(q_3,q_{5,6})\) and \(E_{5,6}\) Partioning of PSX4_35_36_FKSS() with \({\rm min}\big(s_{45}, s_{46}\big)>{\rm min}\big(s_{15},s_{25},s_{35},s_{16},s_{26},s_{36}\big)\).

Parameters:

• ra (8) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

subroutine  PSX4_coP_45_46_FKSS(ra, q1, m1, q2, m2, q3, m3, q4, m4, q5, q6, weight)

The corner piece to PSX4_P_45_46_FKSS()

Parameters:

• ra (5) [real(kind=prec),in] :: the random numbers

• qi (4) [real(kind=prec),out] :: the momenta

• mi [real(kind=prec),in] :: the masses

• weight [real(kind=prec),out] :: the Jacobian

[1]

The user is allowed to further split mixed contributions at NNLO, i.e. contributions with emissions connecting different fermion lines. This is achieved via the optional parameter mx of the auxiliary function combonts. The latter sets the desired flavour combination for ntssoft, and mx=1 allows to choose among different mixed contributions. For example, for \(\ell_1\ell_2\to\ell_1\ell_2\) scattering, if a formal charge \(Q_{1(2)}\) is assigned for each photon emission from \(\ell_{1(2)}\), ntssoft will be able to distinguish among the contributions labelled by \(Q_1^5\,Q_2^3\), \(Q_1^4\,Q_2^4\) and \(Q_1^3\,Q_2^5\).