libmcmule user guide
McMule’s Fortran code has hundreds of functions and subroutine and we will not document all of them here.
However, we will list those exposed by the mcmule
module.
McMule uses a custom type for real numbers that usually defaults to double precision and is used in the following
- 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
Momenta in McMule are of the form (/ px, py, pz, E /)
though it is often advisable to instead use the functions defined in Section Scalar quantities.
Parameters, variables, and functions of libmcmule
The main module of McMule, mcmule
, exposes the following parameters, variables, and functions to the user
Some parameters are denoted as fixed (cannot be changed without recompiling McMule) or protected (can be changed only through the use of functions)
Variables to be set by the user file
- 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.mcmule
will ensure that the size of this array is correct. Exposed asmcmule_pass_cut
in C.
- names (nr_q) [character(len=namesLen)]
The names of the distributions the user wishes the calculate.
mcmule
will ensure that the size of this array is correct Exposed asmcmule_names
in C.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.
- musq [real(kind=prec)]
The renormalisation scale \(\mu^2\). This variable needs to be set by the user, otherwise McMule will fail. Exposed as
mcmule_musq
in C.
- pol1 (4) [real(kind=prec)]
The polarisation of the first polarised particle. Exposed as
mcmule_pol1
in C.
- pol2 (4) [real(kind=prec)]
The polarisation of the second polarised particle Exposed as
mcmule_pol2
in C.
- filenamesuffix [character(len=filenamesuffixLen)]
The observable-specific suffix to the vegas file Exposed as
mcmule_filenamesuffix
in C.
- userweight [real(kind=prec)]
The weight the user wishes to attach to a given event, may be set by
user.userevent()
. Exposed asmcmule_userweight
in C.
Couplings and masses
Couplings are (usually) hard-coded in McMule and cannot be changed.
Masses can be changed using the initflavour()
routine.
- GF [real(kind=prec),fixed]
The Fermi constant. For predominantly historic reasons, this is set to
1._prec
. Exposed asmcmule_Gf
in C.
- alpha [real(kind=prec),fixed]
The fine-structure constant in the OS scheme. For predominantly historic reasons, this is set to
1._prec
. Exposed asmcmule_alpha
in C.
- sw2 [real(kind=prec)]
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. Exposed as
mcmule_sw2
in C.
- Mel [real(kind=prec),fixed]
The numerical value of the electron mass in MeV, irregardless of the
flavour
- Mmu [real(kind=prec),fixed]
The numerical value of the muon mass in MeV, irregardless of the
flavour
- Mtau [real(kind=prec),fixed]
The numerical value of the tau mass in MeV, irregardless of the
flavour
- Mproton [real(kind=prec),fixed]
The numerical value of the proton mass in MeV, irregardless of the
flavour
- MZ [real(kind=prec),fixed]
The numerical value of the Z boson mass in MeV
- Mm [real(kind=prec),protected]
The actual value of the
m
particle, usually the muon mass but ifflavour
is eg.tau-e
the tau mass. Exposed asmcmule_Mm
in C.
- Me [real(kind=prec),protected]
The actual value of the
e
particle, usually the electron mass but ifflavour
is eg.tau-mu
the muon mass. Exposed asmcmule_Me
in C.
- Mt [real(kind=prec),protected]
The actual value of the
t
particle, usually the tau mass Exposed asmcmule_Mtau
in C.
- scms [real(kind=prec),protected]
The numerical value of the centre-of-mass energy Exposed as
mcmule_scms
in C.
- subroutine initflavour(flavour, scms)
This sets the values of the masses and centre-of-mass energy as per Section Flavour.
- Parameters:
flavour [character(len=*)] :: the flavour to use, will be written into
flavour
scms [real(kind=prec),optional] :: the Mandelstam \(s = (p_1+p_2)^2\) to use.
Warning
Calling this function at run-time is possible (eg. during
userevent()
) but requires updating ofmasses
(cf. Issue#57)
- Mj [real(kind=prec)]
The mass of a New Physics particle. Exposed as
mcmule_Mj
in C.Warning
Calling this function at run-time is possible (eg. during
userevent()
) but requires updating ofmasses
(cf. Issue#57)
- CLj, CRj
The values of the left- and right-handed \(J\bar\psi\psi\) coupling Exposed as
mcmule_CLj
andmcmule_CRj
in C.
- the_proton_ff [procedure,pointer]
A function pointer to the proton form factor. This needs to be of the form
subroutine proton_formfactor(Q2, Gel, Gmag) real(kind=prec), intent(in) :: Q2 real(kind=prec), intent(out) :: Gel, Gmag ... end subroutine proton_formfactor
where
Q2
is the momentum transfer \(-(p_e-p_p)^2\) andGel
andGmag
are the Sachs form factors.
- lambda [real(kind=prec)]
When using a dipole or monopole proton form factor, this is the width of the dipole in MeV<sup>2</sup>.
- kappa [real(kind=prec)]
When using a dipole or monopole proton form factor, this the proton’s magnetic dipole moment
Scalar quantities
- function s(p1, p2) [real(kind=prec)]
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 sq(p) [real(kind=prec)]
The Lorentz square \(p^2\)
- Parameters:
p (4) [real(kind=prec)] :: the momentum \(p\)
- function asymtensor(p1, p2, p3, p4) [real(kind=prec)]
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 cos_th(p1, p2) [real(kind=prec)]
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
- function eta(p) [real(kind=prec)]
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 rap(p) [real(kind=prec)]
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 pt(p) [real(kind=prec)]
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 absvec(p) [real(kind=prec)]
The length of the three-vector part \(|\vec p\)
- Parameters:
p (4) [real(kind=prec)] :: the momentum \(p\)
- function phi(p) [real(kind=prec)]
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 rij(p1, p2) [real(kind=prec)]
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\)
Transformations
- function boost_back(rec, mo)
boosts the momentum
mo
from the frame whererec
is at rest to the frame whererec
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 ofrec
, 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
Main routine
If the user wishes to call McMule from their own executable rather than using the default mcmule
, they may use
- subroutine runmcmule(ncall_ad, itmx_ad, ncall, itmx, initial_ran_seed, xi1, xi2, piece, flav, tfilename)
This function initialises McMule and performs the integration. Exposed as
mcmule_runmcmule()
in C.- Parameters:
ncall_ad [integer] :: number of calls per iteration during pre-conditioning
itmxl_ad [integer] :: number of iterations during pre-conditioning. No pre-conditioning is performed if this value is zero
ncall [integer] :: number of calls per iteration during the main run
itmxl [integer] :: number of iterations during the main run (cf. Section sec_stat).
initial_ran_seed [integer] :: the initial random seed (cf. Section Random number generation).
xi1 [real(kind=prec)] :: the \(\xi\) values, either used for subtraction or the integrated eikonal. It is recommended that they are equal.
piece [character(len=*)] :: the
mcmule/which_piece
to calculateflav [character(len=*)] :: the
mcmule/flavour
to usetfilename [character(len=*),optional] :: the filename to use for the output, autogenerated if not present
- subroutine load_quant(filename)
Loads the observable defined in filename using
dl_open
- Parameters:
filename [character(len=*)] :: the library to open
- subroutine set_observable(number_hist, number_bins, lower_bounds, upper_bounds, measurement_function, user_integration, user_initialisation, user_integration_dimension, names_length, filenamesuffix_length)
Set the observable to an already loaded function.
- Parameters:
number_hist [integer] :: The number of distributions to calculate, similar to
user/nq
number_bins [integer] :: The number of bins in the distributions to calculate, similar to
user/nbins
lower_bounds (number_hist) [real(kind=prec),target] :: The lower bounds of the distributions, similar to
user/min_val
upper_bounds (number_hist) [real(kind=prec),target] :: The upper bounds of the distributions, similar to
user/max_val
measurement_function [quantfunc] :: a pointer to the measurement function like
user/quant()
user_integration [usereventfunc,optional] :: a pointer to the user event function like
user/userevent()
user_initialisation [inituserfunc,optional] :: a pointer to the user initialisation
user_integration_dimension [integer,optional] :: the number of integrations the user wishes to carry out to account eg. for beam effects, similar to
user/userdim
names_length [integer] :: the maximally allowed length of the histogram
mcmule/names
, similar touser/namesLen
filenamesuffix_length [integer,optional] :: the maximally allowed length of the observable name as specified in
mcmule/filenamesuffix
, similar touser/filenamesuffixLen
Technical parameters
Warning
The following parameters are exposed for the user to view for finer control over the measurement function. Extreme care must be taken when using these.
- nel [integer]
Set to 1 if electron VP loops are to be included, set to 0 otherwise. More options may be added later. Exposed as
mcmule_Nel
in C.
- nmu [integer]
Set to 1 if muon VP loops are to be included, set to 0 otherwise. More options may be added later. Exposed as
mcmule_Nmu
in C.
- ntau [integer]
Set to 1 if tau VP loops are to be included, set to 0 otherwise. More options may be added later. Exposed as
mcmule_Ntau
in C.
- nhad [integer]
Set to 1 if HVP loops are to be included, set to 0 otherwise. More options may be added later. Exposed as
mcmule_Nhad
in C.
- which_piece [character(len=25),protected]
The piece being integrated, cf. Section Available processes and which_piece and which_piece. This variable is read-only and is set by McMule.
- flavour [character(len=15),protected]
The flavour configuration being used. This can only be changed using
initflavour()
The following parameters should not be modified by the user unless especially advised to do so
- softcut [real(kind=prec),advanced]
The value of \(\xi\) below which the integrand is set to zero without subtraction
- collcut [real(kind=prec),advanced]
The value of \(\cos\theta\) below which the integrand is set to zero
- sSwitch [real(kind=prec),advanced]
The value of \(\xi\) below which the matrix element is approximated at LP. This is only available for some matrix elements
- ntsSwitch [real(kind=prec),advanced]
The value of \(\xi\) below which the matrix element is approximated at NLP. This is only available for some matrix elements
- pcSwitch [real(kind=prec),advanced]
The value of \(1-y\) below which the matrix element is approximated. This is only available for some matrix elements
- xicut1 [real(kind=prec),protected]
The \(\xi_c\) value below which to subtract the first photon
- xicut2 [real(kind=prec),protected]
The \(\xi_c\) value below which to subtract the second photon
- xieik1 [real(kind=prec),protected]
The argument to the first integrated eikonal
- xicut2 [real(kind=prec),protected]
The argument to the second integrated eikonal
- zero [real(kind=prec),fixed]
Real numbers below this value may be truncated
- pi [real(kind=prec),fixed]
- xiout [protected]
The current values of \(\xi\) for single-real corrections.
- xioutA, xioutB
The current values of \(\xi_i\) for double-real corrections.
- yout [protected]
The current values of \(y\) for single-real corrections.
- youtA, youtB
The current values of \(y_i\) for double-real corrections.
- qsoft [protected]
The momentum of the soft photon without \(\xi\) for single-real corrections.
- qsoftA, qsoftB
The momenta of the soft photons without \(\xi_i\) for double-real corrections.
- nr_q [protected]
This variable is obtained by reading
user.nq
. Exposed asmcmule_nr_q
in C;
- function sachs_gel(q2) [real(kind=prec)]
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 sachs_gmag(q2) [real(kind=prec)]
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\)
Symbols exposed to C/C++
McMule is inherently a Fortran code. However, it can be used with a user file written in C/C++ even though many of the quality-of-life functions defined in Sections Scalar quantities and Transformations are not available.
-
type real
The real number type used in McMule
-
char *mcmule_names
Preprocessor alias for
names
. It is recommened to use theMCMULE_SET_NAME
macro to interact with this variable.
-
char *mcmule_filenamesuffix
Preprocessor alias, points to the first element of
filenamesuffix
.
-
proton_ff mcmule_protonff
Preprocessor alias for
the_proton_ff
.
-
int mcmule_Nel
Preprocessor alias for
Nel
.
-
int mcmule_Nmu
Preprocessor alias for
Nmu
.
-
int mcmule_Ntau
Preprocessor alias for
Ntau
.
-
int mcmule_Nhad
Preprocessor alias for
Nhad
.
-
MCMULE_SET_NAME(q, name)
Sets the q-th observable to the name
name
-
void mcmule_initflavour(char *flavour, double *scms)
Interface fo f:subr:initflavour.
- Parameters:
flavour (char*) – the flavour combnination to use
scms (double*) – a pointer to the centre-of-mass energy. Passing a null pointer will is only allowed for decay processes or when using a flavour preset.
-
void mcmule_runmcmule(int ncall_ad, int itmx_ad, int ncall, int itmx, int initial_ran_seed, real xicut1, real xicut2, char *piece, char *flav, char *tfilename);
Interface for
runmcmule()
.- Parameters:
tfilename (char*) – set to
NULL
to use default
-
void mcmule_set_observable(int number_hist, int number_bins, real *lower_bounds, real *upper_bounds, quantfunc measurement_function, usereventfunc user_integration, inituserfunc user_initialisation, int user_integration_dimension, int names_length, int filenamesuffix_length)
Interface for
set_observable()
.
The user file in Fortran
A Fortran user file must be in the form of a module user
that defines the following symbols.
Warning
None of the values may be parameters
in Fortran as they are not exposed to libmcmule
- namesLen [integer,default = 6]
The maximally allowed length of the histogram
mcmule/names
.
- filenamesuffixLen [integer,default = 10]
The maximally allowed length of the observable name as specified in
mcmule/filenamesuffix
.
- nq [integer,default=1]
The number of distributions the user intends to calculate.
- nbins [integer,default=100]
The number of bins in the distributions the user intends to calculate
- userdim [integer,default=0]
The number of integrations the user wishes to carry out to account eg. for beam effects
- bin_kind [integer,default=0]
The binning mechanism being used,
0
for \({\rm d}\sigma/{\rm d}Q\) and1
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
- subroutine inituser()
If defined, this is called without arguments once as soon as McMule starts and has read all other configuration, meaning that it can access
mcmule/which_piece
It needs to at least setmcmule/flavour
throughmcmule/initflavour()
. 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
mcmule/filenamesuffix
variable which is appended to the name of the VEGAS file.Example for reading a cut and setting the flavour:
SUBROUTINE INITUSER integer cut call initflavour("mu-e", 1000.) read*,cut write(filenamesuffix,'(I2)') cut END SUBROUTINE INITUSER
with a global variable
cut
- subroutine userevent(x, ndim)
If defined, the user may use this routine in combination with
mcmule/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 equaluserdim
.
- function quant(q1, q2, q3, q4, q5, q6, q7)
The measurement function the user wishes to calculate. This needs to at least set
mcmule/pass_cut
but also returns the values of the observables that are to be computed. It also fixes fix the renormalisation scalemusq
though this could be done elsewhere. If the user wishes to consider polarised scattering,mcmule/pol1
andmcmule/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 ==== !! musq = me**2 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
The user file in C/C++
A C/C++ user file must define the following symbols. The defaults are identical to Fortran.
-
int mcmule_namelen
The maximally allowed length of the histogram
mcmule_names
.
-
int mcmule_filename_suffix_length
The maximally allowed length of the observable name as specified in
mcmule_filenamesuffix
.
-
int mcmule_number_hist
The number of distributions the user intends to calculate.
-
int mcmule_number_bins
The number of bins in the distributions the user intends to calculate
-
int mcmule_user_integration_dimension
The number of integrations the user wishes to carry out to account eg. for beam effects
-
int mcmule_bin_kind
The binning mechanism being used,
0
for \({\rm d}\sigma/{\rm d}Q\) and1
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
-
void mcmule_user_initialisation()
If defined, this is called without arguments once as soon as McMule starts and has read all other configuration, meaning that it can access
mcmule/which_piece
andmcmule/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
mcmule_filenamesuffix
variable which is appended to the name of the VEGAS file.
-
void mcmule_user_integration(real *x, int ndim)
If defined, the user may use this routine in combination with
mcmule_userweight
to integrate over further parameters.- Parameters:
x (real*) – a pointer to the first element of the values of the integration
ndim (int) – the dimension of
x
, should equalmcmule_user_integration_dimension
.
-
void mcmule_measurement_function(double **res, real *q1, real *q2, real *q3, real *q4, real *q5, real *q6, real *q7)
The measurement function the user wishes to calculate. This needs to at least set
mcmule_pass_cut
but also returns the values of the observables that are to be computed. It also fixes fix the renormalisation scalemcmule_musq
though this could be done elsewhere. If the user wishes to consider polarised scattering,mcmule_pol1
andmcmule_pol2
need to be set. This may use the macroMCMULE_SET_NAME
- Parameters:
void mcmule_measurement_function(double**res,double*p1,double*p2,double*p3,double*p4,double*p5,double*p6,double*p7) { mcmule_musq = 1.; res[0][0] = p3[2] / p3[3]; // Set the second quant to p3z/E3 res[0][1] = p4[2] / p4[3]; // Set the second quant to p4z/E4 mcmule_pass_cut[0] = 1; // Accept event mcmule_pass_cut[1] = 1; // Accept event MCMULE_SET_NAME(0, "p3z/E3") MCMULE_SET_NAME(1, "p4z/E4") }