In [1]:

from pymule import *
from pymule.plot import twopanel
import tarfile
import pymule.mpl_axes_aligner

setup(cachefolder='/tmp/mcmule')

Populating the interactive namespace from numpy and matplotlib

We take $E_\text{beam} = 150\,{\rm GeV}$ negative muons ($\mu^-$) and scatter on atomic electrons. We always require that $E_e > 1\,{\rm GeV}$.

Load Setup 2¶

In [2]:

setup(folder='em2em/out.tar.bz2')

In [3]:

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

Load NLO¶

In [4]:

nloEE = scaleset(mergefks(
    sigma('em2emREE15'), sigma('em2emREE35'),
    sigma('em2emFEE')
), conv*alpha**3)
nloMM = scaleset(mergefks(
    sigma('em2emRMM'), sigma('em2emFMM')
), conv*alpha**3)
nloEM = scaleset(mergefks(
    sigma('em2emREM'), sigma('em2emFEM')
), conv*alpha**3)

In [5]:

nlo = {
    k: combineNplots(addplots, [
        mergebins(nloEE[k], 2), nloMM[k], nloEM[k]
    ])
    for k in ['Ee', 'Em', 'tee', 'tmm', 'thetae', 'thetam']
}
nlo['value'] = plusnumbers(
    nloEE['value'], nloMM['value'], nloEM['value']
)
nlo['time'] = nloEE['time'] + nloMM['time'] + nloEM['time']

Load NNLO¶

In [6]:

nnlo = scaleset(mergefks(
    sigma('em2emRREE1'), sigma('em2emRREE3'),
    sigma('em2emRFEE15', obs='8'), sigma('em2emRFEE35', obs='8'),
    sigma('em2emFFEE', folder='ff2.tar.gz')
), conv*alpha**4)

Load Setup 4¶

Here, we further require that $\Big|\pi-|\phi_e-\phi_\mu|\Big| < 3.5\,{\rm mrad}$

In [7]:

setup(folder='em2em_aco/out.tar.bz2')

Load NLO¶

In [8]:

aconloEE = scaleset(mergefks(
    sigma('em2emREE15'), sigma('em2emREE35'),
    sigma('em2emFEE')
), conv*alpha**3)
aconloMM = scaleset(mergefks(
    sigma('em2emRMM'), sigma('em2emFMM')
), conv*alpha**3)
aconloEM = scaleset(mergefks(
    sigma('em2emREM'), sigma('em2emFEM')
), conv*alpha**3)

Here, we create $\xi_c$ plot

In [9]:

fignlo, aconlo = mergefkswithplot([
    [
        sigma('em2emFEE'), sigma('em2emFMM'), sigma('em2emFEM')
    ],
    [
        sigma('em2emREE15'), sigma('em2emREE35'),
        sigma('em2emRMM'),
        sigma('em2emREM')
    ]
], scale=conv*alpha**3)

mulify(fignlo)
fignlo.axes[0].set_ylim(-55, 90)
fignlo.axes[0].set_xlim(1.5e-3, 1.4)
savefig("plots/xi-nlo.pdf")

Load NNLO¶

Make $\xi_c$ plot

In [10]:

fignnlo, aconnlo = mergefkswithplot([
    [sigma('em2emFFEE', folder='ff2.tar.gz')],
    [sigma('em2emRFEE15'), sigma('em2emRFEE35')],
    [sigma('em2emRREE1'), sigma('em2emRREE3')]
], conv*alpha**4)

mulify(fignnlo)
fignnlo.axes[0].set_xlim(1.5e-3, 1.4)
fignnlo.axes[0].set_ylim(-7, 14)
savefig("plots/xi-nnlo.pdf")

Results¶

cross section¶

In [11]:

def parsestr(s):
    n = s.index("(")
    p = s.index('.')
    y = float(s[:n])
    e = float(s[n+1:-1])
    return np.array([y, e * 10**(p-n+1)])

In [12]:

loPV = parsestr('245.038910(1)')

In [13]:

entries = [
    lo['value'], loPV,

    timesnumbers(loPV, parsestr('0.04289(1)')),
    timesnumbers(loPV, parsestr('-0.08817(1)')),
    nloEE['value'],
    aconloEE['value'],
    timesnumbers(loPV, parsestr('-0.00028(1)')),
    timesnumbers(loPV, parsestr('-0.00256(1)')),
    nloMM['value'],
    [1, 100] * aconloMM['value'],
    timesnumbers(loPV, parsestr('-0.00147(2)')),
    timesnumbers(loPV, parsestr('0.00017(2)')),
    nloEM['value'],
    [1, 10] * aconloEM['value'],
    timesnumbers(loPV, parsestr('0.04114(1)')),
    timesnumbers(loPV, parsestr('-0.09055(1)')),
    nlo['value'],
    aconlo['value'],

    nnlo['value'], aconnlo['value'],
    plusnumbers(lo['value'], nlo['value'], nnlo['value']),
    plusnumbers(lo['value'], aconlo['value'], aconnlo['value'])
]
entries = [printnumber(i) for i in entries]

kfac = [
    nloEE['value'][0] / lo['value'][0], aconloEE['value'][0] / lo['value'][0],
    nloMM['value'][0] / lo['value'][0], aconloMM['value'][0] / lo['value'][0],
    nloEM['value'][0] / lo['value'][0], aconloEM['value'][0] / lo['value'][0],

    nlo['value'][0] / lo['value'][0], aconlo['value'][0] / lo['value'][0],

    nnlo['value'][0] / nlo['value'][0],
    aconnlo['value'][0] / aconlo['value'][0]
]

entries.insert(6, kfac[0])
entries.insert(7, kfac[1])
entries.insert(12, kfac[2])
entries.insert(13, kfac[3])
entries.insert(18, kfac[4])
entries.insert(19, kfac[5])
entries.insert(24, kfac[6])
entries.insert(25, kfac[7])
entries.insert(28, kfac[8])
entries.insert(29, kfac[9])

In [14]:

tex = "\n".join([
    r"%%!TEX root=thesis",
    r"\begin{tabular}{l|r|r||r|r||r|r}",
    r" & \multicolumn{2}{c||}{\cite{Alacevich:2018vez}}  & "
     + r"\multicolumn{2}{c||}{\mcmule{}} & \multicolumn{2}{c}{$K$} "
     + r"\\\cline{2-7}",
    r" & \multicolumn{1}{c|}{Setup 2} & \multicolumn{1}{c||}{Setup 4}",
    r" & \multicolumn{1}{c|}{Setup 2} & \multicolumn{1}{c||}{Setup 4}",
    r" & \multicolumn{1}{c|}{Setup 2} & \multicolumn{1}{c  }{Setup 4}               \\\hline",  # nopep8
    r"$\sigma  ^{(0)}$ & \multicolumn{2}{c||}{\tt%13s}  & \multicolumn{2}{c||}{\tt%13s}   \\",  # nopep8
    "\\hline",
    r"$\sigma_e  ^{(1)}$ & \tt %13s &\tt  %13s  &\tt %13s &\tt %13s &\tt %7.4f &\tt %7.4f \\",  # nopep8
    r"$\sigma_\mu^{(1)}$ & \tt %13s &\tt  %13s  &\tt %13s &\tt %13s &\tt %7.4f &\tt %7.4f \\",  # nopep8
    r"$\sigma_m  ^{(1)}$ & \tt %13s &\tt  %13s  &\tt %13s &\tt %13s &\tt %7.4f &\tt %7.4f \\",  # nopep8
    r"$\sigma    ^{(1)}$ & \tt %13s &\tt  %13s  &\tt %13s &\tt %13s &\tt %7.4f &\tt %7.4f \\",  # nopep8
    "\\hline",
    r"$\sigma_e  ^{(2)}$ &\multicolumn{2}{c||}{}&\tt %13s &\tt %13s &\tt %7.4f &\tt %7.4f \\",  # nopep8
    "\\hline",
    "\\hline",
    r"$\sigma          $ &\multicolumn{2}{c||}{}&\tt %13s &\tt %13s &\multicolumn{2}{c}{} \\",  # nopep8
    r"\end{tabular}"
]) % tuple(entries)

In [15]:

with open("plots/xs.tex", "w") as fp:
    fp.write(tex)

$\theta_e$¶

NLO split into $\sigma_e$, $\sigma_\mu$, and $\sigma_m$¶

define function¶

In [16]:

def dogaugeplot(lo, nloEE, nloEM, nloMM, nlo):
    def prep(x):
        if len(x) == 100:
            return scaleplot(x[:37], 1e-3)
        elif len(x) == 200:
            return scaleplot(mergebins(x, 2)[:37], 1e-3)
        elif len(x) == 202:
            return scaleplot(mergebins(x[1:-1], 2)[:37], 1e-3)
    distLO = prep(lo['thetae'])
    distEE = prep(nloEE['thetae'])
    distEM = prep(nloEM['thetae'])
    distMM = prep(nloMM['thetae'])
    distFU = prep(nlo['thetae'])

    distEE[:, 0] = distLO[:, 0]
    distEM[:, 0] = distLO[:, 0]
    distMM[:, 0] = distLO[:, 0]
    distFU[:, 0] = distLO[:, 0]
    distLO[36, 1:] = 0
    distEE[36, 1:] = 0
    distEM[36, 1:] = 0
    distMM[36, 1:] = 0
    distFU[36, 1:] = 0

    fig, (ax1, ax2, ax3) = twopanel(
        labelx="$\\theta_e\,/\,{\\rm mrad}$",
        upleft=[distLO, addplots(distLO, distFU)],
        colupleft=['C2', 'C0'],
        labupleft="$\\D\\sigma/\\D\\theta_e\ /\ {\\rm\\upmu b}$",

        downleft=[divideplots(distEE, distLO, 1)[4:-1]],
        labdownleft={
            'ylabel': "$K^{(1)}_e$",
            'color': 'C0'
        },
        coldownleft=['C0'],

        downright=[
            divideplots(distMM, distLO, 1)[:-1],
            divideplots(distEM, distLO, 1)[1:-1]
        ],
        labdownright={
            'ylabel': "$K^{(1)}_{\\mu,m}$",
            'color': 'C1'
        },
        downalign=[1, 1],
        coldownright=['C1', 'C1']
    )
    # ax3.collections = []
    ax3.lines[0].set_linestyle(":")
    mulify(fig)
    fig.legend(
        ax1.lines[:1] + ax2.lines + ax3.lines,
        [
            "$\\sigma^{(0)}$",
            "$\\sigma^{(1)}_e$",
            "$\\sigma^{(1)}_\\mu$",
            "$\\sigma^{(1)}_m$"
        ], loc='upper center', ncol=4
    )
    return fig

Make plots¶

Setup 2, i.e. no acoplanarity cut

In [17]:

fig = dogaugeplot(lo, nloEE, nloEM, nloMM, nlo)
savefig("plots/gauge.pdf")

Setup 4, i.e. with acoplanarity cut

In [18]:

acofig = dogaugeplot(lo, aconloEE, aconloEM, aconloMM, aconlo)
savefig("plots/gauge-aco.pdf")

NNLO¶

In [19]:

fig, (ax1, ax2, ax3) = kplot(
    {
        'lo': mergebins(scaleplot(lo['thetae'], 1e-3), 2)[:37],
        'nlo': mergebins(scaleplot(nloEE['thetae'], 1e-3), 2)[:37],
        'nnlo': mergebins(scaleplot(nnlo['thetae'], 1e-3), 2)[:37]
    },
    labelx="$\\theta_e\,/\,{\\rm mrad}$",
    labelsigma="$\\D\\sigma/\\D\\theta_e\ /\ {\\rm\\upmu b}$",
    legend={
        'lo': '$\\sigma^{(0)}$',
        'nlo': '$\\sigma_e^{(1)}$',
        'nnlo': '$\\sigma_e^{(2)}$'
    },
    legendopts={'what': 'u', 'loc': 'lower right'}
)
ax2.set_ylim(0.8, 2.5)
pymule.mpl_axes_aligner.yaxes(ax2, ax3)
savefig("plots/nnlo.pdf")

In [20]:

fig, (ax1, ax2, ax3) = kplot(
    {
        'lo': mergebins(scaleplot(lo['thetae'], 1e-3), 2)[:37],
        'nlo': mergebins(scaleplot(aconloEE['thetae'], 1e-3), 2)[1:38],
        'nnlo': mergebins(scaleplot(aconnlo['thetae'], 1e-3), 2)[1:38]
    },
    labelx="$\\theta_e\,/\,{\\rm mrad}$",
    labelsigma="$\\D\\sigma/\\D\\theta_e\ /\ {\\rm\\upmu b}$",
    legend={
        'lo': '$\\sigma^{(0)}$',
        'nlo': '$\\sigma_e^{(1)}$',
        'nnlo': '$\\sigma_e^{(2)}$'
    },
    legendopts={'what': 'u', 'loc': 'lower right'}
)
savefig("plots/nnlo-aco.pdf")

Compare with Pavia¶

In [21]:

t = tarfile.open('PV.tar.bz2')
def getPV(name):
    for i in t:
        if name == i.name:
            dat = np.loadtxt(t.extractfile(i))
            dat[:,0] += np.diff(dat[:,0])[0]/2
            return dat

$\theta_e$¶

In [22]:

from pymule import mpl_axes_aligner

Pavia provides the total electronic NNLO, i.e. $\sigma_2=\sigma^{(0)} + \sigma^{(1)} + \sigma^{(2)}$ in mrad. Hence, we need to rescale by 1000 and subtract our full NLO

In [23]:

resPV = scaleplot(getPV('PV/nocut/ethlab_PV_200.txt'),1e3)
nnloPV = mergebins(
    addplots(resPV, addplots(lo['thetae'], nloEE['thetae']), sb=-1),
    2
)

fullLO = mergebins(lo['thetae'],2)[:37]
fullNLO = addplots(fullLO, nlo['thetae'])[:37]
fullNNLO = addplots(fullNLO, mergebins(nnlo['thetae'],2))[:37]

In [24]:

fig, axs = plt.subplots(
    3, sharex=True, gridspec_kw={'hspace': 0}
)
axs = list(axs)

sca(axs[0])
errorband(fullLO, col='C2')
errorband(fullNNLO, col='C3')
legend(
    [Line2D([0], [0], color=i, lw=1.5) for i in ['C2', 'C0', 'C3']],
    ['LO', 'NLO', 'NNLO'],
    loc=4
)
ylabel(r"$\D\sigma/\D\theta_e\,/\,{\rm nb}$")

sca(axs[1])
errorband(divideplots(fullNLO, fullLO, offset=-1), col='C0')
gca().tick_params(labelcolor='C0', axis='y')
ylabel(r'$\delta K^{(1)}$', color='C0')

axs.append(gca().twinx())
errorband(divideplots(fullNNLO, fullNLO, offset=-1), col='C3')
gca().tick_params(labelcolor='C3', axis='y')
axhline(0, color='black', linewidth=0.4, zorder=1)
ylabel(r'$\delta K^{(2)}$', color='C3')
mpl_axes_aligner.yaxes(axs[1], gca(), 0,0)

sca(axs[2])
errorband(addplots(nnloPV,mergebins(nnlo['thetae'],2),sb=-1), col='C3')
ylabel(r"$\sigma^{(2)}_{\rm PV}-\sigma^{(2)}_{\textsc{McMule}}\,/\,{\rm nb}$")
axhline(0, color='black', linewidth=0.4, zorder=1)

xlabel(r"$\theta_e\,/\,{\rm mrad}$")
mulify(gcf())
gcf().savefig("plots/comparison.pdf")