Improved comparison

michel-decay/f-and-g/f-and-g.py

@@ -154,29 +154,35 @@ flo, glo, tlo = getfandg(lfunclo, 'lo')
 fnlo, gnlo, tnlo = getfandg(lfuncnlo, 'nlo')
 fnnlo, gnnlo, tnnlo = getfandg(lfuncnnlo, 'nnlo')
 
+Ee = flo[:, 0]
+mask = (Mel < Ee) & (Ee < Mmu/2 + Mel**2/Mmu**2)
+
+flo = flo[mask]
+glo = glo[mask]
+fnlo = fnlo[mask]
+gnlo = gnlo[mask]
+fnnlo = fnnlo[mask]
+gnnlo = gnnlo[mask]
 
 F = addplots(flo, addplots(fnlo, fnnlo))
 G = addplots(glo, addplots(gnlo, gnnlo))
 ###########################################################}}}
 ##########################################################################}}}
-## Compare at LO and NLO{{{
-# Here we compare the results with Arbuzov [hep-ph/0110047]
-Ee = flo[:, 0]
+## Compare at LO, NLO, and NNLO{{{
+# Here we compare the results with Arbuzov
+Ee = Ee[mask]
+
 xe = 2*Ee/Mmu
 beta = sqrt(1-Mel**2/Ee**2)
 
 ### Compare at LO{{{
 fLOref = 2/Mmu * xe**2*beta * (3-2*xe+xe/4*(3*xe-4)*(1-beta**2))
 gLOref = 2/Mmu * xe**2*beta * ((1-2*xe)*beta + 3*xe**2/4 * beta * (1-beta**2))
-fLOref[Ee < Mel] = 0.
-gLOref[Ee < Mel] = 0.
-fLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
-gLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
 
 
 twopanel(
     "$E_e$", labupleft="$f_0$",
-    upleft=[flo, np.column_stack((Ee, fLOref, np.zeros((4000))))],
+    upleft=[flo, np.column_stack((Ee, fLOref, np.zeros(Ee.shape)))],
     downleft=[np.column_stack((Ee, 1-flo[:, 1]/fLOref, flo[:, 2]/fLOref))]
 )
 ylim(-1e-3, 1e-3)
@@ -184,7 +190,7 @@ ylim(-1e-3, 1e-3)
 
 twopanel(
     "$E_e$", labupleft="$g_0$",
-    upleft=[glo, np.column_stack((Ee, gLOref, np.zeros((4000))))],
+    upleft=[glo, np.column_stack((Ee, gLOref, np.zeros(Ee.shape)))],
     downleft=[np.column_stack((Ee, 1-glo[:, 1]/gLOref, glo[:, 2]/gLOref))]
 )
 ylim(-5e-2, 5e-2)
@@ -192,31 +198,39 @@ ylim(-5e-2, 5e-2)
 ### Compare at NLO{{{
 fNLOref = alpha/(2*pi) * 2/Mmu * xe**2*beta * arbuzov.fnlo(xe, Mel/Mmu)
 gNLOref = alpha/(2*pi) * 2/Mmu * xe**2*beta * arbuzov.gnlo(xe, Mel/Mmu)
-fNLOref[Ee < Mel] = 0.
-gNLOref[Ee < Mel] = 0.
-fNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
-gNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
 
 twopanel(
     "$E_e$", labupleft="$f_1$",
-    upleft=[fnlo, np.column_stack((Ee, fNLOref, np.zeros((4000))))],
+    upleft=[fnlo, np.column_stack((Ee, fNLOref, np.zeros(Ee.shape)))],
     downleft=[np.column_stack((Ee, 1-fnlo[:, 1]/fNLOref, fnlo[:, 2]/fNLOref))]
 )
 ylim(-2e-3, 2e-3)
 twopanel(
     "$E_e$", labupleft="$g_1$",
-    upleft=[gnlo, np.column_stack((Ee, gNLOref, np.zeros((4000))))],
+    upleft=[gnlo, np.column_stack((Ee, gNLOref, np.zeros(Ee.shape)))],
     downleft=[np.column_stack((Ee, 1-gnlo[:, 1]/gNLOref, gnlo[:, 2]/gNLOref))]
 )
 ylim(-2e-3, 2e-3)
 ###########################################################}}}
 ### Compare at NNLO{{{
-fNNLOref = (alpha/(2*pi))**2 * 2/Mmu * (np.log(Mmu**2/Mel**2)-1)**2/2 * arbuzov.f2LL(xe)  # nopep8
-gNNLOref = (alpha/(2*pi))**2 * 2/Mmu * (np.log(Mmu**2/Mel**2)-1)**2/2 * arbuzov.g2LL(xe)  # nopep8
-fNNLOref[Ee < Mel] = 0.
-gNNLOref[Ee < Mel] = 0.
-fNNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
-gNNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
+L = np.log(Mmu**2/Mel**2)
+pref2loop = (alpha/(2*pi))**2 * 2/Mmu
+fNNLOref = pref2loop * (0.5*L**2*arbuzov.f2LL(xe) + L*arbuzov.f2NLL(xe))
+gNNLOref = pref2loop * (0.5*L**2*arbuzov.g2LL(xe) + L*arbuzov.g2NLL(xe))
+
+figure()
+plot(Ee, fNNLOref)
+errorband(fnnlo)
+ylim(-2e-5, 5e-5)
+xlabel("$E_e$")
+ylabel("$f_2$")
+
+figure()
+plot(Ee, gNNLOref)
+errorband(gnnlo)
+ylim(-2e-5, 5e-5)
+xlabel("$E_e$")
+ylabel("$g_2$")
 ###########################################################}}}
 
 ##########################################################################}}}