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.mcmulewill ensure that the size of this array is correct. Exposed asmcmule_pass_cutin C.
- names (nr_q) [character(len=namesLen)]
The names of the distributions the user wishes the calculate.
mcmulewill ensure that the size of this array is correct Exposed asmcmule_namesin 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_musqin C.
- pol1 (4) [real(kind=prec)]
The polarisation of the first polarised particle. Exposed as
mcmule_pol1in C.
- pol2 (4) [real(kind=prec)]
The polarisation of the second polarised particle Exposed as
mcmule_pol2in C.
- filenamesuffix [character(len=filenamesuffixLen)]
The observable-specific suffix to the vegas file Exposed as
mcmule_filenamesuffixin C.
- userweight [real(kind=prec)]
The weight the user wishes to attach to a given event, may be set by
user.userevent(). Exposed asmcmule_userweightin 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_Gfin 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_alphain 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_sw2in 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
mparticle, usually the muon mass but ifflavouris eg.tau-ethe tau mass. Exposed asmcmule_Mmin C.
- Me [real(kind=prec),protected]
The actual value of the
eparticle, usually the electron mass but ifflavouris eg.tau-muthe muon mass. Exposed asmcmule_Mein C.
- Mt [real(kind=prec),protected]
The actual value of the
tparticle, usually the tau mass Exposed asmcmule_Mtauin C.
- scms [real(kind=prec),protected]
The numerical value of the centre-of-mass energy Exposed as
mcmule_scmsin 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
flavourscms [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_Mjin 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_CLjandmcmule_CRjin 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
Q2is the momentum transfer \(-(p_e-p_p)^2\) andGelandGmagare 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
mofrom the frame whererecis at rest to the frame whererecis 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
moto (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 Statistics).
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_pieceto calculateflav [character(len=*)] :: the
mcmule/flavourto 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/nqnumber_bins [integer] :: The number of bins in the distributions to calculate, similar to
user/nbinslower_bounds (number_hist) [real(kind=prec),target] :: The lower bounds of the distributions, similar to
user/min_valupper_bounds (number_hist) [real(kind=prec),target] :: The upper bounds of the distributions, similar to
user/max_valmeasurement_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/userdimnames_length [integer] :: the maximally allowed length of the histogram
mcmule/names, similar touser/namesLenfilenamesuffix_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_Nelin 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_Nmuin 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_Ntauin 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_Nhadin 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_qin 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_NAMEmacro 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
NULLto 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,
0for \({\rm d}\sigma/{\rm d}Q\) and1for \(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_pieceIt needs to at least setmcmule/flavourthroughmcmule/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/filenamesuffixvariable 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/userweightto 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_cutbut also returns the values of the observables that are to be computed. It also fixes fix the renormalisation scalemusqthough this could be done elsewhere. If the user wishes to consider polarised scattering,mcmule/pol1andmcmule/pol2need 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,
0for \({\rm d}\sigma/{\rm d}Q\) and1for \(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_pieceandmcmule/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_filenamesuffixvariable 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_userweightto 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_cutbut also returns the values of the observables that are to be computed. It also fixes fix the renormalisation scalemcmule_musqthough this could be done elsewhere. If the user wishes to consider polarised scattering,mcmule_pol1andmcmule_pol2need 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") }