pymule user guide

This section describes all public functions and classes in pymule.

Working with files

pymule.mergefks(*sets, **kwargs)

performs the FKS merge

Parameters:

  • sets – random-seed-merged results (usually from sigma())

  • binwisechi – bool, optional, default False; if set to True, also return extra distributions containing the \(\chi^2\) of the bin-wise FKS merge. This cannot be used together with anyxi and the result should not be passed to scaleset for obvious reasons.

Returns:

the FKS-merged final set containing cross sections, distributions, and run-time information. The chi2a return is a list of the following

  • the \(\chi^2\) of the FKS merge

  • a list of \(\chi^2\) from previous operations, such as random seed merging or the integration.

Note

Optional argument anyxi (or anything starting with anyxi): Sometimes it is necessary to merge \(\xi_c\)-dependent runs (such as a counter term) and \(\xi_c\)-independent runs (such as the one-loop term). Do not use this together with binwisechi

Example:

Load the LO results for the muon decay using sigma()

>>> mergefks(sigma("m2enn0"))

Load the NLO results

>>> mergefks(sigma("m2ennV"), sigma("m2ennR"))

Load the NNLO results where m2ennNF does not depend on \(\xi_c\)

>>> mergefks(sigma("m2ennFF"), sigma("m2ennRF"), sigma("m2ennRR"), anyxi=sigma("m2ennNF"))
pymule.setup(**kwargs)

sets the default arguemnts for sigma().

Parameters:

  • folder – str, optional; file name, optional; folder or tarball to search for vegas files Initialised to current directory (.).

  • flavour – str, optional; the flavour to load, defaults to everything Initialised to everything, i.e. .*.

  • obs – str, optional; the observable to load (the bit after the O), defaults to everything Initialised to everything, i.e. ''.

  • folderp – str, optional; a regular expression to match directory structures of a tar file, defaults to everything Initialised to everything, i.e. .*.

  • filenames – list, optional; list of files to loads, defaults to all files in folder (recurisvely if tar ball) Initialised to None, meaning everything.

  • merge – dict, optional: a dict of histograms {'name': n} to merge n bins in the histogram name. Initialised to to no merging, i.e. {}

  • types – list of callables, optional; functions that convert the groups matched by r into python objects. Common examples would be int or float. Initialised to [int,float, float] as per McMule filename convention.

  • sanitycheck – callable, optional; a function that, given a vegas dict, whether to include the file in the output (return True) or to skip (return False). Initialised to lambda x : True, i.e. include everything.

  • cache – folder name, optional; if existing folder, use as cache for compressed tarballs

Example:

Setup some folders, ensure that /tmp/mcmule exists

>>> setup(folder="path/to/data.tar.bz2", cachefolder="/tmp/mcmule")
Example:

Restrict observable

>>> setup(obs="3")
Example:

Drop runs with a \(\chi^2 > 10\)

>>> setup(sanitycheck=lambda x : x['chi2a'] < 10)
pymule.sigma(piece, **kwargs)

loads a which_piece and statistically combines the random seed.

Parm piece:

str; which_piece to load

Parameters:

  • folder – str, optional; file name, optional; folder or tarball to search for vegas files Initialised to current directory (.).

  • flavour – str, optional; the flavour to load, defaults to everything Initialised to everything, i.e. .*.

  • obs – str, optional; the observable to load (the bit after the O), defaults to everything Initialised to everything, i.e. ''.

  • folderp – str, optional; a regular expression to match directory structures of a tar file, defaults to everything Initialised to everything, i.e. .*.

  • filenames – list, optional; list of files to loads, defaults to all files in folder (recurisvely if tar ball) Initialised to None, meaning everything.

  • merge – dict, optional: a dict of histograms {'name': n} to merge n bins in the histogram name. Initialised to to no merging, i.e. {}

  • types – list of callables, optional; functions that convert the groups matched by r into python objects. Common examples would be int or float. Initialised to [int,float, float] as per McMule filename convention.

  • sanitycheck – callable, optional; a function that, given a vegas dict, whether to include the file in the output (return True) or to skip (return False). Initialised to lambda x : True, i.e. include everything.

  • cache – folder name, optional; if existing folder, use as cache for compressed tarballs

Returns:

a dict with the tuples of FKS parameters as keys and vegas datasets as values.

Note

Use setup() to set the defaults. Arguments provided here override the defaults

Example:

Load the leading order muon decay

>>> sigma("m2enn0")

Load only observable O3

>>> sigma("m2enn0", obs="3")

Working with errors

pymule.addplots(a, b, sa=1.0, sb=1.0)

adds or subtracts two plots

Parameters:

  • a – Nx3 numpy matrix; the first plot

  • b – Nx3 numpy matrix; the second plot

  • sa – float, optional; the coefficient of the first plot

  • sb – float, optional; the coefficient of the second plot

Returns:

a Nx3 numpy matrix with \(s_a\cdot a + s_b\cdot b\)

Note

a and b must share x values, otherwise entries are dropped

Example:

subtract two plots a and b

>>> addplots(a, b, sb=-1)
Example:

Given the LO plots thetaLO and the NLO corrections thetadNLO, we calculate the \(K\) factor as either

>>> thetaNLO = addplots(thetaLO, thetadNLO)
pymule.chisq(values)

calculates the \(\chi^2/\textrm{d.o.f.}\) of numbers

Parameters:

value – Nx2 numpy matrix or list of lists; the values as [[y1, dy1], [y2, dy2], ...]

Returns:

float; the \(\chi^2/\textrm{d.o.f.} = \frac{1}{n} \sum_{n=1}^n(\frac{y_i-\bar y}{\delta y_i})^2\) with the average value \(\bar y\)

Example:

a good example

>>> chisq([[20.0, 0.8],
...        [21.6, 0.9],
...        [18.7, 1.2]])
1.3348808062205872

and a bad example

>>> chisq([[16.2, 0.8],
...        [22.9, 0.9],
...        [8.81, 1.2]])
30.173852184366673
pymule.dividenumbers(a, b)

divides numbers

Parameters:

  • a – list of floats; the numerator with error [a, da]

  • b – list of floats; the denominator with error [b, db]

Returns:

the result of the division a/b [y, dy]

Example:

Divide \((2.3\pm0.1) / (45\pm0.01)\)

>>> dividenumbers([2.3, 0.1], [45., 0.01])
array([0.05111111, 0.00222225])
pymule.divideplots(a, b, offset=0.0)

divides two plots

Parameters:

  • a – Nx3 numpy matrix; the numerator plot

  • b – Nx3 numpy matrix; the denominator plot

  • offset – float, optional; shifts the result

Returns:

a Nx3 numpy matrix with \(a/b + offset\)

Note

a and b must share x values, otherwise entries are dropped

Example:

Given the LO plots thetaLO and the NLO corrections thetadNLO, we calculate the \(K\) factor as either

>>> thetaNLO = addplots(thetaLO, thetadNLO)
>>> thetaK = divideplots(thetaNLO, thetaLO)
>>> thetaK = divideplots(thetadNLO, thetaLO, offset=+1.)
pymule.integratehistogram(hist)

integrates a histogram

Parameters:

hist – Nx3 numpy matrix; the histogram to integrate \(d\sigma/dx\) as np.array([[x1, y1, e1], [x2, y2, e2], ...])

Returns:

float; the integrated histogram \(\int d\sigma/dx dx\) without error estimate

Example:

Integrate a histogram

>>> hist
array([[          -inf, 0.00000000e+00, 0.00000000e+00],
       [5.00000000e-02, 4.77330751e+01, 2.26798977e-01],
       [1.50000000e-01, 7.40641192e+01, 2.36498021e-01],
       ...,
       [8.85000000e+00, 1.67513948e+00, 1.16218116e-01],
       [8.95000000e+00, 0.00000000e+00, 0.00000000e+00],
       [           inf, 0.00000000e+00, 0.00000000e+00]])
>>> integratehistogram(hist)
4188.519369660588
pymule.mergebins(p, n)

merges n adjacent bins into one larger bin, reducing the uncertainty.

Parameters:

  • p – Nx3 numpy matrix; the plot

  • n – int; how many bins to merge

Returns:

a (N/n)x3 numpy matrix

Note

This process loses len(p)%n bins at the end of the histogram

Example:

merge five bins

>>> len(p)
200
>>> len(mergebins(p, 5))
40

Bins may be lost

>>> len(p)
203
>>> len(mergebins(p, 5))
40
pymule.mergenumbers(values, quiet=False)

statistically combines values with uncertainties

Parameters:

  • values – Nx2 numpy matrix or list of lists; the values as [[y1, dy1], [y2, dy2], ...]

  • quiet – bool, optional; whether to print or return the \(\chi^2\) for the combination

Returns:

either answer as numpy array [y, dy] or tuple of \(chi^2\) and answer

Example:

If quiet is not specified this will print the \(chi^2\)

>>> mergenumbers([[20.0, 0.8],
...               [21.6, 0.9],
...               [18.7, 1.2]])
1.3348808062205872
array([20.30718232,  0.53517179])

Otherwise, it will return it

>>> mergenumbers([[20.0, 0.8],
...               [21.6, 0.9],
...               [18.7, 1.2]], quiet=True)
(1.3348808062205872, array([20.30718232,  0.53517179]))
pymule.mergeplots(ps, returnchi=False)

statistically combines a list of plots

Parameters:

  • ps – list of Nx3 numpy matrices; the plots to combine as [np.array([[x1, y1, e1], [x2, y2, e2], ...]), ...]

  • returnchi – bool, optional; if True returns two plots, the requested combination and the bin-wise \(\chi^2\)

Returns:

a Nx3 numpy matrix if returnchi=False

Example:

Load a number of vegas files and merge them

>>> data = [
...     importvegas(i)['thetae']
...     for i in glob.glob('out/em2em0*')
... ]
>>> mergeplots(data)
pymule.plusnumbers(*args)

adds numbers and errors

Parameters:

yi – list of floats; a number with error [yi, dyi]

Returns:

the result of the addition [y, dy]

Example:

Adding \((10\pm1)+(20\pm0.5)+(-5\pm2)\)

>>> plusnumbers([10, 1], [20, 0.5], [-5, 2])
array([25.        ,  2.29128785])
pymule.printnumber(x, prec=0)

returns a string representation of a number with uncertainties to one significant digit

Parameters:

  • x – a list with two floats; the number as [x, dx]

  • prec – int, otpional; number of extra signficant figures

Returns:

str; the formatted string

Example:

printing \(53.2\pm0.1\) to one significant figure

>>> printnumber([53.2, 0.1])
"53.2(1)"
pymule.scaleplot(a, sx, sy=None)

rescales a plot such that the integrated plot remains unchanged, i.e. rescale \(x\to x/s\) and \(y\to y\cdot s\). This is useful to, for example, change units.

Parameters:

  • a – Nx3 numpy matrix; the plot

  • sx – float; the inverse scale factor for the x direction

  • sy – float, optional; if present, sy will be used for the y direction instead of sx

Returns:

a Nx3 numpy matrix

Example:

rescaling units from rad to mrad

>>> scaleplot(data, 1e-3)
pymule.timesnumbers(a, b)

multiplies numbers

Parameters:

  • a – list of floats; the first factor with error [a, da]

  • b – list of floats; the second factor with error [b, db]

Returns:

the result of the multiplication a*b [y, dy]

Example:

Divide \((0.5\pm0.02) * (45\pm0.01)\)

>>> timesnumbers([0.5,0.02], [45, 0.1])
array([22.5       ,  0.90138782])

Plotting

pymule.errorband(p, ax=None, col='default', underflow=False, overflow=False, linestyle='solid')

plots an errorband of a compatible histogram

Parameters:

  • p – Nx3 numpy matrix; the histogram to plot as np.array([[x1, y1, e1], [x2, y2, e2], ...])

  • ax – axes, optional: the axes object to use, defaults to gca() which may create a new axes.

  • col – the colour to be used for the plot. Per default matplotlib decides using the order specified in colours

  • underflow – bool, optional; whether to plot the underflow bin. Either logical or number indicating the how much bigger it shall be

  • overflow – bool, optional; whether to plot the overflow bin. Either logical or number indicating the how much bigger it shall be

  • linestyle – str, optional; which line style to use

Returns:

the artis of the main line but not the one of the errorbars

Example:

Make a simple plot

>>> errorband(dat)

Make a plot in red with dashed lines

>>> errorband(dat, 'red', 'dashed')
pymule.kplot(sigma, labelx='$x_e$', labelsigma=None, labelknlo='$\\delta K^{(1)}$', labelknnlo='$\\delta K^{(2)}$', legend={'lo': '$\\rm LO$', 'nlo': '$\\rm NLO$', 'nnlo': '$\\rm NNLO$'}, legendopts={'loc': 'upper right', 'what': 'l'}, linestyle2=':', show=[0, -1], showk=[1, 2], nomule=False)

produces a K factor plot in line with McMule’s design, i.e. a two-panel plot showing in the upper panel the cross sections and in the lower panel the K factor defined as

\[K^{(i)} = d\sigma^{(i)} / d\sigma^{(i-1)}\]
Parameters:

  • sigma – dict; the data to plot, given as a dict with keys lo, nlo, and possibly nnlo. Only pass the corrections, not the full distribution

  • labelx – str, optional; label for the x axis (supports LaTeX maths)

  • labelsigma – str, optional; label for the upper y axis (supports LaTeX maths)

  • labelknlo – str, optional; the labels for the NLO K factor

  • labelknnlo – str, optional; the labels for the NNLO K factor

  • show – list, optional; a list which cross sections to show, 0 indicates the LO cross section, 1 the NLO etc. -1 indicates the last given cross section

  • showk – list, optional; a list which K factors to show, 0 indicates the LO cross section, 1 the NLO etc. -1 indicates the last given cross section

  • legend – dict, optional; a dict with the legend for lo, nlo, nnlo. The keys nlo2 and nnlo2 are optional and will be drawn dashed in the lower panel.

  • legendopts –

    dict, optional; a kwargs dict of options to be passed to legend(..) as well as the what key indicating whether the legend such be placed in the lower panel (l, default), upper panel (u), or as a figlegend (fig). Notable is the loc-key that places the legend inside the object specified by what. Possible values are (cf. legend)

    • upper right

    • upper left

    • lower left

    • lower right

    • right

    • center left

    • center right

    • lower center

    • upper center

    • center

  • nomule – bool, optional; if set to True, no mule will be printed

Returns:

the figure as well as all axis created

Example:

An NNLO K factor plot

>>> fig, (ax1, ax2, ax3) = kplot(
...   {
...     'lo':  lodata['thetae'],
...     'nlo': nlodata['thetae'],
...     'nnlo':nnlodata['thetae'],
...   },
...   labelx="$\theta_e\,/\,{\rm mrad}$",
...   labelsigma="$\D\sigma/\D\theta_e\ /\ {\rm\upmu b}$",
...   legend={
...     'lo': '$\sigma^{(0)}$',
...     'nlo': '$\sigma^{(1)}$',
...     'nnlo': '$\sigma^{(2)}$'
...   },
...   legendopts={'what': 'u', 'loc': 'lower right'}
... )
pymule.mergefkswithplot(sets, scale=1.0, showfit=[True, True], xlim=[-7, 0])

performs and FKS merge like mergefks() but it also produces a \(\xi_c\) independence plot

Note

In contrast mergefks(), here phase-space partioned results need to be merged first. This is done by grouping those into an array first, sorted by number of particles in the final state, i.e. we start with the n-particle corrections.

Parameters:

  • sets – list of list of random-seed-merged results (usually from sigma()), starting with the lowest particle number and going up

  • scale – float, optional; rescale factor for the plot and result

  • showfit – [bool, bool], optional; whether to show the fit lines in the overview plot (first element) and the zoomed in plot (second element)

  • xlim – tuple of floats, optional; upper and lower bounds for \(\log\xi_c\)

Returns:

a figure and the FKS-merged final set containing cross sections, distributions, and run-time information. The chi2a return is a list of the following

  • the \(\chi^2\) of the FKS merge

  • a list of \(\chi^2\) from previous operations, such as random seed merging or the integration.

Example:

In the partioned muon-electron scattering case

>>> fig, res = mergefkswithplot([
...     [
...         sigma('em2emFEE'), sigma('em2emFMM'), sigma('em2emFEM')
...     ], [
...         sigma('em2emREE15'), sigma('em2emREE35'),
...         sigma('em2emRMM'),
...         sigma('em2emREM')
...     ]
... ])
pymule.mulify(fig, delx=0, dely=0, col='lightgray', realpha=True)

adds the McMule logo to a figure

Parameters:

  • fig – figure to add the logo

  • delx – float, optional; shift the logo in x direction

  • dely – float, optional; shift the logo in x direction

  • col – colour specifier, optional; colour to use for the logo

  • realpha – bool, optional; whether to re-run the alpha channel

pymule.watermark(fig, txt='PRELIMINARY', fontsize=60, rotation=20)

watermarks a figure

Parameters:

  • fig – the figure to watermark

  • txt – str, optional; the watermark text to use

  • fontsize – int, optional; the fontsize of the watermark

  • rotation – int, optional; the angle of the watermark in deg

Example:

Watermark a figure as preliminary

>>> fig = figure()
>>> ...
>>> watermark(fig)

Watermark a figure as incomplete

>>> fig = figure()
>>> ...
>>> watermark(fig, "INCOMPLETE")
pymule.xiresidue(sets, n, xlim=[-7, 0], scale=1)

creates a residue plot for a \(\xi_c\) fit

Parameters:

  • sets – dict or list; random-seed-merged results (usually from sigma()) or list thereof

  • n – int; order of the fit, 1 at NLO, 2 at NNLO

  • xlim – tuple of floats, optional; upper and lower bounds for \(\log\xi_c\)

  • scale – float, optional; rescale factor for the plots

Returns:

a figure and the fit coefficients as a matrix