Getting started

McMule is written in Fortran 95 with helper and analysis tools written in python. This guide will help you to get started using McMule by describing in detail how to calculate the NLO corrections to \(\tau\to[\nu\bar\nu] e\gamma\). Since the neutrinos are not detected, we average over them, indicated by the brackets. Hence, we have to be fully inclusive w.r.t. the neutrinos. Still, the code allows to make any cut on the other final-state particles. As we will see, the BR for this process, as measured by BaBar [14, 18] has a discrepancy of more than \(3\,\sigma\) from the SM value. This will illustrate the importance of fully differential NLO corrections in QED.

Obtaining the code

McMule is distributed multiple ways

• as a precompiled executable for recent-ish GNU Linux distributions. To be precise, your version of glibc needs to be newer than 2.17. The currently supported versions of most popular distributions (CentOS, Debian, Ubuntu, Fedora, RHEL) should be fine.

• as a Docker image that can be used on any platform.

• as the source code on Gitlab that can be compiled by the user. This contains the current release in the default release branch as well as the developer preview (devel).

Here we will focus on the first method as it is by far the easiest. For developers and tinkerers, we refer to Section Manual compilation on how to compile the code yourself.

First, obtain the McMule distribution from our website and extract the tarball

$ tar xzvf mcmule-release.tar.gz
mcmule-release/mcmule
mcmule-release/mcmule.mod

That’s it. You can, if you want, install McMule to use it from any directory with the following commands but this is not required

$ cp mcmule-release/mcmule /usr/local/bin
$ cp mcmule-release/mcmule.mod /usr/local/include

To make use of McMule results, we also require the pymule python package It can be used from the tools/ folder of the McMule repository but it is recommended that the user installs it

$ pip3 install git+https://gitlab.com/mule-tools/pymule.git

Simple runs at LO

Setting McMule up

In this example we want to compute two distributions, the invariant mass of the \(e\gamma\) pair, \(m_{\gamma e}\equiv \sqrt{(p_e+p_\gamma)^2}\), and the energy of the electron, \(E_e\), in the rest frame of the tau. To avoid an IR singularity in the BR, we have to require a minimum energy of the photon. We choose this to be \(E_\gamma \ge 10\,{\rm MeV}\) as used in [14, 18].

At first, we need to find out how the process \(\tau\to\nu\bar\nu e\gamma\) is implemented in McMule. For this, we refer to the table in Section Available processes and which_piece that specifies the pieces (sometimes called which_piece) that is required for a generic processes. The generic process is a prototype for the physical process such as \(\ell \to\nu\bar\nu \ell' \gamma\) where the flavour of the lepton \(\ell\) is left open. In our case, we need to consider the row for \(\mu\to\nu\bar\nu e\gamma\). Since we are only interested in LO, the only which_piece we need is m2enng0. To change from the generic process \(\mu\to\nu\bar\nu e\gamma\) to the process we are actually interested in, \(\tau\to\nu\bar\nu e\gamma\), we pick the flavour tau-e which refers to a \(\tau \to e \cdots\) transition. Other options would be tau-mu for \(\tau\to\mu \cdots`\) or mu-e for \(\mu\to e\cdots\).

Next, we need to find out which particle ordering is used in McMule for this piece, i.e. which variable will contain eg. the electron momentum. This is called the particle identification or PID. We can refer to the table in Section Particle ID to find that for the which_piece m2enng0, we have

\[\mu^-(p_1) \to e^-(p_2) \big[\bar\nu_e\nu_\mu](p_3,p_4) \gamma(p_5)\]

We can now implement our observables. For this, we need to define a user.f95 file in the src folder. An empty template can be found in the file tools/user-empty.f95. We can use this file to the measurement function we want to calculate, i.e. which distributions and cuts we want to apply. We can further add some code that will execute at the beginning of the Monte Carlo run (allowing us eg. to further configure our calculation) and for each event (to simulate beam spread).

We begin by specifying the metadata of our histograms: we want two histograms (nr_q = 2) with 90 bins each (nr_bins = 90). The ranges should be \(0 < m_{\gamma e} < 1800\,{\rm MeV}\) and \(0 \le E_e \le 900\,{\rm MeV}\).

Listing 1 The metadata for our calculation with two histograms (nr_q = 2) with 90 bins each (nr_bins = 90). The ranges should be \(0 < m_{\gamma e} < 1800\,{\rm MeV}\) and \(0 \le E_e \le 900\,{\rm MeV}\).

  real:: &
     min_val(nrq) = (/    0.,   0. /)
  real:: &
     max_val(nrq) = (/ 1800., 900. /)
  integer :: userdim = 0

Note

Finding suitable values for the ranges can be tricky beyond LO and usually requires a few test runs. Since all histograms have the same number of bins, one is often forced to have empty bins to ensure ‘nice’ bin widths.

We can now define the actual measurement function called quant(). We need to

• calculate the invariant mass of the \(e\gamma\) pair. This is done using the momentum-squaring function sq(). The result is store in the first distribution, quant(1).

• store the electron energy in quant(2). Since this is frame-dependent, we need to know that McMule generates the particles in the tau rest frame. However, in general it is better to boost into that frame. Further, McMule stores momenta as (/px,py,pz,E/), meaning the energy is q2(4).

• cut on the photon energy q5(4). The variable pass_cut controls the cuts. Initially it is set to .true., to indicate that the event is kept. Applying a cut amounts to setting pass_cut to .false..

Listing 2 The measurement function at LO

  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

  pass_cut = .true.
  if(q5(4) < 10._prec) pass_cut = .false.

  names(1) = 'minv'
  quant(1) = sqrt(sq(q2+q5))

  names(2) = 'Ee'
  quant(2) = q2(4)

  END FUNCTION QUANT

Additionally to the numeric value in quant(i) we store a human-readable name in names(i).

Warning

The maximal length of these names is defined in the variable namesLen which defaults to 6 characters.

Also note that this measurement function is not IR-safe!

We now need to compile our observable into a shared library so that McMule can load it. To do this, we run

$ gfortran -fPIC --shared -o user.so user.f95

This requires the mcmule.mod file to either be in the current directory or installed somewhere the compiler can find it. Otherwise, one needs to add the corresponding flag

$ gfortran -I/path/to/the/folder/of/mcmule.mod/ -fPIC --shared -o user.so user.f95

We now need to re-compile McMule to ensure that we have the correct version of user.f95.

Warning

The mcmule.mod header file is autogenerated by GFortran during the compilation of McMule. If you are using a copy of GFortran prior to version 8, this means you will have to regenerate the header file manually. To do this, you can use the build-header.sh script.

Running McMule manually

Now the mule is ready to trot. For quick and dirty runs of McMule, the easiest way is to just execute the mcmule binary in the same directory as the user.so file and input the configuration by hand. However, since this is not how the code is meant to be used, it will not prompt the user but just expect the correct input.

We now need to choose the statistics we want. For this example, we pick 10 iterations with \(1000\times 10^3\) points each for pre-conditioning and 50 iterations with \(1000\times 10^3\) points each for the actual numerical evaluation (cf. Section Statistics for some heuristics to determine the statistics needed). We pick a random seed between 0 and \(2^{31}-1\) (cf. Section Random number generation), say 70998, and for the input variable which_piece we enter m2enng0 as discussed above. The flavour variable is now set to tau-e to change from the generic process \(\mu\to\nu\bar\nu e\gamma\) to the process we are actually interested in, \(\tau\to\nu\bar\nu e\gamma\). This system is used for other processes as well. The input variable which_piece determines the generic process and the part of it that is to be computed (i.e. tree level, real, double virtual etc.). In a second step, the input flavour associates actual numbers to the parameters entering the matrix elements and phase-space generation. This means that we need to input the following (the specifications for the input can be found in Table 1):

Warning

When running mcmule outside the normal repository, you need to make sure that an out/ folder exists.

$ gfortran -fPIC --shared -o user.so user.f95
$ ./mcmule
1000
10
10000
50
70998
1.0
1.0
m2enng0
tau-e


         *******************************************
         *              C O L L I E R              *
         *                                         *
         *        Complex One-Loop Library         *
         *      In Extended Regularizations        *
         *                                         *
         *    by A.Denner, S.Dittmaier, L.Hofer    *
         *                                         *
         *              version 1.2.3              *
         *                                         *
         *******************************************

   - * - * - * - * - * - * - * -
        Version information
  Full SHA: 3342511
  Git  SHA: 1fbc291
  Git branch: HEAD
   - * - * - * - * - * - * - * -
 Calculating tau->e nu nu gamma at LO
   - * - * - * - * - * - * - * -

 internal avgi, sd:    31902651645147.434        3242845143300.6875
 internal avgi, sd:    36962119569527.797        1491060763146.8340
 internal avgi, sd:    39908483081760.562        701506532475.22485
 internal avgi, sd:    41908326436302.352        183707731215.21738
 internal avgi, sd:    41771416194096.336        55441877946.459091
 internal avgi, sd:    41871492562379.680        27645368422.638184
 internal avgi, sd:    41870973597620.547        21172712863.774796
 internal avgi, sd:    41881968277900.094        17287639806.820400
 internal avgi, sd:    41894819976244.469        15148087824.181145
 internal avgi, sd:    41892443511666.180        13860145189.905710
 internal avgi, sd:    41883909931737.320        9081654737.2369480
 internal avgi, sd:    41891877400107.203        5996399688.5281315
 internal avgi, sd:    41887401454137.172        4967120009.5028763
 internal avgi, sd:    41894988589984.109        4318086453.2893734
 internal avgi, sd:    41895218930734.938        3855189670.2831044
 internal avgi, sd:    41893628691682.039        3569029881.5161963
 internal avgi, sd:    41895702521658.094        3301046354.0162683
 internal avgi, sd:    41894921420510.164        3068605199.3146548
 internal avgi, sd:    41894380982483.836        2884341089.4262500
 internal avgi, sd:    41894136940077.953        2719744511.4164872
 internal avgi, sd:    41894755045877.328        2575585554.2240872
 internal avgi, sd:    41894180414331.000        2448105950.1048141
 internal avgi, sd:    41892974586371.242        2335940686.3421283
 internal avgi, sd:    41892018243977.422        2237743190.2910728
 internal avgi, sd:    41892128399199.852        2151548541.7080536
 internal avgi, sd:    41891054946079.172        2079987935.2725692
 internal avgi, sd:    41890529496649.336        2009220806.5744867
 internal avgi, sd:    41889627128683.867        1952022430.9618585
 internal avgi, sd:    41889091169697.750        1901209181.0766854
 internal avgi, sd:    41889491086513.711        1851335964.4142988
 internal avgi, sd:    41889024177143.492        1804584253.3775585
 internal avgi, sd:    41888652800094.414        1763879105.2402925
 internal avgi, sd:    41888186242209.695        1734933321.8742726
 internal avgi, sd:    41888838662647.031        1702558920.3787835
 internal avgi, sd:    41888878166048.805        1664687915.8245957
 internal avgi, sd:    41888871786161.102        1628412032.4284451
 internal avgi, sd:    41888673286317.961        1598222188.0324447
 internal avgi, sd:    41888363240043.000        1570566301.1038322
 internal avgi, sd:    41888533287695.047        1541834130.0455377
 internal avgi, sd:    41888087919550.688        1514438031.2038479
 internal avgi, sd:    41887838382975.297        1487390337.2575607
 internal avgi, sd:    41887692329889.953        1465777595.5131192
 internal avgi, sd:    41887528786746.531        1450928637.2665195
 internal avgi, sd:    41887814931451.086        1432276674.4390638
 internal avgi, sd:    41887764015763.508        1409275835.0397925
 internal avgi, sd:    41887871329949.469        1388512123.7455208
 internal avgi, sd:    41887728279057.234        1369450030.9152539
 internal avgi, sd:    41887673022843.117        1350554230.4978600
 internal avgi, sd:    41887824787223.562        1332418851.3856776
 internal avgi, sd:    41887720562604.266        1315515744.3643637
 internal avgi, sd:    41887747627717.641        1297906527.3172038
 internal avgi, sd:    41887385610296.359        1280408319.1259799
 internal avgi, sd:    41887163475026.672        1262887224.9655898
 internal avgi, sd:    41887020587065.422        1248392301.3985555
 internal avgi, sd:    41886965905979.375        1236043197.7830524
 internal avgi, sd:    41887132288349.984        1225259563.1290646
 internal avgi, sd:    41887118281531.000        1211732414.0191844
 internal avgi, sd:    41887256099447.883        1199948076.9753296
 internal avgi, sd:    41887425753145.656        1188759649.9116085
 internal avgi, sd:    41887079359692.539        1176252324.5589268
 points:    10*  1M +  50* 10M         random seed:      70998
 part: m2enng0                      xicut: 1.00000    delcut: 1.00000
 points on phase space  321225751 thereof fucked up          0



 result, error:  {  4.18871E+13,  1.17625E+09 };   chisq:  1.27

   - * - * - * - * - * - * - * -

McMule begins by printing some auto-versioning information (the SHA1 hashes of the source code and the git version) as well as some user-defined information from the subroutine inituser(). Next, the integration begins. After every iteration, McMule prints the current best estimate and error of the total cross section or decay rate. Before exiting, it will also print again the input used as well as the number of points evaluated and the final result. This run took approximately 15 minutes.

Table 1 The options read from stdin by McMule. The calls are multiplied by 1000.

Variable name	Data type	Comment
nenter_ad	integer	calls / iteration during pre-conditioning
itmx_ad	integer	iterations during pre-conditioning
nenter	integer	calls / iteration during main run
itmx	integer	iterations during main run
ran_seed	integer	random seed \(z_1\)
xinormcut	real(prec)	the \(0<\xi_{c}\le1\) parameter
delcut	real(prec)	the \(\delta_{\text{cut}}\) parameter (or at NNLO the second \(\xi_{c}\))
which_piece	char(10)	the part of the calculation to perform
flavour	char(8)	the particles involved
(opt)	unknown	the user can request further input during inituser()

Analysing the output

After running McMule we want to calculate the actual cross section or decay rate and make plots. The McMule output is saved to the out/ folder as a .vegas file that contains the entire state of the integrator (cf. Section Differential distributions and intermediary state files). We can open this file in python and make plots.

While it is possible to open just a single file using importvegas(), this is rarely done as real-world calculations can involve hundreds of .vegas files. Instead, we move the .vegas file into a new directory, say example1 and then use sigma() and mergefks().

from pymule import *

lifetime = 1/(1000*(6.582119e-25)/(2.903e-13))
# define vegas directory
setup(folder="example1/")
dat = scaleset(
    mergefks(sigma("m2enng0")),
    GF**2*alpha*lifetime
)

dat.keys()
# dict_keys(['time', 'chi2a', 'value', 'Ee', 'minv'])

Warning

In McMule the numerical value of the Fermi constant \(G_F\) and the fine-structure constant \(\alpha\) are set to one for predominately historical reasons. This needs to be restored in python, eg. using scaleset()

The variable dat now contains the runtime (time), branching ratio (after multiplication with the lifetime, value), and \(\chi^2\) of the integration (chi2a) as well as our distributions (Ee and minv). Numerical values such as cross sections or branching ratios are stored as numpy arrays with errors as np.array([y, dy]). Distributions are stored as numpy \(N\times 3\) matrices

np.array([[x1, y1, dy1],
          [x2, y2, dy2],
          [x3, y3, dy3],
          [x4, y4, dy4],
          [x5, y5, dy5],
          ...
          [xn, yn, dyn]])

These can be manipulated eg. using the tools of pymule described in Section pymule user guide. For now, we will just make a plot of the \(E_e\) distribution Ee

from matplotlib.pyplot import *
fig = plt.figure()
errorband(dat['Ee'])
plt.ylabel(r'$\D\mathcal{B}/\D E_e$')
plt.xlabel(r'$E_e\,/\,{\rm MeV}$')
mulify(fig)
fig.savefig("dummy.svg")

Figure 1 Result of the LO test run for the \(E_e\) distribution

Running at NLO and beyond

A few things change once we go beyond LO since we can have extra radiation. To account for this, more which_piece need to be ran and then correctly combined. This also increases the number of runs necessary, meaning that the manual approach from above is no longer feasible.

Setting McMule up

Referring back to Section Available processes and which_piece we find that we need the pieces m2enngF and m2enngR for virtual and real corrections respectively. The PID table of Section Particle ID tells us that the real photon can is going to be q6.

We first need to decide whether we want to calculate exclusive or inclusive decays. The details here depend on the exact experimental situation which can be tricky to properly implement. Following the BaBar analysis [14, 18] we will consider the exclusive radiative decay, i.e. we request precisely one photon with energy \(E_\gamma>10\,{\rm MeV}\). The function quant() will have to take this into account with the additional argument q6, the momentum of the second photon.

To ensure IR safety, we define the harder and softer photon gh and gs, respectively, and require that the former (latter) has energy larger (smaller) than \(10\,{\rm MeV}\). This new version of quant() is also suitable for the LO calculation and it is generally advisable to use a single quant() function for all parts of a computation.

Listing 3 The measurement function beyond LO. The changes w.r.t. to LO are highlighted.

  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)
  real (kind=prec) :: gs(4), gh(4)
  !! ==== keep the line below in any case ==== !!
  call fix_mu

  if (q5(4) > q6(4)) then
    gh = q5 ; gs = q6
  else
    gh = q6 ; gs = q5
  endif

  if (gh(4) < 10.) pass_cut = .false.
  if (gs(4) > 10.) pass_cut = .false.

  names(1) = 'minv'
  quant(1) = sqrt(sq(q2+qh))

  names(2) = 'Ee'
  quant(2) = q2(4)

  END FUNCTION QUANT

Running McMule

The FKS scheme used in McMule introduces an unphysical parameter called \(\xi_c\) that can be varied between

\[0<\xi_{c}\le\xi_{\text{max}} = 1-\frac{\big(\sum_i m_i\big)^2}{s}\]

Checking the independence of physical results on the latter serves as a consistency check, both of the implementation of McMule but also of the IR safety of the measurement function. To do this, it can help to disentangle m2enngF into m2enngV and m2enngC though this is not necessary. Only the latter depends on \(\xi_c\) and this part is typically much faster in the numerical evaluation.

A particularly convenient way to run McMule is using menu files [1]. A menu file contains a list of jobs to be computed s.t. the user will only have to vary the random seed and \(\xi_{c}\) by hand as the statistical requirements are defined globally in a config file. This is completed by a submission script, usually called submit.sh. The submit script is what will need to be launched which in turn will take care of the starting of different jobs. It can be run on a normal computer or on a SLURM cluster [27]. To prepare the run in this way we can use pymule

Listing 4 The steps necessary to use pymule to prepare running McMule. Note that numbers listed as seeds are random and hence not reproducible.

$ pymule create -i
What generic process? [m2enn] m2enng
Which flavour combination? [mu-e] tau-e
How many / which seeds? [5]
Which xi cuts? [[0.5, 0.25, 0.125]]
Where to store data? [m2enngtau-e]  example2
Which pieces? [['0', 'V', 'R']] 0, V, C, R
How much statistics for 0 (pc, pi, c, i)? [(10000, 20, 100000, 100)] 1000,10,1000,50
How much statistics for  V (pc, pi, c, i)? [(10000, 20, 100000, 100)] 1000,10,1000,50
How much statistics for  C (pc, pi, c, i)? [(10000, 20, 100000, 100)] 1000,10,1000,50
How much statistics for  R (pc, pi, c, i)? [(10000, 20, 100000, 100)] 5000,50,10000,100
Building files. To rerun this, execute
pymule create\
      --seeds 70998 66707 69184 75845 63937 \
      -xi 0.5 0.25 0.125 \
      --flavour tau-e \
      --genprocess m2enng \
      --output-dir babar-tau-e \
      --prog mcmule \
      --stat  R,5000,50,10000,100 \
      --stat 0,1000,10,1000,50 \
      --stat  V,1000,10,1000,50 \
      --stat  C,1000,10,1000,50
Expect 3750 iterations, 20.250000G calls
Created menu, config and submit script in  example2
Please change the ntasks and time options accordingly

When using the tool, we are asked various questions, most of which have a default answer in square brackets. In the end pymule will create a directory that the user decided to call example2, where all results will be stored. The menu and config files generated by pymule are shown in Listing 5 and Listing 6

Listing 5 menu file for the present calculation

## Generated at 16:00 on February 28 2020 by yannickulrich
# git version: redesign (b558978)

conf example2/m2enng-tau-e.conf

run 19397 1.000000 m2enng0 tau-e 0
run 52088 1.000000 m2enng0 tau-e 0
run 83215 1.000000 m2enng0 tau-e 0
run 93857 1.000000 m2enng0 tau-e 0
run 86361 1.000000 m2enng0 tau-e 0


run 19397 1.000000 m2enngV tau-e 0
run 52088 1.000000 m2enngV tau-e 0
run 83215 1.000000 m2enngV tau-e 0
run 93857 1.000000 m2enngV tau-e 0
run 86361 1.000000 m2enngV tau-e 0


run 19397 0.500000 m2enngC tau-e 0
run 52088 0.500000 m2enngC tau-e 0
run 83215 0.500000 m2enngC tau-e 0
run 93857 0.500000 m2enngC tau-e 0
run 86361 0.500000 m2enngC tau-e 0
run 19397 0.500000 m2enngR tau-e 0
run 52088 0.500000 m2enngR tau-e 0
run 83215 0.500000 m2enngR tau-e 0
run 93857 0.500000 m2enngR tau-e 0
run 86361 0.500000 m2enngR tau-e 0

run 19397 0.250000 m2enngC tau-e 0
run 52088 0.250000 m2enngC tau-e 0
run 83215 0.250000 m2enngC tau-e 0
run 93857 0.250000 m2enngC tau-e 0
run 86361 0.250000 m2enngC tau-e 0
run 19397 0.250000 m2enngR tau-e 0
run 52088 0.250000 m2enngR tau-e 0
run 83215 0.250000 m2enngR tau-e 0
run 93857 0.250000 m2enngR tau-e 0
run 86361 0.250000 m2enngR tau-e 0

run 19397 0.125000 m2enngC tau-e 0
run 52088 0.125000 m2enngC tau-e 0
run 83215 0.125000 m2enngC tau-e 0
run 93857 0.125000 m2enngC tau-e 0
run 86361 0.125000 m2enngC tau-e 0
run 19397 0.125000 m2enngR tau-e 0
run 52088 0.125000 m2enngR tau-e 0
run 83215 0.125000 m2enngR tau-e 0
run 93857 0.125000 m2enngR tau-e 0
run 86361 0.125000 m2enngR tau-e 0

Listing 6 Configuration file for the present calculation

## Generated at 16:00 on February 28 2020 by yannickulrich
# git version: redesign (b558978)

# specify the program to run relative to `pwd`
binary=mcmule

# specify the output folder
folder=example2/

# Specify the variables nenter_ad, itmx_ad, nenter and itmx
# for each piece you want to run.
declare -A STAT=(
  ["m2enng0"]="1000\n10\n1000\n50"
  ["m2enngC"]="1000\n10\n1000\n50"
  ["m2enngR"]="5000\n50\n10000\n100"
  ["m2enngV"]="1000\n10\n1000\n50"
)

To start mcmule, we now just need to execute the created example2/submit.sh after copying the user library user.so into the same folder. Note that per default this will spawn at most as many jobs as the computer pymule ran on had CPU cores. If the user wishes a different number of parallel jobs, change the fifth line of example2/submit.sh to

#SBATCH --ntasks=<number of cores>

To now run McMule, just execute

$ nohup ./example2/submit.sh &

The nohup is not technically necessary but advisable, especially on remote systems. When running on a SLURM system, the other SLURM parameters --partition, --time, and --clusters need to be adapted as well.

Warning

The submission script will call itself on multiple occasions. Therefore, it is not advisable to change its name or the name of the run directory without taking precautions.

Analysing the results

After running the code, we need to combine the various which_piece into physical results that we will want to use to create plots. This is the moment where pymule’s mergefks() shines.

Listing 7 An example code to analyse the results for \(\tau\to\nu\bar\nu e\gamma\) in pymule. Note that, in the Fortran code \(G_F = \alpha = 1\). In pymule they are at their physical values.

from pymule import *

# To normalise branching ratios, we need the tau lifetime
lifetime = 1/(1000*(6.582119e-25)/(2.903e-13))

# The folder where McMule has stored the statefiles
setup(folder='example2/out.tar.bz2')

# Import LO data and re-scale to branching ratio
LO = scaleset(mergefks(sigma('m2enng0')), GF**2*lifetime*alpha)

# Import NLO corrections from the three pieces
NLO = scaleset(mergefks(
    sigma('m2enngR'),      # real corrections
    sigma('m2enngC'),     # counter term
    anyxi=sigma('m2enngV') # virtual corrections
), GF**2*lifetime*alpha**2)

# The branching ratio at NLO = LO + correction
fullNLO = plusnumbers(LO['value'], NLO['value'])

# Print results
print("BR_0 = ", printnumber(LO['value']))
print("dBR  = ", printnumber(NLO['value']))

# Produce energy plot
fig1, (ax1, ax2) = kplot(
    {'lo': LO['Ee'], 'nlo': NLO['Ee']},
    labelx=r"$E_e\,/\,{\rm MeV}$",
    labelsigma=r"$\D\mathcal{B}/\D E_e$"
)
ax2.set_ylim(-0.2,0.01)

# Produce visible mass plot
fig2, (ax1, ax2) = kplot(
    {'lo': LO['minv'], 'nlo': NLO['minv']},
    labelx=r"$m_{e\gamma}\,/\,{\rm MeV}$",
    labelsigma=r"$\D\mathcal{B}/\D m_{e\gamma}$"
)
ax1.set_yscale('log')
ax1.set_xlim(1000,0) ; ax1.set_ylim(5e-9,1e-3)
ax2.set_ylim(-0.2,0.)

Once pymule is imported and setup, we import the LO and NLO which_piece and combine them using two central pymule commands that we have encountered above: sigma() and mergefks(). sigma() takes the which_piece as an argument and imports matching results, already merging different random seeds. pymule`mergefks() takes the results of (multiple) sigma() invocations, adds results with matching \(\xi_{c}\) values and combines the result. In the present case, \(\sigma_n^{(1)}\) is split into multiple contributions, namely m2enngV and m2enngC. This is indicated by the anyxi argument.

Next, we can use some of pymule’s tools (cf. Listing Listing 7) to calculate the full NLO BRs from the corrections and the LO results

\[\begin{split}\mathcal{B}|_{\text{LO}} &= 1.8339(1) \times 10^{-2}\\ \mathcal{B}|_{\text{NLO}} &= 1.6451(1) \times 10^{-2}\end{split}\]

which agree with [10, 21], but \(\mathcal{B}|_{\text{NLO}}\) is in tension with the value \(\mathcal{B}|_{\text{exp}} = 1.847(54)\times 10^{-2}\) reported by BaBar [14, 18]. As discussed in [21, 24] it is very likely that this tension would be removed if a full NLO result was used to take into account the effects of the stringent experimental cuts to extract the signal. This issue has been explained in detail in [21, 24, 25].

As a last step, we can use the matplotlib-backed kplot() command to present the results for the distributions (logarithmic for \(m_{e\gamma}\) and linear for \(E_e\)). The upper panel of Figure 2 shows the results for the invariant mass \(m_{e\gamma}\) at LO (green) and NLO (blue) in the range \(0\le m_{e\gamma} \le 1\,{\rm GeV}\). Note that this, for the purposes of the demonstration, does not correspond to the boundaries given in the run.

Figure 2 Results of the toy run to compute \(m_{e\gamma}\) for \(\tau\to\nu\bar\nu e\gamma\). Upper panels show the LO (green) and NLO (blue) results, the lower panels show the NLO K factor.

The distribution falls sharply for large \(m_{e\gamma}\). Consequently, there are only few events generated in the tail and the statistical error becomes large. This can be seen clearly in the lower panel, where the NLO \(K\) factor is shown. It is defined as

\[K^{(1)} = 1+ \frac{\mathrm{d}\sigma^{(1)}} {\mathrm{d}\sigma^{(0)}}\]

and the band represents the statistical error of the Monte Carlo integration. To obtain a reliable prediction for larger values of \(m_{e\gamma}\), i.e. the tail of the distribution, we would have to perform tailored runs. To this end, we should introduce a cut \(m_{\text{cut}}\ll m_\tau\) on \(m_{e\gamma}\) to eliminate events with larger invariant mass. Due to the adaption in the numerical integration, we then obtain reliable and precise results for values of \(m_{e\gamma} \lesssim m_{\text{cut}}\).

Figure 3 shows the electron energy distribution, again at LO (green) and NLO (blue). As for \(m_{e\gamma}\) the corrections are negative and amount to roughly \(10\%\). Since this plot is linear, they can be clearly seen by comparing LO and NLO. In the lower panel once more the \(K\) factor is depicted. Unsurprisingly, at the very end of the distribution, \(E_e\sim 900\,{\rm MeV}\), the statistics is out of control.

Figure 3 Results of the toy run to compute \(E_e\) for \(\tau\to\nu\bar\nu e\gamma\). Upper panels show the LO (green) and NLO (blue) results, the lower panels show the NLO K factor.

More complicated runs

To demonstrate some of McMule capabilities, we tweak the observable a bit. Since the tau is usually produced in \(e^+e^-\to \tau^+\tau^-\), we instead use the LO cross section

(1)\[ \frac{{\rm d}\sigma}{{\rm d}(\cos\theta)} \propto \Big(1+\frac{4m_\tau^2}{s}\Big) + \Big(1+\frac{4m_\tau^2}{s}\Big) \cos\theta\]

with \(\sqrt s = m_{\Upsilon(4S)} = 10.58\,{\rm GeV}\).

To accurately simulate this situation, we need to

• choose a random value for \(\theta\),

• construct the tau momentum \(p_1\) in the lab frame,

• boost the momenta from McMule into this frame, and

• apply a correction weight from (1)

for every event. We require the following cuts in the lab frame

• the produced electron and hard photon have \(-0.75 \le \cos\theta_{i,e^-} \le +0.95\)

• the hard photon energy is bigger than 220 MeV

Further, we want to have a switch for inclusive and exclusive measurements without having to adapt the user file.

Asking for user input

To be able to switch cuts on and off, we need to read input from the user at runtime. This can be done in the inituser() routine where input can be read. We can store the result in a global variable (exclusiveQ) so we can later use it in quant(). Further, we need modify the name of the vegas file by changing filenamesuffix. It is also good practice to print the configuration chosen for documentation purposes.

  SUBROUTINE INITUSER
  read*, exclusiveQ

  if(exclusiveQ == 1) then
    print*, "Calculating tau->e nu nu gamma in ee->tau tau exclusive"
    filenamesuffix = "e"
  else
    print*, "Calculating tau->e nu nu gamma in ee->tau tau inclusive"
    filenamesuffix = "i"
  endif

  ! Let the tau be unpolarised
  pol1 = (/ 0._prec, 0._prec, 0._prec, 0._prec /)
  END SUBROUTINE

Note

When using the menu file system, this can only be a single integer. To read multiple bits of information, you need to encode the data somehow.

Generation of the tau momentum

We can use the user integration feature of McMule to generate \(\cos\theta\). This allows us to write

\[\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\). Since (1) is sufficiently simple, we will sample \(\cos\theta\) with a uniform distribution and apply a correction weight rather trying to sample it directly. We set the variable userdim to one to indicate that we want to carry out \(m=1\) extra integrations and define the function userevent() that sets global variable cth for \(\cos\theta\)

  SUBROUTINE USEREVENT(X, NDIM)
  integer :: ndim
  real(kind=prec) :: x(ndim)

  cth = 2*x(1) - 1
  userweight = (1+Mm**2/Etau**2) + (1-Mm**2/Etau**2) * cth
  END SUBROUTINE USEREVENT

Warning

This function can be used to change the centre-of-mass energy and masses of the particles. However, one must the re-compute the flux factors and \(\xi_{\text{max}}\) relations.

Boosting into the lab frame

We begin by writing down the momentum of tau in the lab frame as

\[p_1 = \Big( 0, |\vec p| \sqrt{1-\cos\theta^2}, |\vec p| \cos\theta, E\Big)\]

with \(|\vec p| = \sqrt{E^2-m_\tau^2}\). Next, we use the McMule function boost_back() to boost the momenta we are given into the lab frame. From there we can continue applying our cuts as before, utilising the McMule function cos_th to calculate the angle between the particle and the beam axis.

  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) :: ptau, cos_e, cos_g
  real (kind=prec) :: p1Lab(4), p2Lab(4), p5Lab(4), p6Lab(4)
  real (kind=prec) :: quant(nr_q)
  real (kind=prec) :: gs(4), gh(4)
  real (kind=prec), parameter :: ez(4) = (/ 0., 0., 1., 0. /)
  !! ==== keep the line below in any case ==== !!
  call fix_mu
  pass_cut = .true.

  ptau = sqrt(Etau**2-Mtau**2)
  p1Lab = (/ 0., ptau*sqrt(1-cth**2), ptau*cth, Etau /)

  p1Lab = boost_back(p1Lab, q1)
  p2Lab = boost_back(p1Lab, q2)
  p5Lab = boost_back(p1Lab, q5)
  p6Lab = boost_back(p1Lab, q6)

  if (p5Lab(4) > p6Lab(4)) then
    gh = p5Lab ; gs = p6Lab
  else
    gh = p6Lab ; gs = p5Lab
  endif

  cos_e = cos_th(p2Lab, ez)
  cos_g = cos_th(gh   , ez)

  if ( (cos_e > 0.95 .or. cos_e < -0.75) .or. (cos_g > 0.95 .or. cos_g < -0.75) ) then
    pass_cut = .false.
    return
  endif

  if(exclusiveQ == 1) then
    if (gh(4) < 220. .or. gs(4) > 220.) then
      pass_cut = .false.
      return
    endif
  else
    if (gh(4) < 220.) then
      pass_cut = .false.
      return
    endif
  endif

  names(1) = 'minv'
  quant(1) = sqrt(sq(q2+gh))

  names(2) = 'Ee'
  quant(2) = p2Lab(4)

  names(3) = 'cos_e'
  quant(3) = cos_e

  names(4) = 'cos_g'
  quant(4) = cos_g
  END FUNCTION QUANT

Running and analysis

At this point we can run McMule and proceed with the analysis as before. We need to do two runs, one for the exclusive and one for the inclusive. However, only the real corrections differ, therefore we only need 24 runs and not 36. The last argument of the run command in the menu file will be passed as the observable we have defined in inituser(). We need to pass 1 (exclusive) or 0 (inclusive) as shown in Listing 8 [2]

Listing 8 The menu file for the present calculation

image registry.gitlab.com/mule-tools/mcmule:redesign example3/user.f95

conf example3/m2enng-tau-e.conf

run 75217 1.000000 m2enng0 tau-e 0
run 52506 1.000000 m2enng0 tau-e 0
run 22671 1.000000 m2enng0 tau-e 0

run 53796 1.000000 m2enngV tau-e 0
run 15282 1.000000 m2enngV tau-e 0
run 89444 1.000000 m2enngV tau-e 0

run 98870 0.600000 m2enngC tau-e 0
run 91991 0.600000 m2enngC tau-e 0
run 79769 0.600000 m2enngC tau-e 0

run 21175 0.800000 m2enngC tau-e 0
run 57581 0.800000 m2enngC tau-e 0
run 81929 0.800000 m2enngC tau-e 0

run 70604 0.600000 m2enngR tau-e 0
run 33013 0.600000 m2enngR tau-e 0
run 22530 0.600000 m2enngR tau-e 0

run 82222 0.800000 m2enngR tau-e 0
run 30935 0.800000 m2enngR tau-e 0
run 40689 0.800000 m2enngR tau-e 0


run 70604 0.600000 m2enngR tau-e 1
run 33013 0.600000 m2enngR tau-e 1
run 22530 0.600000 m2enngR tau-e 1

run 82222 0.800000 m2enngR tau-e 1
run 30935 0.800000 m2enngR tau-e 1
run 40689 0.800000 m2enngR tau-e 1

We can now run McMule. When analysing the output we need to take care to not mix the different observables which we do by passing the optional argument obs to sigma(). The resulting plot is shown in :numref:fig_Eboost:

Listing 9 The analysis pipeline for this calculation

# Loading the LO is the same as before
setup(folder='example3/out.tar.bz2')
LO = scaleset(mergefks(sigma('m2enng0')), GF**2*lifetime*alpha)

# Import the excl. NLO by specifying the observable
# for the real corrections
NLOexcl = scaleset(mergefks(
    sigma('m2enngR', obs='e'),
    sigma('m2enngC'),
    anyxi=sigma('m2enngV')
), GF**2*lifetime*alpha**2)
fullNLOexcl = addsets([LO, NLOexcl])

NLOincl = scaleset(mergefks(
    sigma('m2enngR', obs='i'),
    sigma('m2enngC'),
    anyxi=sigma('m2enngV')
), GF**2*lifetime*alpha**2)
fullNLOincl = addsets([LO, NLOincl])

print("BR_0 = ", printnumber(LO['value']))
print("BRexcl = ", printnumber(fullNLOexcl['value']))
print("BRincl = ", printnumber(fullNLOincl['value']))

fig3, (ax1, ax2) = kplot(
    {
        'lo': scaleplot(LO['Ee'], 1e3),
        'nlo': scaleplot(NLOexcl['Ee'], 1e3),
        'nlo2': scaleplot(NLOincl['Ee'], 1e3)
    },
    labelx=r"$E_e\,/\,{\rm GeV}$",
    labelsigma=r"$\D\mathcal{B}/\D E_{e}$",
    legend={
        'lo': '$\\rm LO$',
        'nlo': '$\\rm NLO\ exclusive$',
        'nlo2': '$\\rm NLO\ inclusive$'
    },
    legendopts={'what': 'l', 'loc': 'lower left'}
)
ax2.set_ylim(-0.12,0.02)

Figure 4 Results of the toy run to compute \(E_{e}\) in the labframe.

[1]

The name menu was originally used by the cryptanalysts at Bletchley Park to describe a particular set of configurations for the ‘computer’ to try

[2]

Looking at the results from our previous run, we can deduce that \(\xi_c\sim0.7\) is the optimal place for running since \(\sigma_n(\xi_c=0.7)\sim \sigma_{n+1}(\xi_c=0.7)\) which reduces cancellation between the different pieces. This optimisation is not strictly necessary and we still run for two values of \(\xi_c\), 0.6 and 0.8.