.. f:currentmodule:: mcmule .. _sec_getting_started: 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 :math:`\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 :term:`BR` for this process, as measured by BaBar :cite:`Lees:2015gea, Oberhof:2015snl` has a discrepancy of more than :math:`3\,\sigma` from the :term:`SM` value. This will illustrate the importance of fully differential :term:`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 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 :ref:`sec_compilation` on how to compile the code yourself. macOS users will also have to manually compile McMule, cf. Section :ref:`sec_compilation`. First, obtain the McMule distribution from `our website `_ and extract the tarball .. code:: bash $ 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 .. code:: bash $ cp mcmule-release/mcmule /usr/local/bin $ cp mcmule-release/mcmule.mod /usr/local/include .. _sec_getting_started_deps: Installing dependencies ^^^^^^^^^^^^^^^^^^^^^^^ Even if you do not plan to develop McMule, you will need a GFortran compiler to specify your observables. If you have not already done so, you can install it with .. code:: bash sudo apt install gfortran # on Debian derivatives such as Ubuntu sudo dnf install gcc-gfortran # on RHEL derivatives sudo brew install gcc # on macOS, requires homebrew to be installed 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 .. code:: bash $ pip3 install pymule --index-url https://gitlab.com/api/v4/projects/19995364/packages/pypi/simple Please note that most modern linux distributions have implemented `PEP 668 `_ and will refuse to install packages using pip3 outside of a virtual environment. If the above command fails, you should create a `virtual environment `_ to work in .. code:: bash $ python3 -m venv mcmule-environment $ source mcmule-environment/bin/activate $ pip3 install pymule --index-url https://gitlab.com/api/v4/projects/19995364/packages/pypi/simple You will have to run the ``source`` command every time you start working with McMule. We apologise for the inconvenience this causes. .. _sec_getting_started_lo: Simple runs at LO ----------------- Setting McMule up ^^^^^^^^^^^^^^^^^ In this example we want to compute two distributions, the invariant mass of the :math:`e\gamma` pair, :math:`m_{\gamma e}\equiv \sqrt{(p_e+p_\gamma)^2}`, and the energy of the electron, :math:`E_e`, in the rest frame of the tau. To avoid an :term:`IR` singularity in the :term:`BR`, we have to require a minimum energy of the photon. We choose this to be :math:`E_\gamma \ge 10\,{\rm MeV}` as used in :cite:`Lees:2015gea, Oberhof:2015snl`. At first, we need to find out how the process :math:`\tau\to\nu\bar\nu e\gamma` is implemented in McMule. For this, we refer to the table in Section :ref:`sec_pieces` that specifies the pieces (sometimes called :f:var:`which_piece`) that is required for a :term:`generic process`. The generic process is a prototype for the physical process such as :math:`\ell \to\nu\bar\nu \ell' \gamma` where the flavour of the lepton :math:`\ell` is left open. In our case, we need to consider the row for :math:`\mu\to\nu\bar\nu e\gamma`. Since we are only interested in :term:`LO`, the only :f:var:`which_piece` we need is ``m2enng0``. To change from the generic process :math:`\mu\to\nu\bar\nu e\gamma` to the process we are actually interested in, :math:`\tau\to\nu\bar\nu e\gamma`, we pick the :f:var:`flavour` ``tau-e`` which refers to a :math:`\tau \to e \cdots` transition. Other options would be ``tau-mu`` for :math:`\tau\to\mu \cdots` or ``mu-e`` for :math:`\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 :term:`PID`. We can refer to the table in Section :ref:`sec_pid` to find that for the :f:var:`which_piece` ``m2enng0``, we have .. math:: \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. An empty template can be found on `our website `__. We can use this file to implement the :term:`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 :math:`0 < m_{\gamma e} < 1800\,{\rm MeV}` and :math:`0 \le E_e \le 900\,{\rm MeV}`. .. literalinclude:: mcmule/example1/user.f90 :language: fortran :lineno-start: 10 :lines: 10-15 :caption: The metadata for our calculation with two histograms (``nr_q = 2``) with 90 bins each (``nr_bins = 90``). The ranges should be :math:`0 < m_{\gamma e} < 1800\,{\rm MeV}` and :math:`0 \le E_e \le 900\,{\rm MeV}`. .. note:: Finding suitable values for the ranges can be tricky beyond :term:`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 :term:`measurement function` called :f:func:`quant`. We need to * calculate the invariant mass of the :math:`e\gamma` pair. This is done using the momentum-squaring function :f:func:`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 :f:var:`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.``. .. literalinclude:: mcmule/example1/user.f90 :language: fortran :lineno-start: 63 :lines: 63-80 :caption: The :term:`measurement function` at :term:`LO` Additionally to the numeric value in ``quant(i)`` we store a human-readable name in :f:var:`names`. .. warning:: The maximal length of these names is defined in the variable :f:var:`namesLen` which defaults to 6 characters. Also note that this :term:`measurement function` is not :term:`IR`-safe! Finally, we need to specify the masses of the particles involved. This is called the :f:var:`flavour` in McMule. We do this by defining a routine called :f:subr:`inituser` which calls the McMule function :f:subr:`initflavour`. Additionally, we also use this for some self-documentation .. literalinclude:: mcmule/example1/user.f90 :language: fortran :lineno-start: 55 :lines: 55-61 :caption: The :term:`measurement function` at :term:`LO` This sets the :f:var:`flavour` variable is to ``tau-e`` to change from the generic process :math:`\mu\to\nu\bar\nu e\gamma` to the process we are actually interested in, :math:`\tau\to\nu\bar\nu e\gamma`. We now need to compile our observable into a :term:`shared library` so that McMule can load it. To do this, we run .. code:: bash $ gfortran -fPIC --shared -o user.so user.f95 On macOS you need to run .. code:: bash $ gfortran -fPIC -bundle -undefined dynamic_lookup -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 .. code:: bash $ gfortran -I/path/to/the/folder/of/mcmule.mod/ -fPIC --shared -o user.so user.f95 .. warning:: A version of the ``mcmule.mod`` header file is distributed as part 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. .. note:: When compiling our observable we may have to use ``-fdefault-real-8`` in front of ``user.f95``. 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 :math:`1000\times 10^3` points each for pre-conditioning and 50 iterations with :math:`1000\times 10^3` points each for the actual numerical evaluation (cf. Section :ref:`sec_stat` for some heuristics to determine the statistics needed). We pick a :term:`random seed` between 0 and :math:`2^{31}-1` (cf. Section :ref:`sec_rng`), say 70998, and for the input variable :f:var:`which_piece` we enter ``m2enng0`` as discussed above. This system is used for other processes as well. The input variable :f:var:`which_piece` determines the generic process and the part of it that is to be computed (i.e. tree level, real, double virtual etc.). This means that we need to input the following (the specifications for the input can be found in :numref:`tab_mcmuleinput`): .. warning:: When running ``mcmule`` outside the normal repository, you need to make sure that an ``out/`` folder exists. .. code-block:: $ gfortran -fPIC --shared -o user.so user.f95 $ ./mcmule 1000 10 10000 50 70998 1.0 1.0 m2enng0 ******************************************* * 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 :term:`SHA1` hashes of the source code and the git version) as well as some user-defined information from the subroutine :f:subr:`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. .. _tab_mcmuleinput: .. table:: The options read from ``stdin`` by McMule. The calls are multiplied by 1000. :widths: auto ====================== ============== ================================================================================= **Variable name** **Data type** **Comment** ====================== ============== ================================================================================= ``ncall_ad`` ``integer`` calls / iteration during pre-conditioning ``itmx_ad`` ``integer`` iterations during pre-conditioning ``ncall`` ``integer`` calls / iteration during main run ``itmx`` ``integer`` iterations during main run ``ran_seed`` ``integer`` :term:`random seed` :math:`z_1` ``xinormcut`` ``real(prec)`` the :math:`0<\xi_{c}\le1` parameter ``delcut`` ``real(prec)`` the :math:`\delta_{\text{cut}}` parameter (or at :term:`NNLO` the second :math:`\xi_{c}`) :f:var:`which_piece` ``char(10)`` the part of the calculation to perform (opt) unknown the user can request further input during :f:subr:`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 :ref:`sec_vegasff`). We can open this file in python and make plots. While it is possible to open just a single file using :func:`pymule.vegas.FileRecord.read`, 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 :func:`~pymule.loader.sigma` and :func:`~pymule.loader.mergefks`. .. literalinclude:: mcmule/dummy.py :language: python :linenos: :lines: 1-11 .. warning:: In McMule the numerical value of the Fermi constant :math:`G_F` and the fine-structure constant :math:`\alpha` are set to one for predominately historical reasons. This needs to be restored in python The variable ``dat`` now contains the runtime (``time``), branching ratio (after multiplication with the lifetime, ``value``), and :math:`\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 :math:`N\times 3` matrices .. code:: python 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 :ref:`sec_pymule`. For now, we will just make a plot of the :math:`E_e` distribution ``Ee`` .. literalinclude:: mcmule/dummy.py :language: python :lineno-start: 14 :lines: 14-19 .. figure:: mcmule/dummy.svg Result of the :term:`LO` test run for the :math:`E_e` distribution .. _sec_getting_started_nlo: Running at NLO and beyond ------------------------- A few things change once we go beyond :term:`LO` since we can have extra radiation. To account for this, more :f:var:`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 :ref:`sec_pieces` we find that we need the pieces ``m2enngF`` and ``m2enngR`` for virtual and real corrections respectively. The :term:`PID` table of Section :ref:`sec_pid` 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 :cite:`Lees:2015gea, Oberhof:2015snl` we will consider the exclusive radiative decay, i.e. we request precisely one photon with energy :math:`E_\gamma>10\,{\rm MeV}`. The function :f:func:`quant` will have to take this into account with the additional argument ``q6``, the momentum of the second photon. To ensure :term:`IR` safety, we define the harder and softer photon ``gh`` and ``gs``, respectively, and require that the former (latter) has energy larger (smaller) than :math:`10\,{\rm MeV}`. This new version of :f:func:`quant` is also suitable for the :term:`LO` calculation and it is generally advisable to use a single :f:func:`quant` function for all parts of a computation. .. literalinclude:: mcmule/example2/user.f90 :language: fortran :lineno-start: 64 :lines: 64-87 :caption: The measurement function beyond :term:`LO`. The changes w.r.t. to :term:`LO` are highlighted. :emphasize-lines: 5,9-16 Running McMule ^^^^^^^^^^^^^^ The :term:`FKS` scheme used in McMule introduces an unphysical parameter called :math:`\xi_c` that can be varied between .. math:: 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 :term:`IR` safety of the :term:`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 :math:`\xi_c` and this part is typically much faster in the numerical evaluation. A particularly convenient way to run McMule is using :term:`menu files` [1]_. A :term:`menu file` contains a list of jobs to be computed s.t. the user will only have to vary the :term:`random seed` and :math:`\xi_{c}` by hand as the statistical requirements are defined globally in a :term:`config file`. This is completed by a :term:`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 :cite:`Yoo:2003slurm`. To prepare the run in this way we can use pymule .. literalinclude:: listings/pymule-int.txt :caption: The steps necessary to use pymule to prepare running McMule. Note that numbers listed as seeds are random and hence not reproducible. :language: bash :name: lst_pymulein 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 :term:`menu` and :term:`config` files generated by pymule are shown in :numref:`lst_menu` and :numref:`lst_conf` .. literalinclude:: mcmule/example2/menu-m2enng.menu :caption: :term:`menu file` for the present calculation :language: bash :name: lst_menu .. literalinclude:: mcmule/example2/m2enng.toml :caption: Configuration file for the present calculation :language: toml :name: lst_conf We now need to ensure pymule finds the ``mcmule`` executable, either by copying it into the directory pymule just created or by adapting the path in the :term:`config`. After copying the user library ``user.so`` into the same folder, we can start ``mcmule`` by going into the directory ``example2`` that pymule just created and running .. code:: bash cd example2 cp /path/to/mcmule . cp /path/to/user.so . pymule batch shell menu-m2enng.menu 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 .. code:: bash pymule batch shell -np menu-m2enng.menu When running on a SLURM system, we need to use a SLURM :term:`submission script` , adapted to the needs of the cluster ``mcmule`` is running on .. code:: bash #!/bin/bash #SBATCH --partition= #SBATCH --time=