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 as mcmule_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 as mcmule_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 as mcmule_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 as mcmule_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 as mcmule_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 if flavour is eg. tau-e the tau mass. Exposed as mcmule_Mm in C.

Me [real(kind=prec),protected]

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

Mt [real(kind=prec),protected]

The actual value of the t particle, usually the tau mass Exposed as mcmule_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 of masses (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 of masses (cf. Issue#57)

CLj, CRj

The values of the left- and right-handed \(J\bar\psi\psi\) coupling Exposed as mcmule_CLj and mcmule_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\) and Gel and Gmag 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 where rec is at rest to the frame where rec is specified, i.e.

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

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

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

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

Return:

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

function  boost_rf(rec, mo)

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

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

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

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

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

Return:

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

function  euler_mat(a, b, c)

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

\[\begin{split}\begin{pmatrix} c_{\alpha}c_{\beta }c_{\gamma}-s_{\alpha}s_{\gamma}&-c_{\alpha}c_{\beta } s_{\gamma}-c_{\gamma}s_{\alpha}&c_{\alpha}s_{\beta }&0\\ c_{\beta }c_{\gamma}s_{\alpha}+c_{\alpha}s_{\gamma}& c_{\alpha}c_{\gamma}-c_{\beta } s_{\alpha}s_{\gamma}&s_{\alpha}s_{\beta }&0\\ -c_{\gamma}s_{\beta } & s_{\beta }s_{\gamma} &c_{\beta } &0\\ 0 & 0 & 0 &1 \end{pmatrix}\end{split}\]
Parameters:
  • a [real(kind=prec)] :: the angle \(\alpha\)

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

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

Return:

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

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 calculate

  • flav [character(len=*)] :: the mcmule/flavour to use

  • tfilename [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 to user/namesLen

  • filenamesuffix_length [integer,optional] :: the maximally allowed length of the observable name as specified in mcmule/filenamesuffix, similar to user/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 as mcmule_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

bool *mcmule_pass_cut

Preprocessor alias for pass_cut.

char *mcmule_names

Preprocessor alias for names. It is recommened to use the MCMULE_SET_NAME macro to interact with this variable.

real mcmule_musq

Preprocessor alias for musq.

real *mcmule_pol1

Preprocessor alias, points to the first element of pol1.

real *mcmule_pol2

Preprocessor alias, points to the first element of pol2.

char *mcmule_filenamesuffix

Preprocessor alias, points to the first element of filenamesuffix.

real *mcmule_userweight

Preprocessor alias, points to the user weight

real mcmule_Gf

Preprocessor alias for GF.

real mcmule_alpha

Preprocessor alias for alpha.

real mcmule_sw2

Preprocessor alias for sw2.

real mcmule_Mm

Preprocessor alias for Mm.

real mcmule_Me

Preprocessor alias for Me.

real mcmule_Mtau

Preprocessor alias for Mt.

real mcmule_scms

Preprocessor alias for scms.

real mcmule_Mj

Preprocessor alias for Mj.

real mcmule_CLj

Preprocessor alias for CLj.

real mcmule_CRj

Preprocessor alias for CRj.

proton_ff mcmule_protonff

Preprocessor alias for the_proton_ff.

real mcmule_protonff_lambda

Preprocessor alias for lambda.

real mcmule_protonff_kappa

Preprocessor alias for kappa.

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.

int mcmule_nr_q

Preprocessor alias for nr_q.

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\) and 1 for \(Q {\rm d}\sigma/{\rm d}Q\).

Warning

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

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

The lower bounds of the distributions

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

The upper bounds of the distributions

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 set mcmule/flavour through mcmule/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 equal userdim.

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 scale musq though this could be done elsewhere. If the user wishes to consider polarised scattering, mcmule/pol1 and mcmule/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\) and 1 for \(Q {\rm d}\sigma/{\rm d}Q\).

Warning

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

real mcmule_lower_bounds[nr_q]

The lower bounds of the distributions

real mcmule_upper_bounds[nr_q]

The upper bounds of the distributions

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 and mcmule/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:
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 scale mcmule_musq though this could be done elsewhere. If the user wishes to consider polarised scattering, mcmule_pol1 and mcmule_pol2 need to be set. This may use the macro MCMULE_SET_NAME

Parameters:
  • res (real**) – a double pointer to the result.

  • qi (real*) – pointers to the first component of the momenta

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")
}