Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

IC CC #109

Merged
merged 2 commits into from
May 5, 2021
Merged

IC CC #109

Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions benchmarks/runners/sandbox.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,11 +79,11 @@ def generate_observables():
kinematics = []
kinematics.extend(
#[dict(x=0.1,Q2=90, y=0)]
[dict(x=x, Q2=20.0) for x in xgrid[::5]]
#[dict(x=x, Q2=20.0, y=0) for x in np.geomspace(1e-4, 1, 10)]
#[dict(x=x, Q2=20.0) for x in xgrid[:-1:5]]
[dict(x=x, Q2=20.0, y=0) for x in np.geomspace(1e-4, .9, 10)]
)
kinematics.extend(
[dict(x=x, Q2=1.51**2, y=0) for x in np.geomspace(1e-4, 1, 10)]
[dict(x=x, Q2=1.51**2, y=0) for x in np.geomspace(1e-4, .9, 10)]
)
# kinematics.extend([dict(x=x, Q2=90) for x in np.linspace(.8, .99, 10).tolist()])
#kinematics.extend([dict(x=0.10914375746330703, Q2=Q2) for Q2 in np.geomspace(4, 1e3, 10).tolist()])
Expand All @@ -94,12 +94,12 @@ def generate_observables():
kinematics.extend([dict(x=0.001, Q2=Q2,y=0) for Q2 in np.geomspace(4, 20, 10).tolist()])
observable_names = [
#"F2_light",
"F2_charm",
#"F2_charm",
# "F2_bottom",
# "F2_top",
#"F2_total",
#"FL_light",
"FL_charm",
#"FL_charm",
# "FL_bottom",
# "FL_total",
#"F3_light",
Expand All @@ -124,7 +124,7 @@ def generate_observables():
return observables.build(observable_names=observable_names,kinematics=kinematics,update=update)

def doit(self):
self.run([{"PTO": 1, "FNS": "FONLL-A", "NfFF": 4, "mc": 1.51}], self.generate_observables(), ["gonly"])
self.run([{"IC":1,"PTO": 1, "mc": 1.51}], self.generate_observables(), ["conly"])

if __name__ == "__main__":
sand = Sandbox()
Expand Down
12 changes: 12 additions & 0 deletions docs/sphinx/source/refs/ic.bib
Original file line number Diff line number Diff line change
Expand Up @@ -11,3 +11,15 @@ @article{luca-intrinsic
pages = "122",
year = "2015"
}
@article{kretzer-schienbein,
author = "Kretzer, S. and Schienbein, I.",
title = "{Heavy quark initiated contributions to deep inelastic structure functions}",
eprint = "hep-ph/9805233",
archivePrefix = "arXiv",
reportNumber = "DO-TH-98-05",
doi = "10.1103/PhysRevD.58.094035",
journal = "Phys. Rev. D",
volume = "58",
pages = "094035",
year = "1998"
}
2 changes: 1 addition & 1 deletion extras/intrinsic/ic/data.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,7 @@

M2 = """Splus * Del * x / (Q2IC)"""

M3 = """2 * Rplus"""
M3 = """x * Rplus"""

I1 = """dlog( ( Spp + Del ) / ( Spp - Del ) ) / Del"""

Expand Down
41 changes: 41 additions & 0 deletions src/yadism/coefficient_functions/intrinsic/f2_cc.py
Original file line number Diff line number Diff line change
@@ -1 +1,42 @@
# -*- coding: utf-8 -*-

import numpy as np

from eko import constants

from yadism.esf.distribution_vec import rsl_from_distr_coeffs

from .. import partonic_channel as pc

from . import raw_cc

from . import f2_nc


class Splus(f2_nc.Splus):
def __init__(self, ESF, m1sq):
super().__init__(ESF, m1sq, 0.01)


# class Splus(pc.PartonicChannel):
# """
# The convolution point simplifies to :math:`x` when m2=0,
# see :eqref:`6` of :cite:`kretzer-schienbein`.
# """
# def __init__(self, ESF, m1sq):
# super().__init__(ESF)
# self.m1sq = m1sq
# self.y = -ESF.Q2 / m1sq
# self.lo = (1.-1./self.y)

# def LO(self):
# return 0, 0, self.lo

# def NLO(self):
# norm = self.lo * constants.CF
# lnomx = raw_cc.lnomx * norm
# omx = raw_cc.omx(self.y) * norm
# delta = raw_cc.f2sv(self.y) * norm
# def reg(z):
# return (norm * raw_cc.f2r(self.y,z) - (lnomx*np.log(1.-z) + omx)/(1.-z))
# return rsl_from_distr_coeffs(reg, delta, omx, lnomx)
7 changes: 7 additions & 0 deletions src/yadism/coefficient_functions/intrinsic/f3_cc.py
Original file line number Diff line number Diff line change
@@ -1 +1,8 @@
# -*- coding: utf-8 -*-

from . import f3_nc


class Rplus(f3_nc.Rplus):
def __init__(self, ESF, m1sq):
super().__init__(ESF, m1sq, 0.01)
7 changes: 7 additions & 0 deletions src/yadism/coefficient_functions/intrinsic/fl_cc.py
Original file line number Diff line number Diff line change
@@ -1 +1,8 @@
# -*- coding: utf-8 -*-

from . import fl_nc


class Splus(fl_nc.Splus):
def __init__(self, ESF, m1sq):
super().__init__(ESF, m1sq, 0.01)
16 changes: 7 additions & 9 deletions src/yadism/coefficient_functions/intrinsic/kernels.py
Original file line number Diff line number Diff line change
Expand Up @@ -27,15 +27,13 @@ def generate(esf, ihq):
cfs = import_pc_module(kind, esf.process)
m2hq = esf.sf.m2hq[ihq - 4]
if esf.process == "CC":
# w = weights(esf.sf.coupling_constants, esf.Q2, kind, flavors[ihq - 1], ihq)
# wq = {k: v for k, v in w["q"].items() if abs(k) == ihq}
# if kind == "F3":
# return (kernels.Kernel(wq, cfs["intrinsic"]["Rp"](esf, m1sq=m2hq, m2sq=0.0)),)
# return (
# kernels.Kernel(wq, cfs["intrinsic"]["Sp"](esf, m1sq=m2hq, m2sq=0.0)),
# kernels.Kernel(wq, cfs["intrinsic"]["Sm"](esf, m1sq=m2hq, m2sq=0.0)),
# )
return ()
w = kernels.cc_weights(
esf.sf.coupling_constants, esf.Q2, kind, kernels.flavors[ihq - 1], ihq
)
wq = {k: v for k, v in w["ns"].items() if abs(k) == ihq}
if kind == "F3":
return (kernels.Kernel(wq, cfs.Rplus(esf, m1sq=m2hq)),)
return (kernels.Kernel(wq, cfs.Splus(esf, m1sq=m2hq)),)
if kind == "F3":
wVA = esf.sf.coupling_constants.get_weight(ihq, esf.Q2, "VA")
wAV = esf.sf.coupling_constants.get_weight(ihq, esf.Q2, "AV")
Expand Down
167 changes: 167 additions & 0 deletions src/yadism/coefficient_functions/intrinsic/raw_cc.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,167 @@
# -*- coding: utf-8 -*-
import numpy as np
from scipy.special import spence


def li2(x):
return spence(1 - x)


lnomx = -2.0


def omx(y):
return -7.0 / 2.0 + 2.0 * np.log(1.0 - y)


def f2sv(y):
return (
1.0
/ 6.0
* (
-15.0
- 4.0 * np.pi ** 2
+ 9.0 * np.log(1.0 - y)
+ 12.0 * li2(y / (-1.0 + y))
)
)


def f1sv(y):
return f2sv(y) - 1.0 / y * np.log(1.0 - y)


def f3sv(y):
return f1sv(y)


def f1r(y, z):
P1 = 1.0 / (1.0 - z)
P2 = np.log(1.0 - z) * P1
return -(
(
P2
* (
y ** 4
+ y ** 2 * (-13 + (-8 + y) * y) * z ** 2
+ 8 * y * (3 + 5 * y) * z ** 3
- 2 * (5 + 8 * y * (4 + y)) * z ** 4
+ 8 * (3 + 5 * y) * z ** 5
+ (-13 + (-8 + y) * y) * z ** 6
+ z ** 8
)
)
/ ((y - z) * (y - z ** 2) ** 3)
) + (
P1
* (
y ** 4
- 6 * y ** 3 * (2 + y) * z
- 2 * y ** 2 * (-19 + (-13 + y) * y) * z ** 2
+ 2 * y * (-10 + y * (-34 + 7 * y)) * z ** 3
+ 2 * (10 + (34 - 7 * y) * y) * z ** 5
+ 2 * (-19 + (-13 + y) * y) * z ** 6
+ 6 * (2 + y) * z ** 7
- z ** 8
- 2
* (
y ** 4
+ y ** 2 * (-13 + (-8 + y) * y) * z ** 2
+ 8 * y * (3 + 5 * y) * z ** 3
- 2 * (5 + 8 * y * (4 + y)) * z ** 4
+ 8 * (3 + 5 * y) * z ** 5
+ (-13 + (-8 + y) * y) * z ** 6
+ z ** 8
)
* np.log(z / (-y + z))
)
) / (
2.0 * (y - z) * (y - z ** 2) ** 3
)


def f2r(y, z):
P1 = 1.0 / (1.0 - z)
P2 = np.log(1.0 - z) * P1
return -(
(
P2
* (-1 + y)
* z ** 2
* (
-(z ** 6 * (5 + (-4 + z) * z))
+ y ** 5 * (1 + z ** 2)
+ y ** 4 * (-1 + z * (8 - 21 * z + 4 * z ** 3))
+ y ** 2 * z ** 2 * (-5 - (-2 + z) * z * (34 + 7 * z * (-10 + 3 * z)))
+ y ** 3 * z * (4 + z * (-47 + z * (112 - 50 * z + z ** 3)))
+ y * z ** 4 * (-20 + z * (68 + z * (-47 + z * (8 + z))))
)
)
/ ((y - z) * (y - z ** 2) ** 5)
) - (
P1
* (-1 + y)
* z ** 2
* (
-(y ** 4 * (3 + y))
+ 2 * y ** 3 * (5 + y * (22 + y)) * z
+ 2 * y ** 3 * (-54 + y * (-41 + 3 * y)) * z ** 2
+ 6 * y ** 2 * (5 + y * (30 + y)) * z ** 3
- 6 * y * (5 + y * (30 + y)) * z ** 5
+ 2 * y * (54 + (41 - 3 * y) * y) * z ** 6
- 2 * (5 + y * (22 + y)) * z ** 7
+ (3 + y) * z ** 8
+ 2
* (
-(z ** 6 * (5 + (-4 + z) * z))
+ y ** 5 * (1 + z ** 2)
+ y ** 4 * (-1 + z * (8 - 21 * z + 4 * z ** 3))
+ y ** 2 * z ** 2 * (-5 - (-2 + z) * z * (34 + 7 * z * (-10 + 3 * z)))
+ y ** 3 * z * (4 + z * (-47 + z * (112 - 50 * z + z ** 3)))
+ y * z ** 4 * (-20 + z * (68 + z * (-47 + z * (8 + z))))
)
* np.log(z / (-y + z))
)
) / (
2.0 * (y - z) * (y - z ** 2) ** 5
)


def f3r(y, z):
P1 = -1.0 / (1.0 - z)
P2 = -np.log(1.0 - z) * P1
return (
P2
* (-1 + y)
* z
* (
y ** 3
+ y * (-3 + (-7 + y) * y) * z ** 2
+ 16 * y * z ** 3
+ (-3 + (-7 + y) * y) * z ** 4
+ z ** 6
)
) / ((y - z) * (y - z ** 2) ** 3) + (
P1
* (-1 + y)
* z
* (
(y - z ** 2)
* (
z ** 3 * (6 + z)
+ 2 * y * z * (3 + (-11 + z) * z)
+ y ** 2 * (1 + 2 * z + 4 * z ** 2)
)
+ 2
* (
y ** 3
+ y * (-3 + (-7 + y) * y) * z ** 2
+ 16 * y * z ** 3
+ (-3 + (-7 + y) * y) * z ** 4
+ z ** 6
)
* np.log(z / (-y + z))
)
) / (
2.0 * (y - z) * (y - z ** 2) ** 3
)
12 changes: 6 additions & 6 deletions src/yadism/coefficient_functions/intrinsic/raw_nc.py
Original file line number Diff line number Diff line change
Expand Up @@ -64,22 +64,22 @@ def fl_sminus_virt(pc):
return ((np.sqrt(pc.m1sq * pc.m2sq)*((pc.C1m + pc.C1p)*pc.delta**2 + 4*pc.CRm*pc.Q2) + 2*pc.Cplus*(pc.delta**2 - pc.Q2*pc.sigma_pp))*pc.x)/(2.*pc.delta*pc.Q2)

def f3_rplus_raw(pc):
return (2*pc.s1hat*((-2*pc.delta**2*(2*pc.m2sq + pc.s1hat))/pc.s1hat**2 + pc.sigma_mp - 3*pc.sigma_pm - (2*(pc.delta**2 + 2*pc.m2sq*pc.sigma_pm))/pc.s1hat - ((pc.s1hat - pc.sigma_mp)*(pc.s1hat + pc.sigma_pp))/(2.*(pc.m2sq + pc.s1hat))))/(pc.deltap*(pc.m2sq + pc.s1hat)) - (2*pc.L_xi*(pc.s1hat*(pc.s1hat**2 - 4*pc.m1sq*pc.sigma_mp + 3*pc.s1hat*pc.sigma_pm) + 2*pc.delta**2*(2*pc.s1hat + pc.sigma_pp)))/(pc.deltap**2*pc.s1hat)
return (pc.s1hat*((-2*pc.delta**2*(2*pc.m2sq + pc.s1hat))/pc.s1hat**2 + pc.sigma_mp - 3*pc.sigma_pm - (2*(pc.delta**2 + 2*pc.m2sq*pc.sigma_pm))/pc.s1hat - ((pc.s1hat - pc.sigma_mp)*(pc.s1hat + pc.sigma_pp))/(2.*(pc.m2sq + pc.s1hat)))*pc.x)/(pc.deltap*(pc.m2sq + pc.s1hat)) - (pc.L_xi*(pc.s1hat*(pc.s1hat**2 - 4*pc.m1sq*pc.sigma_mp + 3*pc.s1hat*pc.sigma_pm) + 2*pc.delta**2*(2*pc.s1hat + pc.sigma_pp))*pc.x)/(pc.deltap**2*pc.s1hat)

def f3_rplus_soft(pc):
return -8 - (4*pc.L_xisoft*pc.sigma_pp)/pc.delta
return (-2*(2*pc.delta + pc.L_xisoft*pc.sigma_pp)*pc.x)/pc.delta

def f3_rplus_virt(pc):
return 2*pc.CRm
return pc.CRm*pc.x

def f3_rminus_raw(pc):
return (4*np.sqrt(pc.m1sq * pc.m2sq)*(pc.s1hat - pc.sigma_mp))/(pc.deltap*(pc.m2sq + pc.s1hat)) + (4*pc.L_xi*np.sqrt(pc.m1sq * pc.m2sq)*(pc.s1hat + pc.sigma_pm))/pc.deltap**2
return (2*np.sqrt(pc.m1sq * pc.m2sq)*(pc.s1hat - pc.sigma_mp)*pc.x)/(pc.deltap*(pc.m2sq + pc.s1hat)) + (2*pc.L_xi*np.sqrt(pc.m1sq * pc.m2sq)*(pc.s1hat + pc.sigma_pm)*pc.x)/pc.deltap**2

def f3_rminus_soft(pc):
return 0

def f3_rminus_virt(pc):
return 2*pc.Cplus
return pc.Cplus*pc.x

def m1_splus(pc):
return pc.sigma_pp/(2.*pc.delta)
Expand All @@ -100,7 +100,7 @@ def ml_sminus(pc):
return (2*np.sqrt(pc.m1sq * pc.m2sq)*pc.x)/pc.delta

def m3_rplus(pc):
return 2
return pc.x

def m3_rminus(pc):
return 0
Expand Down