In [1]:

from pymule import *

setup(folder='ep2ep_mesa_paper/out.tar.bz2', cachefolder='/tmp/mcmule/')

Populating the interactive namespace from numpy and matplotlib

Load¶

In [2]:

lo = scaleset(mergefks(sigma('mp2mp0')), alpha**2*conv)

In [3]:

fig, nloNoVP = mergefkswithplot([
    [sigma('mp2mpF')],
    [sigma('mp2mpR15'), sigma('mp2mpR35')]
], scale=alpha**3*conv, xlim=[-3, 0])
fig.savefig('plots/xi-nlo.pdf')

In [4]:

fig, nnloNoVP = mergefkswithplot([
    [sigma('mp2mpFF')],
    [sigma('mp2mpRF15'), sigma('mp2mpRF35')],
    [sigma('mp2mpRR1516'), sigma('mp2mpRR3536')],
], scale=alpha**4*conv, xlim=[-3, 0])
fig.savefig('plots/xi-nnlo.pdf')

In [5]:

nlo = scaleset(mergefks(
    sigma('mp2mpF'),
    sigma('mp2mpR15'), sigma('mp2mpR35'),
    anyxiVP=sigma('mp2mpA', obs='1')
), alpha**3*conv)

In [6]:

nnlo = scaleset(mergefks(
    sigma('mp2mpFF'), sigma('mp2mpAF', obs='1'),
    sigma('mp2mpRF15'), sigma('mp2mpRF35'),
    sigma('mp2mpAR15', obs='1'), sigma('mp2mpAR35', obs='1'),
    sigma('mp2mpRR1516'), sigma('mp2mpRR3536'),
    anyxiF=sigma('mp2mpAA', obs='1'),
    anyxiNF=sigma('mp2mpNF', obs='1')
), alpha**4*conv)

In [7]:

nnloNF = scaleset(mergefks(
    sigma('mp2mpFF'),
    sigma('mp2mpRF15'), sigma('mp2mpRF35'),
    sigma('mp2mpRR1516'), sigma('mp2mpRR3536'),
    anyxiNF=sigma('mp2mpNF', obs='2')
), alpha**4*conv)

Print numbers¶

In [8]:

print("run time %f days" % (
    (lo['time'] + nlo['time'] + nnlo['time']) / 3600. / 24.)
)
print("\\sigma^(0) = " + printnumber(lo['value']))
print("\\sigma^(1) = +" + printnumber(nlo['value']))
print("\\sigma^(2) = " + printnumber(nnlo['value']))

run time 306.743920 days
\sigma^(0) = 34.5392485(6)
\sigma^(1) = +1.776362(5)
\sigma^(2) = -0.023728(2)

Make pictures¶

$\theta_e$¶

In [9]:

fig, (ax1, ax2, ax3) = kplot(
    {
        'lo':    mergebins(lo['thetae'], 4),
        'nlo':   mergebins(nlo['thetae'], 4),
        'nnlo':  mergebins(nnlo['thetae'], 4),
        'nlo2':  mergebins(nloNoVP['thetae'], 4),
        'nnlo2': mergebins(nnloNoVP['thetae'], 4)
    },
    labelx="$\\theta_e\,/\,{\\rm deg}$",
    labelsigma="$\\D\\sigma/\\D\\theta_e\ /\ {\\rm\\upmu b}$",
    legend={
        'lo': '$\\sigma^{(0)}$',
        'nlo': '$\\sigma^{(1)}$',
        'nnlo': '$\\sigma^{(2)}$',
        'nlo2': '$\\sigma^{(1)}_\\text{no VP}$',
        'nnlo2': '$\\sigma^{(2)}_\\text{no VP}$'
    },
    legendopts={'what': 'u', 'loc': 'upper right'}
)

fig.savefig('plots/thetae.pdf')

$Q_e^2$ and $Q_p^2$¶

In [10]:

fig, (ax1, ax2, ax3) = kplot(
    {
        'lo':    mergebins(lo['qsqP'][10:189], 4),
        'nlo':   mergebins(nlo['qsqP'][10:189], 4),
        'nnlo':  mergebins(nnlo['qsqP'][10:189], 4),
        'nlo2':  mergebins(nlo['qsqE'][10:189], 4),
        'nnlo2': mergebins(nnlo['qsqE'][10:189], 4),

    },
    labelx="$|t|\,/\,{\\rm MeV}^2$",
    labelsigma="$\\D\\sigma/\\D|t|\ /\ {\\rm\\upmu b}$",
    legend={
        'lo': '$\\sigma^{(0)}$',
        'nlo': '$\\sigma^{(1)}$',  'nlo2': '$\\sigma^{(1)}(|t|)$',
        'nnlo': '$\\sigma^{(2)}$', 'nnlo2': '$\\sigma^{(2)}(|t_e|)$'
    },
    legendopts={'what': 'u', 'loc': 'upper right'}
)
ax2.set_ylim(-0.07, 0.07)
ax3.set_ylim(-0.0042, 0.0042)

fig.savefig('plots/qsq.pdf')