Added muone analysis

mu-e-scattering/muone/muone.py

+# vim: foldmethod=marker
+from pymule import *
+from pymule.plot import twopanel
+import pymule.mpl_axes_aligner
+
+## Load Setup 2{{{
+setup(folder='em2em/out.tar.bz2')
+
+lo = scaleset(mergefks(sigma('em2em0')), conv*alpha**2)
+### Load NLO{{{
+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)
+
+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{{{
+
+nnlo = scaleset(mergefks(
+    sigma('em2emRREE1'), sigma('em2emRREE3'),
+    sigma('em2emRFEE15', obs='8'), sigma('em2emRFEE35', obs='8'),
+    sigma('em2emFF', folder='ff.tar.bz2')
+), conv*alpha**4)
+
+###########################################################}}}
+#####################################################################}}}
+## Load Setup 4{{{
+setup(folder='em2em_aco/out.tar.bz2')
+
+### Load NLO{{{
+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)
+
+
+fignlo, aconlo = mergefkswithplot([
+    [
+        sigma('em2emREE15'), sigma('em2emREE35'),
+        sigma('em2emRMM'),
+        sigma('em2emREM')
+    ],
+    [
+        sigma('em2emFEE'), sigma('em2emFMM'), sigma('em2emFEM')
+    ]
+], scale=conv*alpha**3)
+
+mulify(fignlo)
+fignlo.axes[0].set_ylim(-55,90)
+fignlo.axes[0].set_xlim(1.5e-3,1.4)
+
+
+###########################################################}}}
+### Load NNLO{{{
+
+fignnlo, aconnlo = mergefkswithplot([
+    [sigma('em2emRREE1') , sigma('em2emRREE3' )],
+    [sigma('em2emRFEE15'), sigma('em2emRFEE35')],
+    [sigma('em2emFF', folder='ff.tar.bz2')]
+], conv*alpha**4)
+
+mulify(fignnlo)
+fignnlo.axes[0].set_xlim(1.5e-3,1.4)
+fignnlo.axes[0].set_ylim(-7,14)
+
+###########################################################}}}
+#####################################################################}}}
+## Results{{{
+### cross section{{{
+
+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)])
+
+loPV = parsestr('245.038910(1)')
+
+
+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])
+
+tex = "\n".join([
+    r"%%!TEX root=thesis",
+    r"\begin{tabular}{l|r|r||r|r||r|r}",
+    r" & \multicolumn{2}{c||}{\cite{Alacevich:2018vez}}  & \multicolumn{2}{c||}{\mcmule{}} & \multicolumn{2}{c}{$K$} \\\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",
+    r"$\sigma  ^{(0)}$ & \multicolumn{2}{c||}{\tt%13s}  & \multicolumn{2}{c||}{\tt%13s}\\",
+    "\\hline",
+    r"$\sigma_e  ^{(1)}$ & \tt %13s &\tt  %13s  &\tt  %13s &\tt  %13s &\tt  %7.4f  &\tt %7.4f \\",
+    r"$\sigma_\mu^{(1)}$ & \tt %13s &\tt  %13s  &\tt  %13s &\tt  %13s &\tt  %7.4f  &\tt %7.4f \\",
+    r"$\sigma_m  ^{(1)}$ & \tt %13s &\tt  %13s  &\tt  %13s &\tt  %13s &\tt  %7.4f  &\tt %7.4f \\",
+    r"$\sigma    ^{(1)}$ & \tt %13s &\tt  %13s  &\tt  %13s &\tt  %13s &\tt  %7.4f  &\tt %7.4f \\",
+    "\\hline",
+    r"$\sigma_e  ^{(2)}$ & \multicolumn{2}{c||}{}         &\tt %13s & \tt %13s &\tt %7.4f &\tt %7.4f \\",
+    "\\hline",
+    "\\hline",
+    r"$\sigma          $ & \multicolumn{2}{c||}{}         &\tt %13s & \tt %13s &    \multicolumn{2}{c}{}     \\",
+    r"\end{tabular}"
+]) % tuple(entries)
+
+###########################################################}}}
+### theta e{{{
+### NLO split{{{
+### define function{{{
+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{{{
+fig = dogaugeplot(lo, nloEE, nloEM, nloMM, nlo)
+acofig = dogaugeplot(lo, aconloEE, aconloEM, aconloMM, aconlo)
+#######################################}}}
+#################################################}}}
+### NNLO {{{
+
+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,1,ax3,1)
+
+
+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'}
+)
+#################################################}}}
+###########################################################}}}
+#####################################################################}}}
+