"""Collection of quantum number conservation rules for particle reactions.
This module is the place where the 'expert' defines the rules that verify quantum
numbers of the reaction.
A rule is a function that takes quantum numbers as input and outputs a boolean. There
are three different types of rules:
1. `GraphElementRule` that work on individual graph edges or nodes.
2. `EdgeQNConservationRule` that work on the interaction level, which use ingoing edges,
outgoing edges as arguments. E.g.: `.ChargeConservation`.
3. `ConservationRule` that work on the interaction level, which use ingoing edges,
outgoing edges and a interaction node as arguments. E.g: `.parity_conservation`.
The arguments can be any type of quantum number. However a rule argument resembling
edges only accepts `~.quantum_numbers.EdgeQuantumNumbers`. Similarly arguments that
resemble a node only accept `~.quantum_numbers.NodeQuantumNumbers`. The argument types
do not have to be limited to a single quantum number, but can be a composite (see
`.CParityEdgeInput`).
.. warning::
Besides the rule logic itself, a rule also has the responsibility of
stating its run conditions. These run conditions **must** be stated by
the type annotations of its :code:`__call__` method. The type annotations
therefore are not just there for *static* type checking: they also
carry more information about the rule that is extracted *dynamically*
by the :mod:`.solving` module.
Generally, the conditions can be separated into two categories:
* variable conditions
* toplogical conditions
Currently, only variable conditions are being used. Topological conditions could be
created in the form of `~typing.Tuple` instead of `~typing.List`.
For additive quantum numbers, the decorator `additive_quantum_number_rule` can be used
to automatically generate the appropriate behavior.
The module is therefore strongly typed (both for the reader of the code and for type
checking with :doc:`mypy <mypy:index>`). An example is `.HelicityParityEdgeInput`, which
has been defined to provide type checks on `.parity_conservation_helicity`.
.. seealso:: :doc:`/usage/conservation`
"""
import sys
from copy import deepcopy
from functools import reduce
from textwrap import dedent
from typing import Any, Callable, List, Optional, Set, Tuple, Type, Union
from attrs import define, field, frozen
from attrs.converters import optional
from .quantum_numbers import EdgeQuantumNumbers as EdgeQN
from .quantum_numbers import NodeQuantumNumbers as NodeQN
from .quantum_numbers import arange
if sys.version_info >= (3, 8):
from typing import Protocol
else:
from typing_extensions import Protocol
def _is_boson(spin_magnitude: float) -> bool:
return abs(spin_magnitude % 1) < 0.01
def _is_particle_antiparticle_pair(pid1: int, pid2: int) -> bool:
# we just check if the pid is opposite in sign
# this is a requirement of the pid numbers of course
return pid1 == -pid2
[docs]class GraphElementRule(Protocol):
[docs] def __call__(self, __qns: Any) -> bool:
...
[docs]class EdgeQNConservationRule(Protocol):
[docs] def __call__(
self, __ingoing_edge_qns: List[Any], __outgoing_edge_qns: List[Any]
) -> bool:
...
[docs]class ConservationRule(Protocol):
[docs] def __call__(
self,
__ingoing_edge_qns: List[Any],
__outgoing_edge_qns: List[Any],
__node_qns: Any,
) -> bool:
...
# Note a generic would be more fitting here. However the type annotations of
# __call__ method in a concrete version of the generic are still containing the
# TypeVar types. See https://github.com/python/typing/issues/762
[docs]def additive_quantum_number_rule(
quantum_number: type,
) -> Callable[[Any], EdgeQNConservationRule]:
r"""Class decorator for creating an additive conservation rule.
Use this decorator to create a `EdgeQNConservationRule` for a quantum number to
which an additive conservation rule applies:
.. math:: \sum q_{in} = \sum q_{out}
Args:
quantum_number: Quantum number to which you want to apply the additive
conservation check. An example would be `.EdgeQuantumNumbers.charge`.
"""
def decorator(rule_class: Any) -> EdgeQNConservationRule:
def new_call(
# pylint: disable=unused-argument
self: Type[EdgeQNConservationRule],
ingoing_edge_qns: List[quantum_number], # type: ignore[valid-type]
outgoing_edge_qns: List[quantum_number], # type: ignore[valid-type]
) -> bool:
return sum(ingoing_edge_qns) == sum(outgoing_edge_qns)
rule_class.__call__ = new_call
rule_class.__doc__ = dedent(f"""
Decorated via `{additive_quantum_number_rule.__name__}`.
Check for `~.EdgeQuantumNumbers.{quantum_number.__name__}` conservation.
""")
return rule_class
return decorator
[docs]@additive_quantum_number_rule(EdgeQN.charge)
class ChargeConservation(EdgeQNConservationRule):
pass
[docs]@additive_quantum_number_rule(EdgeQN.baryon_number)
class BaryonNumberConservation(EdgeQNConservationRule):
pass
[docs]@additive_quantum_number_rule(EdgeQN.electron_lepton_number)
class ElectronLNConservation(EdgeQNConservationRule):
pass
[docs]@additive_quantum_number_rule(EdgeQN.muon_lepton_number)
class MuonLNConservation(EdgeQNConservationRule):
pass
[docs]@additive_quantum_number_rule(EdgeQN.tau_lepton_number)
class TauLNConservation(EdgeQNConservationRule):
pass
[docs]@additive_quantum_number_rule(EdgeQN.strangeness)
class StrangenessConservation(EdgeQNConservationRule):
pass
[docs]@additive_quantum_number_rule(EdgeQN.charmness)
class CharmConservation(EdgeQNConservationRule):
pass
[docs]@additive_quantum_number_rule(EdgeQN.bottomness)
class BottomnessConservation(EdgeQNConservationRule):
pass
[docs]def parity_conservation(
ingoing_edge_qns: List[EdgeQN.parity],
outgoing_edge_qns: List[EdgeQN.parity],
l_magnitude: NodeQN.l_magnitude,
) -> bool:
r"""Implement :math:`P_{in} = P_{out} \cdot (-1)^L`."""
if any(p is None for p in [*ingoing_edge_qns, *outgoing_edge_qns]):
return False
if len(ingoing_edge_qns) == 1 and len(outgoing_edge_qns) == 2:
parity_in = reduce(lambda x, y: x * y.value, ingoing_edge_qns, 1)
parity_out = reduce(lambda x, y: x * y.value, outgoing_edge_qns, 1)
return parity_in == (parity_out * (-1) ** l_magnitude)
return True
[docs]def parity_conservation_helicity(
ingoing_edge_qns: List[HelicityParityEdgeInput],
outgoing_edge_qns: List[HelicityParityEdgeInput],
parity_prefactor: NodeQN.parity_prefactor,
) -> bool:
r"""Implements parity conservation for helicity formalism.
Check the following:
.. math:: A_{-\lambda_1-\lambda_2} = P_1 P_2 P_3 (-1)^{S_2+S_3-S_1}
A_{\lambda_1\lambda_2}
.. math:: \mathrm{parity\,prefactor} = P_1 P_2 P_3 (-1)^{S_2+S_3-S_1}
.. note:: Only the special case :math:`\lambda_1=\lambda_2=0` may return `False`
independent on the parity prefactor.
"""
if len(ingoing_edge_qns) == 1 and len(outgoing_edge_qns) == 2:
out_spins = [x.spin_magnitude for x in outgoing_edge_qns]
parity_product = reduce(
lambda x, y: x * y.parity.value if y.parity else x,
ingoing_edge_qns + outgoing_edge_qns,
1,
)
prefactor = parity_product * (-1.0) ** (
sum(out_spins) - ingoing_edge_qns[0].spin_magnitude
)
if all(x.spin_projection == 0.0 for x in outgoing_edge_qns) and prefactor == -1:
return False
return prefactor == parity_prefactor
return True
[docs]def c_parity_conservation(
ingoing_edge_qns: List[CParityEdgeInput],
outgoing_edge_qns: List[CParityEdgeInput],
interaction_node_qns: CParityNodeInput,
) -> bool:
"""Check for :math:`C`-parity conservation.
Implements :math:`C_{in} = C_{out}`.
"""
def _get_c_parity_multiparticle(
part_qns: List[CParityEdgeInput], interaction_qns: CParityNodeInput
) -> Optional[int]:
c_parities_part = [x.c_parity.value for x in part_qns if x.c_parity]
# if all states have C parity defined, then just multiply them
if len(c_parities_part) == len(part_qns):
return reduce(lambda x, y: x * y, c_parities_part, 1)
# two particle case
if len(part_qns) == 2:
if _is_particle_antiparticle_pair(part_qns[0].pid, part_qns[1].pid):
ang_mom = interaction_qns.l_magnitude
# if boson
if _is_boson(part_qns[0].spin_magnitude):
return (-1) ** int(ang_mom)
coupled_spin = interaction_qns.s_magnitude
if isinstance(coupled_spin, int) or coupled_spin.is_integer():
return (-1) ** int(ang_mom + coupled_spin)
return None
c_parity_in = _get_c_parity_multiparticle(ingoing_edge_qns, interaction_node_qns)
if c_parity_in is None:
return True
c_parity_out = _get_c_parity_multiparticle(outgoing_edge_qns, interaction_node_qns)
if c_parity_out is None:
return True
return c_parity_in == c_parity_out
[docs]def g_parity_conservation(
ingoing_edge_qns: List[GParityEdgeInput],
outgoing_edge_qns: List[GParityEdgeInput],
interaction_qns: GParityNodeInput,
) -> bool:
"""Check for :math:`G`-parity conservation.
Implements for :math:`G_{in} = G_{out}`.
"""
def check_multistate_g_parity(
isospin: EdgeQN.isospin_magnitude,
double_state_qns: Tuple[GParityEdgeInput, GParityEdgeInput],
) -> Optional[int]:
if _is_particle_antiparticle_pair(
double_state_qns[0].pid, double_state_qns[1].pid
):
ang_mom = interaction_qns.l_magnitude
if isinstance(isospin, int) or isospin.is_integer():
# if boson
if _is_boson(double_state_qns[0].spin_magnitude):
return (-1) ** int(ang_mom + isospin)
coupled_spin = interaction_qns.s_magnitude
if isinstance(coupled_spin, int) or coupled_spin.is_integer():
return (-1) ** int(ang_mom + coupled_spin + isospin)
return None
def check_g_parity_isobar(
single_state: GParityEdgeInput,
couple_state: Tuple[GParityEdgeInput, GParityEdgeInput],
) -> bool:
couple_state_g_parity = check_multistate_g_parity(
single_state.isospin_magnitude,
(couple_state[0], couple_state[1]),
)
single_state_g_parity = (
ingoing_edge_qns[0].g_parity.value if ingoing_edge_qns[0].g_parity else None
)
if not couple_state_g_parity or not single_state_g_parity:
return True
return couple_state_g_parity == single_state_g_parity
no_g_parity_in_part = [True for x in ingoing_edge_qns if x.g_parity is None]
no_g_parity_out_part = [True for x in outgoing_edge_qns if x.g_parity is None]
# if all states have G parity defined, then just multiply them
if not any(no_g_parity_in_part + no_g_parity_out_part):
in_g_parity = reduce(
lambda x, y: x * y.g_parity.value if y.g_parity else x,
ingoing_edge_qns,
1,
)
out_g_parity = reduce(
lambda x, y: x * y.g_parity.value if y.g_parity else x,
outgoing_edge_qns,
1,
)
return in_g_parity == out_g_parity
# two particle case
particle_counts = (len(ingoing_edge_qns), len(outgoing_edge_qns))
if particle_counts == (1, 2):
return check_g_parity_isobar(
ingoing_edge_qns[0],
(outgoing_edge_qns[0], outgoing_edge_qns[1]),
)
if particle_counts == (2, 1):
return check_g_parity_isobar(
outgoing_edge_qns[0],
(ingoing_edge_qns[0], ingoing_edge_qns[1]),
)
return True
[docs]@frozen
class IdenticalParticleSymmetryOutEdgeInput:
spin_magnitude: EdgeQN.spin_magnitude = field(converter=EdgeQN.spin_magnitude)
spin_projection: EdgeQN.spin_projection = field(converter=EdgeQN.spin_projection)
pid: EdgeQN.pid = field(converter=EdgeQN.pid)
[docs]def identical_particle_symmetrization(
ingoing_parities: List[EdgeQN.parity],
outgoing_edge_qns: List[IdenticalParticleSymmetryOutEdgeInput],
) -> bool:
"""Verifies multi particle state symmetrization for identical particles.
In case of a multi particle state with identical particles, their exchange symmetry
has to follow the spin statistic theorem.
For bosonic systems the total exchange symmetry (parity) has to be even (+1). For
fermionic systems the total exchange symmetry (parity) has to be odd (-1).
In case of a particle decaying into N identical particles (N>1), the decaying
particle has to have the same parity as required by the spin statistic theorem of
the multi body state.
"""
def _check_particles_identical(
particles: List[IdenticalParticleSymmetryOutEdgeInput],
) -> bool:
"""Check if pids and spins match."""
if len(particles) == 1:
return False
reference_pid = particles[0].pid
reference_spin_proj = particles[0].spin_projection
for particle in particles[1:]:
if particle.pid != reference_pid:
return False
if particle.spin_projection != reference_spin_proj:
return False
return True
if len(ingoing_parities) == 1:
if _check_particles_identical(outgoing_edge_qns):
if _is_boson(outgoing_edge_qns[0].spin_magnitude):
# we have a boson, check if parity of mother is even
parity = ingoing_parities[0]
if parity == -1:
# if its odd then return False
return False
else:
# its fermion
parity = ingoing_parities[0]
if parity == 1:
return False
return True
@frozen
class _Spin:
magnitude: float
projection: float
def _is_clebsch_gordan_coefficient_zero(
spin1: _Spin, spin2: _Spin, spin_coupled: _Spin
) -> bool:
m_1 = spin1.projection
j_1 = spin1.magnitude
m_2 = spin2.projection
j_2 = spin2.magnitude
proj = spin_coupled.projection
mag = spin_coupled.magnitude
if (j_1 == j_2 and m_1 == m_2) or (m_1 == 0.0 and m_2 == 0.0):
if abs(mag - j_1 - j_2) % 2 == 1:
return True
if j_1 == mag and m_1 == -proj:
if abs(j_2 - j_1 - mag) % 2 == 1:
return True
if j_2 == mag and m_2 == -proj:
if abs(j_1 - j_2 - mag) % 2 == 1:
return True
return False
[docs]def ls_spin_validity(spin_input: SpinNodeInput) -> bool:
r"""Check for valid isospin magnitude and projection."""
return _check_spin_valid(
float(spin_input.l_magnitude), float(spin_input.l_projection)
) and _check_spin_valid(
float(spin_input.s_magnitude), float(spin_input.s_projection)
)
def _check_magnitude(
in_part: List[float],
out_part: List[float],
interaction_qns: Optional[Union[SpinMagnitudeNodeInput, SpinNodeInput]],
) -> bool:
def couple_mags(j_1: float, j_2: float) -> List[float]:
return [
x / 2.0 for x in range(int(2 * abs(j_1 - j_2)), int(2 * (j_1 + j_2 + 1)), 2)
]
def couple_magnitudes(
magnitudes: List[float],
interaction_qns: Optional[Union[SpinMagnitudeNodeInput, SpinNodeInput]],
) -> Set[float]:
if len(magnitudes) == 1:
return set(magnitudes)
coupled_magnitudes = {magnitudes[0]}
for mag in magnitudes[1:]:
temp_set = coupled_magnitudes
coupled_magnitudes = set()
for ref_mag in temp_set:
coupled_magnitudes.update(couple_mags(mag, ref_mag))
if interaction_qns:
if interaction_qns.s_magnitude in coupled_magnitudes:
return set(
couple_mags(
interaction_qns.s_magnitude,
interaction_qns.l_magnitude,
)
)
return set() # in case there the spin coupling fails
return coupled_magnitudes
in_tot_spins = couple_magnitudes(in_part, interaction_qns)
out_tot_spins = couple_magnitudes(out_part, interaction_qns)
matching_spins = in_tot_spins.intersection(out_tot_spins)
if len(matching_spins) > 0:
return True
return False
def _check_spin_couplings(
in_part: List[_Spin],
out_part: List[_Spin],
interaction_qns: Optional[SpinNodeInput],
) -> bool:
in_tot_spins = __calculate_total_spins(in_part, interaction_qns)
out_tot_spins = __calculate_total_spins(out_part, interaction_qns)
matching_spins = in_tot_spins & out_tot_spins
if len(matching_spins) > 0:
return True
return False
def __calculate_total_spins(
spins: List[_Spin],
interaction_qns: Optional[SpinNodeInput] = None,
) -> Set[_Spin]:
total_spins = set()
if len(spins) == 1:
return set(spins)
total_spins = __create_coupled_spins(spins)
if interaction_qns:
coupled_spin = _Spin(interaction_qns.s_magnitude, interaction_qns.s_projection)
if coupled_spin in total_spins:
return __spin_couplings(
coupled_spin,
_Spin(interaction_qns.l_magnitude, interaction_qns.l_projection),
)
total_spins = set()
return total_spins
def __create_coupled_spins(spins: List[_Spin]) -> Set[_Spin]:
"""Creates all combinations of coupled spins."""
spins_daughters_coupled: Set[_Spin] = set()
spin_list = deepcopy(spins)
while spin_list:
if spins_daughters_coupled:
temp_coupled_spins = set()
tempspin = spin_list.pop()
for spin in spins_daughters_coupled:
coupled_spins = __spin_couplings(spin, tempspin)
temp_coupled_spins.update(coupled_spins)
spins_daughters_coupled = temp_coupled_spins
else:
spins_daughters_coupled.add(spin_list.pop())
return spins_daughters_coupled
def __spin_couplings(spin1: _Spin, spin2: _Spin) -> Set[_Spin]:
r"""Implement the coupling of two spins.
:math:`|S_1 - S_2| \leq S \leq |S_1 + S_2|` and :math:`M_1 + M_2 = M`
"""
s_1 = spin1.magnitude
s_2 = spin2.magnitude
sum_proj = spin1.projection + spin2.projection
return {
_Spin(x, sum_proj)
for x in arange(abs(s_1 - s_2), s_1 + s_2 + 1, 1.0)
if x >= abs(sum_proj)
and not _is_clebsch_gordan_coefficient_zero(spin1, spin2, _Spin(x, sum_proj))
}
def _check_spin_valid(magnitude: float, projection: float) -> bool:
if magnitude % 0.5 != 0.0:
return False
if abs(projection) > magnitude:
return False
return float(projection - magnitude).is_integer()
[docs]def isospin_validity(isospin: IsoSpinEdgeInput) -> bool:
r"""Check for valid isospin magnitude and projection."""
return _check_spin_valid(
float(isospin.isospin_magnitude), float(isospin.isospin_projection)
)
[docs]def isospin_conservation(
ingoing_isospins: List[IsoSpinEdgeInput],
outgoing_isospins: List[IsoSpinEdgeInput],
) -> bool:
r"""Check for isospin conservation.
Implements
.. math::
|I_1 - I_2| \leq I \leq |I_1 + I_2|
Also checks :math:`I_{1,z} + I_{2,z} = I_z` and if Clebsch-Gordan coefficients are
all 0.
"""
if not sum(x.isospin_projection for x in ingoing_isospins) == sum(
x.isospin_projection for x in outgoing_isospins
):
return False
if not all(isospin_validity(x) for x in ingoing_isospins + outgoing_isospins):
return False
return _check_spin_couplings(
[_Spin(x.isospin_magnitude, x.isospin_projection) for x in ingoing_isospins],
[_Spin(x.isospin_magnitude, x.isospin_projection) for x in outgoing_isospins],
None,
)
[docs]def spin_validity(spin: SpinEdgeInput) -> bool:
r"""Check for valid spin magnitude and projection."""
return _check_spin_valid(float(spin.spin_magnitude), float(spin.spin_projection))
[docs]def spin_conservation(
ingoing_spins: List[SpinEdgeInput],
outgoing_spins: List[SpinEdgeInput],
interaction_qns: SpinNodeInput,
) -> bool:
r"""Check for spin conservation.
Implements
.. math::
|S_1 - S_2| \leq S \leq |S_1 + S_2|
and
.. math::
|L - S| \leq J \leq |L + S|
Also checks :math:`M_1 + M_2 = M` and if Clebsch-Gordan coefficients are all 0.
.. seealso:: /docs/usage/ls-coupling
"""
# L and S can only be used if one side is a single state
# and the other side contains of two states (isobar)
# So do a full check if this is the case
if (len(ingoing_spins) == 1 and len(outgoing_spins) == 2) or (
len(ingoing_spins) == 2 and len(outgoing_spins) == 1
):
return _check_spin_couplings(
[_Spin(x.spin_magnitude, x.spin_projection) for x in ingoing_spins],
[_Spin(x.spin_magnitude, x.spin_projection) for x in outgoing_spins],
interaction_qns,
)
# otherwise don't use S and L and just check magnitude
# are integral or non integral on both sides
return (
sum(float(x.spin_magnitude) for x in ingoing_spins).is_integer() # type: ignore[union-attr]
== sum(float(x.spin_magnitude) for x in outgoing_spins).is_integer() # type: ignore[union-attr]
)
[docs]def spin_magnitude_conservation(
ingoing_spins: List[SpinEdgeInput],
outgoing_spins: List[SpinEdgeInput],
interaction_qns: SpinMagnitudeNodeInput,
) -> bool:
r"""Check for spin conservation.
Implements
.. math::
|S_1 - S_2| \leq S \leq |S_1 + S_2|
and
.. math::
|L - S| \leq J \leq |L + S|
"""
# L and S can only be used if one side is a single state
# and the other side contains of two states (isobar)
# So do a full check if this is the case
if (len(ingoing_spins) == 1 and len(outgoing_spins) == 2) or (
len(ingoing_spins) == 2 and len(outgoing_spins) == 1
):
return _check_magnitude(
[x.spin_magnitude for x in ingoing_spins],
[x.spin_magnitude for x in outgoing_spins],
interaction_qns,
)
# otherwise don't use S and L and just check magnitude
# are integral or non integral on both sides
return (
sum(float(x.spin_magnitude) for x in ingoing_spins).is_integer() # type: ignore[union-attr]
== sum(float(x.spin_magnitude) for x in outgoing_spins).is_integer() # type: ignore[union-attr]
)
[docs]def clebsch_gordan_helicity_to_canonical(
ingoing_spins: List[SpinEdgeInput],
outgoing_spins: List[SpinEdgeInput],
interaction_qns: SpinNodeInput,
) -> bool:
"""Implement Clebsch-Gordan checks.
For :math:`S_1, S_2` to :math:`S` and the :math:`L,S` to :math:`J` coupling based on
the conversion of helicity to canonical amplitude sums.
.. note:: This rule does not check that the spin magnitudes couple correctly to
:math:`L` and :math:`S`, as this is already performed by
`.spin_magnitude_conservation`.
"""
if len(ingoing_spins) == 1 and len(outgoing_spins) == 2:
out_spin1 = _Spin(
outgoing_spins[0].spin_magnitude,
outgoing_spins[0].spin_projection,
)
out_spin2 = _Spin(
outgoing_spins[1].spin_magnitude,
-outgoing_spins[1].spin_projection,
)
helicity_diff = out_spin1.projection + out_spin2.projection
if helicity_diff != interaction_qns.s_projection:
return False
ang_mom = _Spin(interaction_qns.l_magnitude, interaction_qns.l_projection)
coupled_spin = _Spin(interaction_qns.s_magnitude, interaction_qns.s_projection)
parent_spin = ingoing_spins[0].spin_magnitude
coupled_spin = _Spin(coupled_spin.magnitude, helicity_diff)
if not _check_spin_valid(coupled_spin.magnitude, coupled_spin.projection):
return False
in_spin = _Spin(parent_spin, helicity_diff)
if not _check_spin_valid(in_spin.magnitude, in_spin.projection):
return False
if _is_clebsch_gordan_coefficient_zero(out_spin1, out_spin2, coupled_spin):
return False
return not _is_clebsch_gordan_coefficient_zero(ang_mom, coupled_spin, in_spin)
return False
[docs]def helicity_conservation(
ingoing_spin_mags: List[EdgeQN.spin_magnitude],
outgoing_helicities: List[EdgeQN.spin_projection],
) -> bool:
r"""Implementation of helicity conservation.
Check for :math:`|\lambda_2-\lambda_3| \leq S_1`.
"""
if len(ingoing_spin_mags) == 1 and len(outgoing_helicities) == 2:
mother_spin = ingoing_spin_mags[0]
if mother_spin >= abs(outgoing_helicities[0] - outgoing_helicities[1]):
return True
return False
[docs]@frozen
class GellMannNishijimaInput:
# pylint: disable=too-many-instance-attributes
charge: EdgeQN.charge = field(converter=EdgeQN.charge)
isospin_projection: Optional[EdgeQN.isospin_projection] = field(
converter=optional(EdgeQN.isospin_projection), default=None
)
strangeness: Optional[EdgeQN.strangeness] = field(
converter=optional(EdgeQN.strangeness), default=None
)
charmness: Optional[EdgeQN.charmness] = field(
converter=optional(EdgeQN.charmness), default=None
)
bottomness: Optional[EdgeQN.bottomness] = field(
converter=optional(EdgeQN.bottomness), default=None
)
topness: Optional[EdgeQN.topness] = field(
converter=optional(EdgeQN.topness), default=None
)
baryon_number: Optional[EdgeQN.baryon_number] = field(
converter=optional(EdgeQN.baryon_number), default=None
)
electron_lepton_number: Optional[EdgeQN.electron_lepton_number] = field(
converter=optional(EdgeQN.electron_lepton_number), default=None
)
muon_lepton_number: Optional[EdgeQN.muon_lepton_number] = field(
converter=optional(EdgeQN.muon_lepton_number), default=None
)
tau_lepton_number: Optional[EdgeQN.tau_lepton_number] = field(
converter=optional(EdgeQN.tau_lepton_number), default=None
)
[docs]def gellmann_nishijima(edge_qns: GellMannNishijimaInput) -> bool:
r"""Check the Gell-Mann–Nishijima formula.
`Gell-Mann–Nishijima formula
<https://en.wikipedia.org/wiki/Gell-Mann%E2%80%93Nishijima_formula>`_:
.. math::
Q = I_3 + \frac{1}{2}(B+S+C+B'+T)
where
:math:`Q` is charge (computed),
:math:`I_3` is `.Spin.projection` of `~.Particle.isospin`,
:math:`B` is `~.Particle.baryon_number`,
:math:`S` is `~.Particle.strangeness`,
:math:`C` is `~.Particle.charmness`,
:math:`B'` is `~.Particle.bottomness`, and
:math:`T` is `~.Particle.topness`.
"""
def calculate_hypercharge(
particle: GellMannNishijimaInput,
) -> float:
"""Calculate the hypercharge :math:`Y=S+C+B+T+B`."""
return sum(
0.0 if x is None else x
for x in [
particle.strangeness,
particle.charmness,
particle.bottomness,
particle.topness,
particle.baryon_number,
]
)
if (
edge_qns.electron_lepton_number
or edge_qns.muon_lepton_number
or edge_qns.tau_lepton_number
):
return True
isospin_3 = 0.0
if edge_qns.isospin_projection:
isospin_3 = edge_qns.isospin_projection
if float(edge_qns.charge) != (isospin_3 + 0.5 * calculate_hypercharge(edge_qns)):
return False
return True
[docs]class MassConservation:
"""Mass conservation rule."""
def __init__(self, width_factor: float):
self.__width_factor = width_factor
[docs] def __call__(
self,
ingoing_edge_qns: List[MassEdgeInput],
outgoing_edge_qns: List[MassEdgeInput],
) -> bool:
r"""Implements mass conservation.
:math:`M_{out} - N \cdot W_{out} < M_{in} + N \cdot W_{in}`
It makes sure that the net mass outgoing state :math:`M_{out}` is smaller than
the net mass of the ingoing state :math:`M_{in}`. Also the width :math:`W` of
the states is taken into account.
"""
mass_in = sum(x.mass for x in ingoing_edge_qns)
width_in = sum(x.width for x in ingoing_edge_qns if x.width)
mass_out = sum(x.mass for x in outgoing_edge_qns)
width_out = sum(x.width for x in outgoing_edge_qns if x.width)
return (mass_out - self.__width_factor * width_out) < (
mass_in + self.__width_factor * width_in
)
