"""Collection of quantum number conservation rules for particle reactions.
Contains:
- Functors for quantum number condition checks.
"""
import logging
from abc import ABC, abstractmethod
from copy import deepcopy
from functools import reduce
from numpy import arange
from .particle import (
InteractionQuantumNumberNames,
ParticleDecayPropertyNames,
ParticlePropertyNames,
QNNameClassMapping,
QuantumNumberClasses,
Spin,
StateQuantumNumberNames,
is_boson,
)
[docs]class AbstractConditionFunctor(ABC):
"""Abstract interface of a condition functor."""
[docs] @abstractmethod
def check(self, qn_names, in_edges, out_edges, int_node):
pass
[docs]class DefinedForAllEdges(AbstractConditionFunctor):
"""Check if a graph has all edges defined."""
[docs] def check(self, qn_names, in_edges, out_edges, int_node):
for qn_name in qn_names:
for edge in in_edges + out_edges:
if qn_name not in edge:
return False
return True
[docs]class DefinedForAllOutgoingEdges(AbstractConditionFunctor):
"""Check if all outgoing edges are defined."""
[docs] def check(self, qn_names, in_edges, out_edges, int_node):
for qn_name in qn_names:
for edge in out_edges:
if qn_name not in edge:
return False
return True
[docs]class DefinedForInteractionNode(AbstractConditionFunctor):
"""Check if all interaction nodes are defined."""
[docs] def check(self, qn_names, in_edges, out_edges, int_node):
for qn_name in qn_names:
if qn_name not in int_node:
return False
return True
[docs]class DefinedIfOtherQnNotDefinedInOutSeparate(AbstractConditionFunctor):
"""Implements logic for..."""
def __init__(self, other_qn_names):
self.other_qn_names = other_qn_names
[docs] def check(self, qn_names, in_edges, out_edges, int_node):
return self.check_edge_set(
qn_names, in_edges, int_node
) and self.check_edge_set(qn_names, out_edges, int_node)
[docs] def check_edge_set(self, qn_names, edges, int_node):
found_for_all = True
for qn_name in self.other_qn_names:
if isinstance(
qn_name, (StateQuantumNumberNames, ParticlePropertyNames)
):
found_for_all = True
for edge_props in edges:
if not self.find_in_dict(qn_name, edge_props):
found_for_all = False
break
if not found_for_all:
break
else:
if not self.find_in_dict(qn_name, int_node):
found_for_all = False
break
if not found_for_all:
for qn_name in qn_names:
if isinstance(
qn_name, (StateQuantumNumberNames, ParticlePropertyNames)
):
for edge_props in edges:
if not self.find_in_dict(qn_name, edge_props):
return False
else:
if not self.find_in_dict(qn_name, int_node):
return False
return True
[docs] @staticmethod
def find_in_dict(name, props):
found = False
for key, val in props.items():
if name == key and val is not None:
found = True
break
return found
[docs]def is_particle_antiparticle_pair(pid1, pid2):
# 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 AbstractRule(ABC):
"""Abstract interface for a conservation rule."""
def __init__(self, rule_name):
self.rule_name = str(rule_name)
self.required_qn_names = []
self.qn_conditions = []
self.specify_required_qns()
def __repr__(self):
return str(self.rule_name)
def __str__(self):
return str(self.rule_name)
[docs] def get_qn_conditions(self):
return self.qn_conditions
[docs] @abstractmethod
def specify_required_qns(self):
pass
[docs] def add_required_qn(self, qn_name, qn_condition_functions=None):
if not (
isinstance(
qn_name,
(
StateQuantumNumberNames,
InteractionQuantumNumberNames,
ParticlePropertyNames,
ParticleDecayPropertyNames,
),
)
):
raise TypeError(
"qn_name has to be of type "
+ "ParticleQuantumNumberNames or "
+ "InteractionQuantumNumberNames or "
+ "ParticlePropertyNames"
)
self.required_qn_names.append(qn_name)
if qn_condition_functions:
for condition in qn_condition_functions:
self.qn_conditions.append(([qn_name], condition))
[docs] def get_required_qn_names(self):
return self.required_qn_names
[docs] @abstractmethod
def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
pass
[docs] def check_requirements(self, in_edges, out_edges, int_node):
for (qn_name_list, condition_functor) in self.get_qn_conditions():
# part_props = [x for x in qn_name_list if isinstance(
# x, ParticlePropertyNames)]
# if part_props:
# return False
if not condition_functor.check(
qn_name_list, in_edges, out_edges, int_node
):
logging.debug(
"condition %s for quantum numbers %s for rule %s not satisfied",
condition_functor.__class__,
qn_name_list,
self.__class__,
)
return False
return True
[docs]class AdditiveQuantumNumberConservation(AbstractRule):
r"""Check for conservation of an additive quantum numbers.
:math:`\sum q_{in} = \sum q_{out}`
Additive quantum numbers are, for example:
- electric charge
- baryon number
- lepton number
"""
def __init__(self, qn_name):
self.qn_name = qn_name
super().__init__(qn_name.name + "Conservation")
[docs] def specify_required_qns(self):
self.add_required_qn(self.qn_name, [DefinedForAllEdges()])
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
in_qn_sum = sum(
[
part[self.qn_name]
for part in ingoing_part_qns
if self.qn_name in part
]
)
out_qn_sum = sum(
[
part[self.qn_name]
for part in outgoing_part_qns
if (self.qn_name in part)
]
)
return in_qn_sum == out_qn_sum
[docs]class ParityConservation(AbstractRule):
"""Check parity conservation."""
def __init__(self):
super().__init__("ParityConservation")
[docs] def specify_required_qns(self):
self.add_required_qn(
StateQuantumNumberNames.Parity, [DefinedForAllEdges()]
)
self.add_required_qn(
InteractionQuantumNumberNames.L, [DefinedForInteractionNode()]
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
r"""Implement :math:`P_{in} = P_{out} \cdot (-1)^L`."""
# is this valid for two outgoing particles only?
parity_label = StateQuantumNumberNames.Parity
parity_in = reduce(
lambda x, y: x * y[parity_label], ingoing_part_qns, 1
)
parity_out = reduce(
lambda x, y: x * y[parity_label], outgoing_part_qns, 1
)
ang_mom = interaction_qns[InteractionQuantumNumberNames.L].magnitude()
if parity_in == (parity_out * (-1) ** ang_mom):
return True
return False
[docs]class ParityConservationHelicity(AbstractRule):
"""Check parity conservation for the helicity formalism."""
def __init__(self):
super().__init__("ParityConservationHelicity")
[docs] def specify_required_qns(self):
self.add_required_qn(
StateQuantumNumberNames.Parity, [DefinedForAllEdges()]
)
self.add_required_qn(
StateQuantumNumberNames.Spin, [DefinedForAllEdges()]
)
self.add_required_qn(
InteractionQuantumNumberNames.ParityPrefactor,
[DefinedForInteractionNode()],
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
r"""Implements the check parity conservation check.
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}
Notice that only the special case :math:`\lambda_1=\lambda_2=0` may
return False.
"""
if len(ingoing_part_qns) == 1 and len(outgoing_part_qns) == 2:
spin_label = StateQuantumNumberNames.Spin
parity_label = StateQuantumNumberNames.Parity
spins = [
x[spin_label].magnitude()
for x in ingoing_part_qns + outgoing_part_qns
]
parity_product = reduce(
lambda x, y: x * y[parity_label],
ingoing_part_qns + outgoing_part_qns,
1,
)
prefactor = parity_product * (-1.0) ** (
spins[1] + spins[2] - spins[0]
)
daughter_hel = [
0
for x in outgoing_part_qns
if x[spin_label].projection() == 0.0
]
if len(daughter_hel) == 2:
if prefactor == -1:
return False
pf_label = InteractionQuantumNumberNames.ParityPrefactor
return prefactor == interaction_qns[pf_label]
return True
[docs]class CParityConservation(AbstractRule):
"""Check for :math:`C`-parity conservation."""
def __init__(self):
super().__init__("CParityConservation")
[docs] def specify_required_qns(self):
self.add_required_qn(StateQuantumNumberNames.CParity)
# the spin quantum number is required to check if the daughter
# particles are fermions or bosons
self.add_required_qn(
StateQuantumNumberNames.Spin,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.CParity]
)
],
)
self.add_required_qn(
InteractionQuantumNumberNames.L,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.CParity]
)
],
)
self.add_required_qn(
InteractionQuantumNumberNames.S,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.CParity]
)
],
)
self.add_required_qn(
ParticlePropertyNames.Pid,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.CParity]
)
],
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
"""Check for :math:`C_{in} = C_{out}`."""
c_parity_in = self.get_c_parity_multiparticle(
ingoing_part_qns, interaction_qns
)
if c_parity_in is None:
return True
c_parity_out = self.get_c_parity_multiparticle(
outgoing_part_qns, interaction_qns
)
if c_parity_out is None:
return True
return c_parity_in == c_parity_out
[docs] @staticmethod
def get_c_parity_multiparticle(part_qns, interaction_qns):
c_parity_label = StateQuantumNumberNames.CParity
pid_label = ParticlePropertyNames.Pid
ang_mom_label = InteractionQuantumNumberNames.L
int_spin_label = InteractionQuantumNumberNames.S
no_c_parity_part = [
part_qns.index(x)
for x in part_qns
if c_parity_label not in x or x[c_parity_label] is None
]
# if all states have C parity defined, then just multiply them
if not no_c_parity_part:
return reduce(lambda x, y: x * y[c_parity_label], part_qns, 1)
# two particle case
if len(part_qns) == 2:
if is_particle_antiparticle_pair(
part_qns[0][pid_label], part_qns[1][pid_label]
):
ang_mom = interaction_qns[ang_mom_label].magnitude()
# if boson
if is_boson(part_qns[0]):
return (-1) ** ang_mom
coupled_spin = interaction_qns[int_spin_label].magnitude()
return (-1) ** (ang_mom + coupled_spin)
return None
[docs]class GParityConservation(AbstractRule):
"""Check for :math:`G`-parity conservation."""
def __init__(self):
super().__init__("GParityConservation")
[docs] def specify_required_qns(self):
self.add_required_qn(StateQuantumNumberNames.GParity)
# the spin quantum number is required to check if the daughter
# particles are fermions or bosons
self.add_required_qn(
StateQuantumNumberNames.Spin,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.GParity]
)
],
)
self.add_required_qn(
InteractionQuantumNumberNames.L,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.GParity]
)
],
)
self.add_required_qn(
InteractionQuantumNumberNames.S,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.GParity]
)
],
)
self.add_required_qn(
StateQuantumNumberNames.IsoSpin,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.GParity]
)
],
)
self.add_required_qn(
ParticlePropertyNames.Pid,
[
DefinedIfOtherQnNotDefinedInOutSeparate(
[StateQuantumNumberNames.GParity]
)
],
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
"""Check for :math:`G_{in} = G_{out}`."""
g_parity_label = StateQuantumNumberNames.GParity
no_g_parity_in_part = [
ingoing_part_qns.index(x)
for x in ingoing_part_qns
if g_parity_label not in x or x[g_parity_label] is None
]
no_g_parity_out_part = [
outgoing_part_qns.index(x)
for x in outgoing_part_qns
if g_parity_label not in x or x[g_parity_label] is None
]
# if all states have G parity defined, then just multiply them
if not no_g_parity_in_part + no_g_parity_out_part:
in_g_parity = reduce(
lambda x, y: x * y[g_parity_label], ingoing_part_qns, 1
)
out_g_parity = reduce(
lambda x, y: x * y[g_parity_label], outgoing_part_qns, 1
)
return in_g_parity == out_g_parity
# two particle case
particle_counts = (len(ingoing_part_qns), len(outgoing_part_qns))
if particle_counts == (1, 2):
if g_parity_label in ingoing_part_qns[0]:
out_g_parity = self.check_multistate_g_parity(
ingoing_part_qns, outgoing_part_qns, interaction_qns
)
in_g_parity = ingoing_part_qns[0][g_parity_label]
if out_g_parity is not None and in_g_parity is not None:
return out_g_parity == in_g_parity
if particle_counts == (2, 1):
if g_parity_label in outgoing_part_qns[0]:
in_g_parity = self.check_multistate_g_parity(
outgoing_part_qns, ingoing_part_qns, interaction_qns
)
out_g_parity = outgoing_part_qns[0][g_parity_label]
if out_g_parity is not None and in_g_parity is not None:
return out_g_parity == in_g_parity
return True
[docs] @staticmethod
def check_multistate_g_parity(
single_state_qns, double_state_qns, interaction_qns
):
isospin_label = StateQuantumNumberNames.IsoSpin
pid_label = ParticlePropertyNames.Pid
ang_mom_label = InteractionQuantumNumberNames.L
int_spin_label = InteractionQuantumNumberNames.S
if is_particle_antiparticle_pair(
double_state_qns[0][pid_label], double_state_qns[1][pid_label]
):
ang_mom = interaction_qns[ang_mom_label].magnitude()
isospin = single_state_qns[0][isospin_label].magnitude()
# if boson
if is_boson(double_state_qns[0]):
return (-1) ** (ang_mom + isospin)
coupled_spin = interaction_qns[int_spin_label].magnitude()
return (-1) ** (ang_mom + coupled_spin + isospin)
return None
[docs]class IdenticalParticleSymmetrization(AbstractRule):
"""Implementation of particle symmetrization."""
def __init__(self):
super().__init__("IdenticalParticleSymmetrization")
[docs] def specify_required_qns(self):
self.add_required_qn(StateQuantumNumberNames.Parity)
self.add_required_qn(
ParticlePropertyNames.Pid, [DefinedForAllOutgoingEdges()]
)
self.add_required_qn(
StateQuantumNumberNames.Spin, [DefinedForAllOutgoingEdges()]
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
if self.check_particles_identical(outgoing_part_qns):
parity_label = StateQuantumNumberNames.Parity
if is_boson(outgoing_part_qns[0]):
# we have a boson, check if parity of mother is even
parity = ingoing_part_qns[0][parity_label]
if parity == -1:
# if its odd then return False
return False
else:
# its fermion
parity = ingoing_part_qns[0][parity_label]
if parity == 1:
return False
return True
[docs] @staticmethod
def check_particles_identical(particles):
"""Check if pids and spins match."""
pid_label = ParticlePropertyNames.Pid
spin_label = StateQuantumNumberNames.Spin
reference_pid = particles[0][pid_label]
reference_spin = particles[0][spin_label]
for particle in particles[1:]:
if particle[pid_label] != reference_pid:
return False
if particle[spin_label] != reference_spin:
return False
return True
[docs]class SpinConservation(AbstractRule):
"""Implementation of conservation of a spin-like quantum number.
That is, for a two body decay (coupling of two particle states). See
:meth:`~.SpinConservation.check` for details.
"""
def __init__(self, spinlike_qn, use_projection=True):
if not isinstance(spinlike_qn, StateQuantumNumberNames):
raise TypeError(
"Expecting Enum of the type "
"ParticleQuantumNumberNames for spinlike_qn"
)
if spinlike_qn not in QNNameClassMapping:
raise ValueError("spinlike_qn is not associated with a QN class")
if QNNameClassMapping[spinlike_qn] is not QuantumNumberClasses.Spin:
raise ValueError("spinlike_qn is not of class Spin")
self.spinlike_qn = spinlike_qn
self.use_projection = use_projection
super().__init__(spinlike_qn.name + "Conservation")
[docs] def specify_required_qns(self):
self.add_required_qn(self.spinlike_qn, [DefinedForAllEdges()])
# for actual spins we include the angular momentum
if self.spinlike_qn is StateQuantumNumberNames.Spin:
self.add_required_qn(
InteractionQuantumNumberNames.L, [DefinedForInteractionNode()]
)
self.add_required_qn(
InteractionQuantumNumberNames.S, [DefinedForInteractionNode()]
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
r"""Check for spin conservation.
Implements
.. math::
|S_1 - S_2| \leq S \leq |S_1 + S_2|
and optionally
.. math::
|L - S| \leq J \leq |L + S|
Also checks :math:`M_1 + M_2 = M` and if Clebsch-Gordan coefficients
are all 0.
"""
spin_label = self.spinlike_qn
in_spins = [x[spin_label] for x in ingoing_part_qns]
out_spins = [x[spin_label] for x in outgoing_part_qns]
if self.use_projection and not self.check_projections(
in_spins, out_spins
):
return False
return self.check_magnitude(in_spins, out_spins, interaction_qns)
[docs] @staticmethod
def check_projections(in_part, out_part):
in_proj = [x.projection() for x in in_part]
out_proj = [x.projection() for x in out_part]
return sum(in_proj) == sum(out_proj)
[docs] def check_magnitude(self, in_part, out_part, interaction_qns):
in_tot_spins = self.calculate_total_spins(in_part, interaction_qns)
out_tot_spins = self.calculate_total_spins(out_part, interaction_qns)
matching_spins = in_tot_spins.intersection(out_tot_spins)
if len(matching_spins) > 0:
return True
return False
[docs] def calculate_total_spins(self, part_list, interaction_qns):
# pylint: disable=too-many-branches
total_spins = set()
if len(part_list) == 1:
if self.use_projection:
total_spins.add(part_list[0])
else:
total_spins.add(Spin(part_list[0].magnitude(), 0))
else:
# first couple all spins together
spins_daughters_coupled = set()
spin_list = deepcopy(part_list)
while spin_list:
if spins_daughters_coupled:
temp_coupled_spins = set()
tempspin = spin_list.pop()
for spin in spins_daughters_coupled:
coupled_spins = self.spin_couplings(spin, tempspin)
temp_coupled_spins.update(coupled_spins)
spins_daughters_coupled = temp_coupled_spins
else:
spins_daughters_coupled.add(spin_list.pop())
if InteractionQuantumNumberNames.L in interaction_qns:
ang_mom = interaction_qns[InteractionQuantumNumberNames.L]
spin = interaction_qns[InteractionQuantumNumberNames.S]
if self.use_projection:
if spin in spins_daughters_coupled:
total_spins.update(self.spin_couplings(spin, ang_mom))
else:
if spin.magnitude() in [
x.magnitude() for x in spins_daughters_coupled
]:
total_spins.update(self.spin_couplings(spin, ang_mom))
else:
if self.use_projection:
total_spins = spins_daughters_coupled
else:
total_spins = [
Spin(x.magnitude(), 0.0)
for x in spins_daughters_coupled
]
return total_spins
[docs] def spin_couplings(self, spin1, spin2):
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`
"""
j_1 = spin1.magnitude()
j_2 = spin2.magnitude()
if self.use_projection:
sum_proj = spin1.projection() + spin2.projection()
possible_spins = [
Spin(x, sum_proj)
for x in arange(abs(j_1 - j_2), j_1 + j_2 + 1, 1).tolist()
if x >= abs(sum_proj)
]
return [
x
for x in possible_spins
if not is_clebsch_gordan_coefficient_zero(spin1, spin2, x)
]
return [
Spin(x, 0)
for x in arange(abs(j_1 - j_2), j_1 + j_2 + 1, 1).tolist()
]
[docs]def is_clebsch_gordan_coefficient_zero(spin1, spin2, spin_coupled):
m_1 = spin1.projection()
j_1 = spin1.magnitude()
m_2 = spin2.projection()
j_2 = spin2.magnitude()
proj = spin_coupled.projection()
mag = spin_coupled.magnitude()
is_zero = False
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:
is_zero = True
elif j_1 == mag and m_1 == -proj:
if abs(j_2 - j_1 - mag) % 2 == 1:
is_zero = True
elif j_2 == mag and m_2 == -proj:
if abs(j_1 - j_2 - mag) % 2 == 1:
is_zero = True
return is_zero
[docs]class ClebschGordanCheckHelicityToCanonical(AbstractRule):
"""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.
"""
def __init__(self):
super().__init__("ClebschGordanCheckHelicityToCanonical")
[docs] def specify_required_qns(self):
self.add_required_qn(
StateQuantumNumberNames.Spin, [DefinedForAllEdges()]
)
self.add_required_qn(
InteractionQuantumNumberNames.L, [DefinedForInteractionNode()]
)
self.add_required_qn(
InteractionQuantumNumberNames.S, [DefinedForInteractionNode()]
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
if len(ingoing_part_qns) == 1 and len(outgoing_part_qns) == 2:
spin_label = StateQuantumNumberNames.Spin
in_spins = [x[spin_label] for x in ingoing_part_qns]
out_spins = [x[spin_label] for x in outgoing_part_qns]
out_spins[1] = Spin(
out_spins[1].magnitude(), -out_spins[1].projection()
)
helicity_diff = sum([x.projection() for x in out_spins])
ang_mom = interaction_qns[InteractionQuantumNumberNames.L]
spin = interaction_qns[InteractionQuantumNumberNames.S]
if spin.magnitude() < abs(helicity_diff) or in_spins[
0
].magnitude() < abs(helicity_diff):
return False
spin = Spin(spin.magnitude(), helicity_diff)
if is_clebsch_gordan_coefficient_zero(
out_spins[0], out_spins[1], spin
):
return False
in_spins[0] = Spin(in_spins[0].magnitude(), helicity_diff)
return not is_clebsch_gordan_coefficient_zero(
ang_mom, spin, in_spins[0]
)
return False
[docs]class HelicityConservation(AbstractRule):
"""Implementation of helicity conservation."""
def __init__(self):
super().__init__("HelicityConservation")
[docs] def specify_required_qns(self):
self.add_required_qn(
StateQuantumNumberNames.Spin, [DefinedForAllEdges()]
)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
r"""Check for :math:`|\lambda_2-\lambda_3| \leq S_1`."""
if len(ingoing_part_qns) == 1 and len(outgoing_part_qns) == 2:
spin_label = StateQuantumNumberNames.Spin
mother_spin = ingoing_part_qns[0][spin_label].magnitude()
daughter_hel = [
x[spin_label].projection() for x in outgoing_part_qns
]
if mother_spin >= abs(daughter_hel[0] - daughter_hel[1]):
return True
return False
[docs]class GellMannNishijimaRule(AbstractRule):
"""Conservation rule for Gell-Mann-Nishijima."""
def __init__(self):
super().__init__("GellMannNishijimaRule")
[docs] def specify_required_qns(self):
self.add_required_qn(
StateQuantumNumberNames.Charge, [DefinedForAllEdges()]
)
self.add_required_qn(
StateQuantumNumberNames.IsoSpin, [DefinedForAllEdges()]
)
self.add_required_qn(StateQuantumNumberNames.Strangeness)
self.add_required_qn(StateQuantumNumberNames.Charm)
self.add_required_qn(StateQuantumNumberNames.Bottomness)
self.add_required_qn(StateQuantumNumberNames.Topness)
self.add_required_qn(StateQuantumNumberNames.BaryonNumber)
self.add_required_qn(StateQuantumNumberNames.ElectronLN)
self.add_required_qn(StateQuantumNumberNames.MuonLN)
self.add_required_qn(StateQuantumNumberNames.TauLN)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
r"""Check the Gell-Mann–Nishijima formula.
:math:`Q=I_3+\frac{Y}{2}` for each particle.
"""
charge_label = StateQuantumNumberNames.Charge
isospin_label = StateQuantumNumberNames.IsoSpin
eln = StateQuantumNumberNames.ElectronLN
mln = StateQuantumNumberNames.MuonLN
tln = StateQuantumNumberNames.TauLN
for particle in ingoing_part_qns + outgoing_part_qns:
if (
sum(
[
abs(particle[x])
for x in [eln, mln, tln]
if x in particle
]
)
> 0.0
):
# if particle is a lepton then skip the check
continue
isospin_3 = 0
if isospin_label in particle:
isospin_3 = particle[isospin_label].projection()
if float(particle[charge_label]) != (
isospin_3 + 0.5 * self.calculate_hypercharge(particle)
):
return False
return True
[docs] @staticmethod
def calculate_hypercharge(particle):
"""Calculate the hypercharge :math:`Y=S+C+B+T+B`."""
qn_labels = [
StateQuantumNumberNames.Strangeness,
StateQuantumNumberNames.Charm,
StateQuantumNumberNames.Bottomness,
StateQuantumNumberNames.Topness,
StateQuantumNumberNames.BaryonNumber,
]
qn_values = [particle[x] for x in qn_labels if x in particle]
return sum(qn_values)
[docs]class MassConservation(AbstractRule):
"""Mass conservation rule."""
def __init__(self, width_factor=3):
self.width_factor = width_factor
super().__init__("MassConservation")
[docs] def specify_required_qns(self):
self.add_required_qn(
ParticlePropertyNames.Mass, [DefinedForAllEdges()]
)
self.add_required_qn(ParticleDecayPropertyNames.Width)
[docs] def check(self, ingoing_part_qns, outgoing_part_qns, interaction_qns):
r"""Implement the mass check.
: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_label = ParticlePropertyNames.Mass
width_label = ParticleDecayPropertyNames.Width
mass_in = sum([x[mass_label] for x in ingoing_part_qns])
width_in = sum(
[x[width_label] for x in ingoing_part_qns if width_label in x]
)
mass_out = sum([x[mass_label] for x in outgoing_part_qns])
width_out = sum(
[x[width_label] for x in outgoing_part_qns if width_label in x]
)
return (mass_out - self.width_factor * width_out) < (
mass_in + self.width_factor * width_in
)