Skip to content

Commit

Permalink
implement the fl11 partonic coupling
Browse files Browse the repository at this point in the history
  • Loading branch information
@giacomomagni
giacomomagni committed May 1, 2024
1 parent 74befd5 commit afae238
Show file tree
Hide file tree
Showing 2 changed files with 146 additions and 75 deletions.
116 changes: 72 additions & 44 deletions src/yadism/coefficient_functions/coupling_constants.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,9 +132,7 @@ def nc_partonic_coupling(self, pid):
qZ = {"V": self.vectorial_coupling(pid), "A": self.weak_isospin_3[pid]}
return qph, qZ

def partonic_coupling_fl2(
self, mode, pid, quark_coupling_type, cc_mask=None, nf=None
):
def partonic_coupling(self, mode, pid, quark_coupling_type, cc_mask=None):
"""Compute the coupling of the boson to the parton.
Parameters
Expand Down Expand Up @@ -176,14 +174,22 @@ def partonic_coupling_fl11(self, mode, pid, nf, quark_coupling_type):
"""Compute the coupling of the boson to the parton for the flavor class :math:`fl_{11}`.
This is a generalization of :cite:`Larin:1996wd` (Table 2) for NC.
See also pag. 27. This coupling is given by :math:`Tr(Q_f) Q_f`.
See also pag. 27. This coupling is given by :
.. :math:
W_{q,bb'} = \frac{Tr(Q_b)}{n_f} Q_b'
the gluon and pure singlet coupling are then build after summing on all the
different electroweak channels.
Parameters
----------
mode : str
scattered bosons
pid : int
parton identifier
nf : int
number of active flavors
quark_coupling_type : str
flag to distinguish for heavy quarks between vectorial and axial-vectorial
coupling
Expand All @@ -193,30 +199,28 @@ def partonic_coupling_fl11(self, mode, pid, nf, quark_coupling_type):
float
partonic coupling
"""
# select quark coupling
# In CC there is not such flavor class
if mode == "WW":
return 0

first = quark_coupling_type[0]
second = quark_coupling_type[1]

def switch_mode(quark):
def switch_mode(quark, coupling_type):
qph, qZ = self.nc_partonic_coupling(quark)
if mode == "phph":
return qph[first], qph[second]
if mode == "phZ":
return qph[first], qZ[second]
if mode == "ZZ":
return qZ[first], qZ[second]
return ValueError(f"Unknown mode {mode}")

# load linear coupling
g1, _ = switch_mode(abs(pid))

# compute trace
if mode in ["phph", "Zph"]:
return qph[coupling_type]
if mode in ["phZ", "ZZ"]:
return qZ[coupling_type]
raise ValueError(f"Unknown mode: {mode}")

# first coupling contribute with a mean
pids = range(1, nf + 1)
trace = 0
for quark in pids:
_, g2 = switch_mode(quark)
trace += g2
return g1 * trace
g1 = np.mean([switch_mode(quark, first) for quark in pids])
# the second coupling is just the linear coupling
g2 = switch_mode(abs(pid), second)

return g1 * g2

def propagator_factor(self, mode, Q2):
r"""Compute propagator correction to account for different bosons (:math:`\\eta` in PDG).
Expand Down Expand Up @@ -254,9 +258,12 @@ def propagator_factor(self, mode, Q2):
return eta_W
raise ValueError(f"Unknown mode: {mode}")

def get_weight(self, pid, Q2, quark_coupling_type, cc_mask=None, nf=None):
def get_weight(self, pid, Q2, quark_coupling_type, cc_mask=None):
"""Compute the weight for the pid contributions to the structure function.
Combine the charges, both on the leptonic side and the hadronic side, as well
as propagator changes and/or corrections.
Parameters
----------
pid : int
Expand Down Expand Up @@ -297,10 +304,33 @@ def get_weight(self, pid, Q2, quark_coupling_type, cc_mask=None, nf=None):
if abs(pid) != pos_pid:
return 0.0

partonic_coupling = self.partonic_coupling_fl2
return self.get_nc_weight(partonic_coupling, pid, Q2, quark_coupling_type)
w_phph = (
self.leptonic_coupling("phph", quark_coupling_type)
* self.propagator_factor("phph", Q2)
* self.partonic_coupling("phph", pid, quark_coupling_type)
)
# pure photon exchange
if self.obs_config["process"] == "EM":
return w_phph
# allow Z to be mixed in
if self.obs_config["process"] == "NC":
# photon-Z interference
w_phZ = (
2
* self.leptonic_coupling("phZ", quark_coupling_type)
* self.propagator_factor("phZ", Q2)
* self.partonic_coupling("phZ", pid, quark_coupling_type)
)
# true Z contributions
w_ZZ = (
self.leptonic_coupling("ZZ", quark_coupling_type)
* self.propagator_factor("ZZ", Q2)
* self.partonic_coupling("ZZ", pid, quark_coupling_type)
)
return w_phph + w_phZ + w_ZZ
raise ValueError(f"Unknown process: {self.obs_config['process']}")

def get_fl11_weight(self, pid, Q2, quark_coupling_type, nf):
def get_fl11_weight(self, pid, Q2, nf, quark_coupling_type):
"""Same as :func:`get_weight`but now for the NC flavor class :math:`fl_{11}`.
Combine the charges, both on the leptonic side and the hadronic side,
Expand All @@ -312,11 +342,11 @@ def get_fl11_weight(self, pid, Q2, quark_coupling_type, nf):
particle identifier
Q2 : float
DIS virtuality
nf : int
number of active flavors
quark_coupling_type : str
flag to distinguish for heavy quarks between vectorial and axial-vectorial
coupling
nf : float
number of active flavors
Returns
-------
Expand All @@ -327,43 +357,41 @@ def get_fl11_weight(self, pid, Q2, quark_coupling_type, nf):
if self.obs_config["process"] == "CC":
return 0.0

# if we are are we in positivity mode, we fallback to the standard coupling
# TODO: what does positivity mode means for this coupling ?
if (
self.obs_config["nc_pos_charge"] is not None
and self.obs_config["nc_pos_charge"] != "all"
):
return self.get_weight(pid, Q2, quark_coupling_type)

partonic_coupling = self.partonic_coupling_fl11
return self.get_nc_weight(partonic_coupling, pid, Q2, quark_coupling_type, nf)
return 0.0

def get_nc_weight(self, partonic_coupling, pid, Q2, quark_coupling_type, nf=None):
"""Combine the NC couplings."""
# NC or EM
w_phph = (
self.leptonic_coupling("phph", quark_coupling_type)
* self.propagator_factor("phph", Q2)
* partonic_coupling("phph", pid, quark_coupling_type, nf=nf)
* self.partonic_coupling_fl11("phph", pid, nf, quark_coupling_type)
)
# pure photon exchange
if self.obs_config["process"] == "EM":
return w_phph
# allow Z to be mixed in
if self.obs_config["process"] == "NC":
# photon-Z interference
# photon-Z interference, this class is not symmetric
w_phZ = (
2
* self.leptonic_coupling("phZ", quark_coupling_type)
self.leptonic_coupling("phZ", quark_coupling_type)
* self.propagator_factor("phZ", Q2)
* partonic_coupling("phZ", pid, quark_coupling_type, nf=nf)
* self.partonic_coupling_fl11("phZ", pid, nf, quark_coupling_type)
)
w_Zph = (
self.leptonic_coupling("phZ", quark_coupling_type)
* self.propagator_factor("phZ", Q2)
* self.partonic_coupling_fl11("Zph", pid, nf, quark_coupling_type)
)
# true Z contributions
w_ZZ = (
self.leptonic_coupling("ZZ", quark_coupling_type)
* self.propagator_factor("ZZ", Q2)
* partonic_coupling("ZZ", pid, nf, quark_coupling_type, nf=nf)
* self.partonic_coupling_fl11("ZZ", pid, nf, quark_coupling_type)
)
return w_phph + w_phZ + w_ZZ
return w_phph + w_phZ + w_Zph + w_ZZ
raise ValueError(f"Unknown process: {self.obs_config['process']}")

@classmethod
Expand Down
105 changes: 74 additions & 31 deletions src/yadism/coefficient_functions/light/kernels.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,49 +88,46 @@ def generate(esf, nf):
kernels_list = [ns, g, s]

# at N3LO we need to add also the fl11 diagrams
if 3 in esf.orders:
weights_fl11 = nc_weights(
esf.info.coupling_constants, esf.Q2, nf, is_pv=False, is_fl11=True
)
if esf.info.configs.theory["pto"] == 3:

gluon_fl11 = pcs.GluonFL11(esf, nf)
quark_fl11 = pcs.QuarkFL11(esf, nf)
weights_fl11 = nc_fl11_weights(esf.info.coupling_constants, esf.Q2, nf)
ns_fl11 = kernels.Kernel(
weights_fl11["ns"],
pcs.NonSingletFL11(esf, nf),
)
g_fl11 = kernels.Kernel(
weights_fl11["g"],
pcs.GluonFL11(esf, nf),
quark_fl11,
)
g_fl11 = kernels.Kernel(weights_fl11["g"], gluon_fl11)
s_fl11 = kernels.Kernel(
weights_fl11["s"],
pcs.SingletFL11(esf, nf),
quark_fl11,
)
kernels_list.extend([ns_fl11, g_fl11, s_fl11])
return kernels_list


def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False, is_fl11=False):
def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False):
"""
Compute light NC weights.
Parameters
----------
coupling_constants : CouplingConstants
manager for coupling constants
Q2 : float
W virtuality
nf : int
number of light flavors
is_pv: bool
True if observable violates parity conservation
skip_heavylight : bool
prevent the last quark to couple to the boson
is_fl11 : bool
if True compute the coupling of the flavor class :math:`fl_{11}`
coupling_constants : CouplingConstants
manager for coupling constants
Q2 : float
W virtuality
nf : int
number of light flavors
is_pv: bool
True if observable violates parity conservation
skip_heavylight : bool
prevent the last quark to couple to the boson
Returns
-------
weights : dict
mapping pid -> weight for ns, g and s channel
weights : dict
mapping pid -> weight for ns, g and s channel
"""
# quark couplings
tot_ch_sq = 0
Expand All @@ -147,12 +144,9 @@ def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False, is_fl11
q, Q2, "VA"
) + coupling_constants.get_weight(q, Q2, "AV")
else:
coupling = coupling_constants.get_weight

# TODO: do AA contribute to fl11 ??
if is_fl11:
coupling = coupling_constants.get_fl11_weight
w = coupling(q, Q2, "VV", nf=nf) + coupling(q, Q2, "AA", nf=nf)
w = coupling_constants.get_weight(
q, Q2, "VV", nf=nf
) + coupling_constants.get_weight(q, Q2, "AA", nf=nf)
ns_partons[q] = w
ns_partons[-q] = w if not is_pv else -w
tot_ch_sq += w
Expand All @@ -166,3 +160,52 @@ def nc_weights(coupling_constants, Q2, nf, is_pv, skip_heavylight=False, is_fl11
# gluon and singlet coupling = charge average (omitting the *2/2)
s_partons = {q: ch_av for q in [*pids, *(-q for q in pids)]}
return {"ns": ns_partons, "g": {21: ch_av}, "s": s_partons}


def nc_fl11_weights(coupling_constants, Q2, nf, skip_heavylight=False):
"""Compute the NC weights for the flavor class :math:`fl_{11}`.
For the time being we don't have such diagrams for parity violating
structure functions.
Parameters
----------
coupling_constants : CouplingConstants
manager for coupling constants
Q2 : float
W virtuality
nf : int
number of light flavors
skip_heavylight : bool
prevent the last quark to couple to the boson
Returns
-------
weights : dict
mapping pid -> weight for ns, g and s channel
"""
# quark couplings
ns_partons = {}
tot_ch_sq = 0
pids = range(1, nf + 1)

for q in pids:
# we don't want this quark to couple to the photon (because it is treated separately),
# but still let it take part in the average
if skip_heavylight and q == nf:
continue
w = coupling_constants.get_fl11_weight(
q, Q2, nf, "VV"
) + coupling_constants.get_fl11_weight(q, Q2, nf, "AA")
ns_partons[q] = w
ns_partons[-q] = w
tot_ch_sq += w

# compute gluon, singlet (S) (omitting the *2/2)
# NOTE: what we call singlet (S) is actually the pure singlet (PS) coefficient.
# Here the contribution of the singlet PDF to the NS coefficient
# and the contribution of the singlet PDF to the S coefficient
# do not cancel completely.
ch_av = tot_ch_sq / len(pids)
s_partons = {q: ch_av - ns_partons[q] for q in [*pids, *(-q for q in pids)]}
return {"ns": ns_partons, "g": {21: ch_av}, "s": s_partons}

0 comments on commit afae238

Please sign in to comment.