In [1]:

from pymule import *
from pymule.maths import Li2, Li3, zeta2, zeta3
import pymule.compress
import arbuzov
from pymule.plot import threepanel

gamma0 = Mmu**5/192/pi**3

Populating the interactive namespace from numpy and matplotlib

Validation¶

Load data¶

In [2]:

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

In [3]:

openLeptons = mergefks(sigma('m2ennee0', merge={'xh': 2, 'xs': 2}))
openLeptons = addplots(openLeptons['xh'], openLeptons['xs'], sa=0.5, sb=0.5)

In [4]:

nnlo = mergefks(
    sigma('m2ennRF', merge={'xe': 5}),
    sigma('m2ennRR', merge={'xe': 2}),
    sigma('m2ennFF', merge={'xe': 2}, folder='../ff.tar.bz2', flavour='mu-e')
)
xev = np.array([1, 1 / gamma0, 1 / gamma0])*addplots(nnlo['xe'], openLeptons)

Load logarithms from [hep-ph/0205172] (Arbuzov et al.)¶

In [5]:

def f2loop(xe):
    fLL = arbuzov.f2LLns(xe)/3 + arbuzov.f2LLs(xe)/2 + arbuzov.f2LL(xe)/2
    fNLL = (
        arbuzov.f2NLLns(xe) + arbuzov.f2NLLs(xe)
        + arbuzov.f2NLL(xe) + arbuzov.f2NLLint(xe)
    )
    z = Mel / Mmu
    f0 = -np.sqrt(xe**2 - 4*z**2)*(2*xe**2 + 4*z**2 - 3*xe*(1 + z**2))
    f00 = -xe**2*(-3 + 2*xe)
    L = np.log(Mmu**2/Mel**2)
    return f0 / f00 * (L**2 * fLL + L * fNLL) / 2 / pi / pi

Load data from [hep-ph/0505069] (Anastasiou et al.)¶

In [6]:

const = np.array(pymule.compress.uncompress(
    "1:eJx11WtQU0cUAOAQ8AFqLRGFgAqCxAfYagBFCnNT5SEWTGtIHB21QuFCGST0BsSqgHU"
    "oBWsrAuWRUARtaUu0gkAtFguCRYQEgwgBEiCIQsBoRBRx1LQn/ursTs6Pndn5ZnfP2d3Z"
    "XREh3ElyaTSayAKakM9FiST9/z3KoCrSmbO7CKrGpWwm22qUQLwjz578+zz4BKtiuV4wg"
    "rrTa3mc0vYuQWWwrkZrL2hQpwmZDvNzwX1D8sKbZ6lRP7I94lumfTdBzfdTpDn39qD+jp"
    "gT5VADbueyYg7vZifqiqvM5o6oHoJa/XjN0SXyW6i7p7Hv/zQIPndVv1ow1Ix6Tmmlo/d"
    "eJUFxqtZbSizrUc+Qt9DHNOCVL+4nhJdVYvXNTvquKq6XoBo/qpmqrcpHPVOSyd78xujP"
    "VNvN+bkNiJ8rMkYfQZXyAm3cbX5H/V27FP24Qz/4e7p5pw1/oJ6afyDApwTcybWzKJHZg"
    "Pr0F7/VBrmqYP+WNnzTKGpGPSb0qSv/gtFrvty4x/Mm6rGeXX6pPmqC+pjHOnttuA31Q4"
    "eZCZFycM5CiT5oSwfq4dvavAWhAwQVR1fMyU5XoN7ecvH78tvgd463Pzsw2In61nPThmr"
    "BIEHFzK6TvFrZhbqvO/XaUQOuyD/Td2PHXdT9FvFsguOHCGrj5qCYjG3dqN/4laxTGMBL"
    "H/UFcRg9qC9vUtrot2gIil4TfPipFPPe91/LkjLAS76SRsSzlaj3PKtrd2sHv7Nz9lruJ"
    "czv1CW3DS8cJihvN570pV0v6h2pXrcKwsBlVjmjB4WYy7ZO3QwtAhdPNqmklzFvnVvVYj"
    "4ITlPrLu7QYZ5+qDv15Ip7BEW2LDk+6dqH+vEPclM4JLiskrMz71PMj9E+OTZdAa4o+cy"
    "s+DTm3vFJnimPjX4qo/J8A+aeHsYYgfxOSJhJasw3zDzZMJ4MvjG58kTgS8x1UYsnhH+B"
    "k3H/PLRd0o/6+NourSvtPkEVh/eHjbMxH9Nnj6n8wWP3aMWOgZiX7as9G5EJLhPor13lY"
    "l7iLCqxk4OLuW/MJyMwLx5b/2OH9QOTLuAbA5zNm9bwd2MeZr8obJ74gcn1F2Rump8zAG"
    "6xfzg53Rdzq9DnVtucR036HEa1pYEcNVl/U5ryenrFqMn8G/1/aPTRgyv8GQeb0jCvt+Q"
    "1PPEYI6hCP6fqmRzMk98GuILttbTtEuaJvl6H1tWDG1Zt/1oiwzzBbCpxhKaF82WS1nOn"
    "MHcT2q6L9QeXvXP0SN9CFeqrPbvdnLLA6RbZIxXumI+UvRwpl4EXPpcKbXmYa8gr9/Zaj"
    "0P+E429WiHmA25JwwwBOKlWR/LzMS94+/6DF3dOylnVmOft31XIHQSPbPt30b4WzHNcFh"
    "fMcpkw6SFnNnNPkeAG+cyT6/WYBwtehG6RTphcP9ChNmRGb3p+s6w+izSPh7B/9Q6lt3S"
    "YG0ILzL0Og8eUB1BxLDXqrxh8+sN68MKslJMBIsxfKBlmZWY6uH9RzdqwLswvB2yqWxUA"
    "XvihZVCkxwDql6yeXxnIMo5fxv1ZlIt57Nv/C5w2LWEtf4N5tJ8y3p7xCOrv0P7SSg6i7"
    "ix9xSrjg7eWe60RKTB3TKhz3SUG90w7UbHMZwj1pV7JKxcMgYt33XZrPY95zzmJOtf5MU"
    "F55FtrM+ZpUO+K3q0KjgaPpP+5ei9H0/Affiqv0w=="
))

const[:, 1] = (
    const[:, 1] * (
        -np.sqrt(const[:, 0]**2 - 4*(Mel/Mmu)**2)*(
            2*const[:, 0]**2 + 4*(Mel/Mmu)**2
            - 3*const[:, 0]*(1 + (Mel/Mmu)**2)
        )
    ) / alpha**2 * 2e-4
)

Make plot¶

In [7]:

xe = np.linspace(0, 1, 200)[2:-1]

In [8]:

ax = gca()
errorband(xev, ax=ax, col='C3')
ax.plot(xe, f2loop(xe), 'C2', linestyle='dashed')
ax.plot(const[:, 0], const[:, 1] + f2loop(const[:, 0]), 'C1')
ax.set_ylim(-15, 25)
ax.set_xlabel('$x_e$')
ax.set_ylabel("$\\D\\sigma^{(2)} / \\D x_e / \\sigma_0$")
mulify(ax.figure)

ax.legend([
    '$\\textsc{McMule}$',
    '\\cite{Arbuzov:2002pp,Arbuzov:2002cn}',
    '\\cite{Anastasiou2005The-electron}'
], loc='upper center')
savefig("comparison.pdf")

kplot¶

In [9]:

def f0loop(xe):
    beta = sqrt(1-Mel**2/(xe*Mmu/2)**2)
    return 2* xe**2*beta * (3-2*xe+xe/4*(3*xe-4)*(1-beta**2))
def f1loop(xe):
    beta = sqrt(1-Mel**2/(xe*Mmu/2)**2)
    return f0loop(xe) + alpha/(2*pi) * 2 * xe**2*beta * arbuzov.fnlo(xe, Mel/Mmu)

setup(folder='nlo.tar.bz2')
lo = mergefks(sigma('m2enn0', merge={'xe': 2}))
nlo = mergefks(
    sigma('m2ennF', merge={'xe': 2}),
    sigma('m2ennR', merge={'xe': 2})
)
fulllo = [1,1/gamma0,1/gamma0]*lo['xe']
fullnlo = addplots(fulllo, nlo['xe'], sb=alpha/gamma0)
fullnnlo = addplots(fullnlo, xev, sb=alpha**2)

fig, (ax1,ax2,ax3) = threepanel(r'$x = 2E_e/m_\mu$',
    upleft=[fulllo, fullnnlo],
    labupleft=r'$\D\mathcal{B}/\D x$',
    colupleft=[colours.orderscheme['lo'], colours.orderscheme['nnlo']],
    middleleft=divideplots(fullnlo, fulllo,-1),
    labmiddleleft=r'$\delta K_1 = \delta\mathcal{B}_1/\mathcal{B}_0$',
    colmiddleleft=[colours.orderscheme['nlo']],
    downleft=divideplots([1,1e4*alpha**2,1e4*alpha**2]*xev, fullnlo),
    labdownleft=r'$10^4\times(\delta K_2 = \delta\mathcal{B}_1/\mathcal{B}_1)$',
    coldownleft=[colours.orderscheme['nnlo']]
)

sca(ax2)
plot(xe, f1loop(xe) / f0loop(xe) - 1, 'C0', linestyle='dashed')
ylim(-0.12,1.22)

sca(ax3)
plot(xe, 1e4*alpha**2*f2loop(xe) / f1loop(xe), 'C2', linestyle='dashed')
plot(const[:, 0], 1e4*alpha**2*(const[:, 1] + f2loop(const[:, 0]))/f1loop(const[:,0]), 'C1')
ylim(-0.0005e4,0.0008e4)

ax1.legend(['LO','NNLO'],loc=4)
ax2.legend([r'\textsc{McMule}', 'Arbuzov 01'])
ax3.legend([r'\textsc{McMule}', 'Arbuzov 02', 'Anastasiou et al 05'],loc=2)
mulify(fig)
fig.savefig("full-plot.pdf")