pymule user guide

This section describes all public functions and classes in pymule.

Working with files

class pymule.Chi2

An object to represent a \(\chi^2\) value.

This can be created using Chi2.from_single() or Chi2.from_merge(). This class supports addition and float.

Attribute type:

an enum to indicate the type of \(\chi^2\) this is

Attribute value:

the value of \(\chi^2\) for single and mergers

Attribute children:

the children Chi2 objects for sums and mergers

class Type(value, names=None, *, module=None, qualname=None, type=None, start=1, boundary=None)

The \(\chi^2\) type enum

Attribute NONE:

empty entry

Attribute SINGLE:

a single entry

Attribute SUM:

a sum of entries

Attribute MERGE:

a statistical merge of entries

classmethod from_merge(chi: float, *s: Chi2) → Chi2

combines multiple \(\chi^2\) into a single one

Parameters:

  • chi – the value of \(\chi^2\) of the merge, must be positive

  • s – a sequence of Chi2 objects

Returns:

a new Chi2 object

Return type:

Chi2

Raises:

  • ValueError – the parameter chi is not positive or no Chi2 is be provided in s

  • TypeError – the parameter chi is not a number or the objects provided are not Chi2

Example:

Provide a merged set

>>> c1 = Chi2.from_single(1.02)
>>> c2 = Chi2.from_single(1.41)
>>> c3 = Chi2.from_single(2.02)
>>> tot = Chi2.from_merge(1.3, c1, c2, c3)
classmethod from_single(v: float) → Chi2

creates a single Chi2 object

Parameters:

v – the value of \(\chi^2\), must be positive

Returns:

a new Chi2 object

Return type:

Chi2

Raises:

  • ValueError – the parameter v is not positive

  • TypeError – the parameter v is not a number

Example:

Create a \(\chi^2 = 1.02\)

>>> Chi2.from_single(1.02)
<Chi^2: single 1.02>
class pymule.FileRecord

A subclass of Record to track McMule output files.

Additionally to the attributes of Record this contains all relevant properties of a McMule run

Attribute version:

the file version

Attribute inttype:

the integer type used in McMule in the format for struct

Attribute sha:

the sha hash provided by the file

Attribute sampler:

the sampler used, currently either vegas or foam

Attribute iteration:

the current iteration of McMule

Attribute randy:

the current random seed

Attribute msg:

the integration message

Example:

Load a record stored at /path/to/record.mcmule and print its value

>>> rec = FileRecord.read("/path/to/record.mcmule")
>>> print(rec.value)
[1.00646123e-04 1.06589201e-08]
classmethod read(filename: str = '', fp: IO[bytes] | None = None) → FileRecord

opens a McMule file, infers the correct file version, and creates a subclass of FileRecord

Parameters:

  • filename – file name to open, optional

  • fp – file pointer to read, optional

Returns:

a recording corresponding to the file

Return type:

subclass of FileRecord

Note

Either filename xor fp need to be specified

Note

If less than two iterations have been completed, no histograms will be returned

Example:

Load a file for the muon decay

>>> FileRecord.read('m2ennRR_mu-e_S0000068031X0.50000D0.50000_ITMX080x150M_O12.vegas')
<Record m2ennRR_mu-e_S0000068031X0.50000D0.50000_ITMX080x150M_O12.vegas xsec = 0.00010065(1), 12 histograms>
pymule.mergefks(*sets: XiRecord, **kwargs: XiRecord) → Record

performs the FKS merge

Parameters:

sets – random-seed-merged results of type XiRecord (usually from sigma())

Returns:

the FKS-merged final set containing cross sections, distributions, and run-time information.

Return type:

Record

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).

Example:

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

>>> mergefks(sigma("m2enn0"))
>>> <Record m2enn0 xsec = 2211500(2), 1 histograms>

Load the NLO results

>>> mergefks(sigma("m2ennV"), sigma("m2ennR"))
<Record xsec = -4143205.8581450013, 1 histograms>

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

>>> mergefks(
...     sigma("m2ennFF"), sigma("m2ennRF"), sigma("m2ennRR"),
...     anyxi=sigma("m2ennNF")
... )
<Record xsec = -131.13(9), 1 histograms>
pymule.setup(folder: str | None = None, **kwargs: Any) → None

sets the default arguments 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. .*.

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

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

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: str, **kwargs: Any) → XiRecord

loads a which_piece and returns a XiRecord

Parameters:

  • piece – str; which_piece to load

  • 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. .*.

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

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

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")
>>> <XiRecord set m2enn0, xic = [1.0], 1 histograms>

Load only observable O3

>>> sigma("m2enn0", obs="3")
>>> <XiRecord set m2enn0, xic = [1.0], 1 histograms>

Working with errors

pymule.addplots(a: <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, b: <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, sa: float = 1.0, sb: float = 1.0) → <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>

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: ~typing.List[~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>) → float

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: ~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, b: ~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>) → __getitem__()' id='127546590996240'>

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: <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, b: <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, offset: float = 0.0) → <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>

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: <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>) → float

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: <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, n: int) → <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>

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: ~typing.List[~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>) → __getitem__()' id='127546590996240'>

statistically combines values with uncertainties

Parameters:

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

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: ~typing.List[<MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>) → <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>

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], ...]), ...]

Returns:

a Nx3 numpy matrix

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: ~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>) → __getitem__()' id='127546590996240'>

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: ~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, prec: int = 0) → str

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: <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, sx: float, sy: float | None = None) → <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>

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: ~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>, b: ~typing.List[float] | <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>) → __getitem__()' id='127546590996240'>

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(*psets: XiRecord, showfit: tuple[bool, bool] = (True, True), xlim: tuple[float, float] = (-7.0, 0.0)) → tuple[Any, Record]

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 summed first.

Parameters:

  • sets – a sequence of XiRecord (or sums thereof, usually from sigma()), starting with the lowest particle number and going up

  • 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 Record containing cross sections, distributions, and run-time information.

Example:

In the partioned muon-electron scattering case

>>> fig, res = mergefkswithplot(
...     alpha*(sigma('em2emFEE')+sigma('em2emFMM')+sigma('em2emFEM')),
...     alpha*(
...         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(pset: XiRecord, n: int, xlim: tuple[float, float] = (-7.0, 0.0)) → tuple[~typing.Any, <MagicMock name='mock.typing.NDArray.__getitem__()' id='127546590996240'>]

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

Parameters:

  • sets – a XiRecord

  • 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\)

Returns:

a figure and the fit coefficients as a matrix