Initialisation¶

In [1]:

from pymule import *
from pymule.plot import mpl_axes_aligner
from pymule.plot import twopanel, threepanel

Populating the interactive namespace from numpy and matplotlib

In [2]:

setup(folder='ee2ee_prad_prod/out.tar.bz2',cache='/tmp/mcmule')

In [3]:

def delete_zeros(plot):
    index = []
    for i in range(0,len(plot)):
        if(plot[i,2]==0):
            index.append(i)
    return delete(plot,index,axis=0)

In [4]:

def mergebins_part(plot,n_start,n_end,nbins):
    plot_merged = plot[0:n_start]
    plot_merged = concatenate((plot_merged, mergebins(plot[n_start:n_end],nbins)))
    plot_merged = concatenate((plot_merged, plot[n_end:]))
    return plot_merged

settings¶

In [5]:

SMALL_SIZE = 10
MEDIUM_SIZE = 12
BIGGER_SIZE = 14

rc('font', size=SMALL_SIZE)          # controls default text sizes
rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
rc('xtick', labelsize=MEDIUM_SIZE)   # fontsize of the tick labels
rc('ytick', labelsize=MEDIUM_SIZE)   # fontsize of the tick labels
rc('legend', fontsize=MEDIUM_SIZE)   # legend fontsize
rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title

import complete dataset¶

observable '14' corresponds to the fiducial cuts and observable '24' to the simplified cuts. the sanity checks filter out random seeds that cannot be trusted.

In [6]:

setup(obs='14')
loF = scaleset(mergefks(sigma('ee2ee0')), alpha**2*conv)
nloF_g = scaleset(mergefks(
    sigma('ee2eeF'), sigma('ee2eeR125'), sigma('ee2eeR345')
), alpha**3*conv)
nloF_vp = scaleset(
    sigma('ee2eeA', sanitycheck=lambda d: d['chi2a'] < 10)[(1.0,1.0)]
 , alpha**3*conv)
nnloF_g = scaleset(mergefks(
    sigma('ee2eeFF'), sigma('ee2eeRF125'), sigma('ee2eeRF345'),
    sigma('ee2eeRR15162526'), sigma('ee2eeRR35364546'),
), alpha**4*conv)
nnloF_vp = scaleset(mergefks(
    sigma('ee2eeAF'), sigma('ee2eeAR125'), sigma('ee2eeAR345'),
    anyxiAA=sigma('ee2eeAA'), anyxiNF=sigma('ee2eeNF')
), alpha**4*conv)

In [7]:

setup(obs='24')
loS = scaleset(mergefks(sigma('ee2ee0')), alpha**2*conv)
nloS_g = scaleset(mergefks(
    sigma('ee2eeF'), sigma('ee2eeR125'), sigma('ee2eeR345')
), alpha**3*conv)
nloS_vp = scaleset(
    sigma('ee2eeA', sanitycheck=lambda d: d['chi2a'] < 10)[(1.0,1.0)]
 , alpha**3*conv)
nnloS_g = scaleset(mergefks(
    sigma('ee2eeFF'), sigma('ee2eeRF125'), sigma('ee2eeRF345'),
    sigma('ee2eeRR15162526'), sigma('ee2eeRR35364546'),
), alpha**4*conv)
nnloS_vp = scaleset(mergefks(
    sigma('ee2eeAF'), sigma('ee2eeAR125'), sigma('ee2eeAR345'),
    anyxiAA=sigma('ee2eeAA'), anyxiNF=sigma('ee2eeNF')
), alpha**4*conv)

total cross section & $K$-factors¶

corrections fiducial:

In [8]:

printnumber(loF['value'])

Out[8]:

'2291.022(2)'

In [9]:

printnumber(nloF_g['value'])

Out[9]:

'-148.361(1)'

In [10]:

printnumber(nloF_vp['value'])

Out[10]:

'19.32910(2)'

In [11]:

printnumber(nnloF_g['value'])

Out[11]:

'2.8243(6)'

In [12]:

printnumber(nnloF_vp['value'])

Out[12]:

'-1.31228(1)'

total cross section fiducial:

In [13]:

sigma2F = plusnumbers(loF['value'],nloF_g['value'],nloF_vp['value'],nnloF_g['value'],nnloF_vp['value'])

In [14]:

printnumber(sigma2F)

Out[14]:

'2163.502(2)'

corrections simple:

In [15]:

printnumber(loS['value'])

Out[15]:

'2291.022(2)'

In [16]:

printnumber(nloS_g['value'])

Out[16]:

'78.226(6)'

In [17]:

printnumber(nloS_vp['value'])

Out[17]:

'19.32910(2)'

In [18]:

printnumber(nnloS_g['value'])

Out[18]:

'-0.172(3)'

In [19]:

printnumber(nnloS_vp['value'])

Out[19]:

'-0.05701(2)'

total cross section simple:

In [20]:

sigma2S = plusnumbers(loS['value'],nloS_g['value'],nloS_vp['value'],nnloS_g['value'],nnloS_vp['value'])

In [21]:

printnumber(sigma2S)

Out[21]:

'2388.349(6)'

$K$-factors fiducial:

In [22]:

sigma0F = loF['value']
sigma1F_g = plusnumbers([sigma0F,nloF_g['value']])
sigma1F_vp = plusnumbers([sigma0F,nloF_vp['value']])
sigma1F = plusnumbers([sigma1F_vp,nloF_g['value']])
sigma2F_g = plusnumbers([sigma1F,nnloF_g['value']])
sigma2F_vp = plusnumbers([sigma1F,nnloF_vp['value']])

In [23]:

dk1F_g = 100*dividenumbers(nloF_g['value'], sigma0F)
dk1F_vp = 100*dividenumbers(nloF_vp['value'], sigma0F)
dk2F_g = 100*dividenumbers(nnloF_g['value'], sigma1F)
dk2F_vp = 100*dividenumbers(nnloF_vp['value'], sigma1F)

In [24]:

printnumber(dk1F_g)

Out[24]:

'-6.47575(5)'

In [25]:

printnumber(dk1F_vp)

Out[25]:

'0.8436891(9)'

In [26]:

printnumber(dk2F_g)

Out[26]:

'0.13063(3)'

In [27]:

printnumber(dk2F_vp)

Out[27]:

'-0.0606977(6)'

$K$-factors simple:

In [28]:

sigma0S = loS['value']
sigma1S_g = plusnumbers([sigma0S,nloS_g['value']])
sigma1S_vp = plusnumbers([sigma0S,nloS_vp['value']])
sigma1S = plusnumbers([sigma1S_vp,nloS_g['value']])
sigma2S_g = plusnumbers([sigma1S,nnloS_g['value']])
sigma2S_vp = plusnumbers([sigma1S,nnloS_vp['value']])

In [29]:

dk1S_g = 100*dividenumbers(nloS_g['value'], sigma0S)
dk1S_vp = 100*dividenumbers(nloS_vp['value'], sigma0S)
dk2S_g = 100*dividenumbers(nnloS_g['value'], sigma1S)
dk2S_vp = 100*dividenumbers(nnloS_vp['value'], sigma1S)

In [30]:

printnumber(dk1S_g)

Out[30]:

'3.4145(2)'

In [31]:

printnumber(dk1S_vp)

Out[31]:

'0.8436891(9)'

In [32]:

printnumber(dk2S_g)

Out[32]:

'-0.0072(1)'

In [33]:

printnumber(dk2S_vp)

Out[33]:

'-0.0023868(6)'

fiducial energy distribution¶

extract distribution:

In [34]:

loFden = loF['e_n']
nloFden = addplots(nloF_g['e_n'],nloF_vp['e_n'])
nnloFden = addplots(nnloF_g['e_n'],nnloF_vp['e_n'])

merge bins:

In [35]:

nbins = 3
n1 = 0
n2 = 70
loFdenM = mergebins_part(loFden, n1, n2, nbins)
nloFdenM = mergebins_part(nloFden, n1, n2, nbins)
nnloFdenM = mergebins_part(nnloFden, n1, n2, nbins)

calculate $K$-factors:

In [36]:

sigma0Fden = loFdenM
sigma1Fden = addplots(sigma0Fden,nloFdenM)
sigma2Fden = addplots(sigma1Fden,nnloFdenM)
dk1den = divideplots(nloFdenM, sigma0Fden)
dk2den = divideplots(nnloFdenM, sigma1Fden)

remove undefined $K$-factor regions:

In [37]:

dk1den = delete_zeros(dk1den)
dk2den = delete_zeros(dk2den)

plot energy distribution:

In [38]:

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

sca(axs[0])
ylabel('$\\D\\sigma/\\D E_n\ /\ {\\rm\\upmu b}$')
errorband(sigma0Fden,col='C2')
errorband(sigma2Fden,col='C3')
yscale('log')

lines = [Line2D([0], [0], color='C2', lw=1.5),
                Line2D([0], [0], color='C0', lw=1.5),
                Line2D([0], [0], color='C3', lw=1.5)]

sca(axs[1])
ylabel('$\\delta K$')
errorband(dk1den,col='C0')
errorband(dk2den,col='C3')
axhline(0, color='black', linewidth=0.4, zorder=1)

sca(axs[2])
ylabel('$\\delta K$')
errorband(dk2den,col='C3')
axhline(0, color='black', linewidth=0.4, zorder=1)
ylim(-0.003,0.007)
xlabel('$E_n\,/\,{\\rm MeV}$')
xlim(100,1300)
legend(lines,['$\\mathrm{LO}$','$\\mathrm{NLO}$','$\\mathrm{NNLO}$'],loc='lower left')

mulify(fig)
fig.savefig('energy.pdf')

fiducial angular distribution¶

extract distributions:

In [39]:

loFdthn = loF['th_n']
nloFdthn = addplots(nloF_g['th_n'],nloF_vp['th_n'])
nnloFdthn = addplots(nnloF_g['th_n'],nnloF_vp['th_n'])

In [40]:

loFdthw = loF['th_w']
nloFdthw = addplots(nloF_g['th_w'],nloF_vp['th_w'])
nnloFdthw = addplots(nnloF_g['th_w'],nnloF_vp['th_w'])

merge bins:

In [41]:

loFdthnM = mergebins(loFdthn,1)
nloFdthnM = mergebins(nloFdthn,1)
nnloFdthnM = mergebins(nnloFdthn,1)

In [42]:

loFdthwM = mergebins(loFdthw,2)
nloFdthwM = mergebins(nloFdthw,2)
nnloFdthwM = mergebins(nnloFdthw,2)

calculate $K$-factors:

In [43]:

sigma0Fdthn = loFdthnM
sigma1Fdthn = addplots(sigma0Fdthn,nloFdthnM)
sigma2Fdthn = addplots(sigma1Fdthn,nnloFdthnM)
dk1dthn = divideplots(nloFdthnM, sigma0Fdthn)
dk2dthn = divideplots(nnloFdthnM, sigma1Fdthn)

In [44]:

sigma0Fdthw = loFdthwM
sigma1Fdthw = addplots(sigma0Fdthw,nloFdthwM)
sigma2Fdthw = addplots(sigma1Fdthw,nnloFdthwM)
dk1dthw = divideplots(nloFdthwM, sigma0Fdthw)
dk2dthw = divideplots(nnloFdthwM, sigma1Fdthw)

remove undefined $K$-factor regions:

In [45]:

dk1dthn = delete_zeros(dk1dthn)
dk2dthn = delete_zeros(dk2dthn)
dk1dthw = delete_zeros(dk1dthw)
dk2dthw = delete_zeros(dk2dthw)

plots angular distributions:

In [46]:

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

sca(axs[0])
ylabel('$\\D\\sigma/\\D\\theta\ /\ {\\rm\\upmu b}$')
errorband(sigma2Fdthn,col='C3')
errorband(sigma2Fdthw,col='C3',linestyle=(0,(0.4,0.4)))
yscale('log')

lines = [Line2D([0], [0], color='C0', lw=1.5),
         Line2D([0], [0], color='C3', lw=1.5),
         Line2D([0], [0], color='C0', lw=1.5,linestyle=(0,(0.4,0.4))),
         Line2D([0], [0], color='C3', lw=1.5,linestyle=(0,(0.4,0.4)))]

sca(axs[1])
ylabel('$\\delta K$')
axhline(0, color='black', linewidth=0.4, zorder=1)
errorband(dk1dthn,col='C0')
errorband(dk2dthn,col='C3')
errorband(dk1dthw[1:],col='C0',linestyle=(0,(0.4,0.4)))
errorband(dk2dthw,col='C3',linestyle=(0,(0.4,0.4)))

sca(axs[2])
ylabel('$\\delta K$')
axhline(0, color='black', linewidth=0.4, zorder=1)
errorband(dk2dthn,col='C3')
errorband(dk2dthw,col='C3',linestyle=(0,(0.4,0.4)))
ylim(-0.003,0.009)
xlim(0.5,7.2)
legend(lines,['$\\mathrm{NLO \\ narrow}$','$\\mathrm{NNLO \\ narrow}$',
              '$\\mathrm{NLO \\ wide}$','$\\mathrm{NNLO \\ wide}$'])

xlabel('$\\theta\,/\,{\\rm deg}$')
mulify(fig)
fig.savefig('angle.pdf')