$\mu-e$ scattering at NNLO

This work was done in collaboration with Alessandro Broggio, William J. Torres Bobadilla, Andrea Ferroglia, Manoj Kumar Mandal, Pierpaolo Mastrolia, Jonathan Ronca, Max Zoller

The results for muon-electron scattering are tailored to the characteristics of the MUonE experiment. The kinematics of the process is defined by

\begin{equation} e\,(p_1, m) + \mu\,(p_2, M, E) \to e\,(p_3, E_e, \theta_e) + \mu\,(p_4, E_\mu, \theta_\mu) \end{equation}

The muon beam energy is set to $E=160$ GeV, consistent with the M2 beam line at CERN North Area.

Cuts are imposed on the outgoing electron energy, $E_e > 1$ GeV, and on the muon scattering angle, $\theta_\mu > 0.3$ mrad.

An elasticity cut can be also enforced, \begin{equation} 0.9 < \frac{\theta_\mu}{\theta^{\rm el}_\mu} < 1.1 \,, \end{equation} using the relation between $\theta_e$ and $\theta_\mu$ in the absence of radiation \begin{equation} \tan\theta^{\rm el}_\mu = \frac{2\tan\theta_e}{(1+\gamma^2 \tan^2\theta_e) (1+g^*_\mu)-2} \,, \end{equation} where \begin{equation} \gamma = \frac{E+m}{\sqrt{s}}\,, \qquad g^*_\mu = \frac{E m+M^2}{E m+m^2}\,, \end{equation} and $s$ is the centre-of-mass energy.

Based on the constraints above, we can define four scenarios:

• S1 : $E_e > 1$ GeV, $\theta_\mu > 0.3$ mrad
• S2 : $E_e > 1$ GeV, $\theta_\mu > 0.3$ mrad, $0.9 < \theta_\mu/\theta^{\rm el}_\mu < 1.1$
• S1' : $E_e > 1$ GeV
• S2' : $E_e > 1$ GeV, $0.9 < \theta_\mu/\theta^{\rm el}_\mu < 1.1$

NLO and NNLO corrections can be split into gauge-invariant subsets according to their nominal leptonic charge. If $q$ resp. $Q$ is assigned to each photon emission from an electron resp. muon line, we will have, on top of the LO charge $q^2\,Q^2$, three subsets at NLO, i.e. $\{q^2,\, qQ,\, Q^2\}$, and five subsets at NNLO, i.e. $\{q^4,\,q^3 Q,\, q^2 Q^2,\,q Q^3,\,Q^4\}$.

The corrections can be generically split into photonic, hadronic and leptonic. The last two, when summed together, will be referred to as fermionic and include all the contributions from diagrams with VP insertions. Only photonic corrections are presented split into the subsets described above. Thus we have:

• NLO photonic (@PH) contributions, where @ can be (E, X, M)
• NLO hadronic resp. leptonic (HVP resp. LVP) contribution
• NNLO photonic (@PH) contributions, where @ can be (E, X31, X22, X13, M)
• NNLO hadronic resp. leptonic (HVP resp LVP) contribution

In addition, the last two characters in the names of the different contributions, i.e. m or p, and 0(2) or 1(3), refer to the polarisation of the muons (minus or plus) and to the scenarios where the elasticity cut is enforced, 1(3), or not, 0(2).

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

Populating the interactive namespace from numpy and matplotlib

Load McMule data

setup(cachefolder='/tmp/mcmule/')

S1: $\,\mu^- e\to \mu^- e$

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic-.tar.bz2')
lo          = scaleset(mergefks(sigma('em2em0')), alpha**2*conv)
nloEPHm0    = scaleset(mergefks(sigma('em2emFEE'),sigma('em2emREE15'),sigma('em2emREE35')), alpha**3*conv)
nloXPHm0    = scaleset(mergefks(sigma('em2emFEM'),sigma('em2emREM')), alpha**3*conv)
nloMPHm0    = scaleset(mergefks(sigma('em2emFMM'),sigma('em2emRMM')), alpha**3*conv)
nnloEPHm0   = scaleset(mergefks(sigma('em2emFFEEEE'),sigma('em2emRFEEEE15'),sigma('em2emRFEEEE35'),
                                sigma('em2emRREEEE1516'),sigma('em2emRREEEE3536')), alpha**4*conv)
nnloX31PHm0 = scaleset(mergefks(sigma('em2emFF31z'),sigma('em2emRF3115'),sigma('em2emRF3135'),
                                sigma('em2emRR311516'),sigma('em2emRR313536')), alpha**4*conv)
nnloX22PHm0 = scaleset(mergefks(sigma('em2emFF22z'),sigma('em2emRF2215'),sigma('em2emRF2235'),
                                sigma('em2emRR221516'),sigma('em2emRR223536')), alpha**4*conv)
nnloX13PHm0 = scaleset(mergefks(sigma('em2emFF13z'),sigma('em2emRF1315'),sigma('em2emRF1335'),
                                sigma('em2emRR131516'),sigma('em2emRR133536')), alpha**4*conv)
nnloMPHm0   = scaleset(mergefks(sigma('em2emFFMMMM'),sigma('em2emRFMMMM'),sigma('em2emRRMMMM')), alpha**4*conv)
setup(obs='2')
setup(folder='em2em_muone_paper/out-hadronic-.tar.bz2')
nloHVPm0    = scaleset(mergefks(sigma('em2emA')), alpha**3*conv)
nnloHVPm0   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic-.tar.bz2')
nloLVPm0    = scaleset(mergefks(sigma('em2emA')), alpha**3*conv)
nnloLVPm0   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)

S2: $\,\mu^- e\to \mu^- e$

setup(obs='1')
setup(folder='em2em_muone_paper/out-photonic-.tar.bz2')
nloEPHm1    = scaleset(mergefks(sigma('em2emFEE',obs='0'),sigma('em2emREE15'),sigma('em2emREE35')), alpha**3*conv)
nloXPHm1    = scaleset(mergefks(sigma('em2emFEM',obs='0'),sigma('em2emREM')), alpha**3*conv)
nloMPHm1    = scaleset(mergefks(sigma('em2emFMM',obs='0'),sigma('em2emRMM')), alpha**3*conv)
nnloEPHm1   = scaleset(mergefks(sigma('em2emFFEEEE',obs='0'),sigma('em2emRFEEEE15'),sigma('em2emRFEEEE35'),
                                sigma('em2emRREEEE1516'),sigma('em2emRREEEE3536')), alpha**4*conv)
nnloX31PHm1 = scaleset(mergefks(sigma('em2emFF31z',obs='0'),sigma('em2emRF3115'),sigma('em2emRF3135'),
                                sigma('em2emRR311516'),sigma('em2emRR313536')),alpha**4*conv)
nnloX22PHm1 = scaleset(mergefks(sigma('em2emFF22z',obs='0'),sigma('em2emRF2215'),sigma('em2emRF2235'),
                                sigma('em2emRR221516'),sigma('em2emRR223536')),alpha**4*conv)
nnloX13PHm1 = scaleset(mergefks(sigma('em2emFF13z',obs='0'),sigma('em2emRF1315'),sigma('em2emRF1335'),
                                sigma('em2emRR131516'),sigma('em2emRR133536')),alpha**4*conv)
nnloMPHm1   = scaleset(mergefks(sigma('em2emFFMMMM',obs='0'),sigma('em2emRFMMMM'),sigma('em2emRRMMMM')), alpha**4*conv)
setup(obs='3')
setup(folder='em2em_muone_paper/out-hadronic-.tar.bz2')
nnloHVPm1   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF',obs='2'),anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic-.tar.bz2')
nnloLVPm1   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF',obs='2'),anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)

S1: $\,\mu^+ e\to \mu^+ e$

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic+.tar.bz2')
nloXPHp0    = scaleset(mergefks(sigma('em2emFEM'),sigma('em2emREM')), alpha**3*conv)
setup(obs='2')
setup(folder='em2em_muone_paper/out-hadronic+.tar.bz2')
nnloHVPp0   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic+.tar.bz2')
nnloLVPp0   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)

S2: $\,\mu^+ e\to \mu^+ e$

setup(obs='1')
setup(folder='em2em_muone_paper/out-photonic+.tar.bz2')
nloXPHp1    = scaleset(mergefks(sigma('em2emFEM',obs='0'),sigma('em2emREM')), alpha**3*conv)
setup(obs='3')
setup(folder='em2em_muone_paper/out-hadronic+.tar.bz2')
nnloHVPp1   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF',obs='2'),anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic+.tar.bz2')
nnloLVPp1   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'),sigma('em2emAR35'),
                                anyxi1=sigma('em2emNF',obs='2'),anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)

S1': $\,\mu^- e\to \mu^- e$

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic-_thm.tar.bz2')
loh          = scaleset(mergefks(sigma('em2em0')), alpha**2*conv)
nloEPHm0h    = scaleset(mergefks(sigma('em2emFEE'),sigma('em2emREE15'), sigma('em2emREE35')), alpha**3*conv)
nloXPHm0h    = scaleset(mergefks(sigma('em2emFEM'),sigma('em2emREM')), alpha**3*conv)
nloMPHm0h    = scaleset(mergefks(sigma('em2emFMM'),sigma('em2emRMM')), alpha**3*conv)
nnloEPHm0h   = scaleset(mergefks(sigma('em2emFFEEEE'),sigma('em2emRFEEEE15'),sigma('em2emRFEEEE35'),
                                 sigma('em2emRREEEE1516'),sigma('em2emRREEEE3536')), alpha**4*conv)
nnloX31PHm0h = scaleset(mergefks(sigma('em2emFF31z'),sigma('em2emRF3115'),sigma('em2emRF3135'),
                                 sigma('em2emRR311516'),sigma('em2emRR313536')), alpha**4*conv)
nnloX22PHm0h = scaleset(mergefks(sigma('em2emFF22z'),sigma('em2emRF2215'),sigma('em2emRF2235'),
                                 sigma('em2emRR221516'),sigma('em2emRR223536')), alpha**4*conv)
nnloX13PHm0h = scaleset(mergefks(sigma('em2emFF13z'),sigma('em2emRF1315'),sigma('em2emRF1335'),
                                 sigma('em2emRR131516'),sigma('em2emRR133536')), alpha**4*conv)
nnloMPHm0h   = scaleset(mergefks(sigma('em2emFFMMMM'),sigma('em2emRFMMMM'),sigma('em2emRRMMMM')), alpha**4*conv)
setup(obs='2')
setup(folder='em2em_muone_paper/out-hadronic-_thm.tar.bz2')
nloHVPm0h    = scaleset(mergefks(sigma('em2emA')), alpha**3*conv)
nnloHVPm0h   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic-_thm.tar.bz2')
nloLVPm0h    = scaleset(mergefks(sigma('em2emA')), alpha**3*conv)
nnloLVPm0h   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)

S2': $\,\mu^- e\to \mu^- e$

setup(obs='1')
setup(folder='em2em_muone_paper/out-photonic-_thm.tar.bz2')
nloEPHm1h    = scaleset(mergefks(sigma('em2emFEE',obs='0'),sigma('em2emREE15'),sigma('em2emREE35')), alpha**3*conv)
nloXPHm1h    = scaleset(mergefks(sigma('em2emFEM',obs='0'),sigma('em2emREM')), alpha**3*conv)
nloMPHm1h    = scaleset(mergefks(sigma('em2emFMM',obs='0'),sigma('em2emRMM')), alpha**3*conv)
nnloEPHm1h   = scaleset(mergefks(sigma('em2emFFEEEE',obs='0'),sigma('em2emRFEEEE15'),sigma('em2emRFEEEE35'),
                                 sigma('em2emRREEEE1516'),sigma('em2emRREEEE3536')), alpha**4*conv)
nnloX31PHm1h = scaleset(mergefks(sigma('em2emFF31z',obs='0'),sigma('em2emRF3115'),sigma('em2emRF3135'),
                                 sigma('em2emRR311516'),sigma('em2emRR313536')),alpha**4*conv)
nnloX22PHm1h = scaleset(mergefks(sigma('em2emFF22z',obs='0'),sigma('em2emRF2215'),sigma('em2emRF2235'),
                                 sigma('em2emRR221516'),sigma('em2emRR223536')),alpha**4*conv)
nnloX13PHm1h = scaleset(mergefks(sigma('em2emFF13z',obs='0'),sigma('em2emRF1315'),sigma('em2emRF1335'),
                                 sigma('em2emRR131516'),sigma('em2emRR133536')),alpha**4*conv)
nnloMPHm1h   = scaleset(mergefks(sigma('em2emFFMMMM',obs='0'),sigma('em2emRFMMMM'),sigma('em2emRRMMMM')), alpha**4*conv)
setup(obs='3')
setup(folder='em2em_muone_paper/out-hadronic-_thm.tar.bz2')
nnloHVPm1h   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'),sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF',obs='2'), anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic-_thm.tar.bz2')
nnloLVPm1h   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'),sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF',obs='2'), anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)

S1': $\,\mu^+ e\to \mu^+ e$

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic+_thm.tar.bz2')
nloXPHp0h    = scaleset(mergefks(sigma('em2emFEM'),sigma('em2emREM')), alpha**3*conv)
setup(obs='2')
setup(folder='em2em_muone_paper/out-hadronic+_thm.tar.bz2')
nnloHVPp0h   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic+_thm.tar.bz2')
nnloLVPp0h   = scaleset(mergefks(sigma('em2emAF'),sigma('em2emAR15'),sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF'),anyxi2=sigma('em2emAA')), alpha**4*conv)

S2': $\,\mu^+ e\to \mu^+ e$

setup(obs='1')
setup(folder='em2em_muone_paper/out-photonic+_thm.tar.bz2')
nloXPHp1h    = scaleset(mergefks(sigma('em2emFEM',obs='0'),sigma('em2emREM')), alpha**3*conv)
setup(obs='3')
setup(folder='em2em_muone_paper/out-hadronic+_thm.tar.bz2')
nnloHVPp1h   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'),sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF',obs='2'),anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-leptonic+_thm.tar.bz2')
nnloLVPp1h   = scaleset(mergefks(sigma('em2emAF',obs='2'),sigma('em2emAR15'), sigma('em2emAR35'),
                                 anyxi1=sigma('em2emNF',obs='2'),anyxi2=sigma('em2emAA',obs='2')), alpha**4*conv)

Cross sections ($\mu$b) - S1 and S2

def pn(a):
    return printnumber(a)
def dn(a,b):
    return dividenumbers(a,b)

lo0       = lo['value']
nlo0m     = plusnumbers(nloEPHm0['value'],  nloXPHm0['value'],  nloMPHm0['value'],  nloHVPm0['value'],  nloLVPm0['value'],  lo0)
nlo1m     = plusnumbers(nloEPHm1['value'],  nloXPHm1['value'],  nloMPHm1['value'],  nloHVPm0['value'],  nloLVPm0['value'],  lo0)
nnloYPHm0 = plusnumbers( nnloX31PHm0['value'], nnloX22PHm0['value'],  nnloX13PHm0['value'])
nnloYPHm1 = plusnumbers( nnloX31PHm1['value'], nnloX22PHm1['value'],  nnloX13PHm1['value'])
nnlo0m    = plusnumbers(nnloEPHm0['value'], nnloYPHm0, nnloMPHm0['value'], nnloHVPm0['value'], nnloLVPm0['value'], nlo0m)
nnlo1m    = plusnumbers(nnloEPHm1['value'], nnloYPHm1, nnloMPHm1['value'], nnloHVPm1['value'], nnloLVPm1['value'], nlo1m)

print("MUONE- [LO]                                       :: ",
      pn(lo0))
print("MUONE- [NLOm0, NLOm1, dK[NLOm0]%, dK[NLOm1]%]     :: ",
      pn(nlo0m),  pn(nlo1m),  pn(dn(plusnumbers(nlo0m, -lo0),  lo0)  *100.), pn(dn(plusnumbers(nlo1m, -lo0),  lo0)  *100.))
print("MUONE- [NNLOm0, NNLOm1, dK[NNLOm0]%, dK[NNLOm1]%] :: ",
      pn(nnlo0m), pn(nnlo1m), pn(dn(plusnumbers(nnlo0m,-nlo0m),nlo0m)*100.), pn(dn(plusnumbers(nnlo1m,-nlo1m),nlo1m)*100.))

MUONE- [LO]                                       ::  106.443563(5)
MUONE- [NLOm0, NLOm1, dK[NLOm0]%, dK[NLOm1]%]     ::  106.99038(3) 102.86304(3) 0.51372(3) -3.36377(3)
MUONE- [NNLOm0, NNLOm1, dK[NNLOm0]%, dK[NNLOm1]%] ::  106.97977(3) 102.88154(3) -0.00992(4) 0.01799(4)

nlo0p     = plusnumbers(nloEPHm0['value'],  nloXPHp0['value'],  nloMPHm0['value'],  nloHVPm0['value'],  nloLVPm0['value'],  lo0)
nlo1p     = plusnumbers(nloEPHm1['value'],  nloXPHp1['value'],  nloMPHm1['value'],  nloHVPm0['value'],  nloLVPm0['value'],  lo0)
nnloYPHp0 = plusnumbers(-nnloX31PHm0['value'], nnloX22PHm0['value'], -nnloX13PHm0['value'])
nnloYPHp1 = plusnumbers(-nnloX31PHm1['value'], nnloX22PHm1['value'], -nnloX13PHm1['value'])
nnlo0p    = plusnumbers(nnloEPHm0['value'], nnloYPHp0, nnloMPHm0['value'], nnloHVPp0['value'], nnloLVPp0['value'], nlo0p)
nnlo1p    = plusnumbers(nnloEPHm1['value'], nnloYPHp1, nnloMPHm1['value'], nnloHVPp1['value'], nnloLVPp1['value'], nlo1p)

print("MUONE+ [LO]                                       :: ",
      pn(lo0))
print("MUONE+ [NLOp0, NLOp1, dK[NLOp0]%, dK[NLOp1]%]     :: ",
      pn(nlo0p),  pn(nlo1p),  pn(dn(plusnumbers(nlo0p, -lo0),  lo0)  *100.), pn(dn(plusnumbers(nlo1p, -lo0),  lo0)  *100.))
print("MUONE+ [NNLOp0, NNLOp1, dK[NNLOp0]%, dK[NNLOp1]%] :: ",
      pn(nnlo0p), pn(nnlo1p), pn(dn(plusnumbers(nnlo0p,-nlo0p),nlo0p)*100.), pn(dn(plusnumbers(nnlo1p,-nlo1p),nlo1p)*100.))

MUONE+ [LO]                                       ::  106.443563(5)
MUONE+ [NLOp0, NLOp1, dK[NLOp0]%, dK[NLOp1]%]     ::  107.41847(3) 103.18338(3) 0.91589(3) -3.06283(3)
MUONE+ [NNLOp0, NNLOp1, dK[NNLOp0]%, dK[NNLOp1]%] ::  107.41832(3) 103.19386(3) -0.00013(4) 0.01016(4)

print("                  |    S1    S2   |   dK1%    dK2%  ")
print("--------------------------------------------------------------------------------------------")
print("\sigma^{(0)}      | ", pn(lo['value']), " ", pn(lo['value']))
print("--------------------------------------------------------------------------------------------")
print("\sigma^{(1)}_e    | ", pn(nloEPHm0['value']), " ", pn(nloEPHm1['value'])," | ",   pn(dn(nloEPHm0['value'],lo0)*100.), " ", pn(dn(nloEPHm1['value'],lo0)*100.))
print("\sigma^{(1)}_emM  | ", pn(nloXPHm0['value']), " ", pn(nloXPHm1['value'])," | ",   pn(dn(nloXPHm0['value'],lo0)*100.), " ", pn(dn(nloXPHm1['value'],lo0)*100.))
print("\sigma^{(1)}_emP  | ", pn(nloXPHp0['value']), " ", pn(nloXPHp1['value'])," | ",   pn(dn(nloXPHp0['value'],lo0)*100.), " ", pn(dn(nloXPHp1['value'],lo0)*100.))
print("\sigma^{(1)}_m    | ", pn(nloMPHm0['value']), " ", pn(nloMPHm1['value'])," | ",   pn(dn(nloMPHm0['value'],lo0)*100.), " ", pn(dn(nloMPHm1['value'],lo0)*100.))
print("\sigma^{(1)}_lep  | ", pn(nloLVPm0['value']), " ", pn(nloLVPm0['value'])," | ",   pn(dn(nloLVPm0['value'],lo0)*100.), " ", pn(dn(nloLVPm0['value'],lo0)*100.))
print("\sigma^{(1)}_had  | ", pn(nloHVPm0['value']), " ", pn(nloHVPm0['value'])," | ",   pn(dn(nloHVPm0['value'],lo0)*100.), " ", pn(dn(nloHVPm0['value'],lo0)*100.))
print("--------------------------------------------------------------------------------------------")
print("\sigma^{(2)}_e    | ", pn(nnloEPHm0['value']), " ", pn(nnloEPHm1['value'])," | ", pn(dn(nnloEPHm0['value'],nlo0m)*100.), " ", pn(dn(nnloEPHm1['value'],nlo1m)*100.))
print("\sigma^{(2)}_emM  | ", pn(nnloYPHm0), " ", pn(nnloYPHm1)," | ", pn(dn(nnloYPHm0,nlo0m)*100.), " ", pn(dn(nnloYPHm1,nlo1m)*100.))
print("\sigma^{(2)}_emP  | ", pn(nnloYPHp0), " ", pn(nnloYPHp1)," | ", pn(dn(nnloYPHp0,nlo0p)*100.), " ", pn(dn(nnloYPHp1,nlo1p)*100.))
print("\sigma^{(2)}_m    | ", pn(nnloMPHm0['value']), " ", pn(nnloMPHm1['value'])," | ", pn(dn(nnloMPHm0['value'],nlo0m)*100.), " ", pn(dn(nnloMPHm1['value'],nlo1m)*100.))
print("\sigma^{(2)}_lepM | ", pn(nnloLVPm0['value']), " ", pn(nnloLVPm1['value'])," | ", pn(dn(nnloLVPm0['value'],nlo0m)*100.), " ", pn(dn(nnloLVPm1['value'],nlo1m)*100.))
print("\sigma^{(2)}_lepP | ", pn(nnloLVPp0['value']), " ", pn(nnloLVPp1['value'])," | ", pn(dn(nnloLVPp0['value'],nlo0p)*100.), " ", pn(dn(nnloLVPp1['value'],nlo1p)*100.))
print("\sigma^{(2)}_hadM | ", pn(nnloHVPm0['value']), " ", pn(nnloHVPm1['value'])," | ", pn(dn(nnloHVPm0['value'],nlo0m)*100.), " ", pn(dn(nnloHVPm1['value'],nlo1m)*100.))
print("\sigma^{(2)}_hadP | ", pn(nnloHVPp0['value']), " ", pn(nnloHVPp1['value'])," | ", pn(dn(nnloHVPp0['value'],nlo0p)*100.), " ", pn(dn(nnloHVPp1['value'],nlo1p)*100.))
print("--------------------------------------------------------------------------------------------")

                  |    S1    S2   |   dK1%    dK2%  
--------------------------------------------------------------------------------------------
\sigma^{(0)}      |  106.443563(5)   106.443563(5)
--------------------------------------------------------------------------------------------
\sigma^{(1)}_e    |  -0.61211(3)   -4.66042(3)  |  -0.57506(3)   -4.37830(3)
\sigma^{(1)}_emM  |  -0.2140417(10)   -0.160169(1)  |  -0.2010846(9)   -0.150473(1)
\sigma^{(1)}_emP  |  0.2140420(10)   0.160169(1)  |  0.2010850(9)   0.150474(1)
\sigma^{(1)}_m    |  -0.0284268(4)   -0.1613361(5)  |  -0.0267060(4)   -0.1515696(4)
\sigma^{(1)}_lep  |  1.3857484(1)   1.3857484(1)  |  1.3018621(1)   1.3018621(1)
\sigma^{(1)}_had  |  0.0156539754(9)   0.0156539754(9)  |  0.014706362(1)   0.014706362(1)
--------------------------------------------------------------------------------------------
\sigma^{(2)}_e    |  0.000896(8)   0.065945(5)  |  0.000838(7)   0.064109(5)
\sigma^{(2)}_emM  |  0.0009469(3)   0.0192568(3)  |  0.0008850(3)   0.0187208(3)
\sigma^{(2)}_emP  |  0.0032930(3)   0.0047869(3)  |  0.0030655(3)   0.0046392(3)
\sigma^{(2)}_m    |  -0.000050708(3)   0.000018293(1)  |  -0.000047395(2)   0.000017783(1)
\sigma^{(2)}_lepM |  -0.0119494(9)   -0.0656753(8)  |  -0.0111687(8)   -0.0638474(8)
\sigma^{(2)}_lepP |  -0.0042408(10)   -0.0595863(9)  |  -0.0039479(9)   -0.0577479(9)
\sigma^{(2)}_hadM |  -0.000453315(6)   -0.001042362(6)  |  -0.000423697(5)   -0.001013350(6)
\sigma^{(2)}_hadP |  -0.000042309(7)   -0.000683733(7)  |  -0.000039387(7)   -0.000662639(7)
--------------------------------------------------------------------------------------------

Cross sections ($\mu$b) - S1' and S2'

lo0h       = loh['value']
nlo0mh     = plusnumbers(nloEPHm0h['value'],  nloXPHm0h['value'],  nloMPHm0h['value'],  nloHVPm0h['value'],  nloLVPm0h['value'],  lo0h)
nlo1mh     = plusnumbers(nloEPHm1h['value'],  nloXPHm1h['value'],  nloMPHm1h['value'],  nloHVPm0h['value'],  nloLVPm0h['value'],  lo0h)
nnloYPHm0h = plusnumbers( nnloX31PHm0h['value'], nnloX22PHm0h['value'],  nnloX13PHm0h['value'])
nnloYPHm1h = plusnumbers( nnloX31PHm1h['value'], nnloX22PHm1h['value'],  nnloX13PHm1h['value'])
nnlo0mh    = plusnumbers(nnloEPHm0h['value'], nnloYPHm0h, nnloMPHm0h['value'], nnloHVPm0h['value'], nnloLVPm0h['value'], nlo0mh)
nnlo1mh    = plusnumbers(nnloEPHm1h['value'], nnloYPHm1h, nnloMPHm1h['value'], nnloHVPm1h['value'], nnloLVPm1h['value'], nlo1mh)

print("MUONEh- [LO]                                       :: ",
      pn(lo0h))
print("MUONEh- [NLOm0, NLOm1, dK[NLOm0]%, dK[NLOm1]%]     :: ",
      pn(nlo0mh),  pn(nlo1mh),  pn(dn(plusnumbers(nlo0mh, -lo0h),  lo0h)  *100.), pn(dn(plusnumbers(nlo1mh, -lo0h),  lo0h)  *100.))
print("MUONEh- [NNLOm0, NNLOm1, dK[NNLOm0]%, dK[NNLOm1]%] :: ",
      pn(nnlo0mh), pn(nnlo1mh), pn(dn(plusnumbers(nnlo0mh,-nlo0mh),nlo0mh)*100.), pn(dn(plusnumbers(nnlo1mh,-nlo1mh),nlo1mh)*100.))

MUONEh- [LO]                                       ::  245.59625(1)
MUONEh- [NLOm0, NLOm1, dK[NLOm0]%, dK[NLOm1]%]     ::  258.9373(2) 236.81343(7) 5.43211(7) -3.57612(3)
MUONEh- [NNLOm0, NNLOm1, dK[NNLOm0]%, dK[NNLOm1]%] ::  258.9977(2) 236.88109(7) 0.02332(9) 0.02857(4)

nlo0ph     = plusnumbers(nloEPHm0h['value'],  nloXPHp0h['value'],  nloMPHm0h['value'],  nloHVPm0h['value'],  nloLVPm0h['value'],  lo0h)
nlo1ph     = plusnumbers(nloEPHm1h['value'],  nloXPHp1h['value'],  nloMPHm1h['value'],  nloHVPm0h['value'],  nloLVPm0h['value'],  lo0h)
nnloYPHp0h = plusnumbers(-nnloX31PHm0h['value'], nnloX22PHm0h['value'], -nnloX13PHm0h['value'])
nnloYPHp1h = plusnumbers(-nnloX31PHm1h['value'], nnloX22PHm1h['value'], -nnloX13PHm1h['value'])
nnlo0ph    = plusnumbers(nnloEPHm0h['value'], nnloYPHp0h, nnloMPHm0h['value'], nnloHVPp0h['value'], nnloLVPp0h['value'], nlo0ph)
nnlo1ph    = plusnumbers(nnloEPHm1h['value'], nnloYPHp1h, nnloMPHm1h['value'], nnloHVPp1h['value'], nnloLVPp1h['value'], nlo1ph)

print("MUONEh+ [LO]                                        :: ",
      pn(lo0h))
print("MUONEh+ [NLOp0, NLOp1, dK[NLOp0]%, dK[NLOp1]%]      :: ",
      pn(nlo0ph),  pn(nlo1ph),  pn(dn(plusnumbers(nlo0ph, -lo0h),  lo0h)  *100.), pn(dn(plusnumbers(nlo1ph, -lo0h),  lo0h)  *100.))
print("MUONEh+ [NNLOp0, NNLOp1, dK[NNLOp0]%, dK[NNLOp1]% ] :: ",
      pn(nnlo0ph), pn(nnlo1ph), pn(dn(plusnumbers(nnlo0ph,-nlo0ph),nlo0ph)*100.), pn(dn(plusnumbers(nnlo1ph,-nlo1ph),nlo1ph)*100.))

MUONEh+ [LO]                                        ::  245.59625(1)
MUONEh+ [NLOp0, NLOp1, dK[NLOp0]%, dK[NLOp1]%]      ::  259.6377(2) 237.38539(7) 5.71731(7) -3.34324(3)
MUONEh+ [NNLOp0, NNLOp1, dK[NNLOp0]%, dK[NNLOp1]% ] ::  259.7680(2) 237.44002(7) 0.05018(9) 0.02301(4)

print("                  |    S1'    S2'   |   dK1'%    dK2'%  ")
print("--------------------------------------------------------------------------------------------")
print("\sigma^{(0)}      | ", pn(loh['value']), " ", pn(loh['value']))
print("--------------------------------------------------------------------------------------------")
print("\sigma^{(1)}_e    | ", pn(nloEPHm0h['value']), " ", pn(nloEPHm1h['value'])," | ",   pn(dn(nloEPHm0h['value'],lo0h)*100.), " ", pn(dn(nloEPHm1h['value'],lo0h)*100.))
print("\sigma^{(1)}_emM  | ", pn(nloXPHm0h['value']), " ", pn(nloXPHm1h['value'])," | ",   pn(dn(nloXPHm0h['value'],lo0h)*100.), " ", pn(dn(nloXPHm1h['value'],lo0h)*100.))
print("\sigma^{(1)}_emP  | ", pn(nloXPHp0h['value']), " ", pn(nloXPHp1h['value'])," | ",   pn(dn(nloXPHp0h['value'],lo0h)*100.), " ", pn(dn(nloXPHp1h['value'],lo0h)*100.))
print("\sigma^{(1)}_m    | ", pn(nloMPHm0h['value']), " ", pn(nloMPHm1h['value'])," | ",   pn(dn(nloMPHm0h['value'],lo0h)*100.), " ", pn(dn(nloMPHm1h['value'],lo0h)*100.))
print("\sigma^{(1)}_lep  | ", pn(nloLVPm0h['value']), " ", pn(nloLVPm0h['value'])," | ",   pn(dn(nloLVPm0h['value'],lo0h)*100.), " ", pn(dn(nloLVPm0h['value'],lo0h)*100.))
print("\sigma^{(1)}_had  | ", pn(nloHVPm0h['value']), " ", pn(nloHVPm0h['value'])," | ",   pn(dn(nloHVPm0h['value'],lo0h)*100.), " ", pn(dn(nloHVPm0h['value'],lo0h)*100.))
print("--------------------------------------------------------------------------------------------")
print("\sigma^{(2)}_e    | ", pn(nnloEPHm0h['value']), " ", pn(nnloEPHm1h['value'])," | ", pn(dn(nnloEPHm0h['value'],nlo0mh)*100.), " ", pn(dn(nnloEPHm1h['value'],nlo1mh)*100.))
print("\sigma^{(2)}_emM  | ", pn(nnloYPHm0h), " ", pn(nnloYPHm1h)," | ", pn(dn(nnloYPHm0h,nlo0mh)*100.), " ", pn(dn(nnloYPHm1h,nlo1mh)*100.))
print("\sigma^{(2)}_emP  | ", pn(nnloYPHp0h), " ", pn(nnloYPHp1h)," | ", pn(dn(nnloYPHp0h,nlo0ph)*100.), " ", pn(dn(nnloYPHp1h,nlo1ph)*100.))
print("\sigma^{(2)}_m    | ", pn(nnloMPHm0h['value']), " ", pn(nnloMPHm1h['value'])," | ", pn(dn(nnloMPHm0h['value'],nlo0mh)*100.), " ", pn(dn(nnloMPHm1h['value'],nlo1mh)*100.))
print("\sigma^{(2)}_lepM | ", pn(nnloLVPm0h['value']), " ", pn(nnloLVPm1h['value'])," | ", pn(dn(nnloLVPm0h['value'],nlo0mh)*100.), " ", pn(dn(nnloLVPm1h['value'],nlo1mh)*100.))
print("\sigma^{(2)}_lepP | ", pn(nnloLVPp0h['value']), " ", pn(nnloLVPp1h['value'])," | ", pn(dn(nnloLVPp0h['value'],nlo0ph)*100.), " ", pn(dn(nnloLVPp1h['value'],nlo1ph)*100.))
print("\sigma^{(2)}_hadM | ", pn(nnloHVPm0h['value']), " ", pn(nnloHVPm1h['value'])," | ", pn(dn(nnloHVPm0h['value'],nlo0mh)*100.), " ", pn(dn(nnloHVPm1h['value'],nlo1mh)*100.))
print("\sigma^{(2)}_hadP | ", pn(nnloHVPp0h['value']), " ", pn(nnloHVPp1h['value'])," | ", pn(dn(nnloHVPp0h['value'],nlo0ph)*100.), " ", pn(dn(nnloHVPp1h['value'],nlo1ph)*100.))
print("--------------------------------------------------------------------------------------------")

                  |    S1'    S2'   |   dK1'%    dK2'%  
--------------------------------------------------------------------------------------------
\sigma^{(0)}      |  245.59625(1)   245.59625(1)
--------------------------------------------------------------------------------------------
\sigma^{(1)}_e    |  10.8488(2)   -11.18183(7)  |  4.41735(7)   -4.55293(3)
\sigma^{(1)}_emM  |  -0.350212(1)   -0.285981(2)  |  -0.1425967(6)   -0.1164435(7)
\sigma^{(1)}_emP  |  0.350212(1)   0.285984(2)  |  0.1425966(6)   0.1164449(7)
\sigma^{(1)}_m    |  -0.0666769(5)   -0.2241327(6)  |  -0.0271490(2)   -0.0912606(2)
\sigma^{(1)}_lep  |  2.8896899(2)   2.8896899(2)  |  1.1766018(1)   1.1766018(1)
\sigma^{(1)}_had  |  0.019430463(1)   0.019430463(1)  |  0.0079115470(7)   0.0079115470(7)
--------------------------------------------------------------------------------------------
\sigma^{(2)}_e    |  0.02713(4)   0.174906(10)  |  0.01048(2)   0.073858(4)
\sigma^{(2)}_emM  |  -0.025255(1)   0.0285224(5)  |  -0.0097532(5)   0.0120442(2)
\sigma^{(2)}_emP  |  0.031647(1)   0.0043945(5)  |  0.0121889(5)   0.0018512(2)
\sigma^{(2)}_m    |  -0.000097445(3)   -0.000038927(2)  |  -0.0000376326(10)   -0.0000164376(6)
\sigma^{(2)}_lepM |  0.059108(2)   -0.134451(2)  |  0.0228273(8)   -0.0567751(8)
\sigma^{(2)}_lepP |  0.071609(2)   -0.123800(2)  |  0.0275803(9)   -0.0521514(8)
\sigma^{(2)}_hadM |  -0.000493546(7)   -0.001277840(7)  |  -0.000190604(3)   -0.000539598(3)
\sigma^{(2)}_hadP |  0.000002888(9)   -0.000834486(9)  |  0.000001112(3)   -0.000351532(4)
--------------------------------------------------------------------------------------------

Differential distributions

[see the attached user file for different observables]

def kvalue(ph0, ph1e,ph1xM,ph1xP,ph1m,  hvp1,lvp1,
                ph2e,ph2xM31,ph2xM22,ph2xM13,ph2m,  hvp2m,lvp2m,hvp2p,lvp2p,
           obs,sc,msc,aa,bb):

    LO = mergebins(scaleplot(ph0[obs],sc)[aa:bb],msc)

    NLOe   = mergebins(scaleplot(ph1e[obs],sc)[aa:bb],msc)
    NLOxM  = mergebins(scaleplot(ph1xM[obs],sc)[aa:bb],msc)
    NLOxP  = mergebins(scaleplot(ph1xP[obs],sc)[aa:bb],msc)
    NLOm   = mergebins(scaleplot(ph1m[obs],sc)[aa:bb],msc)
    NLOhvp = mergebins(scaleplot(hvp1[obs],sc)[aa:bb],msc)
    NLOlvp = mergebins(scaleplot(lvp1[obs],sc)[aa:bb],msc)

    NLO  = addplots(scaleplot(ph1e[obs],sc)[aa:bb],LO)
    NLO  = addplots(scaleplot(ph1m[obs],sc)[aa:bb],NLO)
    NLO  = addplots(scaleplot(hvp1[obs],sc)[aa:bb],NLO)
    NLO  = addplots(scaleplot(lvp1[obs],sc)[aa:bb],NLO)
    NLOM = addplots(scaleplot(ph1xM[obs],sc)[aa:bb],NLO)
    NLOP = addplots(scaleplot(ph1xP[obs],sc)[aa:bb],NLO)

    NNLOe    = mergebins(scaleplot(ph2e[obs],sc)[aa:bb],msc)
    NNLOxM31 = mergebins(scaleplot(ph2xM31[obs],sc)[aa:bb],msc)
    NNLOxM22 = mergebins(scaleplot(ph2xM22[obs],sc)[aa:bb],msc)
    NNLOxM13 = mergebins(scaleplot(ph2xM13[obs],sc)[aa:bb],msc)
    NNLOyM   = addplots(NNLOxM22,NNLOxM31)
    NNLOyM   = addplots(NNLOxM13,NNLOyM)
    NNLOyP   = addplots(NNLOxM22,scaleplot(NNLOxM31,1.,-1.))
    NNLOyP   = addplots(scaleplot(NNLOxM13,1.,-1.),NNLOyP)
    NNLOm    = mergebins(scaleplot(ph2m[obs],sc)[aa:bb],msc)
    NNLOhvpM = mergebins(scaleplot(hvp2m[obs],sc)[aa:bb],msc)
    NNLOlvpM = mergebins(scaleplot(lvp2m[obs],sc)[aa:bb],msc)
    NNLOhvpP = mergebins(scaleplot(hvp2p[obs],sc)[aa:bb],msc)
    NNLOlvpP = mergebins(scaleplot(lvp2p[obs],sc)[aa:bb],msc)

    NNLOM = addplots(scaleplot(ph2e[obs],sc)[aa:bb],NLO)
    NNLOM = addplots(scaleplot(ph2m[obs],sc)[aa:bb],NNLOM)
    NNLOM = addplots(scaleplot(hvp2m[obs],sc)[aa:bb],NNLOM)
    NNLOM = addplots(scaleplot(lvp2m[obs],sc)[aa:bb],NNLOM)
    NNLOM = addplots(scaleplot(ph2xM31[obs],sc)[aa:bb],NNLOM)
    NNLOM = addplots(scaleplot(ph2xM22[obs],sc)[aa:bb],NNLOM)
    NNLOM = addplots(scaleplot(ph2xM13[obs],sc)[aa:bb],NNLOM)

    return ((NLOe,NLOxM,NLOxP,NLOm,NNLOe,NNLOm,NNLOxM31,NNLOxM22,NNLOxM13,NNLOyM,NNLOyP),
            (NLOhvp,NNLOhvpM,NNLOhvpP),
            (NLOlvp,NNLOlvpM,NNLOlvpP),
            (LO,NLOM,NLOP,NNLOM))

def diff_plot(obs, ecut, sc, msc, aa, bb,
              xlab, ylab,
              x1,y1, x2,y2, x3,y3, x4,y4, mulex,muley,
              l1, l2, l3,
              fname):

    if ecut == 0:
        K = kvalue(lo,nloEPHm0,nloXPHm0,nloXPHp0,nloMPHm0,nloHVPm0,nloLVPm0,
                     nnloEPHm0,nnloX31PHm0,nnloX22PHm0,nnloX13PHm0,nnloMPHm0,
                     nnloHVPm0,nnloLVPm0,nnloHVPp0,nnloLVPp0,
                     obs,sc,msc,aa,bb)
    elif ecut == 1:
        K = kvalue(lo,nloEPHm1,nloXPHm1,nloXPHp1,nloMPHm1,nloHVPm0,nloLVPm0,
                 nnloEPHm1,nnloX31PHm1,nnloX22PHm1,nnloX13PHm1,nnloMPHm1,
                 nnloHVPm1,nnloLVPm1,nnloHVPp1,nnloLVPp1,
                 obs,sc,msc,aa,bb)
    elif ecut == 2:
        K = kvalue(loh,nloEPHm0h,nloXPHm0h,nloXPHp0h,nloMPHm0h, nloHVPm0h,nloLVPm0h,
                   nnloEPHm0h,nnloX31PHm0h,nnloX22PHm0h,nnloX13PHm0h,nnloMPHm0h,
                   nnloHVPm0h,nnloLVPm0h,nnloHVPp0h,nnloLVPp0h,
                   obs,sc,msc,aa,bb)

    fig, axs = subplots(
    5, sharex=True, gridspec_kw={'hspace': 0, 'height_ratios':[.9,1,1,1,1]},
    figsize=(5.9, 7.84)
    )

    axs = list(axs)

    sca(axs[0])
    errorband(K[3][0],col='purple')
    errorband(K[3][3],col='orange')
    axs[0].legend([r'{\rm LO}',r'{\rm NNLO$-$}'],ncol=2, loc=l1)
    ylabel(ylab)

    sca(axs[1])
    errorband(divideplots(K[0][0],K[3][0]))
    errorband(divideplots(K[0][1],K[3][0]))
    errorband(divideplots(K[0][3],K[3][0]))
    errorband(divideplots(addplots(K[1][0],K[2][0]),K[3][0]))
    axs[1].legend([r'{\rm electronic}',r'{\rm mixed}',r'{\rm muonic}',r'{\rm fermionic}'], ncol=2,loc=l2)
    axs[1].axhline(0, color='black', linewidth=0.4, zorder=1)
    ylabel("$\delta K_{\\rm NLO}$")
    axs.append(gca().twinx())
    errorband(scaleplot(divideplots(K[1][0],K[3][0]),1.,1.e+2),col='mediumorchid')
    gca().tick_params(labelcolor='mediumorchid', axis='y')
    ylabel(r'$10^2 \times \delta K_{\rm NLO}^{\rm had}$', color='mediumorchid')
    mpl_axes_aligner.yaxes(axs[1], gca(), 0,0)

    sca(axs[2])
    errorband(divideplots(K[0][4], K[3][1]))
    errorband(divideplots(K[0][9],K[3][1]))
    errorband(divideplots(K[0][5], K[3][1]))
    errorband(divideplots(addplots(K[1][1],K[2][1]),K[3][1]))
    axs[2].axhline(0, color='black', linewidth=0.4, zorder=1)
    ylabel("$\delta K_{\\rm NNLO}$")

    sca(axs[3])
    errorband(divideplots(K[0][4], K[3][2]))
    errorband(divideplots(K[0][10],K[3][2]))
    errorband(divideplots(K[0][5], K[3][2]))
    errorband(divideplots(addplots(K[1][2],K[2][2]),K[3][2]))
    axs[3].axhline(0, color='black', linewidth=0.4, zorder=1)
    ylabel("$\delta K_{\\rm NNLO}$")

    sca(axs[4])
    errorband(scaleplot(divideplots(K[0][6], K[3][1]),1.,1.e+1),col='dodgerblue')
    errorband(scaleplot(divideplots(K[0][7], K[3][1]),1.,1.e+1),col='lightseagreen')
    errorband(scaleplot(divideplots(K[0][8],K[3][1]),1.,1.e+1),col='olivedrab')
    ylabel(r'$10 \times \delta K_{\rm NNLO}$')
    axs[4].legend(["$q^3\,Q$","$q^2\,Q^2$","$q\,Q^3$"],ncol=3,loc=l3)
    axs[4].axhline(0, color='black', linewidth=0.4, zorder=1)

    xlabel(xlab)
    figtext(x1,y1,r'$\boxed{e^-\,\mu^-\to e^-\,\mu^-}$')
    figtext(x2,y2,r'$\boxed{e^-\,\mu^-\to e^-\,\mu^-}$')
    figtext(x3,y3,r'$\boxed{e^-\,\mu^+\to e^-\,\mu^+}$')
    figtext(x4,y4,r'$\boxed{e^-\,\mu^-\to e^-\,\mu^-}$')
    draw()
    fig.tight_layout()
    mulify(fig, delx=mulex, dely=muley)
    fig.savefig(fname)

$\theta_e$ in S1

diff_plot('the',0,1e-3,3,5,170,
          '$\\theta_e\,/\,{\\rm mrad}$','$\\D\\sigma/\\D\\theta_e\ /\ {\\rm\\upmu b}$',
          .70,.710, .70,.585, .70,.400, .70,.120, .1,-.4,
          'lower right', 'upper right', 'lower right',
          'plots/the-s1.pdf')

$t_\mu = (p_2-p_4)^2$ in S1

diff_plot('tmm',0,1e+6,5,65,295,
          '$t_\\mu\,/\,{\\rm GeV^2}$','$\\D\\sigma/\\D t_\\mu\ /\ {\\rm\\upmu b}$',
          .68,.635, .68,.585, .68,.400, .68,.220, .1,-.4,
          'upper center', 'upper right', 'lower right',
          'plots/tmu-s1.pdf')

$\theta_\mu$ in S1

diff_plot('thm',0,1e-3,5,15,200,
          '$\\theta_\\mu\,/\,{\\rm mrad}$','$\\D\\sigma/\\D\\theta_\\mu\ /\ {\\rm\\upmu b}$',
          .18,.770, .18,.585, .18,.330, .18,.215, .6,-.4,
          'upper right', 'upper right', 'lower left',
          'plots/thm-s1.pdf')

$\theta_\mu$ in S1'

diff_plot('thm',2,1e-3,5,2,200,
          '$\\theta_\\mu\,/\,{\\rm mrad}$','$\\D\\sigma/\\D\\theta_\\mu\ /\ {\\rm\\upmu b}$',
          .24,.710, .24,.585, .24,.400, .24,.210, .6,-.4,
          'upper right', 'upper left', 'lower right',
          'plots/thm-s1p.pdf')

$\theta_e$ in S2

diff_plot('the',1,1e-3,3,5,170,
          '$\\theta_e\,/\,{\\rm mrad}$','$\\D\\sigma/\\D\\theta_e\ /\ {\\rm\\upmu b}$',
          .68,.710, .68,.450, .68,.265, .68,.120, .1,-.4,
          'lower right', 'upper right', 'lower right',
          'plots/the-s2.pdf')

$t_\mu$ in S2

diff_plot('tmm',1,1e+6,5,65,295,
          '$t_\\mu\,/\,{\\rm GeV^2}$','$\\D\\sigma/\\D t_\\mu\ /\ {\\rm\\upmu b}$',
          .18,.725, .18,.490, .18,.340, .68,.220, .1,-.4,
          'upper center', 'upper right', 'lower right',
          'plots/tmu-s2.pdf')

Validation plots

Loading xic data

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic-.tar.bz2')

def dat_xic(ff,rf1,rf2,rr1,rr2,obs):

    ffa  = np.array([scaleset(sigma(ff,obs='0')[(.1,.1)], alpha**4*conv),
                     scaleset(sigma(ff,obs='0')[(.5,.5)], alpha**4*conv),
                     scaleset(sigma(ff,obs='0')[(1.,1.)], alpha**4*conv)])
    rf1a = np.array([scaleset(sigma(rf1)[(.1,.1)], alpha**4*conv),
                     scaleset(sigma(rf1)[(.5,.5)], alpha**4*conv),
                     scaleset(sigma(rf1)[(1.,1.)], alpha**4*conv)])
    rf2a = np.array([scaleset(sigma(rf2)[(.1,.1)], alpha**4*conv),
                     scaleset(sigma(rf2)[(.5,.5)], alpha**4*conv),
                     scaleset(sigma(rf2)[(1.,1.)], alpha**4*conv)])
    rr1a = np.array([scaleset(sigma(rr1)[(.1,.1)], alpha**4*conv),
                     scaleset(sigma(rr1)[(.5,.5)], alpha**4*conv),
                     scaleset(sigma(rr1)[(1.,1.)], alpha**4*conv)])
    rr2a = np.array([scaleset(sigma(rr2)[(.1,.1)], alpha**4*conv),
                     scaleset(sigma(rr2)[(.5,.5)], alpha**4*conv),
                     scaleset(sigma(rr2)[(1.,1.)], alpha**4*conv)])

    nnlo = []
    for i in range(len(ffa)):
        nnlo.append(ffa[i][obs])
        nnlo[i] = addplots(rf1a[i][obs],nnlo[i])
        nnlo[i] = addplots(rf2a[i][obs],nnlo[i])
        nnlo[i] = addplots(rr1a[i][obs],nnlo[i])
        nnlo[i] = addplots(rr2a[i][obs],nnlo[i])

    return nnlo

obs = 'the'
aa = 5
bb = 175
sc = 1e-3
msc = 1

def val_xic(ff,rf1,rf2,rr1,rr2,obs,sc,msc,aa,bb):

    re  = dat_xic(ff,rf1,rf2,rr1,rr2,obs)
    re0 = scaleset(mergefks(sigma(ff),sigma(rf1),sigma(rf2),sigma(rr1),sigma(rr2)),alpha**4*conv)

    xi0 = mergebins(scaleplot(re0[obs][aa:bb],sc),msc)
    xi1 = mergebins(scaleplot(re[0][aa:bb]   ,sc),msc)
    xi2 = mergebins(scaleplot(re[1][aa:bb]   ,sc),msc)
    xi3 = mergebins(scaleplot(re[2][aa:bb]   ,sc),msc)

    return [xi0,xi1,xi2,xi3]

xicae = val_xic('FFEEEE','RFEEEE15','RFEEEE35','RREEEE1516','RREEEE3536', 'the',sc,msc,aa,bb)
xicax = val_xic('FFMIXDz','RFMIXD15','RFMIXD35','RRMIXD1516','RRMIXD3536','the',sc,msc,aa,bb)

Loading massification data

obs = 'the'
aa = 5
bb = 170
sc = 1e-3
msc = 5

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic-.tar.bz2')
nnloEPHm0_fu = scaleset(mergefks(
    sigma('em2emFFEEEE'),
    sigma('em2emRFEEEE15'), sigma('em2emRFEEEE35'),
    sigma('em2emRREEEE1516'), sigma('em2emRREEEE3536')), alpha**4*conv)
nnloEPHm0_fu_ff = scaleset(mergefks(sigma('em2emFFEEEE')), alpha**4*conv)
nnloEPHm0_ma = scaleset(mergefks(
    sigma('em2emFFEEEEz'),
    sigma('em2emRFEEEE15'), sigma('em2emRFEEEE35'),
    sigma('em2emRREEEE1516'), sigma('em2emRREEEE3536')), alpha**4*conv)
nnloEPHm0_ma_ff = scaleset(mergefks(sigma('em2emFFEEEEz')), alpha**4*conv)
ma = mergebins(scaleplot(nnloEPHm0_ma[obs][aa:bb],sc),msc)
fu = mergebins(scaleplot(nnloEPHm0_fu[obs][aa:bb],sc),msc)
ma_ff = mergebins(scaleplot(nnloEPHm0_ma_ff[obs][aa:bb],sc),msc)
fu_ff = mergebins(scaleplot(nnloEPHm0_fu_ff[obs][aa:bb],sc),msc)

the error for the ratio ma/fu is only due to the double virtual

ma[:,2] = ma_ff[:,2]
fu[:,2] = fu_ff[:,2]

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic-_thm.tar.bz2')
nnloEPHm0h_fu = scaleset(mergefks(
    sigma('em2emFFEEEE'),
    sigma('em2emRFEEEE15'), sigma('em2emRFEEEE35'),
    sigma('em2emRREEEE1516'), sigma('em2emRREEEE3536')), alpha**4*conv)
nnloEPHm0h_fu_ff = scaleset(mergefks(sigma('em2emFFEEEE')), alpha**4*conv)
nnloEPHm0h_ma = scaleset(mergefks(
    sigma('em2emFFEEEEz'),
    sigma('em2emRFEEEE15'), sigma('em2emRFEEEE35'),
    sigma('em2emRREEEE1516'), sigma('em2emRREEEE3536')), alpha**4*conv)
nnloEPHm0h_ma_ff = scaleset(mergefks(sigma('em2emFFEEEEz')), alpha**4*conv)
mah = mergebins(scaleplot(nnloEPHm0h_ma[obs][aa:bb],sc),msc)
fuh = mergebins(scaleplot(nnloEPHm0h_fu[obs][aa:bb],sc),msc)
mah_ff = mergebins(scaleplot(nnloEPHm0h_ma_ff[obs][aa:bb],sc),msc)
fuh_ff = mergebins(scaleplot(nnloEPHm0h_fu_ff[obs][aa:bb],sc),msc)

the error for the ratio mah/fuh is only due to the double virtual

mah[:,2] = mah_ff[:,2]
fuh[:,2] = fuh_ff[:,2]

Loading NTS stabilisation data

obs = 'the'
aa = 5
bb = 170
sc = 1e-3
msc = 5

setup(obs='0')
setup(folder='em2em_muone_paper/out-photonic-.tar.bz2')
nnloXPHm0_o2 = scaleset(mergefks(
    sigma('em2emFFMIXDz'),
    sigma('em2emRFMIXD15'), sigma('em2emRFMIXD35'),
    sigma('em2emRRMIXD1516'), sigma('em2emRRMIXD3536')), alpha**4*conv)
nnloXPHm0_o2_rf = scaleset(mergefks(sigma('em2emRFMIXD15'),sigma('em2emRFMIXD35')), alpha**4*conv)
setup(folder='em2em_muone_paper/out-photonic-_ol4.tar.bz2')
nnloXPHm0_o4 = scaleset(mergefks(
    sigma('em2emFFMIXDz'),
    sigma('em2emRFMIXD15'), sigma('em2emRFMIXD35'),
    sigma('em2emRRMIXD1516'), sigma('em2emRRMIXD3536')), alpha**4*conv)
nnloXPHm0_o4_rf = scaleset(mergefks(sigma('em2emRFMIXD15'),sigma('em2emRFMIXD35')), alpha**4*conv)
o2 = mergebins(scaleplot(nnloXPHm0_o2[obs][aa:bb],sc),msc)
o4 = mergebins(scaleplot(nnloXPHm0_o4[obs][aa:bb],sc),msc)
o2_rf = mergebins(scaleplot(nnloXPHm0_o2_rf[obs][aa:bb],sc),msc)
o4_rf = mergebins(scaleplot(nnloXPHm0_o4_rf[obs][aa:bb],sc),msc)

the error for the ratio o2/o4 is only due to the real virtual

o2[:,2] = o2_rf[:,2]
o4[:,2] = o4_rf[:,2]

Plotting

From top to bottom:

• electronic (S1 and S1') and mixed (S1) NNLO correction
• $\xi_c$ (in)dependence plot for nnloEPHm0
• $\xi_c$ (in)dependence plot for nnloXPHm0
• massification study
• NTS stabilisation study

fig, axs = subplots(
    5, sharex=True, gridspec_kw={'hspace': 0, 'height_ratios':[.9,1,1,1,1]},
    figsize=(5.9, 7.84)
)

axs = list(axs)
sca(axs[0])
errorband(fu,col='darkred')
errorband(fuh,col='mediumorchid')
axs[0].legend(["{\\tt S1}", "{\\tt S1}'"],loc='upper center',fontsize=8,framealpha=0)
axs[0].axhline(0, color='black', linewidth=0.4, zorder=1)
ylabel('$\\D\\sigma/\\D\\theta_e\ /{\\rm\\upmu b} \ {\\tt [ELECTRONIC]}$', color='darkred')
axs.append(gca().twinx())
errorband(scaleplot(o2,1.,1.e+2),col='green')
gca().tick_params(axis='y')
ylabel('$10^2 \\times \\D\\sigma/\\D\\theta_e\ /{\\rm\\upmu b} \ {\\tt [MIXED]}$', color='green')
mpl_axes_aligner.yaxes(axs[0], gca(), 0,0)

sca(axs[1])
dat = np.array([mergebins(i, 15) for i in xicae])
norm = dat[0,:,1]
col = colours.alpha_composite('white', 'darkred', 0.2)
fill_between(
     dat[0,:,0],
    (dat[0,:,1] + dat[0,:,2]/2) / norm - 1,
    (dat[0,:,1] - dat[0,:,2]/2) / norm - 1,
    step='mid',
    facecolor=col, edgecolor=col
)
for i, d in enumerate(dat[1:]):
    errorbar(
        0.4*(i-1)+d[:,0],
        d[:,1] / norm - 1,
        abs(d[:,2] / norm),
        linestyle='none', marker='.'
    )
axs[1].legend(['$\\xi_c\\pm 1\\sigma$','$\\xi_i=0.1$','$\\xi_i=0.5$','$\\xi_i=1.0$'],
              ncol=2,loc='upper right',fontsize=8)
axs[1].axhline(0, color='black', linewidth=0.4, zorder=1)
ylim(-0.06,0.07)
ylabel('$1-\\xi_i/\\xi_c$', color='darkred')
for a,b in list(zip(dat[0,1:,0], dat[0,:-1,0])):
    axvline((a+b)/2, color='black', linewidth=0.4, zorder=1)

sca(axs[2])
dat = np.array([mergebins(i, 15) for i in xicax])
norm = dat[0,:,1]
col = colours.alpha_composite('white', 'green', 0.2)
fill_between(
     dat[0,:,0],
    (dat[0,:,1] + dat[0,:,2]/2) / norm - 1,
    (dat[0,:,1] - dat[0,:,2]/2) / norm - 1,
    step='mid',
    facecolor=col, edgecolor=col
)
for i, d in enumerate(dat[1:]):
    errorbar(
        0.4*(i-1)+d[:,0],
        d[:,1] / norm - 1,
        abs(d[:,2] / norm),
        linestyle='none', marker='.'
    )
axs[2].legend(['$\\xi_c\\pm 1\\sigma$','$\\xi_i=0.1$','$\\xi_i=0.5$','$\\xi_i=1.0$'],
              ncol=2,loc='lower right',fontsize=8)
axs[2].axhline(0, color='black', linewidth=0.4, zorder=1)
ylim(-0.3,0.15)
ylabel('$1-\\xi_i/\\xi_c$', color='green')
for a,b in list(zip(dat[0,1:,0], dat[0,:-1,0])):
    axvline((a+b)/2, color='black', linewidth=0.4, zorder=1)

sca(axs[3])
dat=divideplots(ma,fu,offset=-1.)
col='darkred'
al=0.2
fill_between(
  dat[:,0], dat[:,1]+dat[:,2]/2, dat[:,1]-dat[:,2]/2,
  step='mid',
  facecolor=colours.alpha_composite('white', col, al),
  edgecolor=colours.alpha_composite('white', col, al)
)
step(dat[:,0], dat[:,1], col, where='mid')
dath=divideplots(mah,fuh,offset=-1.)
col='mediumorchid'
al=0.2
fill_between(
  dath[:,0], dath[:,1]+dath[:,2]/2, dath[:,1]-dath[:,2]/2,
  step='mid',
  facecolor=colours.alpha_composite('white', col, al),
  edgecolor=colours.alpha_composite('white', col, al)
)
step(dath[:,0], dath[:,1], col, where='mid')
axs[3].axhline(0, color='black', linewidth=0.4, zorder=1)
ylabel("$\\delta_Z$", color='darkred')
ylim(-0.004,0.0039)

sca(axs[4])
dat=divideplots(o2,o4,offset=-1.)
col='green'
al=0.2
fill_between(
  dat[:,0], dat[:,1]+dat[:,2]/2, dat[:,1]-dat[:,2]/2,
  step='mid',
  facecolor=colours.alpha_composite('white', col, al),
  edgecolor=colours.alpha_composite('white', col, al)
)
step(dat[:,0], dat[:,1], col, where='mid')
axs[4].axhline(0, color='black', linewidth=0.4, zorder=1)
ylabel("$\\delta_{\\rm NTS}$", color='green')
ylim(-0.05,0.05)

xlabel(r'$\theta_e\,/\,{\rm mrad}$')
figtext(.16,.620,"$\\boxed{\\tt ELECTRONIC}$", color='darkred')
figtext(.16,.445,"$\\boxed{\\tt MIXED}$", color='green')
figtext(.16,.265,"$\\boxed{\\tt ELECTRONIC}$", color='darkred')
figtext(.16,.085,"$\\boxed{\\tt MIXED}$", color='green')
draw()
fig.tight_layout()
mulify(fig, delx=.15, dely=-.25)
fig.savefig('plots/val-plot.pdf')

Splitting validation plots

$\xi_c$ (in)dependence plot for nnloXPHm0

fig, axs = subplots(
    2, sharex=True, gridspec_kw={'hspace': 0, 'height_ratios':[.9,1]},
    figsize=(5.9, 3.84)
)

axs = list(axs)
sca(axs[0])
errorband(fu,col='darkred')
axs[0].axhline(0, color='black', linewidth=0.4, zorder=1)
ylabel('$\\D\\sigma/\\D\\theta_e\ /{\\rm\\upmu b}$', color='darkred')

dat = np.array([mergebins(i, 15) for i in xicae])
norm = dat[0,:,1]
sca(axs[1])
col = colours.alpha_composite('white', 'darkred', 0.2)
fill_between(
     dat[0,:,0],
    (dat[0,:,1] + dat[0,:,2]/2) / norm - 1,
    (dat[0,:,1] - dat[0,:,2]/2) / norm - 1,
    step='mid',
    facecolor=col, edgecolor=col
)
for i, d in enumerate(dat[1:]):
    errorbar(
        0.4*(i-1)+d[:,0],
        d[:,1] / norm - 1,
        abs(d[:,2] / norm),
        linestyle='none', marker='.'
    )
axs[1].legend(['$\\xi_c\\pm 1\\sigma$','$\\xi_i=0.1$','$\\xi_i=0.5$','$\\xi_i=1.0$'],
              ncol=2,loc='upper right',fontsize=8)
axs[1].axhline(0, color='black', linewidth=0.4, zorder=1)

plt.ylim(-0.06,0.07)
plt.ylabel('$1-\\xi_i/\\xi_c$', color='darkred')
for a,b in list(zip(dat[0,1:,0], dat[0,:-1,0])):
    plt.axvline((a+b)/2, color='black', linewidth=0.4, zorder=1)

xlabel(r'$\theta_e\,/\,{\rm mrad}$')
figtext(.16,.600,"$\\boxed{\\tt NNLO \,\, ELECTRONIC}$", color='darkred')
draw()
fig.tight_layout()
mulify(fig, delx=.0, dely=-.4)
fig.savefig('plots/val-plot-xicE.pdf')

$\xi_c$ (in)dependence plot for nnloXPHm0

fig, axs = plt.subplots(
    2, sharex=True, gridspec_kw={'hspace': 0, 'height_ratios':[.9,1]},
    figsize=(5.9, 3.84)
)

axs = list(axs)
sca(axs[0])
axs[0].axhline(0, color='black', linewidth=0.4, zorder=1)
errorband(scaleplot(o2,1.,1.e+2),col='green')
ylabel('$10^2 \\times \\D\\sigma/\\D\\theta_e\ /{\\rm\\upmu b}$', color='green')

dat = np.array([mergebins(i, 15) for i in xicax])
norm = dat[0,:,1]
sca(axs[1])
col = colours.alpha_composite('white', 'green', 0.2)
fill_between(
     dat[0,:,0],
    (dat[0,:,1] + dat[0,:,2]/2) / norm - 1,
    (dat[0,:,1] - dat[0,:,2]/2) / norm - 1,
    step='mid',
    facecolor=col, edgecolor=col
)
for i, d in enumerate(dat[1:]):
    errorbar(
        0.4*(i-1)+d[:,0],
        d[:,1] / norm - 1,
        abs(d[:,2] / norm),
        linestyle='none', marker='.'
    )
axs[1].legend(['$\\xi_c\\pm 1\\sigma$','$\\xi_i=0.1$','$\\xi_i=0.5$','$\\xi_i=1.0$'],
              ncol=2,loc='lower right',fontsize=8)
axs[1].axhline(0, color='black', linewidth=0.4, zorder=1)
ylim(-0.3,0.15)
ylabel('$1-\\xi_i/\\xi_c$', color='green')
for a,b in list(zip(dat[0,1:,0], dat[0,:-1,0])):
    axvline((a+b)/2, color='black', linewidth=0.4, zorder=1)

xlabel(r'$\theta_e\,/\,{\rm mrad}$')
figtext(.815,.600,"$\\boxed{\\tt NNLO\,\,MIXED}$", color='green')
draw()
fig.tight_layout()
mulify(fig, delx=.0, dely=-.4)
fig.savefig('plots/val-plot-xicX.pdf')

Massification and NTS stabilisation

fig, axs = plt.subplots(
    3, sharex=True, gridspec_kw={'hspace': 0, 'height_ratios':[.9,1,1]},
    figsize=(6.0, 4.0)
)

axs = list(axs)
sca(axs[0])
errorband(fu,col='darkred')
errorband(fuh,col='mediumorchid')
axs[0].legend(["{\\tt S1}", "{\\tt S1}'"],loc='upper center',fontsize=8,framealpha=0)
axs[0].axhline(0, color='black', linewidth=0.4, zorder=1)
ylabel('$\\D\\sigma/\\D\\theta_e\ /{\\rm\\upmu b} \ {\\tt [E]}$', color='darkred')
axs.append(plt.gca().twinx())
errorband(scaleplot(o2,1.,1.e+2),col='green')
gca().tick_params(axis='y')
ylabel('$10^2 \\times \\D\\sigma/\\D\\theta_e\ /{\\rm\\upmu b} \ {\\tt [X]}$', color='green')
mpl_axes_aligner.yaxes(axs[0], plt.gca(), 0,0)

sca(axs[1])
dat=divideplots(ma,fu,offset=-1.)
col='darkred'
al=0.2
fill_between(
  dat[:,0], dat[:,1]+dat[:,2]/2, dat[:,1]-dat[:,2]/2,
  step='mid',
  facecolor=colours.alpha_composite('white', col, al),
  edgecolor=colours.alpha_composite('white', col, al)
)
step(dat[:,0], dat[:,1], col, where='mid')
dath=divideplots(mah,fuh,offset=-1.)
col='mediumorchid'
al=0.2
fill_between(
  dath[:,0], dath[:,1]+dath[:,2]/2, dath[:,1]-dath[:,2]/2,
  step='mid',
  facecolor=colours.alpha_composite('white', col, al),
  edgecolor=colours.alpha_composite('white', col, al)
)
step(dath[:,0], dath[:,1], col, where='mid')
axs[1].axhline(0, color='black', linewidth=0.4, zorder=1)
ylim(-0.004,0.004)
ylabel("$\\delta_Z$", color='darkred')

sca(axs[2])
dat=divideplots(o2,o4,offset=-1.)
col='green'
al=0.2
plt.fill_between(
  dat[:,0], dat[:,1]+dat[:,2]/2, dat[:,1]-dat[:,2]/2,
  step='mid',
  facecolor=colours.alpha_composite('white', col, al),
  edgecolor=colours.alpha_composite('white', col, al)
)
plt.step(dat[:,0], dat[:,1], col, where='mid')
axs[2].axhline(0, color='black', linewidth=0.4, zorder=1)
ylim(-0.05,0.05)
ylabel("$\\delta_{\\rm NTS}$", color='green')

xlabel(r'$\theta_e\,/\,{\rm mrad}$')
figtext(.17,.430,"$\\boxed{\\tt ELECTRONIC}$", color='darkred')
figtext(.17,.160,"$\\boxed{\\tt MIXED}$", color='green')
draw()
fig.tight_layout()
mulify(fig, delx=.2, dely=-.2)
fig.savefig('plots/val-plot-mass-nts.pdf')