diff --git a/.cspell.json b/.cspell.json index cd1d3817..16f77f36 100644 --- a/.cspell.json +++ b/.cspell.json @@ -197,6 +197,7 @@ "pyplot", "pytest", "qrules", + "regge", "setuptools", "spflueger", "struct", diff --git a/docs/conf.py b/docs/conf.py index 9cec0db1..806e48a8 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -63,15 +63,12 @@ def __get_newtypes(some_type: type) -> list: "NodeQuantumNumber": ("obj", "qrules.quantum_numbers.NodeQuantumNumber"), "NodeQuantumNumberTypes": ("obj", "qrules.quantum_numbers.NodeQuantumNumberTypes"), "NodeType": "typing.TypeVar", - "ParticleWithSpin": ("obj", "qrules.particle.ParticleWithSpin"), "Path": "pathlib.Path", "QNTransition": ("obj", "qrules.workflow.QNTransition"), "qrules.topology.EdgeType": "typing.TypeVar", "qrules.topology.NodeType": "typing.TypeVar", "Rule": ("obj", "qrules.argument_handling.Rule"), "SpinFormalism": ("obj", "qrules.transition.SpinFormalism"), - "StateDefinition": ("obj", "qrules.combinatorics.StateDefinition"), - "StateDefinitionInput": ("obj", "qrules.combinatorics.StateDefinitionInput"), "StateTransition": ("obj", "qrules.transition.StateTransition"), "TypeAliasForwardRef": ("obj", "typing.TypeAlias"), "typing.Literal[-1, 1]": "typing.Literal", diff --git a/docs/index.md b/docs/index.md index d7c5827e..c5c717ad 100644 --- a/docs/index.md +++ b/docs/index.md @@ -27,7 +27,7 @@ quantum numbers it has found. The resulting state transition objects are particularly useful for **amplitude analysis / Partial Wave Analysis** as they contain all information (such as expected masses, -widths, and spin projections) that is needed to formulate an amplitude model. +widths, and spin quantum numbers) that is needed to formulate an amplitude model. The {doc}`/usage` pages illustrate several features of {mod}`qrules`. You can run each of them as Jupyter notebooks with the {fa}`rocket` launch button in the top-right diff --git a/docs/usage.ipynb b/docs/usage.ipynb index cee4d0d2..5e49e7cf 100644 --- a/docs/usage.ipynb +++ b/docs/usage.ipynb @@ -204,6 +204,7 @@ "---\n", "usage/reaction\n", "usage/qn-transitions\n", + "usage/spin-projections\n", "usage/production\n", "usage/particle\n", "usage/visualize\n", diff --git a/docs/usage/conservation.ipynb b/docs/usage/conservation.ipynb index 1d26f56a..efdec5cb 100644 --- a/docs/usage/conservation.ipynb +++ b/docs/usage/conservation.ipynb @@ -191,7 +191,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next, have a look at the edge and node properties, and use the underlying {class}`.Topology` to extract one of the node {class}`.InteractionProperties` with the surrounding states (these are {obj}`tuple`s of a {class}`.Particle` and a {obj}`float` spin projection)." + "Next, have a look at the edge and node properties, and use the underlying {class}`.Topology` to extract one of the node {class}`.InteractionProperties` with the surrounding states (these are {class}`.Particle` instances)." ] }, { @@ -244,10 +244,10 @@ "interaction = transition.interactions[node_id]\n", "\n", "spin_magnitude_conservation(\n", - " ingoing_spin_magnitudes=[{\"spin_magnitude\": incoming_state.particle.spin}],\n", + " ingoing_spin_magnitudes=[{\"spin_magnitude\": incoming_state.spin}],\n", " outgoing_spin_magnitudes=[\n", - " {\"spin_magnitude\": outgoing_state1.particle.spin},\n", - " {\"spin_magnitude\": outgoing_state2.particle.spin},\n", + " {\"spin_magnitude\": outgoing_state1.spin},\n", + " {\"spin_magnitude\": outgoing_state2.spin},\n", " ],\n", " interaction_qns={\n", " \"l_magnitude\": interaction.l_magnitude,\n", @@ -260,41 +260,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Contrary to expectations, this transition does not conserve spin **projection** and therefore {func}`.spin_conservation` returns {obj}`False`:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "spin_conservation(\n", - " ingoing_spins=[\n", - " SpinEdgeInput(\n", - " spin_magnitude=incoming_state.particle.spin,\n", - " spin_projection=incoming_state.spin_projection,\n", - " )\n", - " ],\n", - " outgoing_spins=[\n", - " SpinEdgeInput(\n", - " spin_magnitude=outgoing_state1.particle.spin,\n", - " spin_projection=outgoing_state1.spin_projection,\n", - " ),\n", - " SpinEdgeInput(\n", - " spin_magnitude=outgoing_state2.particle.spin,\n", - " spin_projection=outgoing_state2.spin_projection,\n", - " ),\n", - " ],\n", - " interaction_qns=interaction,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The reason is that AmpForm formulates the {class}`~ampform.helicity.HelicityModel` with the helicity formalism first and then uses a transformation to get the model in the canonical basis (see {func}`~ampform.helicity.formulate_isobar_cg_coefficients`). The canonical basis does not conserve helicity (taken to be {attr}`.State.spin_projection`)." + "{func}`.spin_conservation` additionally checks the spin **projections**, as illustrated in the {ref}`examples above `. The states and interactions generated by the workflow carry no spin projections, so that rule can only be evaluated with manually supplied projections — or by re-enabling spin projections in the problem sets, as demonstrated in {doc}`/usage/spin-projections`." ] }, { diff --git a/docs/usage/production.ipynb b/docs/usage/production.ipynb index af0a8e4d..3b641c05 100644 --- a/docs/usage/production.ipynb +++ b/docs/usage/production.ipynb @@ -239,7 +239,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Production reactions also work through the classic {func}`.generate_transitions` interface, which matches all intermediate edges — including the exchange edges — to particles and generates the spin projections required for a helicity amplitude model. {meth}`.ReactionInfo.group_by_channel` shows which resonances and exchange particles appear in each channel:" + "Production reactions also work through the classic {func}`.generate_transitions` interface, which matches all intermediate edges — including the exchange edges — to particles from a database. {meth}`.ReactionInfo.group_by_channel` shows which resonances and exchange particles appear in each channel:" ] }, { @@ -264,7 +264,7 @@ ")\n", "{\n", " channel: sorted({\n", - " state.particle.name\n", + " state.name\n", " for transition in transitions\n", " for state in transition.intermediate_states.values()\n", " })\n", @@ -288,6 +288,120 @@ "Markdown(qrules.io.asmermaid(particle_reaction, collapse=\"topology\", markdown=True))" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Three-body production and double exchange" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For three or more final-state particles, the topology set contains more than the single-exchange $s$-, $t$-, and $u$-channels, in which the exchange edge fuses into a resonance that continues as a sequence of two-body decays. In the additional **double-exchange** topologies, two exchange edges fuse into a final-state particle — the central-production and ladder diagrams known from double-Regge phenomenology. {func}`.determine_reaction_channel` labels them with one letter per exchange edge:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "topologies_3body = create_isobar_topologies(\n", + " number_of_final_states=3,\n", + " number_of_initial_states=2,\n", + ")\n", + "[determine_reaction_channel(topology) for topology in topologies_3body]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Markdown(qrules.io.asmermaid(topologies_3body, markdown=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As an example with particles, add a pion to the $\\eta$ photoproduction final state: $\\gamma p \\to p\\pi^0\\eta$, with a small selection of allowed intermediate states. All channel types contribute. In the two-letter channels, the final-state assignment determines whether each exchange edge is of $t$- or $u$-type:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "reaction_3body = qrules.generate_transitions(\n", + " initial_state=[\"gamma\", \"p\"],\n", + " final_state=[\"p\", \"pi0\", \"eta\"],\n", + " allowed_intermediate_particles=[\"N(1535)\", \"a(2)(1320)\", \"omega(782)\"],\n", + " allowed_interaction_types=[\"strong\", \"em\"],\n", + " formalism=\"helicity\",\n", + " particle_db=PDG,\n", + ")\n", + "{\n", + " channel: sorted({\n", + " state.name\n", + " for transition in transitions\n", + " for state in transition.intermediate_states.values()\n", + " })\n", + " for channel, transitions in reaction_3body.group_by_channel().items()\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The single-letter channels are the familiar diagrams from the two-body case, now continued with a two-body decay of the resonance. In the double-exchange channels, the collapsed diagrams show how a final-state particle emerges from the two exchange edges. The :code:`\"uu\"` channel proceeds through double $\\omega$ exchange in a ladder configuration:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "channels_3body = reaction_3body.group_by_channel()\n", + "Markdown(qrules.io.asmermaid(channels_3body[\"uu\"], collapse=\"topology\", markdown=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "while the :code:`\"tt\"` channel combines central-production diagrams — both initial states emit a final-state particle and the two baryon exchanges fuse — with double baryon-exchange ladders:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Markdown(qrules.io.asmermaid(channels_3body[\"tt\"], collapse=\"topology\", markdown=True))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{tip}\n", + "If the double-exchange contributions are not of interest, pass e.g. :code:`allowed_channels=[\"s\", \"t\", \"u\"]` to drop those topologies before any solving happens.\n", + ":::\n", + "\n", + ":::{warning}\n", + "For production reactions, always restrict :code:`allowed_intermediate_particles`. The problem-set count grows rapidly with the number of final-state particles (number of topologies × final-state permutations × interaction-type combinations — already 240 problem sets for this $2 \\to 3$ example), and each problem set contains *two* internal edges. Without a particle selection, the quantum-number domains of those edges are unconstrained, so the constraint solver has to backtrack through orders of magnitude more combinations per problem set. A :code:`allowed_channels` selection reduces the number of problem sets further.\n", + ":::" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -366,7 +480,7 @@ "metadata": {}, "source": [ ":::{seealso}\n", - "{doc}`/usage/qn-transitions` for solving reactions at the quantum-number level — without spin projections — which keeps many-body final states tractable, and {doc}`/usage/reaction` for the general workflow.\n", + "{doc}`/usage/qn-transitions` for solving reactions at the quantum-number level, which keeps many-body final states tractable, and {doc}`/usage/reaction` for the general workflow.\n", ":::" ] } diff --git a/docs/usage/qn-transitions.ipynb b/docs/usage/qn-transitions.ipynb index 950ab4e4..1b1e3799 100644 --- a/docs/usage/qn-transitions.ipynb +++ b/docs/usage/qn-transitions.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Transitions without spin projections\n", + "# Quantum-number transitions\n", "\n", ":::{autolink-concat}\n", ":::" @@ -14,7 +14,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The workflows described in {doc}`/usage/reaction` generate a transition for every allowed combination of **spin projections** of the initial and final state, because a helicity amplitude model requires each of those combinations. If you are only interested in which intermediate states and quantum numbers are allowed — for instance, which $J^{PC}$ resonances can appear in a Dalitz-plot decomposition — the spin projections merely multiply the number of {class}`.QNProblemSet`s that have to be solved. This page shows how to generate transitions directly at the $J^{P(C)}$ level with the {mod}`.workflow` module, which is considerably faster." + "The workflows described in {doc}`/usage/reaction` match every intermediate edge of the solved transitions to particles from a database. If you are only interested in which quantum numbers are allowed — for instance, which $J^{PC}$ resonances can appear in a Dalitz-plot decomposition — that matching step is unnecessary. This page shows how to generate transitions directly at the $J^{P(C)}$ level with the {mod}`.workflow` module: the intermediate states then remain *sets* of quantum numbers and are never matched against a particle database." ] }, { @@ -47,36 +47,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Problem sets without spin projections" + "## Quantum-number problem sets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "As an example, take the reaction $J/\\psi \\to \\gamma\\pi^0\\pi^0$ with two $f_0$ resonances as allowed intermediate states. By default, {func}`.create_qn_problem_sets` expands the initial and final state over all combinations of their allowed spin projections:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "expanded = create_qn_problem_sets(\n", - " initial_state=[\"J/psi(1S)\"],\n", - " final_state=[\"gamma\", \"pi0\", \"pi0\"],\n", - " particle_db=PDG,\n", - " allowed_intermediate_particles=[\"f(0)(980)\", \"f(0)(1500)\"],\n", - ")\n", - "sum(len(problems) for problems in expanded.problem_sets.values())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The $J/\\psi$ has three allowed spin projections and the photon has two, so every problem set appears in six spin-projection combinations. With :code:`spin_projections=False`, this Cartesian expansion is skipped altogether and the problem sets only differ in decay topology and interaction types:" + "As an example, take the reaction $J/\\psi \\to \\gamma\\pi^0\\pi^0$ with two $f_0$ resonances as allowed intermediate states. {func}`.create_qn_problem_sets` fans the initial and final state out into a {class}`.QNProblemSet` for every decay topology, kinematic permutation, and combination of allowed interaction types:" ] }, { @@ -90,18 +68,10 @@ " final_state=[\"gamma\", \"pi0\", \"pi0\"],\n", " particle_db=PDG,\n", " allowed_intermediate_particles=[\"f(0)(980)\", \"f(0)(1500)\"],\n", - " spin_projections=False,\n", ")\n", "sum(len(problems) for problems in qn_problem_sets.problem_sets.values())" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The resulting problem sets contain no {attr}`~.EdgeQuantumNumbers.spin_projection`, {attr}`~.NodeQuantumNumbers.l_projection`, or {attr}`~.NodeQuantumNumbers.s_projection` quantum numbers at all. They are identical to what {func}`.strip_spin_projections` produces from the expanded collection, but without ever generating the expansion." - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -113,7 +83,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "{func}`.find_qn_transitions` solves these problem sets purely at the quantum-number level: no particle database is consulted for the intermediate states. It returns {obj}`.QNTransition`s, whose states and interactions are property maps of quantum numbers instead of {class}`.State` and {class}`.InteractionProperties` objects." + "{func}`.find_qn_transitions` solves these problem sets purely at the quantum-number level: no particle database is consulted for the intermediate states. It returns {obj}`.QNTransition`s, whose states and interactions are property maps of quantum numbers instead of {class}`.Particle` and {class}`.InteractionProperties` objects." ] }, { @@ -213,71 +183,6 @@ "Markdown(source)" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Fewer combinatorics for larger reactions" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The spin-projection expansion multiplies the number of problem sets by $\\prod_i(2s_i+1)$ over all initial and final states $i$ (with the exception of massless states, which have no $0$ projection). Skipping it therefore matters most for many-body final states with spin. Take $J/\\psi \\to p\\bar p\\pi^0\\pi^0$, where the expansion factor is $3 \\times 2 \\times 2 = 12$:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%time\n", - "expanded = create_qn_problem_sets(\n", - " initial_state=[\"J/psi(1S)\"],\n", - " final_state=[\"p\", \"p~\", \"pi0\", \"pi0\"],\n", - " particle_db=PDG,\n", - " allowed_intermediate_particles=[\"N(1440)\"],\n", - ")\n", - "sum(len(problems) for problems in expanded.problem_sets.values())" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%time\n", - "qn_problem_sets = create_qn_problem_sets(\n", - " initial_state=[\"J/psi(1S)\"],\n", - " final_state=[\"p\", \"p~\", \"pi0\", \"pi0\"],\n", - " particle_db=PDG,\n", - " allowed_intermediate_particles=[\"N(1440)\"],\n", - " spin_projections=False,\n", - ")\n", - "sum(len(problems) for problems in qn_problem_sets.problem_sets.values())" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The projection-free problem sets are not only fewer, but each of them is also cheaper to solve, because the spin projections do not appear as variables in the constraint problem:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "%%time\n", - "qn_transitions = find_qn_transitions(qn_problem_sets)\n", - "len(qn_transitions)" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -289,7 +194,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Solving at the quantum-number level makes reactions with four- and even five-body final states feasible, for which the workflows *with* spin projections take impractically long (see [ComPWA/qrules#27](https://github.com/ComPWA/qrules/issues/27)). Take $J/\\psi \\to K^+K^-\\pi^+\\pi^-$ with $\\phi(1020)$ and $\\rho(770)$ resonances, where {func}`.create_qn_problem_sets`' :code:`final_state_groupings` argument limits the subsystems to $\\phi\\to K^+K^-$ and $\\rho^0\\to\\pi^+\\pi^-$:" + "Since the intermediate states never have to be matched against a particle database, solving at the quantum-number level keeps reactions with four- and even five-body final states feasible (see [ComPWA/qrules#27](https://github.com/ComPWA/qrules/issues/27)). Take $J/\\psi \\to K^+K^-\\pi^+\\pi^-$ with $\\phi(1020)$ and $f_0(980)$ resonances, where {func}`.create_qn_problem_sets`' :code:`final_state_groupings` argument limits the subsystems to $\\phi\\to K^+K^-$ and $f_0\\to\\pi^+\\pi^-$:" ] }, { @@ -303,7 +208,7 @@ " initial_state=\"J/psi(1S)\",\n", " final_state=[\"K+\", \"K-\", \"pi+\", \"pi-\"],\n", " particle_db=PDG,\n", - " allowed_intermediate_particles=[\"phi(1020)\", \"rho(770)\"],\n", + " allowed_intermediate_particles=[\"phi(1020)\", \"f(0)(980)\"],\n", " allowed_interaction_types=\"strong\",\n", " final_state_groupings=[[[\"K+\", \"K-\"], [\"pi+\", \"pi-\"]]],\n", ")\n", @@ -355,7 +260,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Even a five-body final state remains tractable — here with an additional $\\pi^0$, which lets the number of allowed decay topologies and quantum-number combinations grow substantially:" + "Even a five-body final state remains tractable. Here, an additional $\\pi^0$ and a $\\rho(770)$ resonance let the number of allowed decay topologies and quantum-number combinations grow substantially:" ] }, { @@ -369,7 +274,7 @@ " initial_state=\"J/psi(1S)\",\n", " final_state=[\"K+\", \"K-\", \"pi+\", \"pi-\", \"pi0\"],\n", " particle_db=PDG,\n", - " allowed_intermediate_particles=[\"phi(1020)\", \"rho(770)\"],\n", + " allowed_intermediate_particles=[\"phi(1020)\", \"f(0)(980)\", \"rho(770)\"],\n", " allowed_interaction_types=\"strong\",\n", " final_state_groupings=[[[\"K+\", \"K-\"], [\"pi+\", \"pi-\"]]],\n", ")\n", @@ -381,7 +286,7 @@ "metadata": {}, "source": [ ":::{seealso}\n", - "To generate particle-level transitions *with* spin projections — as required for a helicity amplitude model — use {func}`.find_solutions` or {func}`.generate_transitions` as described in {doc}`/usage/reaction`. If the problem sets have already been created with spin projections, they can still be reduced afterwards with {func}`.strip_spin_projections`. Reactions with two initial states are described in {doc}`/usage/production`.\n", + "To generate particle-level transitions — with the intermediate states matched to a particle database, as required for an amplitude model — use {func}`.find_solutions` or {func}`.generate_transitions` as described in {doc}`/usage/reaction`. The problem sets can be extended with additional quantum numbers and rules before solving, as demonstrated in {doc}`/usage/spin-projections`. Reactions with two initial states are described in {doc}`/usage/production`.\n", ":::" ] } diff --git a/docs/usage/reaction.ipynb b/docs/usage/reaction.ipynb index 3bf31c7c..ef8b3527 100644 --- a/docs/usage/reaction.ipynb +++ b/docs/usage/reaction.ipynb @@ -180,7 +180,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Each {class}`.ProblemSet` provides a mapping of {attr}`~.ProblemSet.initial_facts` that represent the initial and final states with spin projections. The nodes and edges in between these {attr}`~.ProblemSet.initial_facts` are still to be generated. This will be done from the provided {attr}`~.ProblemSet.solving_settings`. There are two mechanisms there:\n", + "Each {class}`.ProblemSet` provides a mapping of {attr}`~.ProblemSet.initial_facts` that represent the known initial and final states. The nodes and edges in between these {attr}`~.ProblemSet.initial_facts` are still to be generated. This will be done from the provided {attr}`~.ProblemSet.solving_settings`. There are two mechanisms there:\n", "\n", "1. One the one hand, the {attr}`.EdgeSettings.qn_domains` and {attr}`.NodeSettings.qn_domains` contained in the {attr}`~.ProblemSet.solving_settings` define the **domain** over which quantum number sets can be generated.\n", "2. On the other, the {attr}`.EdgeSettings.conservation_rules` and {attr}`.NodeSettings.conservation_rules` in {attr}`~.ProblemSet.solving_settings` define which **{mod}`.conservation_rules`** are used to determine which of the sets of generated quantum numbers are valid (rules with a higher priority are evaluated first).\n", @@ -199,11 +199,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "If you are happy with the default settings generated by the {class}`.StateTransitionManager`, just start with solving directly!\n", - "\n", - "```{toggle}\n", - "This step takes about 23 sec on an Intel(R) Core(TM) i7-6820HQ CPU of 2.70GHz running, multi-threaded.\n", - "```" + "If you are happy with the default settings generated by the {class}`.StateTransitionManager`, just start with solving directly!" ] }, { @@ -245,9 +241,9 @@ "class: dropdown\n", "----\n", "\n", - "The \"number of {attr}`~.ReactionInfo.transitions`\" is the total number of allowed {obj}`.MutableTransition` instances that the {class}`.StateTransitionManager` has found. This also includes all allowed **spin projection combinations**. In this channel, we for example consider a $J/\\psi$ with spin projection $\\pm1$ that decays into a $\\gamma$ with spin projection $\\pm1$, which already gives us four possibilities.\n", + "The \"number of {attr}`~.ReactionInfo.transitions`\" is the total number of allowed {obj}`.MutableTransition` instances that the {class}`.StateTransitionManager` has found. Transitions with the same intermediate particle can still differ in, for instance, the interaction types of their nodes or (in the canonical bases) their $LS$-couplings.\n", "\n", - "On the other hand, the intermediate state names that was extracted with {meth}`.ReactionInfo.get_intermediate_particles`, is just a {obj}`set` of the state names on the intermediate edges of the list of {attr}`~.ReactionInfo.transitions`, regardless of spin projection.\n", + "On the other hand, the intermediate state names that was extracted with {meth}`.ReactionInfo.get_intermediate_particles`, is just a {obj}`set` of the state names on the intermediate edges of the list of {attr}`~.ReactionInfo.transitions`.\n", "```" ] }, diff --git a/docs/usage/reaction/intermediate-results.ipynb b/docs/usage/reaction/intermediate-results.ipynb index f5f4e8a3..a7904351 100644 --- a/docs/usage/reaction/intermediate-results.ipynb +++ b/docs/usage/reaction/intermediate-results.ipynb @@ -88,7 +88,6 @@ "particle_db = load_pdg()\n", "interaction_config = InteractionConfig(\n", " create_interaction_settings(\n", - " formalism=\"canonical-helicity\",\n", " particle_db=particle_db,\n", " )\n", ")\n", diff --git a/docs/usage/spin-projections.ipynb b/docs/usage/spin-projections.ipynb new file mode 100644 index 00000000..ed370733 --- /dev/null +++ b/docs/usage/spin-projections.ipynb @@ -0,0 +1,385 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Adding spin projections\n", + "\n", + ":::{autolink-concat}\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The default workflow solves reactions at the level of spin **magnitudes**: it does not track the spin projections of the states, because those only multiply the number of transitions without affecting which intermediate states are allowed. The underlying constraint solver, however, is agnostic about which quantum numbers it propagates — it solves for whatever **facts**, **domains**, and **rules** the {class}`.QNProblemSet`s declare. This page demonstrates that extensibility by re-introducing spin projections for the reaction $J/\\psi \\to \\gamma f_2(1270) \\to \\gamma\\pi^0\\pi^0$ and pruning the helicity combinations with {func}`.helicity_conservation`.\n", + "\n", + ":::{seealso}\n", + "The level of detail can be reduced further as well: with `ls_couplings=False`, {func}`.generate_transitions` and {func}`.create_qn_problem_sets` do not enumerate $LS$-combinations either. The existence rules {class}`.SpinCoupling`, {class}`.SpinParityCoupling`, {class}`.CParityCoupling`, and {class}`.GParityCoupling` then only check that *some* coupling with $L \\leq L_\\mathrm{max}$ can produce the spins and parities, and the allowed $(L, S)$ combinations can be reconstructed from the solutions afterwards (see {func}`.create_interaction_settings`).\n", + ":::" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "tags": [ + "hide-cell" + ] + }, + "outputs": [], + "source": [ + "from collections import defaultdict\n", + "from fractions import Fraction\n", + "\n", + "import attrs\n", + "from IPython.display import Markdown\n", + "\n", + "import qrules\n", + "import qrules.io\n", + "from qrules.argument_handling import GraphEdgePropertyMap\n", + "from qrules.conservation_rules import helicity_conservation, spin_validity\n", + "from qrules.quantum_numbers import EdgeQuantumNumbers\n", + "from qrules.settings import CONSERVATION_LAW_PRIORITIES, EDGE_RULE_PRIORITIES\n", + "from qrules.solving import QNProblemSet\n", + "from qrules.topology import MutableTransition\n", + "from qrules.workflow import create_qn_problem_sets, find_qn_transitions\n", + "\n", + "PDG = qrules.load_pdg()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A projection-free reaction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, create the default problem sets with {func}`.create_qn_problem_sets` (see {doc}`/usage/qn-transitions`). The final-state grouping selects the decay topology in which the two pions form the $f_2(1270)$ resonance:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "collection = create_qn_problem_sets(\n", + " initial_state=[\"J/psi(1S)\"],\n", + " final_state=[\"gamma\", \"pi0\", \"pi0\"],\n", + " particle_db=PDG,\n", + " allowed_intermediate_particles=[\"f(2)(1270)\"],\n", + " allowed_interaction_types=[\"strong\", \"EM\"],\n", + " max_angular_momentum=2,\n", + " final_state_groupings=[[[\"pi0\", \"pi0\"]]],\n", + ")\n", + "sum(len(problem_sets) for problem_sets in collection.problem_sets.values())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Each {class}`.QNProblemSet` declares the known quantum numbers of the initial and final states (its **facts**, {attr}`~.QNProblemSet.initial_facts`) and the value ranges to solve for on the intermediate edges and interaction nodes (its **domains**, part of the {attr}`~.QNProblemSet.solving_settings`). Neither mentions a {attr}`~.EdgeQuantumNumbers.spin_projection`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sorted({\n", + " qn_type.__name__\n", + " for problem_sets in collection.problem_sets.values()\n", + " for problem_set in problem_sets\n", + " for prop_map in problem_set.initial_facts.states.values()\n", + " for qn_type in prop_map\n", + "})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding facts and domains" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Spin projections enter the problem sets through the same two channels as any other quantum number. The external edges get {attr}`~.EdgeQuantumNumbers.spin_projection` **facts**: a single value fixes the projection, while a `list` of values is registered by the {class}`.CSPSolver` as a variable with that list as its domain, so that several projection cases are solved within a single problem set. Here, the $J/\\psi$ is restricted to $\\lambda=\\pm1$ (as it is produced from $e^+e^-$ collisions), the photon to its two helicity states, and the pions are fixed to $0$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spin_projections = {\n", + " \"J/psi(1S)\": [Fraction(-1), Fraction(+1)],\n", + " \"gamma\": [Fraction(-1), Fraction(+1)],\n", + " \"pi0\": Fraction(0),\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The intermediate edge gets a projection **domain**, derived from the spin magnitudes that its {class}`.EdgeSettings` allow, together with the {func}`.spin_validity` rule, which requires $|m| \\leq s$ on a single edge. Every rule carries a priority (higher priorities are evaluated first); take the defaults from {obj}`.EDGE_RULE_PRIORITIES` and {obj}`.CONSERVATION_LAW_PRIORITIES`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def get_particle_name(prop_map: GraphEdgePropertyMap) -> str:\n", + " pid = prop_map[EdgeQuantumNumbers.pid]\n", + " assert pid is not None\n", + " return PDG.find(int(pid)).name\n", + "\n", + "\n", + "def add_spin_projections(problem_set: QNProblemSet) -> QNProblemSet:\n", + " facts = problem_set.initial_facts\n", + " new_states = {\n", + " edge_id: {\n", + " **prop_map,\n", + " EdgeQuantumNumbers.spin_projection: spin_projections[\n", + " get_particle_name(prop_map)\n", + " ],\n", + " }\n", + " for edge_id, prop_map in facts.states.items()\n", + " }\n", + " new_facts = MutableTransition(facts.topology, new_states, dict(facts.interactions))\n", + " settings = problem_set.solving_settings\n", + " new_edge_settings = {}\n", + " for edge_id, edge_settings in settings.states.items():\n", + " if edge_id not in facts.topology.intermediate_edge_ids:\n", + " new_edge_settings[edge_id] = edge_settings\n", + " continue\n", + " max_spin = max(edge_settings.qn_domains[EdgeQuantumNumbers.spin_magnitude])\n", + " projection_domain = [\n", + " Fraction(x, 2) for x in range(-int(2 * max_spin), int(2 * max_spin) + 1)\n", + " ]\n", + " new_edge_settings[edge_id] = attrs.evolve(\n", + " edge_settings,\n", + " conservation_rules={\n", + " **edge_settings.conservation_rules,\n", + " spin_validity: EDGE_RULE_PRIORITIES[spin_validity],\n", + " },\n", + " qn_domains={\n", + " **edge_settings.qn_domains,\n", + " EdgeQuantumNumbers.spin_projection: projection_domain,\n", + " },\n", + " )\n", + " new_settings = MutableTransition(\n", + " settings.topology, new_edge_settings, dict(settings.interactions)\n", + " )\n", + " return QNProblemSet(initial_facts=new_facts, solving_settings=new_settings)\n", + "\n", + "\n", + "collection.problem_sets = {\n", + " strength: [add_spin_projections(problem_set) for problem_set in problem_sets]\n", + " for strength, problem_sets in collection.problem_sets.items()\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The solver now propagates the projections just like any other quantum number:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "qn_transitions = find_qn_transitions(collection)\n", + "len(qn_transitions)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Every transition carries a spin projection for every state, but nothing constrains the projections yet. All combinations of the photon helicity $\\lambda_\\gamma$ and the $f_2(1270)$ helicity $\\lambda_{f_2}$ appear:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "jupyter": { + "source_hidden": true + }, + "mystnb": { + "code_prompt_show": "Helper to collect helicity combinations" + }, + "tags": [ + "hide-input" + ] + }, + "outputs": [], + "source": [ + "def get_helicity_combinations(transitions) -> dict[int, list[int]]:\n", + " combinations = defaultdict(set)\n", + " for transition in transitions:\n", + " topology = transition.topology\n", + " resonance_edge = next(iter(topology.intermediate_edge_ids))\n", + " gamma_edge = next(\n", + " i\n", + " for i in topology.outgoing_edge_ids\n", + " if transition.states[i][EdgeQuantumNumbers.pid] == 22\n", + " )\n", + " lambda_gamma = transition.states[gamma_edge][EdgeQuantumNumbers.spin_projection]\n", + " lambda_f2 = transition.states[resonance_edge][\n", + " EdgeQuantumNumbers.spin_projection\n", + " ]\n", + " combinations[int(lambda_gamma)].add(int(lambda_f2))\n", + " return {k: sorted(v) for k, v in sorted(combinations.items())}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "get_helicity_combinations(qn_transitions)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Adding a conservation rule" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The projections become meaningful once a **rule** constrains them. {func}`.helicity_conservation` checks $|\\lambda_2-\\lambda_3| \\leq S_1$ at each decay node $1 \\to 2\\,3$. Adding it to the conservation rules of the interaction nodes works the same way as adding {func}`.spin_validity` to the edges:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def add_helicity_conservation(problem_set: QNProblemSet) -> QNProblemSet:\n", + " settings = problem_set.solving_settings\n", + " new_node_settings = {\n", + " node_id: attrs.evolve(\n", + " node_settings,\n", + " conservation_rules={\n", + " **node_settings.conservation_rules,\n", + " helicity_conservation: CONSERVATION_LAW_PRIORITIES[\n", + " helicity_conservation\n", + " ],\n", + " },\n", + " )\n", + " for node_id, node_settings in settings.interactions.items()\n", + " }\n", + " new_settings = MutableTransition(\n", + " settings.topology, dict(settings.states), new_node_settings\n", + " )\n", + " return QNProblemSet(\n", + " initial_facts=problem_set.initial_facts, solving_settings=new_settings\n", + " )\n", + "\n", + "\n", + "collection.problem_sets = {\n", + " strength: [add_helicity_conservation(problem_set) for problem_set in problem_sets]\n", + " for strength, problem_sets in collection.problem_sets.items()\n", + "}\n", + "qn_transitions = find_qn_transitions(collection)\n", + "len(qn_transitions)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "At the production node, the rule requires $|\\lambda_\\gamma - \\lambda_{f_2}| \\leq S_{J/\\psi} = 1$, which prunes the helicity combinations to:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "get_helicity_combinations(qn_transitions)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The resulting transitions render like any other quantum-number transition, now with the solved spin projections included:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "source = qrules.io.asmermaid(qn_transitions[0], render_node=True, markdown=True)\n", + "Markdown(source)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + ":::{tip}\n", + "The same mechanism works for any quantum number in {class}`.EdgeQuantumNumbers` and {class}`.NodeQuantumNumbers` and for any rule that consumes them — including rules of your own, as described in {doc}`/usage/conservation`.\n", + ":::" + ] + } + ], + "metadata": { + "colab": { + "toc_visible": true + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.14" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/docs/usage/visualize.ipynb b/docs/usage/visualize.ipynb index 2c0c3283..f5409388 100644 --- a/docs/usage/visualize.ipynb +++ b/docs/usage/visualize.ipynb @@ -65,9 +65,8 @@ "from IPython.display import Markdown\n", "\n", "import qrules\n", - "from qrules.particle import Parity, Spin\n", - "from qrules.topology import create_isobar_topologies, create_n_body_topology\n", - "from qrules.transition import State" + "from qrules.particle import Parity, Particle\n", + "from qrules.topology import create_isobar_topologies, create_n_body_topology" ] }, { @@ -195,7 +194,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "As noted in {ref}`usage/reaction:3. Find solutions`, the {attr}`~.ReactionInfo.transitions` contain all spin projection combinations (which is necessary for the {mod}`ampform` package). It is possible to convert all these solutions to Mermaid language with {func}`~.asmermaid`. To avoid visualizing all solutions, we just take a subset of the {attr}`~.ReactionInfo.transitions`:" + "As noted in {ref}`usage/reaction:3. Find solutions`, the {attr}`~.ReactionInfo.transitions` can contain quite a number of solutions, which differ for instance in their intermediate state or $LS$-couplings. It is possible to convert all these solutions to Mermaid language with {func}`~.asmermaid`. To avoid visualizing all solutions, we just take a subset of the {attr}`~.ReactionInfo.transitions`:" ] }, { @@ -204,7 +203,7 @@ "metadata": {}, "outputs": [], "source": [ - "some_transitions = reaction.transitions[::50][:3]\n", + "some_transitions = reaction.transitions[::30][:3]\n", "source = qrules.io.asmermaid(some_transitions)" ] }, @@ -238,7 +237,7 @@ "metadata": {}, "outputs": [], "source": [ - "qrules.io.write(reaction, \"decay_topologies_with_spin.mmd\")" + "qrules.io.write(reaction, \"decay_topologies.mmd\")" ] }, { @@ -254,7 +253,7 @@ "metadata": {}, "outputs": [], "source": [ - "qrules.io.write(reaction, \"decay_topologies_with_spin.gv\")" + "qrules.io.write(reaction, \"decay_topologies.gv\")" ] }, { @@ -268,7 +267,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Since this list of all possible spin projections {attr}`~.ReactionInfo.transitions` is rather long, it is often useful to collapse comparable graphs. First, {code}`collapse=\"spin\"` allows one to collapse (ignore) the spin projections (we again show a selection only):" + "If the list of {attr}`~.ReactionInfo.transitions` is rather long, it is often useful to collapse comparable graphs. First, {code}`collapse=\"spin\"` hides the interaction properties, such as the $LS$-couplings, and deduplicates the transitions (we again show a selection only):" ] }, { @@ -360,7 +359,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The {meth}`~.FrozenTransition.convert` method makes it possible to convert the types of its {attr}`~.FrozenTransition.states`. This for instance allows us to only render the spin states on in a {class}`.Transition`:\n", + "The {meth}`~.FrozenTransition.convert` method makes it possible to convert the types of its {attr}`~.FrozenTransition.states`. This for instance allows us to only render the spin magnitudes of the states in a {class}`.Transition`:\n", "\n", "::::{margin}\n", "\n", @@ -380,10 +379,9 @@ "outputs": [], "source": [ "spin_transitions = sorted({\n", - " t.convert(lambda s: Spin(s.particle.spin, s.spin_projection))\n", - " for t in reaction.transitions\n", + " t.convert(lambda particle: particle.spin) for t in reaction.transitions\n", "})\n", - "some_spin_transitions = spin_transitions[::67][:3]\n", + "some_spin_transitions = spin_transitions[:3]\n", "source = qrules.io.asmermaid(some_spin_transitions, markdown=True, render_node=True)\n", "Markdown(source)" ] @@ -392,7 +390,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Or any other properties of a {class}`.State`, such as masses or $J^{PC}(I^G)$ numbers:" + "Or any other properties of a {class}`.Particle`, such as masses or $J^{PC}(I^G)$ numbers:" ] }, { @@ -411,9 +409,9 @@ }, "outputs": [], "source": [ - "def render_mass(state: State, digits: int = 3) -> str:\n", - " mass = round(state.particle.mass, digits)\n", - " width = round(state.particle.width, digits)\n", + "def render_mass(particle: Particle, digits: int = 3) -> str:\n", + " mass = round(particle.mass, digits)\n", + " width = round(particle.width, digits)\n", " if width == 0:\n", " return str(mass)\n", " return f\"{mass}±{width}\"\n", @@ -451,8 +449,7 @@ "from fractions import Fraction\n", "\n", "\n", - "def render_jpc_ig(state: State) -> str:\n", - " particle = state.particle\n", + "def render_jpc_ig(particle: Particle) -> str:\n", " text = render_fraction(particle.spin)\n", " superscript = \"\"\n", " if particle.parity is not None:\n", diff --git a/src/qrules/__init__.py b/src/qrules/__init__.py index 1f2c7db3..b14b9d97 100644 --- a/src/qrules/__init__.py +++ b/src/qrules/__init__.py @@ -24,13 +24,7 @@ from qrules import io from qrules import system_control as system_control -from qrules.combinatorics import ( - InitialFacts, - StateDefinitionInput, - create_initial_facts, - to_state_definitions, -) -from qrules.combinatorics import StateDefinition as StateDefinition +from qrules.combinatorics import InitialFacts, create_initial_facts from qrules.conservation_rules import ( BaryonNumberConservation, BottomnessConservation, @@ -43,10 +37,8 @@ StrangenessConservation, TauLNConservation, c_parity_conservation, - clebsch_gordan_helicity_to_canonical, g_parity_conservation, gellmann_nishijima, - identical_particle_symmetrization, isospin_conservation, isospin_validity, parity_conservation, @@ -75,8 +67,8 @@ def check_reaction_violations( # ruff: ignore[complex-structure, too-many-positional-arguments] - initial_state: StateDefinitionInput | Sequence[StateDefinitionInput], - final_state: Sequence[StateDefinitionInput], + initial_state: str | Sequence[str], + final_state: Sequence[str], mass_conservation_factor: float | None = 3.0, particle_db: ParticleCollection | None = None, max_angular_momentum: int = 1, @@ -90,10 +82,8 @@ def check_reaction_violations( # ruff: ignore[complex-structure, too-many-posit correctly. Args: - initial_state: Shortform description of the initial state w/o spin - projections. - final_state: Shortform description of the final state w/o spin - projections. + initial_state: Particle names of the initial state. + final_state: Particle names of the final state. mass_conservation_factor: Factor with which the width is multiplied when checking for `.MassConservation`. Set to `None` in order to deactivate mass conservation. @@ -119,7 +109,8 @@ def check_reaction_violations( # ruff: ignore[complex-structure, too-many-posit .. seealso:: :ref:`usage:Check allowed reactions` """ - initial_state_definitions = to_state_definitions(initial_state) + if isinstance(initial_state, str): + initial_state = [initial_state] if particle_db is None: particle_db = load_pdg() @@ -153,7 +144,7 @@ def check_pure_edge_rules() -> None: } edge_check_result = _check_violations( - initial_facts[0], + initial_facts, node_rules={}, edge_rules=dict.fromkeys( topology.incoming_edge_ids | topology.outgoing_edge_ids, pure_edge_rules @@ -182,13 +173,13 @@ def check_edge_qn_conservation() -> set[frozenset[str]]: TauLNConservation(), isospin_conservation, } - if len(initial_state_definitions) == 1 and mass_conservation_factor is not None: + if len(initial_state) == 1 and mass_conservation_factor is not None: edge_qn_conservation_rules.add(MassConservation(mass_conservation_factor)) return { frozenset((x,)) for x in _check_violations( - initial_facts[0], + initial_facts, node_rules=dict.fromkeys(topology.nodes, edge_qn_conservation_rules), edge_rules={}, ).violated_node_rules[node_id] @@ -198,13 +189,13 @@ def check_edge_qn_conservation() -> set[frozenset[str]]: # since only certain spin rules require the isobar model. These spin rules # are not required here though. topology = create_n_body_topology( - number_of_initial_states=len(initial_state_definitions), + number_of_initial_states=len(initial_state), number_of_final_states=len(final_state), ) node_id = next(iter(topology.nodes)) initial_facts = create_initial_facts( - topology, initial_state_definitions, final_state, particle_db + topology, initial_state, final_state, particle_db ) check_pure_edge_rules() @@ -220,14 +211,13 @@ def check_edge_qn_conservation() -> set[frozenset[str]]: ) ] - initial_facts_list = [] - for ls_combi in ls_combinations: - for facts_combination in initial_facts: - new_facts = attrs.evolve( - facts_combination, - interactions={node_id: ls_combi}, - ) - initial_facts_list.append(new_facts) + initial_facts_list = [ + attrs.evolve( + initial_facts, + interactions={node_id: ls_combi}, + ) + for ls_combi in ls_combinations + ] # Verify each graph with the interaction rules. # Spin projection rules are skipped as they can only be checked reliably @@ -235,11 +225,9 @@ def check_edge_qn_conservation() -> set[frozenset[str]]: conservation_rules: dict[int, set[Rule]] = { node_id: { c_parity_conservation, - clebsch_gordan_helicity_to_canonical, g_parity_conservation, parity_conservation, spin_magnitude_conservation, - identical_particle_symmetrization, } } @@ -273,8 +261,8 @@ def check_edge_qn_conservation() -> set[frozenset[str]]: def generate_transitions( # ruff: ignore[too-many-positional-arguments] - initial_state: StateDefinitionInput | Sequence[StateDefinitionInput], - final_state: Sequence[StateDefinitionInput], + initial_state: str | Sequence[str], + final_state: Sequence[str], allowed_intermediate_particles: list[str] | None = None, allowed_interaction_types: str | Iterable[str] | None = None, formalism: SpinFormalism = "canonical-helicity", @@ -282,6 +270,7 @@ def generate_transitions( # ruff: ignore[too-many-positional-arguments] mass_conservation_factor: float | None = 3.0, max_angular_momentum: int = 2, max_spin_magnitude: float = 2, + ls_couplings: bool = True, topology_building: str = "isobar", number_of_threads: int | None = None, ) -> ReactionInfo: @@ -290,11 +279,7 @@ def generate_transitions( # ruff: ignore[too-many-positional-arguments] Serves as a facade to the `.StateTransitionManager` (see :doc:`/usage/reaction`). Arguments: - initial_state (list): A list of particle names in the initial - state. You can specify spin projections for these particles with a `tuple`, - e.g. :code:`("J/psi(1S)", [-1, 0, +1])`. If spin projections are not - specified, all projections are taken, so the example here would be - equivalent to :code:`"J/psi(1S)"`. + initial_state (list): A list of particle names in the initial state. final_state (list): Same as :code:`initial_state`, but for final state particles. @@ -326,6 +311,11 @@ def generate_transitions( # ruff: ignore[too-many-positional-arguments] max_spin_magnitude: Maximum spin magnitude over which to generate spins. + ls_couplings: If `False`, the solver does not enumerate + :math:`LS`-combinations; the couplings can be reconstructed from the spins + and parities of the solutions afterwards (see + `.create_interaction_settings`). + topology_building (str): Technique with which to build the `.Topology` instances. Allowed values are: @@ -355,9 +345,10 @@ def generate_transitions( # ruff: ignore[too-many-positional-arguments] >>> len(reaction.group_by_topology()) 3 """ - initial_state_definitions = to_state_definitions(initial_state) + if isinstance(initial_state, str): + initial_state = [initial_state] stm = StateTransitionManager( - initial_state=initial_state_definitions, + initial_state=initial_state, final_state=final_state, particle_db=particle_db, allowed_intermediate_particles=allowed_intermediate_particles, @@ -365,6 +356,7 @@ def generate_transitions( # ruff: ignore[too-many-positional-arguments] mass_conservation_factor=mass_conservation_factor, max_angular_momentum=max_angular_momentum, max_spin_magnitude=max_spin_magnitude, + ls_couplings=ls_couplings, topology_building=topology_building, number_of_threads=number_of_threads, ) diff --git a/src/qrules/argument_handling.py b/src/qrules/argument_handling.py index 2694de00..e8fedb0b 100644 --- a/src/qrules/argument_handling.py +++ b/src/qrules/argument_handling.py @@ -233,7 +233,9 @@ def wrapper(states_list: Sequence[Any]) -> list[Any]: def _build_all_arguments(checks: list[Callable]) -> Callable: def wrapper(*args: Any) -> list[Any]: - return [check(arg) for check, arg in zip(checks, args, strict=False) if arg] + # no filtering on falsy args: an empty property map (e.g. a node without LS + # couplings) must still be built into a rule argument + return [check(arg) for check, arg in zip(checks, args, strict=False)] return wrapper diff --git a/src/qrules/combinatorics.py b/src/qrules/combinatorics.py index 24f02680..d39fb580 100644 --- a/src/qrules/combinatorics.py +++ b/src/qrules/combinatorics.py @@ -9,58 +9,25 @@ import itertools from collections import OrderedDict -from collections.abc import Callable, Iterable, Mapping, Sequence from copy import deepcopy -from fractions import Fraction from typing import TYPE_CHECKING, Any -from qrules.argument_handling import Scalar -from qrules.particle import ParticleWithSpin -from qrules.quantum_numbers import InteractionProperties, arange +from qrules.particle import Particle +from qrules.quantum_numbers import InteractionProperties from qrules.topology import MutableTransition, Topology, get_originating_node_list if TYPE_CHECKING: + from collections.abc import Callable, Iterable, Mapping, Sequence + from typing_extensions import TypeIs from qrules.particle import ParticleCollection -StateWithSpins = tuple[str, Sequence[Fraction]] -StateDefinition = str | StateWithSpins -"""Particle name, optionally with a list of spin projections.""" -StateDefinitionInput = str | tuple[str, Sequence[Scalar]] -"""Input type for `StateDefinition` permitting also `int` and `float`""" -InitialFacts = MutableTransition[ParticleWithSpin, InteractionProperties] +InitialFacts = MutableTransition[Particle, InteractionProperties] """A `.Transition` with only initial and final state information.""" -def as_state_definition( - definition: StateDefinitionInput, -) -> StateDefinition: - if type(definition) is str: - return definition - if type(definition) is tuple: - name, state = definition - return name, [Fraction(x) for x in state] - msg = f"value has to be of type {StateDefinitionInput}, got {type(definition)}" - raise ValueError(msg) - - -def to_state_definitions( - state: StateDefinitionInput | Sequence[StateDefinitionInput], / -) -> list[StateDefinitionInput]: - """Wrap a single state definition in a `list`, or copy a sequence of them.""" - if _is_state_definition_input(state): - return [state] - return list(state) - - -def _is_state_definition_input(value: object, /) -> TypeIs[StateDefinitionInput]: - return isinstance(value, str) or ( - isinstance(value, tuple) and len(value) == 2 and isinstance(value[0], str) - ) - - class _KinematicRepresentation: # ruff: ignore[eq-without-hash] def __init__( self, @@ -222,114 +189,37 @@ def fill_groupings( def create_initial_facts( topology: Topology, - initial_state: Sequence[StateDefinitionInput], - final_state: Sequence[StateDefinitionInput], - particle_db: ParticleCollection, - expand_spin_projections: bool = True, -) -> list[InitialFacts]: - """Attach the initial and final states to the external edges of a `.Topology`. - - By default, one `.InitialFacts` is created for every combination of allowed spin - projections of the initial and final state. With - :code:`expand_spin_projections=False`, this Cartesian expansion is skipped and a - single `.InitialFacts` without spin projections is returned, for solving at the - :math:`J^{P(C)}` level (see `.strip_spin_projections`). - """ - states = __create_states_with_spin_projections( - sorted(topology.incoming_edge_ids) + sorted(topology.outgoing_edge_ids), - list(map(as_state_definition, initial_state)) - + list(map(as_state_definition, final_state)), - particle_db, - ) - if not expand_spin_projections: - projection_free_states = { - edge_id: (particle_db[name], None) for edge_id, (name, _) in states.items() - } - return [MutableTransition(topology, projection_free_states)] # type: ignore[arg-type] - spin_states = __generate_spin_combinations(states, particle_db) - return [MutableTransition(topology, states) for states in spin_states] - - -def __create_states_with_spin_projections( - edge_ids: Sequence[int], - state_definitions: Sequence[StateDefinition], + initial_state: Sequence[str], + final_state: Sequence[str], particle_db: ParticleCollection, -) -> dict[int, StateWithSpins]: - if len(edge_ids) != len(state_definitions): +) -> InitialFacts: + """Attach the initial and final state particles to the external edges of a `.Topology`.""" + edge_ids = sorted(topology.incoming_edge_ids) + sorted(topology.outgoing_edge_ids) + particle_names = list(initial_state) + list(final_state) + if len(edge_ids) != len(particle_names): msg = "Number of state definitions is not same as number of edge IDs" raise ValueError(msg) - states = __safe_set_spin_projections(state_definitions, particle_db) - return dict(zip(edge_ids, states, strict=True)) - - -def __safe_set_spin_projections( - state_definitions: Sequence[StateDefinition], - particle_db: ParticleCollection, -) -> Sequence[StateWithSpins]: - def fill_spin_projections(state: StateDefinition) -> StateWithSpins: - if isinstance(state, str): - particle_name = state - particle = particle_db[particle_name] - spin_projections = set(arange(-particle.spin, particle.spin + 1)) - if particle.mass == 0.0 and Fraction(0) in spin_projections: # ruff: ignore[float-equality-comparison] - spin_projections.remove(Fraction(0)) - return particle_name, sorted(spin_projections) - return state - - return [fill_spin_projections(state) for state in state_definitions] - - -def __generate_spin_combinations( - states_with_spin_projections: dict[int, StateWithSpins], - particle_db: ParticleCollection, -) -> list[dict[int, ParticleWithSpin]]: - def populate_edge_with_spin_projections( - permutation: dict[int, ParticleWithSpin], - edge_id: int, - state: StateWithSpins, - ) -> list[dict[int, ParticleWithSpin]]: - particle_name, spin_projections = state - particle = particle_db[particle_name] - new_permutations = [] - for projection in spin_projections: - temp_permutation = deepcopy(permutation) - temp_permutation.update({edge_id: (particle, projection)}) - new_permutations.append(temp_permutation) - return new_permutations - - initial_facts_permutations: list[dict[int, ParticleWithSpin]] = [{}] - for edge_id, state in states_with_spin_projections.items(): - temp_permutations = initial_facts_permutations - initial_facts_permutations = [] - for temp_permutation in temp_permutations: - initial_facts_permutations.extend( - populate_edge_with_spin_projections(temp_permutation, edge_id, state) - ) - - return initial_facts_permutations + states = { + edge_id: particle_db[name] + for edge_id, name in zip(edge_ids, particle_names, strict=True) + } + return MutableTransition(topology, states) def permutate_topology_kinematically( topology: Topology, - initial_state: Sequence[StateDefinitionInput] | Sequence[StateDefinition], - final_state: Sequence[StateDefinitionInput] | Sequence[StateDefinition], + initial_state: Sequence[str], + final_state: Sequence[str], final_state_groupings: list[list[list[str]]] | list[list[str]] | list[str] | None = None, ) -> list[Topology]: - def strip_spin(state: StateDefinitionInput) -> str: - if isinstance(state, tuple): - return state[0] - return state - edge_ids = sorted(topology.incoming_edge_ids) + sorted(topology.outgoing_edge_ids) states = list(initial_state) + list(final_state) return _generate_kinematic_permutations( topology, - particle_names={ - i: strip_spin(s) for i, s in zip(edge_ids, states, strict=True) - }, + particle_names=dict(zip(edge_ids, states, strict=True)), allowed_kinematic_groupings=__get_kinematic_groupings(final_state_groupings), ) @@ -395,7 +285,7 @@ def embed_in_list(some_list: list[Any]) -> list[list[Any]]: def match_external_edges( - graphs: list[MutableTransition[ParticleWithSpin, InteractionProperties]], + graphs: list[MutableTransition[Particle, InteractionProperties]], ) -> None: if not isinstance(graphs, list): msg = "graphs argument is not of type list" @@ -408,7 +298,7 @@ def match_external_edges( def _match_external_edge_ids( - graphs: list[MutableTransition[ParticleWithSpin, InteractionProperties]], + graphs: list[MutableTransition[Particle, InteractionProperties]], ref_graph_id: int, external_edge_getter_function: Callable[[MutableTransition], Iterable[int]], ) -> None: @@ -443,13 +333,13 @@ def _match_external_edge_ids( def __get_initial_state_edge_ids( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], + graph: MutableTransition[Particle, InteractionProperties], ) -> Iterable[int]: return graph.topology.incoming_edge_ids def __get_final_state_edge_ids( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], + graph: MutableTransition[Particle, InteractionProperties], ) -> Iterable[int]: return graph.topology.outgoing_edge_ids @@ -479,7 +369,7 @@ def perform_external_edge_identical_particle_combinatorics( def _external_edge_identical_particle_combinatorics( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], + graph: MutableTransition[Particle, InteractionProperties], external_edge_getter_function: Callable[[MutableTransition], Iterable[int]], ) -> list[MutableTransition]: new_graphs = [graph] @@ -536,7 +426,7 @@ def _calculate_swappings(id_mapping: dict[int, int]) -> OrderedDict: def _create_edge_id_particle_mapping( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], + graph: MutableTransition[Particle, InteractionProperties], edge_ids: Iterable[int], ) -> dict[int, str]: - return {i: graph.states[i][0].name for i in edge_ids} + return {i: graph.states[i].name for i in edge_ids} diff --git a/src/qrules/conservation_rules.py b/src/qrules/conservation_rules.py index 8d11cdf3..7d7555ea 100644 --- a/src/qrules/conservation_rules.py +++ b/src/qrules/conservation_rules.py @@ -45,7 +45,7 @@ """ import operator -from collections.abc import Callable +from collections.abc import Callable, Iterator from copy import deepcopy from fractions import Fraction from functools import reduce @@ -95,9 +95,11 @@ def _is_boson(spin_magnitude: Fraction) -> bool: return abs(spin_magnitude % 1) < 0.01 -def _is_particle_antiparticle_pair(pid1: int, pid2: int) -> bool: +def _is_particle_antiparticle_pair(pid1: int | None, pid2: int | None) -> bool: # we just check if the pid is opposite in sign # this is a requirement of the pid numbers of course + if pid1 is None or pid2 is None: + return False return pid1 == -pid2 @@ -300,15 +302,21 @@ def parity_conservation_helicity( @frozen class CParityEdgeInput: spin_magnitude: EdgeSpinMagnitude = field(converter=to_fraction) - pid: EdgePid = field(converter=int) + pid: EdgePid | None = field(converter=optional(int), default=None) c_parity: EdgeCParity | None = field(converter=optional(to_parity), default=None) @frozen class CParityNodeInput: - # These converters currently do not do anything, as "NewType"s do not have constructors - l_magnitude: NodeLMagnitude = field(converter=to_fraction) - s_magnitude: NodeSMagnitude = field(converter=to_fraction) + # optional, so that the rule remains executable for problem sets without LS + # couplings; the particle-antiparticle branch is then undetermined and is covered + # by the LS-free existence rules instead + l_magnitude: NodeLMagnitude | None = field( + converter=optional(to_fraction), default=None + ) + s_magnitude: NodeSMagnitude | None = field( + converter=optional(to_fraction), default=None + ) def c_parity_conservation( @@ -333,11 +341,13 @@ def _get_c_parity_multiparticle( if len(part_qns) == 2: # ruff: ignore[collapsible-if] if _is_particle_antiparticle_pair(part_qns[0].pid, part_qns[1].pid): ang_mom = interaction_qns.l_magnitude + if ang_mom is None: + return None # if boson if _is_boson(part_qns[0].spin_magnitude): return (-1) ** int(ang_mom) - coupled_spin = Fraction(interaction_qns.s_magnitude) - if isinstance(coupled_spin, int) or coupled_spin.denominator == 1: + coupled_spin = interaction_qns.s_magnitude + if coupled_spin is not None and coupled_spin.denominator == 1: return (-1) ** int(ang_mom + coupled_spin) return None @@ -356,14 +366,21 @@ def _get_c_parity_multiparticle( class GParityEdgeInput: isospin_magnitude: EdgeIsospinMagnitude = field(converter=to_fraction) spin_magnitude: EdgeSpinMagnitude = field(converter=to_fraction) - pid: EdgePid = field(converter=int) + pid: EdgePid | None = field(converter=optional(int), default=None) g_parity: EdgeGParity | None = field(converter=optional(to_parity), default=None) @frozen class GParityNodeInput: - l_magnitude: NodeLMagnitude = field(converter=to_fraction) - s_magnitude: NodeSMagnitude = field(converter=to_fraction) + # optional, so that the rule remains executable for problem sets without LS + # couplings; the particle-antiparticle branch is then undetermined and is covered + # by the LS-free existence rules instead + l_magnitude: NodeLMagnitude | None = field( + converter=optional(to_fraction), default=None + ) + s_magnitude: NodeSMagnitude | None = field( + converter=optional(to_fraction), default=None + ) @frozen @@ -440,12 +457,14 @@ def check_multistate_g_parity( double_state_qns[0].pid, double_state_qns[1].pid ): ang_mom = interaction_qns.l_magnitude + if ang_mom is None: + return None if isinstance(isospin, int) or isospin.denominator == 1: # 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.denominator == 1: + if coupled_spin is not None and coupled_spin.denominator == 1: return (-1) ** int(ang_mom + coupled_spin + isospin) return None @@ -605,6 +624,13 @@ class SpinMagnitudeNodeFacts(TypedDict): s_magnitude: Fraction +class SpinParityFacts(TypedDict): + """Facts required by `SpinParityCoupling`; a subset of `.EdgeFacts`.""" + + spin_magnitude: Fraction + parity: int + + def ls_spin_validity(spin_input: SpinNodeInput) -> bool: r"""Check for valid isospin magnitude and projection.""" return _check_spin_valid( @@ -850,6 +876,337 @@ def spin_magnitude_conservation( ) +def _couple_spins(j1: Fraction, j2: Fraction) -> list[Fraction]: + """List the magnitudes that two spins can couple to. + + >>> _couple_spins(Fraction(1, 2), Fraction(1)) + [Fraction(1, 2), Fraction(3, 2)] + """ + return [ + Fraction(x, 2) for x in range(int(2 * abs(j1 - j2)), int(2 * (j1 + j2)) + 1, 2) + ] + + +def _split_isobar_node( + ingoing: list, outgoing: list +) -> tuple[Any, list[Any]] | tuple[None, None]: + """Split node facts into the single-state side and the two-state side.""" + if (len(ingoing), len(outgoing)) == (1, 2): + return ingoing[0], outgoing + if (len(ingoing), len(outgoing)) == (2, 1): + return outgoing[0], ingoing + return None, None + + +def _iter_ls_couplings( + single_spin: Fraction, + pair_spins: list[Fraction], + max_angular_momentum: int, +) -> Iterator[tuple[int, Fraction]]: + for coupled_spin in _couple_spins(*pair_spins): + for ang_mom in range(max_angular_momentum + 1): + # membership in the coupling series, not just the triangle rule: J must + # differ from both L and S by an integer + if single_spin in _couple_spins(Fraction(ang_mom), coupled_spin): + yield ang_mom, coupled_spin + + +class SpinCoupling(EdgeQNConservationRule): + r"""Check that some :math:`(L, S)` combination couples the spin magnitudes. + + The :math:`LS`-free counterpart of `spin_magnitude_conservation` for problem sets + without `~.NodeQuantumNumbers.l_magnitude` and `~.NodeQuantumNumbers.s_magnitude` + domains: instead of checking specific :math:`(L, S)` values, check whether *any* + combination with :math:`L \leq L_\mathrm{max}` satisfies + :math:`S \in \{|S_1-S_2|, \dots, S_1+S_2\}` and :math:`J \in \{|L-S|, \dots, + L+S\}`. + + >>> from fractions import Fraction + >>> rule = SpinCoupling(max_angular_momentum=0) + >>> vector = {"spin_magnitude": Fraction(1)} + >>> tensor = {"spin_magnitude": Fraction(2)} + >>> rule([vector], [vector, tensor]) + True + >>> scalar = {"spin_magnitude": Fraction(0)} + >>> rule([scalar], [scalar, vector]) + False + >>> SpinCoupling(max_angular_momentum=1)([scalar], [scalar, vector]) + True + + A half-integer pair can never couple to an integer spin, no matter the angular + momentum: + + >>> fermion = {"spin_magnitude": Fraction(1, 2)} + >>> SpinCoupling(max_angular_momentum=10)([vector], [vector, fermion]) + False + """ + + def __init__(self, max_angular_momentum: int) -> None: + self.__max_angular_momentum = max_angular_momentum + + def __call__( + self, + ingoing_spin_magnitudes: list[SpinMagnitudeFacts], + outgoing_spin_magnitudes: list[SpinMagnitudeFacts], + ) -> bool: + single, pair = _split_isobar_node( + ingoing_spin_magnitudes, outgoing_spin_magnitudes + ) + if single is None or pair is None: + return ( + sum( + (float(x["spin_magnitude"]) for x in ingoing_spin_magnitudes), 0.0 + ).is_integer() + == sum( + (float(x["spin_magnitude"]) for x in outgoing_spin_magnitudes), 0.0 + ).is_integer() + ) + couplings = _iter_ls_couplings( + single["spin_magnitude"], + [x["spin_magnitude"] for x in pair], + self.__max_angular_momentum, + ) + return any(couplings) + + +class SpinParityCoupling(EdgeQNConservationRule): + r"""Check that some :math:`(L, S)` combination couples the spins and parities. + + The :math:`LS`-free counterpart of `spin_magnitude_conservation` plus + `parity_conservation`: in addition to the couplings of `SpinCoupling`, the + angular momentum must satisfy :math:`P_{in} = P_{out} \cdot (-1)^L`. + + >>> from fractions import Fraction + >>> rule = SpinParityCoupling(max_angular_momentum=2) + >>> jpsi = {"spin_magnitude": Fraction(1), "parity": -1} + >>> gamma = {"spin_magnitude": Fraction(1), "parity": -1} + >>> f2 = {"spin_magnitude": Fraction(2), "parity": +1} + >>> rule([jpsi], [gamma, f2]) + True + >>> eta = {"spin_magnitude": Fraction(0), "parity": -1} + >>> pion = {"spin_magnitude": Fraction(0), "parity": -1} + >>> rule([eta], [pion, pion]) + False + """ + + def __init__(self, max_angular_momentum: int) -> None: + self.__max_angular_momentum = max_angular_momentum + + def __call__( + self, + ingoing_spin_parities: list[SpinParityFacts], + outgoing_spin_parities: list[SpinParityFacts], + ) -> bool: + single, pair = _split_isobar_node(ingoing_spin_parities, outgoing_spin_parities) + if single is None or pair is None: + return True + parities = [single["parity"], *(x["parity"] for x in pair)] + if any(parity is None for parity in parities): + return False + single_parity, *pair_parities = parities + pair_parity = pair_parities[0] * pair_parities[1] + return any( + single_parity == pair_parity * (-1) ** ang_mom + for ang_mom, _ in _iter_ls_couplings( + single["spin_magnitude"], + [x["spin_magnitude"] for x in pair], + self.__max_angular_momentum, + ) + ) + + +@frozen +class CParityCouplingEdgeInput: + spin_magnitude: EdgeSpinMagnitude = field(converter=to_fraction) + parity: EdgeParity = field(converter=to_parity) + pid: EdgePid | None = field(converter=optional(int), default=None) + c_parity: EdgeCParity | None = field(converter=optional(to_parity), default=None) + + +@frozen +class GParityCouplingEdgeInput: + isospin_magnitude: EdgeIsospinMagnitude = field(converter=to_fraction) + spin_magnitude: EdgeSpinMagnitude = field(converter=to_fraction) + parity: EdgeParity = field(converter=to_parity) + pid: EdgePid | None = field(converter=optional(int), default=None) + g_parity: EdgeGParity | None = field(converter=optional(to_parity), default=None) + + +def _iter_parity_conserving_couplings( + single: CParityCouplingEdgeInput | GParityCouplingEdgeInput, + pair: list[CParityCouplingEdgeInput] | list[GParityCouplingEdgeInput], + max_angular_momentum: int, +) -> Iterator[tuple[int, Fraction]]: + """Iterate over the :math:`(L, S)` couplings that also conserve parity.""" + pair_parity = pair[0].parity.value * pair[1].parity.value + for ang_mom, coupled_spin in _iter_ls_couplings( + single.spin_magnitude, + [x.spin_magnitude for x in pair], + max_angular_momentum, + ): + if single.parity.value == pair_parity * (-1) ** ang_mom: + yield ang_mom, coupled_spin + + +def _get_pair_c_parity( + pair: list[CParityCouplingEdgeInput] | list[GParityCouplingEdgeInput], + ang_mom: int, + coupled_spin: Fraction, +) -> int | None: + """Composite :math:`C`-parity of a particle-antiparticle pair, if defined.""" + if not _is_particle_antiparticle_pair(pair[0].pid, pair[1].pid): + return None + if _is_boson(pair[0].spin_magnitude): + return (-1) ** ang_mom + if coupled_spin.denominator == 1: + return (-1) ** int(ang_mom + coupled_spin) + return None + + +class CParityCoupling(EdgeQNConservationRule): + r"""Check that some :math:`(L, S)` combination conserves :math:`C`-parity. + + The :math:`LS`-free counterpart of `c_parity_conservation`. If all states carry a + :math:`C`-parity, the rule checks :math:`C_{in} = C_{out}` directly. For a + particle-antiparticle pair without individual :math:`C`-parities, the composite + :math:`C = (-1)^L` (bosons) or :math:`C = (-1)^{L+S}` (fermions) depends on the + coupling, so the rule checks whether *some* :math:`(L, S)` combination with + :math:`L \leq L_\mathrm{max}` that conserves parity also conserves + :math:`C`-parity. + + >>> from fractions import Fraction + >>> rule = CParityCoupling(max_angular_momentum=1) + >>> pi_plus = CParityCouplingEdgeInput(spin_magnitude=0, parity=-1, pid=211) + >>> pi_minus = CParityCouplingEdgeInput(spin_magnitude=0, parity=-1, pid=-211) + >>> rho = CParityCouplingEdgeInput( + ... spin_magnitude=1, parity=-1, pid=113, c_parity=-1 + ... ) + >>> rule([rho], [pi_plus, pi_minus]) + True + >>> CParityCoupling(max_angular_momentum=0)([rho], [pi_plus, pi_minus]) + False + + Parity conservation requires an odd :math:`L` here, so a :math:`1^{-+}` state + cannot decay to :math:`\pi^+\pi^-`, even though parity and :math:`C`-parity each + allow some coupling on their own: + + >>> exotic = CParityCouplingEdgeInput( + ... spin_magnitude=1, parity=-1, pid=1, c_parity=+1 + ... ) + >>> CParityCoupling(max_angular_momentum=3)([exotic], [pi_plus, pi_minus]) + False + """ + + def __init__(self, max_angular_momentum: int) -> None: + self.__max_angular_momentum = max_angular_momentum + + def __call__( + self, + ingoing_edge_qns: list[CParityCouplingEdgeInput], + outgoing_edge_qns: list[CParityCouplingEdgeInput], + ) -> bool: + c_parity_in = _multiply_c_parities(ingoing_edge_qns) + c_parity_out = _multiply_c_parities(outgoing_edge_qns) + if c_parity_in is not None and c_parity_out is not None: + return c_parity_in == c_parity_out + single, pair = _split_isobar_node(ingoing_edge_qns, outgoing_edge_qns) + if single is None or pair is None or single.c_parity is None: + return True + if not _is_particle_antiparticle_pair(pair[0].pid, pair[1].pid): + return True + return any( + _get_pair_c_parity(pair, ang_mom, coupled_spin) == single.c_parity.value + for ang_mom, coupled_spin in _iter_parity_conserving_couplings( + single, pair, self.__max_angular_momentum + ) + ) + + +def _multiply_c_parities(states: list[CParityCouplingEdgeInput]) -> int | None: + if any(x.c_parity is None for x in states): + return None + return reduce(operator.mul, (x.c_parity.value for x in states if x.c_parity), 1) + + +class GParityCoupling(EdgeQNConservationRule): + r"""Check that some :math:`(L, S)` combination conserves :math:`G`-parity. + + The :math:`LS`-free counterpart of `g_parity_conservation`. If all states carry a + :math:`G`-parity, the rule checks :math:`G_{in} = G_{out}` directly. For a + particle-antiparticle pair without individual :math:`G`-parities, the composite + :math:`G = C \cdot (-1)^I` depends on the coupling through :math:`C`, so the rule + checks whether *some* :math:`(L, S)` combination with :math:`L \leq + L_\mathrm{max}` that conserves parity also conserves :math:`G`-parity. + + >>> from fractions import Fraction + >>> rule = GParityCoupling(max_angular_momentum=1) + >>> pi_plus = GParityCouplingEdgeInput( + ... isospin_magnitude=1, spin_magnitude=0, parity=-1, pid=211 + ... ) + >>> pi_minus = GParityCouplingEdgeInput( + ... isospin_magnitude=1, spin_magnitude=0, parity=-1, pid=-211 + ... ) + >>> rho = GParityCouplingEdgeInput( + ... isospin_magnitude=1, spin_magnitude=1, parity=-1, pid=113, g_parity=+1 + ... ) + >>> rule([rho], [pi_plus, pi_minus]) + True + >>> GParityCoupling(max_angular_momentum=0)([rho], [pi_plus, pi_minus]) + False + >>> f0 = GParityCouplingEdgeInput( + ... isospin_magnitude=0, spin_magnitude=0, parity=+1, pid=9010221, g_parity=+1 + ... ) + >>> GParityCoupling(max_angular_momentum=0)([f0], [pi_plus, pi_minus]) + True + """ + + def __init__(self, max_angular_momentum: int) -> None: + self.__max_angular_momentum = max_angular_momentum + + def __call__( + self, + ingoing_edge_qns: list[GParityCouplingEdgeInput], + outgoing_edge_qns: list[GParityCouplingEdgeInput], + ) -> bool: + g_parity_in = _multiply_g_parities(ingoing_edge_qns) + g_parity_out = _multiply_g_parities(outgoing_edge_qns) + if g_parity_in is not None and g_parity_out is not None: + return g_parity_in == g_parity_out + single, pair = _split_isobar_node(ingoing_edge_qns, outgoing_edge_qns) + if single is None or pair is None or single.g_parity is None: + return True + if not _is_particle_antiparticle_pair(pair[0].pid, pair[1].pid): + return True + isospin = single.isospin_magnitude + if isospin.denominator != 1: + return True + return any( + _get_pair_g_parity(pair, ang_mom, coupled_spin, isospin) + == single.g_parity.value + for ang_mom, coupled_spin in _iter_parity_conserving_couplings( + single, pair, self.__max_angular_momentum + ) + ) + + +def _get_pair_g_parity( + pair: list[GParityCouplingEdgeInput], + ang_mom: int, + coupled_spin: Fraction, + isospin: Fraction, +) -> int | None: + c_parity = _get_pair_c_parity(pair, ang_mom, coupled_spin) + if c_parity is None: + return None + return c_parity * (-1) ** int(isospin) + + +def _multiply_g_parities(states: list[GParityCouplingEdgeInput]) -> int | None: + if any(x.g_parity is None for x in states): + return None + return reduce(operator.mul, (x.g_parity.value for x in states if x.g_parity), 1) + + def clebsch_gordan_helicity_to_canonical( ingoing_spins: list[SpinEdgeInput], outgoing_spins: list[SpinEdgeInput], diff --git a/src/qrules/io/__init__.py b/src/qrules/io/__init__.py index 6e176d33..de70a709 100644 --- a/src/qrules/io/__init__.py +++ b/src/qrules/io/__init__.py @@ -78,11 +78,12 @@ def asdot( Args: instance: the input `object` that is to be rendered as DOT (graphviz) language. - collapse: Optionally combine transitions. With ``"spin"``, transitions that - differ only in their spin projections are combined. With ``"topology"``, - transitions with equivalent kinematic topologies are grouped and all allowed - particles on each edge are collected. Node properties cannot be rendered - when collapsing by topology. + collapse: Optionally combine transitions. With ``"spin"``, the interaction + properties (such as :math:`LS`-couplings) of the transitions are hidden and + the transitions are deduplicated, leaving only unique graphs with `.Particle` + names on their edges. With ``"topology"``, transitions with equivalent + kinematic topologies are grouped and all allowed particles on each edge are + collected. Node properties cannot be rendered when collapsing by topology. render_node: Whether or not to render node ID (in the case of a `.Topology`) and/or node properties (in the case of a `.MutableTransition`). Meaning of @@ -146,12 +147,14 @@ def asmermaid( instance: the input `object` that is to be rendered as Mermaid flowchart source. - collapse: Optionally combine transitions. With ``"spin"``, transitions that - differ only in their spin projections are combined. With ``"topology"``, - transitions with equivalent kinematic topologies are grouped and all allowed - particles on each edge are collected. Node properties cannot be rendered - when collapsing by topology. With LaTeX rendering, particle lists longer - than six entries are arranged in columns of at most six rows. + collapse: Optionally combine transitions. With ``"spin"``, the interaction + properties (such as :math:`LS`-couplings) of the transitions are hidden and + the transitions are deduplicated, leaving only unique graphs with `.Particle` + names on their edges. With ``"topology"``, transitions with equivalent + kinematic topologies are grouped and all allowed particles on each edge are + collected. Node properties cannot be rendered when collapsing by topology. + With LaTeX rendering, particle lists longer than six entries are arranged in + columns of at most six rows. render_node: Whether or not to render node ID (in the case of a `.Topology`) and/or node properties (in the case of a `.MutableTransition`). Meaning of diff --git a/src/qrules/io/_dict.py b/src/qrules/io/_dict.py index 205a58a0..1d5ff76b 100644 --- a/src/qrules/io/_dict.py +++ b/src/qrules/io/_dict.py @@ -13,7 +13,7 @@ from qrules.particle import Parity, Particle, ParticleCollection, Spin from qrules.quantum_numbers import InteractionProperties from qrules.topology import Edge, FrozenTransition, Topology -from qrules.transition import ReactionInfo, State +from qrules.transition import ReactionInfo def from_particle_collection(particles: ParticleCollection) -> dict: @@ -44,7 +44,9 @@ def _value_serializer(inst: type, field: attrs.Attribute, value: Any) -> Any: # if isinstance(value, (frozenset, set)): return sorted(value) if isinstance(value, abc.Mapping): - if all(isinstance(p, Particle) for p in value.values()): + if not isinstance(inst, FrozenTransition) and all( + isinstance(p, Particle) for p in value.values() + ): return {k: v.name for k, v in value.items()} return { (key.__name__ if callable(key) else key): _value_serializer( @@ -52,7 +54,7 @@ def _value_serializer(inst: type, field: attrs.Attribute, value: Any) -> Any: # ) for key, item in value.items() } - if not isinstance(inst, (ReactionInfo, State, FrozenTransition)): # ruff: ignore[collapsible-if] + if not isinstance(inst, (ReactionInfo, FrozenTransition)): # ruff: ignore[collapsible-if] if isinstance(value, Particle): return value.name if isinstance(value, Parity): @@ -95,12 +97,12 @@ def build_reaction_info(definition: dict) -> ReactionInfo: def build_transition( definition: dict, -) -> FrozenTransition[State, InteractionProperties]: +) -> FrozenTransition[Particle, InteractionProperties]: topology = build_topology(definition["topology"]) states_def: dict[int, dict] = definition["states"] - states: dict[int, State] = {} + states: dict[int, Particle] = {} for i, edge_def in states_def.items(): - states[int(i)] = build_state(edge_def) + states[int(i)] = build_particle(edge_def) interactions_def: dict[int, dict] = definition["interactions"] interactions = { int(i): InteractionProperties(**node_def) @@ -109,18 +111,6 @@ def build_transition( return FrozenTransition(topology, states, interactions) -def build_state(definition: Any) -> State: - if isinstance(definition, (list, tuple)) and len(definition) == 2: - particle = build_particle(definition[0]) - spin_projection = float(definition[1]) - return State(particle, spin_projection) - if isinstance(definition, dict): - particle = build_particle(definition["particle"]) - spin_projection = float(definition["spin_projection"]) - return State(particle, spin_projection) - raise NotImplementedError - - def build_topology(definition: dict) -> Topology: nodes = definition["nodes"] edges_def: dict[int, dict] = definition["edges"] diff --git a/src/qrules/io/_labels.py b/src/qrules/io/_labels.py index 2e306d66..6f9e739e 100644 --- a/src/qrules/io/_labels.py +++ b/src/qrules/io/_labels.py @@ -11,7 +11,7 @@ import attrs from attrs import frozen -from qrules.particle import Particle, ParticleWithSpin, Spin, _render_fraction +from qrules.particle import Particle, Spin, _render_fraction from qrules.quantum_numbers import EdgeQuantumNumbers, InteractionProperties from qrules.solving import ( EdgeSettings, @@ -26,7 +26,7 @@ Topology, Transition, ) -from qrules.transition import ProblemSet, State +from qrules.transition import ProblemSet if TYPE_CHECKING: from collections.abc import Callable, Iterable, Mapping @@ -99,9 +99,10 @@ def select_transitions( ) -> list[Any]: """Reduce a collection of transitions to the graphs that are worth rendering. - The ``collapse`` mode is the printer attribute :code:`collapse`. Spin projections - are only stripped from the interaction nodes if those nodes are not rendered, and - topologies cannot be collapsed at all while node properties are rendered. + The ``collapse`` mode is the printer attribute :code:`collapse`. Interaction + properties are only stripped from the interaction nodes if those nodes are not + rendered, and topologies cannot be collapsed at all while node properties are + rendered. """ if collapse is None: return list(graphs) @@ -134,9 +135,7 @@ def create_edge_label( if isinstance(graph, (ProblemSet, QNProblemSet)): edge_setting = graph.solving_settings.states.get(edge_id) initial_fact = graph.initial_facts.states.get(edge_id) - edge_property: EdgeSettings | GraphEdgePropertyMap | ParticleWithSpin | None = ( - None - ) + edge_property: EdgeSettings | GraphEdgePropertyMap | Particle | None = None if edge_setting: edge_property = edge_setting if initial_fact: @@ -210,7 +209,6 @@ def assignment(self, key: str, value: str, *, compact: bool = False) -> str: ... def membership(self, key: str, domain: str) -> str: ... def particle(self, name: str, latex: str | None) -> str: ... def spin(self, magnitude: str, projection: str) -> str: ... - def state(self, particle: str, projection: str) -> str: ... def superscript(self, base: str, exponent: str) -> str: ... @@ -253,10 +251,6 @@ def particle(name: str, latex: str | None) -> str: def spin(magnitude: str, projection: str) -> str: return f"|{magnitude},{projection}⟩" - @staticmethod - def state(particle: str, projection: str) -> str: - return f"{particle}[{projection}]" - @staticmethod def superscript(base: str, exponent: str) -> str: return base + exponent.translate(_SUPERSCRIPT_SIGNS) @@ -302,10 +296,6 @@ def particle(name: str, latex: str | None) -> str: def spin(magnitude: str, projection: str) -> str: return Rf"\left|{magnitude},{projection}\right\rangle" - @staticmethod - def state(particle: str, projection: str) -> str: - return Rf"{particle}\left[{projection}\right]" - @staticmethod def superscript(base: str, exponent: str) -> str: if not exponent: @@ -546,22 +536,6 @@ def __render_spin(spin: Spin, formatter: _LabelFormatter) -> str: return formatter.spin(spin_magnitude, spin_projection) -@as_string.register(State) -def _(state: State) -> str: - return __render_state(state, _PLAIN_FORMATTER) - - -@as_latex.register(State) -def _(state: State) -> str: - return __render_state(state, _LATEX_FORMATTER) - - -def __render_state(state: State, formatter: _LabelFormatter) -> str: - particle = formatter.render(state.particle) - spin_projection = formatter.fraction(state.spin_projection, plusminus=True) - return formatter.state(particle, spin_projection) - - @frozen class QuantumNumberSignature: """PDG-style :math:`I^G(J^{PC})` summary of a quantum-number property map. @@ -684,11 +658,8 @@ def _(obj: tuple) -> str: def __render_tuple(obj: tuple, formatter: _LabelFormatter) -> str: - if len(obj) == 2: - if isinstance(obj[0], Particle) and isinstance(obj[1], (Fraction, float, int)): - return __render_state(State(*obj), formatter) - if all(isinstance(o, (Fraction, float, int)) for o in obj): - return __render_spin(Spin(*obj), formatter) + if len(obj) == 2 and all(isinstance(o, (Fraction, float, int)) for o in obj): + return __render_spin(Spin(*obj), formatter) rendered_items = [formatter.render(item) for item in obj] if ( formatter is _LATEX_FORMATTER @@ -715,12 +686,12 @@ def _render_latex_columns(items: list[str]) -> str: def get_particle_graphs( - graphs: Iterable[Transition[ParticleWithSpin, InteractionProperties]], + graphs: Iterable[Transition[Particle, InteractionProperties]], ) -> list[FrozenTransition[Particle, None]]: - """Strip `list` of `.Transition` s of the spin projections. + """Strip `list` of `.Transition` s of their interaction properties. - Extract a `list` of `.Transition` instances with only `.Particle` instances on the - edges. + Extract a `list` of unique `.Transition` instances with only `.Particle` instances + on the edges and no interaction properties. .. seealso:: :doc:`/usage/visualize` """ @@ -753,10 +724,8 @@ def strip_projections( def __to_particle(state: Any) -> Particle: - if isinstance(state, State): - return state.particle - if isinstance(state, tuple) and len(state) == 2: - return state[0] + if isinstance(state, Particle): + return state msg = f"Cannot extract a particle from type {type(state).__name__}" raise NotImplementedError(msg) @@ -795,8 +764,6 @@ def collapse_graphs( def _strip_properties(state: Any) -> Any: - if isinstance(state, State): - return state.particle if isinstance(state, abc.Mapping): return FrozenDict(state) return state @@ -825,8 +792,6 @@ def _summarize_property_maps(states: Iterable[Any]) -> set[Any]: def _sorting_key(obj: Any) -> Any: - if isinstance(obj, State): - return obj.particle.name if isinstance(obj, QuantumNumberSignature): return as_string(obj) if isinstance(obj, str): diff --git a/src/qrules/particle/__init__.py b/src/qrules/particle/__init__.py index 228fd932..9b7a0d0d 100644 --- a/src/qrules/particle/__init__.py +++ b/src/qrules/particle/__init__.py @@ -243,10 +243,6 @@ def _get_name_root(name: str) -> str: return re.sub(r"[\*\+\-~\d']", "", name_root) -ParticleWithSpin = tuple[Particle, Fraction] -"""A particle and its spin projection.""" - - class ParticleCollection(abc.MutableSet): # ruff: ignore[eq-without-hash] """Searchable collection of immutable `.Particle` instances.""" diff --git a/src/qrules/settings.py b/src/qrules/settings.py index be252a30..d98f7cd8 100644 --- a/src/qrules/settings.py +++ b/src/qrules/settings.py @@ -22,13 +22,16 @@ BottomnessConservation, ChargeConservation, CharmConservation, + CParityCoupling, ElectronLNConservation, + GParityCoupling, MassConservation, MuonLNConservation, + SpinCoupling, + SpinParityCoupling, StrangenessConservation, TauLNConservation, c_parity_conservation, - clebsch_gordan_helicity_to_canonical, g_parity_conservation, g_parity_validity, gellmann_nishijima, @@ -53,7 +56,6 @@ from qrules.argument_handling import Rule, RuleKey from qrules.particle import Particle, ParticleCollection - from qrules.transition import SpinFormalism __QRULES_PATH = dirname(realpath(__file__)) ADDITIONAL_PARTICLES_DEFINITIONS_PATH: str = join( @@ -72,6 +74,10 @@ ChargeConservation: 100, spin_conservation: 8, spin_magnitude_conservation: 8, + SpinCoupling: 8, + SpinParityCoupling: 6, + CParityCoupling: 5, + GParityCoupling: 3, parity_conservation: 6, c_parity_conservation: 5, g_parity_conservation: 3, @@ -139,74 +145,64 @@ def from_str(description: str) -> InteractionType: def create_interaction_settings( # ruff: ignore[too-many-positional-arguments] - formalism: SpinFormalism, particle_db: ParticleCollection, nbody_topology: bool = False, mass_conservation_factor: float | None = 3.0, max_angular_momentum: int = 2, max_spin_magnitude: float = 2, + ls_couplings: bool = True, ) -> dict[InteractionType, tuple[EdgeSettings, NodeSettings]]: - """Create a container that holds the settings for `.InteractionType`.""" - formalism_edge_settings = EdgeSettings( + """Create a container that holds the settings for `.InteractionType`. + + With :code:`ls_couplings=False`, the settings declare no + `~.NodeQuantumNumbers.l_magnitude` and `~.NodeQuantumNumbers.s_magnitude` + domains, so the solver does not enumerate :math:`LS`-combinations. The + corresponding constraints are imposed by the existence rules `.SpinCoupling`, + `.SpinParityCoupling`, `.CParityCoupling`, and `.GParityCoupling` instead, which + only check whether *some* :math:`(L, S)` combination up to + :code:`max_angular_momentum` exists. The allowed combinations can be reconstructed + from the spins and parities of the solutions afterwards. + """ + default_edge_settings = EdgeSettings( conservation_rules=_with_priorities( { isospin_validity, g_parity_validity, gellmann_nishijima, - spin_validity, }, EDGE_RULE_PRIORITIES, ), qn_domains=_create_domains(particle_db), ) - formalism_node_settings = NodeSettings() - - angular_momentum_domain = __get_ang_mom_magnitudes( - nbody_topology, max_angular_momentum - ) - spin_magnitude_domain = __get_spin_magnitudes(nbody_topology, max_spin_magnitude) - if "helicity" in formalism: - formalism_node_settings.conservation_rules = _with_priorities({ - spin_magnitude_conservation, - helicity_conservation, - }) - formalism_node_settings.qn_domains = { - NodeQN.l_magnitude: angular_momentum_domain, - NodeQN.s_magnitude: spin_magnitude_domain, - } - elif formalism == "canonical": - formalism_node_settings.conservation_rules = _with_priorities({ - spin_magnitude_conservation - }) - if nbody_topology: - formalism_node_settings.conservation_rules = _with_priorities({ - spin_conservation, - ls_spin_validity, - }) - formalism_node_settings.qn_domains = { - NodeQN.l_magnitude: angular_momentum_domain, - NodeQN.l_projection: __extend_negative(angular_momentum_domain), - NodeQN.s_magnitude: spin_magnitude_domain, - NodeQN.s_projection: __extend_negative(spin_magnitude_domain), - } - if formalism == "canonical-helicity": - formalism_node_settings.conservation_rules.update( - _with_priorities({ - clebsch_gordan_helicity_to_canonical, - ls_spin_validity, - }) + if ls_couplings: + default_node_settings = NodeSettings( + conservation_rules=_with_priorities({spin_magnitude_conservation}), + qn_domains={ + NodeQN.l_magnitude: __get_ang_mom_magnitudes( + nbody_topology, max_angular_momentum + ), + NodeQN.s_magnitude: __get_spin_magnitudes( + nbody_topology, max_spin_magnitude + ), + }, ) - formalism_node_settings.qn_domains.update({ - NodeQN.l_projection: [0], - NodeQN.s_projection: __extend_negative(spin_magnitude_domain), - }) + parity_rule: Rule = parity_conservation + c_parity_rule: Rule = c_parity_conservation + g_parity_rule: Rule = g_parity_conservation + else: + default_node_settings = NodeSettings( + conservation_rules=_with_priorities({SpinCoupling(max_angular_momentum)}), + ) + parity_rule = SpinParityCoupling(max_angular_momentum) + c_parity_rule = CParityCoupling(max_angular_momentum) + g_parity_rule = GParityCoupling(max_angular_momentum) if mass_conservation_factor is not None: - formalism_node_settings.conservation_rules.update( + default_node_settings.conservation_rules.update( _with_priorities({MassConservation(mass_conservation_factor)}) ) interaction_type_settings = {} - weak_node_settings = deepcopy(formalism_node_settings) + weak_node_settings = deepcopy(default_node_settings) weak_node_settings.conservation_rules.update( _with_priorities([ ChargeConservation(), @@ -214,11 +210,10 @@ def create_interaction_settings( # ruff: ignore[too-many-positional-arguments] MuonLNConservation(), TauLNConservation(), BaryonNumberConservation(), - identical_particle_symmetrization, ]) ) weak_node_settings.interaction_strength = 10 ** (-4) - weak_edge_settings = deepcopy(formalism_edge_settings) + weak_edge_settings = deepcopy(default_edge_settings) interaction_type_settings[InteractionType.WEAK] = ( weak_edge_settings, @@ -231,15 +226,10 @@ def create_interaction_settings( # ruff: ignore[too-many-positional-arguments] CharmConservation(), StrangenessConservation(), BottomnessConservation(), - parity_conservation, - c_parity_conservation, + parity_rule, + c_parity_rule, }) ) - if "helicity" in formalism: - em_node_settings.conservation_rules.update( - _with_priorities({parity_conservation_helicity}) - ) - em_node_settings.qn_domains.update({NodeQN.parity_prefactor: [-1, 1]}) em_node_settings.interaction_strength = 1 em_edge_settings = deepcopy(weak_edge_settings) @@ -252,7 +242,7 @@ def create_interaction_settings( # ruff: ignore[too-many-positional-arguments] strong_node_settings.conservation_rules.update( _with_priorities({ isospin_conservation, - g_parity_conservation, + g_parity_rule, }) ) @@ -300,7 +290,6 @@ def _create_domains(particle_db: ParticleCollection) -> dict[Any, list]: domains[EdgeQN.spin_magnitude] = __positive_halves_domain( particle_db, lambda p: p.spin ) - domains[EdgeQN.spin_projection] = __extend_negative(domains[EdgeQN.spin_magnitude]) domains[EdgeQN.isospin_magnitude] = __positive_halves_domain( particle_db, lambda p: 0 if p.isospin is None else p.isospin.magnitude, diff --git a/src/qrules/solving.py b/src/qrules/solving.py index eab830d1..44b704d5 100644 --- a/src/qrules/solving.py +++ b/src/qrules/solving.py @@ -38,7 +38,6 @@ ) from qrules.quantum_numbers import ( EdgeQuantumNumber, - EdgeQuantumNumbers, EdgeQuantumNumberTypes, NodeQuantumNumber, NodeQuantumNumberTypes, @@ -289,7 +288,7 @@ def strip_quantum_numbers( def merge_qn_problem_sets( qn_problem_sets: Iterable[QNProblemSet], - merge_qns: Iterable[EdgeQuantumNumberTypes] | None = None, + merge_qns: Iterable[EdgeQuantumNumberTypes], ) -> list[QNProblemSet]: """Merge problem sets whose initial facts differ only in the given quantum numbers. @@ -297,12 +296,10 @@ def merge_qn_problem_sets( grouped. If the differing fact values of a group form a full Cartesian product over the graph edges, the group is merged into a single `QNProblemSet` whose facts carry the values as a `list`. The `CSPSolver` registers such value ranges as regular - variables, so that all cases — by default, all spin-projection combinations of the - initial and final state — are solved within a single constraint problem instead of - one problem set per combination. + variables, so that all cases — for instance, an expansion over spin-projection + combinations of the initial and final state — are solved within a single + constraint problem instead of one problem set per combination. """ - if merge_qns is None: - merge_qns = {EdgeQuantumNumbers.spin_projection} merge_qn_set = set(merge_qns) groups: dict[tuple, list[QNProblemSet]] = defaultdict(list) for problem_set in qn_problem_sets: @@ -313,16 +310,122 @@ def merge_qn_problem_sets( return merged_problem_sets +def remove_dominated_qn_problem_sets( + qn_problem_sets: dict[float, list[QNProblemSet]], +) -> dict[float, list[QNProblemSet]]: + """Remove problem sets whose solutions are a subset of another problem set's. + + Problem sets with equal initial facts and equal quantum-number domains define + constraint problems over the same variables, so a problem set whose + conservation-rule sets contain another's (for every edge and node) can only + produce a subset of that other's solutions. When solving for the *union* of all + solutions — as `.find_qn_transitions` does — such problem sets are redundant. + This typically removes the interaction-type combinations that add rules on top + of a weaker interaction type, e.g. the ``strong`` node settings that extend the + ``EM`` settings by isospin and :math:`G`-parity conservation. Rule priorities + and interaction strengths do not affect the solution set and are ignored in the + comparison. + """ + groups: dict[tuple, list[tuple[float, QNProblemSet]]] = defaultdict(list) + for strength, problem_sets in qn_problem_sets.items(): + for problem_set in problem_sets: + key = ( + problem_set.topology, + _create_facts_key(problem_set.initial_facts), + _create_domains_key(problem_set.solving_settings), + ) + groups[key].append((strength, problem_set)) + surviving: dict[float, list[QNProblemSet]] = defaultdict(list) + for group in groups.values(): + rule_maps = [ + _create_rule_map(problem_set.solving_settings) for _, problem_set in group + ] + for i, (strength, problem_set) in enumerate(group): + is_dominated = any( + _is_sub_rule_map(rule_maps[j], rule_maps[i]) + and (rule_maps[j] != rule_maps[i] or j < i) + for j in range(len(group)) + if j != i + ) + if not is_dominated: + surviving[strength].append(problem_set) + return dict(surviving) + + +def _create_domains_key(settings: GraphSettings) -> tuple: + def element_domains(element_settings: EdgeSettings | NodeSettings) -> tuple: + return tuple( + sorted( + (qn_type.__name__, repr(domain)) + for qn_type, domain in element_settings.qn_domains.items() + ) + ) + + return ( + tuple(sorted((i, element_domains(s)) for i, s in settings.states.items())), + tuple( + sorted((i, element_domains(s)) for i, s in settings.interactions.items()) + ), + ) + + +def _create_rule_map(settings: GraphSettings) -> dict[tuple[str, int], frozenset[str]]: + """Map each graph element to its conservation rules, ignoring priorities. + + Rules are represented by name plus their instance attributes, so that (deep) + copies of a rule compare equal, but parametrized rule instances such as + `.MassConservation` with different parameters do not. + """ + + def rule_id(rule: Any) -> str: + if inspect.isfunction(rule): + return rule.__name__ + if isinstance(rule, str): + return rule + parameters = sorted(getattr(rule, "__dict__", {}).items()) + return f"{type(rule).__name__}{parameters}" + + return { + **{ + ("edge", i): frozenset(map(rule_id, s.conservation_rules)) + for i, s in settings.states.items() + }, + **{ + ("node", i): frozenset(map(rule_id, s.conservation_rules)) + for i, s in settings.interactions.items() + }, + } + + +def _is_sub_rule_map( + weaker: dict[tuple[str, int], frozenset[str]], + stronger: dict[tuple[str, int], frozenset[str]], +) -> bool: + return all(rules <= stronger[element] for element, rules in weaker.items()) + + def _create_merge_key( problem_set: QNProblemSet, merge_qns: set[EdgeQuantumNumberTypes] ) -> tuple: - facts = problem_set.initial_facts + return ( + problem_set.topology, + _create_facts_key(problem_set.initial_facts, exclude_qns=merge_qns), + _create_settings_key(problem_set.solving_settings), + ) + + +def _create_facts_key( + facts: GraphElementProperties, + exclude_qns: set[EdgeQuantumNumberTypes] | None = None, +) -> tuple: + if exclude_qns is None: + exclude_qns = set() states_key = tuple( sorted( (edge_id, qn_type.__name__, repr(value)) for edge_id, prop_map in facts.states.items() for qn_type, value in prop_map.items() - if qn_type not in merge_qns + if qn_type not in exclude_qns ) ) interactions_key = tuple( @@ -332,12 +435,7 @@ def _create_merge_key( for qn_type, value in prop_map.items() ) ) - return ( - problem_set.topology, - states_key, - interactions_key, - _create_settings_key(problem_set.solving_settings), - ) + return (states_key, interactions_key) def _create_settings_key(settings: GraphSettings) -> tuple: @@ -628,7 +726,7 @@ def _insert_allowed_states( for candidate in candidate_states: # need "shallow" copy of the nested dicts new_states = {i: copy(s) for i, s in new_solution.states.items()} - new_states[edge_id].update(candidate) # keep spin_projection + new_states[edge_id].update(candidate) graph = attrs.evolve(new_solution, states=new_states) graphs_with_candidates.append(graph) current_substituted_graphs = graphs_with_candidates @@ -649,8 +747,6 @@ def __is_sub_mapping( state: GraphEdgePropertyMap, reference_state: GraphEdgePropertyMap ) -> bool: for qn_type, qn_value in state.items(): - if qn_type is EdgeQuantumNumbers.spin_projection: - continue if qn_type not in reference_state: return False if qn_value != reference_state[qn_type]: @@ -1133,8 +1229,12 @@ def __convert_solution_keys( MutableTransition[GraphEdgePropertyMap, GraphNodePropertyMap] ] = [] for solution in solutions: - states: dict[int, GraphEdgePropertyMap] = defaultdict(dict) - interactions: dict[int, GraphNodePropertyMap] = defaultdict(dict) + # initialize all graph elements, so that nodes without solved quantum + # numbers (e.g. without LS couplings) still appear in the solution + states: dict[int, GraphEdgePropertyMap] = {i: {} for i in topology.edges} + interactions: dict[int, GraphNodePropertyMap] = { + i: {} for i in topology.nodes + } for var_string, value in solution.items(): ele_id, qn_type = self.__var_string_to_data[var_string] diff --git a/src/qrules/system_control.py b/src/qrules/system_control.py index 6a56e0b5..6cbefe19 100644 --- a/src/qrules/system_control.py +++ b/src/qrules/system_control.py @@ -5,7 +5,6 @@ import logging import operator from abc import ABC, abstractmethod -from fractions import Fraction from typing import TYPE_CHECKING import attrs @@ -25,7 +24,7 @@ if TYPE_CHECKING: from collections.abc import Callable, Iterable - from qrules.particle import Particle, ParticleCollection, ParticleWithSpin + from qrules.particle import Particle, ParticleCollection _LOGGER = logging.getLogger(__name__) @@ -34,10 +33,7 @@ GraphSettingsGroups = dict[Strength, list[tuple[MutableTransition, GraphSettings]]] -def create_edge_properties( - particle: Particle, - spin_projection: float | Fraction | None = None, -) -> GraphEdgePropertyMap: +def create_edge_properties(particle: Particle) -> GraphEdgePropertyMap: edge_qn_mapping: dict[str, type[EdgeQuantumNumber]] = { qn_name: qn_type for qn_name, qn_type in EdgeQuantumNumbers.__dict__.items() @@ -56,8 +52,6 @@ def create_edge_properties( elif "spin" in qn_name: property_map[EdgeQuantumNumbers.spin_magnitude] = value - if spin_projection is not None: - property_map[EdgeQuantumNumbers.spin_projection] = spin_projection if isospin is not None: property_map[EdgeQuantumNumbers.isospin_magnitude] = isospin.magnitude property_map[EdgeQuantumNumbers.isospin_projection] = isospin.projection @@ -85,35 +79,28 @@ def create_node_properties(interactions: InteractionProperties) -> GraphNodeProp return property_map -def find_particle( # ruff: ignore[undocumented-param] +def find_particle( state: GraphEdgePropertyMap, particle_db: ParticleCollection -) -> ParticleWithSpin: - """Create a Particle with spin projection from a qn dictionary. +) -> Particle: + """Find the `.Particle` that matches a quantum number dictionary by PID. The implementation assumes the edge properties match the attributes of a particle inside the `.ParticleCollection`. Args: - states: The quantum number dictionary. particle_db: A `.ParticleCollection` - which is used to retrieve a reference :code:`state` to lower the memory - footprint. + state: The quantum number dictionary. + particle_db: A `.ParticleCollection` which is used to retrieve a reference + `.Particle` to lower the memory footprint. Raises: KeyError: If the edge properties do not contain the pid information or no particle with the same pid is found in the `.ParticleCollection`. - - ValueError: If the edge properties do not contain spin projection info. """ pid = state[EdgeQuantumNumbers.pid] if pid is None: msg = f"{GraphEdgePropertyMap.__name__} does not contain a PID" raise ValueError(msg) - particle = particle_db.find(int(pid)) - spin_projection = state.get(EdgeQuantumNumbers.spin_projection) - if spin_projection is None: - msg = f"{GraphEdgePropertyMap.__name__} does not contain a spin projection" - raise ValueError(msg) - return particle, Fraction(spin_projection) + return particle_db.find(int(pid)) def create_interaction_properties( @@ -154,8 +141,8 @@ class InteractionDeterminator(ABC): @abstractmethod def check( self, - in_states: list[ParticleWithSpin], - out_states: list[ParticleWithSpin], + in_states: list[Particle], + out_states: list[Particle], interactions: InteractionProperties, ) -> list[InteractionType]: pass @@ -166,12 +153,12 @@ class GammaCheck(InteractionDeterminator): def check( # ruff: ignore[no-self-use] self, - in_states: list[ParticleWithSpin], - out_states: list[ParticleWithSpin], + in_states: list[Particle], + out_states: list[Particle], interactions: InteractionProperties, # ruff: ignore[unused-method-argument] ) -> list[InteractionType]: int_types = list(InteractionType) - for particle, _ in in_states + out_states: + for particle in in_states + out_states: if "gamma" in particle.name: int_types = [InteractionType.EM] break @@ -183,12 +170,12 @@ class LeptonCheck(InteractionDeterminator): def check( # ruff: ignore[no-self-use] self, - in_states: list[ParticleWithSpin], - out_states: list[ParticleWithSpin], + in_states: list[Particle], + out_states: list[Particle], interactions: InteractionProperties, # ruff: ignore[unused-method-argument] ) -> list[InteractionType]: node_interaction_types = list(InteractionType) - for particle, _ in in_states + out_states: + for particle in in_states + out_states: if particle.is_lepton(): if particle.name.startswith("nu("): node_interaction_types = [InteractionType.WEAK] @@ -201,10 +188,10 @@ def check( # ruff: ignore[no-self-use] def remove_duplicate_solutions( - solutions: list[MutableTransition[ParticleWithSpin, InteractionProperties]], + solutions: list[MutableTransition[Particle, InteractionProperties]], remove_qns_list: set[type[NodeQuantumNumber]] | None = None, ignore_qns_list: set[type[NodeQuantumNumber]] | None = None, -) -> list[MutableTransition[ParticleWithSpin, InteractionProperties]]: +) -> list[MutableTransition[Particle, InteractionProperties]]: if remove_qns_list is None: remove_qns_list = set() if ignore_qns_list is None: @@ -213,9 +200,7 @@ def remove_duplicate_solutions( _LOGGER.info(f"removing these qns from graphs: {remove_qns_list}") _LOGGER.info(f"ignoring qns in graph comparison: {ignore_qns_list}") - filtered_solutions: list[ - MutableTransition[ParticleWithSpin, InteractionProperties] - ] = [] + filtered_solutions: list[MutableTransition[Particle, InteractionProperties]] = [] remove_counter = 0 for sol_graph in solutions: sol_graph = _remove_qns_from_graph(sol_graph, remove_qns_list) @@ -234,9 +219,9 @@ def remove_duplicate_solutions( def _remove_qns_from_graph( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], + graph: MutableTransition[Particle, InteractionProperties], qn_list: set[type[NodeQuantumNumber]], -) -> MutableTransition[ParticleWithSpin, InteractionProperties]: +) -> MutableTransition[Particle, InteractionProperties]: new_interactions = {} for node_id in graph.topology.nodes: interactions = graph.interactions[node_id] @@ -327,7 +312,7 @@ def require_interaction_property( ingoing_particle_name: str, interaction_qn: type[NodeQuantumNumber], allowed_values: list, -) -> Callable[[MutableTransition[ParticleWithSpin, InteractionProperties]], bool]: +) -> Callable[[MutableTransition[Particle, InteractionProperties]], bool]: """Filter function. Closure, which can be used as a filter function in :func:`.filter_graphs`. @@ -352,7 +337,7 @@ def require_interaction_property( """ def check( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], + graph: MutableTransition[Particle, InteractionProperties], ) -> bool: node_ids = _find_node_ids_with_ingoing_particle_name( graph, ingoing_particle_name @@ -371,15 +356,14 @@ def check( def _find_node_ids_with_ingoing_particle_name( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], + graph: MutableTransition[Particle, InteractionProperties], ingoing_particle_name: str, ) -> list[int]: topology = graph.topology found_node_ids = [] for node_id in topology.nodes: for edge_id in topology.get_edge_ids_ingoing_to_node(node_id): - states = graph.states[edge_id] - edge_particle_name = states[0].name + edge_particle_name = graph.states[edge_id].name if str(ingoing_particle_name) in str(edge_particle_name): found_node_ids.append(node_id) break diff --git a/src/qrules/transition.py b/src/qrules/transition.py index 6a592c12..6e99d125 100644 --- a/src/qrules/transition.py +++ b/src/qrules/transition.py @@ -9,17 +9,11 @@ import attrs from attrs import define, field, frozen -from attrs.validators import in_, instance_of +from attrs.validators import in_ -from qrules._attrs import to_fraction from qrules._implementers import implement_pretty_repr -from qrules.combinatorics import ( - InitialFacts, - StateDefinitionInput, - as_state_definition, - ensure_nested_list, -) -from qrules.particle import Particle, ParticleCollection, ParticleWithSpin, load_pdg +from qrules.combinatorics import InitialFacts, ensure_nested_list +from qrules.particle import Particle, ParticleCollection, load_pdg from qrules.quantum_numbers import InteractionProperties, NodeQuantumNumber from qrules.settings import ( DEFAULT_INTERACTION_TYPES, @@ -51,7 +45,6 @@ if TYPE_CHECKING: from collections.abc import Iterable, Sequence - from fractions import Fraction from qrules.workflow import InteractionConfig @@ -65,10 +58,10 @@ ] """Name for the spin formalism to be used. -The options :code:`"helicity"`, :code:`"canonical-helicity"`, and :code:`"canonical"` -are all used for the helicity formalism, but :code:`"canonical-helicity"` and -:code:`"canonical"` generate angular momentum and coupled spins as well on the -interaction nodes. +The formalism only affects which interaction properties are kept in the output: with +:code:`"helicity"`, the angular momentum and coupled spin magnitudes are filtered from +the interaction nodes, whereas :code:`"canonical-helicity"` and :code:`"canonical"` +keep these :math:`LS`-combinations. """ @@ -129,7 +122,7 @@ class SolutionContainer: by `.collect_reaction_info`. """ - solutions: list[MutableTransition[ParticleWithSpin, InteractionProperties]] = field( + solutions: list[MutableTransition[Particle, InteractionProperties]] = field( factory=list ) """Transitions that satisfy all conservation rules.""" @@ -179,8 +172,7 @@ def to_qn_problem_set(self) -> QNProblemSet: for k, v in self.initial_facts.interactions.items() } states = { - k: create_edge_properties(v[0], v[1]) - for k, v in self.initial_facts.states.items() + k: create_edge_properties(v) for k, v in self.initial_facts.states.items() } return QNProblemSet( initial_facts=MutableTransition(self.topology, states, interactions), @@ -196,8 +188,8 @@ class StateTransitionManager: def __init__( # ruff: ignore[too-many-positional-arguments] self, - initial_state: Sequence[StateDefinitionInput], - final_state: Sequence[StateDefinitionInput], + initial_state: Sequence[str], + final_state: Sequence[str], particle_db: ParticleCollection | None = None, allowed_intermediate_particles: list[str] | None = None, interaction_type_settings: dict[ @@ -211,6 +203,7 @@ def __init__( # ruff: ignore[too-many-positional-arguments] mass_conservation_factor: float | None = 3.0, max_angular_momentum: int = 1, max_spin_magnitude: float = 2, + ls_couplings: bool = True, number_of_threads: int | None = None, ) -> None: if number_of_threads is not None: @@ -232,8 +225,8 @@ def __init__( # ruff: ignore[too-many-positional-arguments] self.__particles = particle_db self.solving_mode = solving_mode """Whether to search for all solutions or stop at the strongest interaction.""" - self.initial_state = list(map(as_state_definition, initial_state)) - self.final_state = list(map(as_state_definition, final_state)) + self.initial_state = list(initial_state) + self.final_state = list(final_state) self.interaction_type_settings = interaction_type_settings self.interaction_determinators: list[InteractionDeterminator] = [ @@ -266,12 +259,12 @@ def __init__( # ruff: ignore[too-many-positional-arguments] if not self.interaction_type_settings: self.interaction_type_settings = create_interaction_settings( - formalism, particle_db=self.__particles, nbody_topology=use_nbody_topology, mass_conservation_factor=mass_conservation_factor, max_angular_momentum=max_angular_momentum, max_spin_magnitude=max_spin_magnitude, + ls_couplings=ls_couplings, ) if allowed_intermediate_particles is None: @@ -387,15 +380,8 @@ def find_quantum_number_transitions( ) -@implement_pretty_repr -@frozen(order=True) -class State: - particle: Particle = field(validator=instance_of(Particle)) - spin_projection: Fraction = field(converter=to_fraction) - - -StateTransition = FrozenTransition[State, InteractionProperties] -"""Transition of some initial `.State` to a final `.State`.""" +StateTransition = FrozenTransition[Particle, InteractionProperties] +"""Transition of initial state `.Particle` instances to final state particles.""" def _sort_tuple( @@ -417,15 +403,13 @@ class ReactionInfo: def __attrs_post_init__(self) -> None: transition = self.transitions[0] - initial = {i: s.particle for i, s in transition.initial_states.items()} - final = {i: s.particle for i, s in transition.final_states.items()} - object.__setattr__(self, "final_state", final) - object.__setattr__(self, "initial_state", initial) + object.__setattr__(self, "final_state", dict(transition.final_states)) + object.__setattr__(self, "initial_state", dict(transition.initial_states)) def get_intermediate_particles(self) -> ParticleCollection: """Extract the names of the intermediate state particles.""" particles = { - state.particle + state for transition in self.transitions for state in transition.intermediate_states.values() } diff --git a/src/qrules/workflow.py b/src/qrules/workflow.py index b39efbb7..4d810c96 100644 --- a/src/qrules/workflow.py +++ b/src/qrules/workflow.py @@ -33,11 +33,9 @@ from qrules._implementers import implement_pretty_repr from qrules.argument_handling import get_required_qns from qrules.combinatorics import ( - as_state_definition, create_initial_facts, match_external_edges, permutate_topology_kinematically, - to_state_definitions, ) from qrules.particle import Particle, ParticleCollection, load_pdg from qrules.quantum_numbers import ( @@ -53,10 +51,8 @@ ) from qrules.solving import ( CSPSolver, - _create_merge_key, complete_intermediate_states, - merge_qn_problem_sets, - strip_quantum_numbers, + remove_dominated_qn_problem_sets, ) from qrules.system_control import ( GammaCheck, @@ -83,14 +79,12 @@ SolutionContainer, SolvingMode, SpinFormalism, - State, ) if TYPE_CHECKING: from collections.abc import Iterable, Sequence - from qrules.combinatorics import InitialFacts, StateDefinition, StateDefinitionInput - from qrules.particle import ParticleWithSpin + from qrules.combinatorics import InitialFacts from qrules.quantum_numbers import EdgeQuantumNumber, NodeQuantumNumber from qrules.solving import ( EdgeSettings, @@ -303,9 +297,6 @@ def create_intermediate_edge_qn_domains() -> dict: # if a list of intermediate states is given by user, # built a domain based on these states intermediate_edge_domains: dict[type[EdgeQuantumNumber], set] = defaultdict(set) - intermediate_edge_domains[EdgeQuantumNumbers.spin_projection].update( - weak_edge_settings.qn_domains[EdgeQuantumNumbers.spin_projection] - ) for particle_props in intermediate_particles.particles: for edge_qn, qn_value in particle_props.items(): if edge_qn in { @@ -431,14 +422,13 @@ def _restrict_domains_to_rules(graph_settings: GraphSettings) -> None: def create_problem_sets( # ruff: ignore[too-many-positional-arguments] - initial_state: Sequence[StateDefinitionInput], - final_state: Sequence[StateDefinitionInput], + initial_state: Sequence[str], + final_state: Sequence[str], particle_db: ParticleCollection, interaction_config: InteractionConfig, intermediate_particles: AllowedIntermediateParticles, topologies: Iterable[Topology], final_state_groupings: list[list[list[str]]] | None = None, - expand_spin_projections: bool = True, allowed_channels: Iterable[str] | None = None, ) -> dict[float, list[ProblemSet]]: """Create a `.ProblemSet` collection over all topologies, grouped by strength. @@ -448,8 +438,6 @@ def create_problem_sets( # ruff: ignore[too-many-positional-arguments] :code:`["s"]` or :code:`["t", "u"]`) are turned into problem sets. """ allowed_channels = _validate_channels(allowed_channels) - initial_state = list(map(as_state_definition, initial_state)) - final_state = list(map(as_state_definition, final_state)) problem_sets = [ ProblemSet(permutation, initial_facts, settings) for topology in topologies @@ -461,13 +449,9 @@ def create_problem_sets( # ruff: ignore[too-many-positional-arguments] ) if allowed_channels is None or determine_reaction_channel(permutation) in allowed_channels - for initial_facts in create_initial_facts( - permutation, - initial_state, - final_state, - particle_db, - expand_spin_projections, - ) + for initial_facts in [ + create_initial_facts(permutation, initial_state, final_state, particle_db) + ] for settings in create_graph_settings( permutation, initial_facts, interaction_config, intermediate_particles ) @@ -647,7 +631,7 @@ def _convert_to_particle_definitions( def collect_reaction_info( # ruff: ignore[complex-structure, too-many-branches] results: dict[float, SolutionContainer], - final_state: Sequence[StateDefinitionInput] | None = None, + final_state: Sequence[str] | None = None, formalism: SpinFormalism = "helicity", filter_remove_qns: set[type[NodeQuantumNumber]] | None = None, filter_ignore_qns: set[type[NodeQuantumNumber]] | None = None, @@ -718,14 +702,10 @@ def collect_reaction_info( # ruff: ignore[complex-structure, too-many-branches] match_external_edges(final_solutions) if final_state is not None: - state_definitions = list(map(as_state_definition, final_state)) final_solutions = [ - _match_final_state_ids(graph, state_definitions) - for graph in final_solutions + _match_final_state_ids(graph, final_state) for graph in final_solutions ] - transitions = [ - graph.freeze().convert(lambda s: State(*s)) for graph in final_solutions - ] + transitions = [graph.freeze() for graph in final_solutions] return ReactionInfo(transitions, formalism) @@ -744,15 +724,15 @@ class QNProblemSetCollection: """`.QNProblemSet` collections, grouped by interaction strength.""" intermediate_particles: AllowedIntermediateParticles """Selection that built the intermediate edge domains; reused by `solve`.""" - final_state: list[StateDefinition] + final_state: list[str] """Final state with which `collect_reaction_info` orders the edge IDs.""" formalism: SpinFormalism = "helicity" """Spin formalism that determines the quantum-number filters.""" def create_qn_problem_sets( # ruff: ignore[too-many-positional-arguments] - initial_state: Sequence[StateDefinitionInput], - final_state: Sequence[StateDefinitionInput], + initial_state: Sequence[str], + final_state: Sequence[str], particle_db: ParticleCollection | None = None, allowed_intermediate_particles: AllowedIntermediateParticles | Iterable[str] @@ -765,10 +745,9 @@ def create_qn_problem_sets( # ruff: ignore[too-many-positional-arguments] mass_conservation_factor: float | None = 3.0, max_angular_momentum: int = 1, max_spin_magnitude: float = 2, + ls_couplings: bool = True, final_state_groupings: list[list[list[str]]] | None = None, allowed_channels: Iterable[str] | None = None, - merge_spin_projections: bool = False, - spin_projections: bool = True, ) -> QNProblemSetCollection: """Create a `.QNProblemSet` collection for a reaction, grouped by strength. @@ -778,27 +757,14 @@ def create_qn_problem_sets( # ruff: ignore[too-many-positional-arguments] with `find_solutions`, optionally after inspecting or modifying its problem sets (see e.g. `.filter_quantum_number_problem_set`). - With :code:`merge_spin_projections`, the spin-projection combinations of the - initial and final state are merged into value ranges on a single `.QNProblemSet` - (see `.merge_qn_problem_sets`), which reduces the number of problem sets and speeds - up solving. - - With :code:`spin_projections=False`, the problem sets contain no spin projections - at all: the Cartesian expansion over spin-projection combinations is skipped - entirely, so the problem sets can only be solved at the :math:`J^{P(C)}` level - with `find_qn_transitions`. This is equivalent to passing the collection through - `strip_spin_projections` afterwards, but much cheaper. - The :code:`allowed_interaction_types` (e.g. :code:`"strong"` or :code:`["em", "weak"]`) restrict the interaction types of the default or given :code:`interaction_config`. For reactions with more than one initial state, :code:`allowed_channels` (e.g. :code:`["s"]` or :code:`["t", "u"]`) restricts the problem sets to specific Mandelstam channels (see - `.determine_reaction_channel`). + `.determine_reaction_channel`). With :code:`ls_couplings=False`, the solver does + not enumerate :math:`LS`-combinations (see `.create_interaction_settings`). """ - if not spin_projections and merge_spin_projections: - msg = "merge_spin_projections has no effect when spin_projections=False" - raise ValueError(msg) _validate_formalism(formalism) if particle_db is None: particle_db = load_pdg() @@ -812,12 +778,12 @@ def create_qn_problem_sets( # ruff: ignore[too-many-positional-arguments] if interaction_config is None: interaction_config = InteractionConfig( type_settings=create_interaction_settings( - formalism, particle_db=particle_db, nbody_topology=use_nbody_topology, mass_conservation_factor=mass_conservation_factor, max_angular_momentum=max_angular_momentum, max_spin_magnitude=max_spin_magnitude, + ls_couplings=ls_couplings, ) ) if allowed_interaction_types is not None: @@ -836,24 +802,14 @@ def create_qn_problem_sets( # ruff: ignore[too-many-positional-arguments] intermediate_particles, topologies, final_state_groupings, - expand_spin_projections=spin_projections, allowed_channels=allowed_channels, ) - qn_problem_sets = _to_qn_problem_sets(problem_sets) - if merge_spin_projections: - qn_problem_sets = { - strength: merge_qn_problem_sets(problems) - for strength, problems in qn_problem_sets.items() - } - collection = QNProblemSetCollection( - problem_sets=qn_problem_sets, + return QNProblemSetCollection( + problem_sets=_to_qn_problem_sets(problem_sets), intermediate_particles=intermediate_particles, - final_state=list(map(as_state_definition, final_state)), + final_state=list(final_state), formalism=formalism, ) - if not spin_projections: - return strip_spin_projections(collection) - return collection def _to_qn_problem_sets( @@ -868,7 +824,7 @@ def _to_qn_problem_sets( def find_solutions( # ruff: ignore[too-many-positional-arguments] qn_problem_sets: QNProblemSetCollection | dict[float, list[QNProblemSet]], particle_db: ParticleCollection, - final_state: Sequence[StateDefinitionInput] | None = None, + final_state: Sequence[str] | None = None, formalism: SpinFormalism | None = None, allowed_intermediate_particles: AllowedIntermediateParticles | Iterable[str] @@ -1004,8 +960,8 @@ def group_by_channel(self) -> dict[str, list[QNTransition]]: def generate_qn_transitions( # ruff: ignore[too-many-positional-arguments] - initial_state: StateDefinitionInput | Sequence[StateDefinitionInput], - final_state: Sequence[StateDefinitionInput], + initial_state: str | Sequence[str], + final_state: Sequence[str], particle_db: ParticleCollection | None = None, allowed_intermediate_particles: Iterable[str] | str | None = None, allowed_interaction_types: str | Iterable[str] | None = None, @@ -1013,24 +969,26 @@ def generate_qn_transitions( # ruff: ignore[too-many-positional-arguments] mass_conservation_factor: float | None = 3.0, max_angular_momentum: int = 1, max_spin_magnitude: float = 2, + ls_couplings: bool = True, final_state_groupings: list[list[list[str]]] | None = None, allowed_channels: Iterable[str] | None = None, topology_building: str = "isobar", ) -> QNReactionInfo: - """Generate allowed transitions without spin projections. + """Generate allowed transitions purely at the quantum-number level. The quantum-number-level counterpart of `.generate_transitions`: chains - `create_qn_problem_sets` with :code:`spin_projections=False` and - `find_qn_transitions`, so that the reaction is solved at the :math:`J^{P(C)}` - level. Since the spin-projection combinatorics is avoided both in the problem sets - and in the constraint problem, this is feasible for reactions with many-body final - states, for which `.generate_transitions` would take impractically long. The - arguments mirror those of `.generate_transitions`. + `create_qn_problem_sets` and `find_qn_transitions`, so that the intermediate + states of the reaction are not matched against a particle database, but remain + the quantum-number sets that were solved for. The arguments mirror those of + `.generate_transitions`; :code:`ls_couplings=False` additionally switches off the + enumeration of :math:`LS`-combinations (see `.create_interaction_settings`). """ + if isinstance(initial_state, str): + initial_state = [initial_state] if particle_db is None: particle_db = load_pdg() qn_problem_sets = create_qn_problem_sets( - to_state_definitions(initial_state), + initial_state, final_state, particle_db=particle_db, allowed_intermediate_particles=allowed_intermediate_particles, @@ -1040,9 +998,9 @@ def generate_qn_transitions( # ruff: ignore[too-many-positional-arguments] mass_conservation_factor=mass_conservation_factor, max_angular_momentum=max_angular_momentum, max_spin_magnitude=max_spin_magnitude, + ls_couplings=ls_couplings, final_state_groupings=final_state_groupings, allowed_channels=allowed_channels, - spin_projections=False, ) transitions = find_qn_transitions(qn_problem_sets, particle_db) if not transitions: @@ -1051,83 +1009,44 @@ def generate_qn_transitions( # ruff: ignore[too-many-positional-arguments] return QNReactionInfo(transitions) -@overload -def strip_spin_projections( - qn_problem_sets: QNProblemSetCollection, -) -> QNProblemSetCollection: ... -@overload -def strip_spin_projections( - qn_problem_sets: dict[float, list[QNProblemSet]], -) -> dict[float, list[QNProblemSet]]: ... -def strip_spin_projections( - qn_problem_sets: QNProblemSetCollection | dict[float, list[QNProblemSet]], -) -> QNProblemSetCollection | dict[float, list[QNProblemSet]]: - """Remove spin projections from the problem sets, for :math:`J^{P(C)}`-level solving. - - Removes the spin projections of the initial and final state from the facts, the - spin-projection domains of the intermediate edges, and the :math:`L`/:math:`S` - projection domains of the interaction nodes. Problem sets that thereby become - identical — such as the expansion over all spin-projection combinations — are - deduplicated. Rules that require spin projections are skipped and reported by the - solver through the not-executed-rules mechanism. - - .. tip:: `create_qn_problem_sets` with :code:`spin_projections=False` produces the - same problem sets without generating the spin-projection expansion in the - first place, which is considerably faster for reactions with many spin-carrying - states. - """ - if isinstance(qn_problem_sets, QNProblemSetCollection): - stripped_collection = copy(qn_problem_sets) - stripped_collection.problem_sets = strip_spin_projections( - qn_problem_sets.problem_sets - ) - return stripped_collection - return { - strength: _unique_problem_sets( - strip_quantum_numbers( - problem_set, - edge_qns={EdgeQuantumNumbers.spin_projection}, - node_qns={ - NodeQuantumNumbers.l_projection, - NodeQuantumNumbers.s_projection, - }, - ) - for problem_set in problem_sets - ) - for strength, problem_sets in qn_problem_sets.items() - } - - -def _unique_problem_sets(problem_sets: Iterable[QNProblemSet]) -> list[QNProblemSet]: - unique: dict[tuple, QNProblemSet] = {} - for problem_set in problem_sets: - unique.setdefault(_create_merge_key(problem_set, set()), problem_set) - return list(unique.values()) +def _solve_qn_problem(qn_problem_set: QNProblemSet) -> tuple[QNProblemSet, QNResult]: + solver = CSPSolver() + return qn_problem_set, solver.find_solutions(qn_problem_set) def find_qn_transitions( qn_problem_sets: QNProblemSetCollection | dict[float, list[QNProblemSet]], particle_db: ParticleCollection | None = None, + number_of_threads: int | None = None, ) -> tuple[QNTransition, ...]: """Find allowed transitions purely at the quantum-number level. Solves the problem sets with the `.CSPSolver` without completing the intermediate states from a particle database: the intermediate states of the returned transitions carry exactly the quantum numbers that the problem sets declare as - domains. Combine with `strip_spin_projections` to obtain :math:`J^{P(C)}`-level - transitions for e.g. a Dalitz-plot decomposition. A :code:`particle_db` is only - used to resolve the initial and final states to `.Particle` instances (see - `collect_qn_transitions`). + domains, e.g. for a Dalitz-plot decomposition at the :math:`J^{P(C)}` level. A + :code:`particle_db` is only used to resolve the initial and final states to + `.Particle` instances (see `collect_qn_transitions`). Problem sets that cannot + contribute additional solutions are dropped + (`.remove_dominated_qn_problem_sets`) and the remaining ones are solved over + :code:`number_of_threads` processes (default: `.NumberOfThreads`). """ if isinstance(qn_problem_sets, QNProblemSetCollection): qn_problem_sets = qn_problem_sets.problem_sets - solver = CSPSolver() + qn_problem_sets = remove_dominated_qn_problem_sets(qn_problem_sets) + if number_of_threads is None: + number_of_threads = NumberOfThreads.get() qn_results: dict[float, list[tuple[QNProblemSet, QNResult]]] = defaultdict(list) - for strength, qn_problems in sorted(qn_problem_sets.items(), reverse=True): - qn_results[strength].extend( - (qn_problem_set, solver.find_solutions(qn_problem_set)) - for qn_problem_set in qn_problems - ) + sorted_problem_sets = sorted(qn_problem_sets.items(), reverse=True) + if number_of_threads > 1: + with Pool(number_of_threads) as pool: + for strength, qn_problems in sorted_problem_sets: + qn_results[strength].extend( + pool.map(_solve_qn_problem, qn_problems, chunksize=1) + ) + else: + for strength, qn_problems in sorted_problem_sets: + qn_results[strength].extend(map(_solve_qn_problem, qn_problems)) return collect_qn_transitions(qn_results, particle_db) @@ -1137,12 +1056,12 @@ def collect_qn_transitions( ) -> tuple[QNTransition, ...]: """Summarize solver results as unique quantum-number-level transitions. - The counterpart of `collect_reaction_info` for a quantum-number-only workflow: no - `.State` objects are created — each transition carries exactly the quantum numbers - that are known from the initial facts or were solved for. If a :code:`particle_db` - is given, the initial and final states — which are fully determined by their - `~.EdgeQuantumNumbers.pid` — are resolved to `.Particle` instances; - intermediate states always remain quantum-number property maps. + The counterpart of `collect_reaction_info` for a quantum-number-only workflow: + each transition carries exactly the quantum numbers that are known from the + initial facts or were solved for. If a :code:`particle_db` is given, the initial + and final states — which are fully determined by their `~.EdgeQuantumNumbers.pid` + — are resolved to `.Particle` instances; intermediate states always remain + quantum-number property maps. """ transitions: dict[QNTransition, None] = {} for qn_result_pairs in qn_results.values(): @@ -1158,12 +1077,16 @@ def collect_qn_transitions( for edge_id in external_edge_ids: pid = states[edge_id][EdgeQuantumNumbers.pid] states[edge_id] = particle_db.find(int(pid)) - interactions = dict(solution.interactions) - for node_id, node_facts in facts.interactions.items(): - interactions[node_id] = { - **node_facts, - **interactions.get(node_id, {}), + # a node without facts and solved quantum numbers (e.g. with + # ls_couplings=False) gets an empty property map + solved_interactions = dict(solution.interactions) + interactions: dict[int, dict] = { + node_id: { + **facts.interactions.get(node_id, {}), + **solved_interactions.get(node_id, {}), } + for node_id in topology.nodes + } transition: QNTransition = FrozenTransition( solution.topology, states={ @@ -1177,14 +1100,13 @@ def collect_qn_transitions( def _match_final_state_ids( - graph: MutableTransition[ParticleWithSpin, InteractionProperties], - state_definition: Sequence[StateDefinition], -) -> MutableTransition[ParticleWithSpin, InteractionProperties]: + graph: MutableTransition[Particle, InteractionProperties], + final_state: Sequence[str], +) -> MutableTransition[Particle, InteractionProperties]: """Temporary fix to https://github.com/ComPWA/qrules/issues/143.""" - particle_names = _strip_spin(state_definition) - name_to_id = {name: i for i, name in enumerate(particle_names)} + name_to_id = {name: i for i, name in enumerate(final_state)} id_remapping = { - name_to_id[graph.states[i][0].name]: i for i in graph.topology.outgoing_edge_ids + name_to_id[graph.states[i].name]: i for i in graph.topology.outgoing_edge_ids } new_topology = graph.topology.relabel_edges(id_remapping) return MutableTransition( @@ -1194,16 +1116,6 @@ def _match_final_state_ids( ) -def _strip_spin(state_definition: Sequence[StateDefinition]) -> list[str]: - particle_names = [] - for state in state_definition: - if isinstance(state, str): - particle_names.append(state) - else: - particle_names.append(state[0]) - return particle_names - - def _parse_interaction_types( description: str | Iterable[str], ) -> list[InteractionType]: diff --git a/tests/channels/test_gammap_to_p_pi0.py b/tests/channels/test_gammap_to_p_pi0.py index d195119a..70d57479 100644 --- a/tests/channels/test_gammap_to_p_pi0.py +++ b/tests/channels/test_gammap_to_p_pi0.py @@ -48,7 +48,7 @@ def test_group_by_channel(particle_database): ) channels = { channel: { - state.particle.name + state.name for transition in transitions for state in transition.intermediate_states.values() } diff --git a/tests/channels/test_jpsi_to_gamma_pi0_pi0.py b/tests/channels/test_jpsi_to_gamma_pi0_pi0.py index 8125d2c9..4a557d0f 100644 --- a/tests/channels/test_jpsi_to_gamma_pi0_pi0.py +++ b/tests/channels/test_jpsi_to_gamma_pi0_pi0.py @@ -8,10 +8,10 @@ def describe_reaction_generation(): @pytest.mark.parametrize( ("allowed_intermediate_particles", "n_topologies", "number_of_solutions"), [ - (["f(0)(1500)"], 1, 4), - (["f(0)(980)", "f(0)(1500)"], 1, 8), - (["f(2)(1270)"], 1, 12), - (["omega(782)"], 1, 8), + (["f(0)(1500)"], 1, 1), + (["f(0)(980)", "f(0)(1500)"], 1, 2), + (["f(2)(1270)"], 1, 1), + (["omega(782)"], 1, 1), ( [ "f(0)(980)", @@ -21,7 +21,7 @@ def describe_reaction_generation(): "omega(782)", ], 2, - 40, + 5, ), ], ) @@ -33,7 +33,7 @@ def it_number_of_solutions( number_of_solutions, ): reaction = qrules.generate_transitions( - initial_state=("J/psi(1S)", [-1, +1]), + initial_state="J/psi(1S)", final_state=["gamma", "pi0", "pi0"], particle_db=particle_database, allowed_interaction_types=["strong", "EM"], @@ -49,7 +49,7 @@ def it_number_of_solutions( def it_id_to_particle_mappings(particle_database): reaction = qrules.generate_transitions( - initial_state=("J/psi(1S)", [-1, +1]), + initial_state="J/psi(1S)", final_state=["gamma", "pi0", "pi0"], particle_db=particle_database, allowed_interaction_types="strong", @@ -57,12 +57,10 @@ def it_id_to_particle_mappings(particle_database): formalism="helicity", ) assert len(reaction.group_by_topology()) == 1 - assert len(reaction.transitions) == 4 + assert len(reaction.transitions) == 1 iter_transitions = iter(reaction.transitions) first_transition = next(iter_transitions) - graph = first_transition.convert( - lambda s: (s.particle, s.spin_projection) - ).unfreeze() + graph = first_transition.unfreeze() ref_mapping_fs = _create_edge_id_particle_mapping( graph, graph.topology.outgoing_edge_ids ) @@ -70,9 +68,7 @@ def it_id_to_particle_mappings(particle_database): graph, graph.topology.incoming_edge_ids ) for transition in iter_transitions: - graph = transition.convert( - lambda s: (s.particle, s.spin_projection) - ).unfreeze() + graph = transition.unfreeze() assert ref_mapping_fs == _create_edge_id_particle_mapping( graph, graph.topology.outgoing_edge_ids ) diff --git a/tests/channels/test_nbody_reactions.py b/tests/channels/test_nbody_reactions.py index e854bfe9..dfd74a6c 100644 --- a/tests/channels/test_nbody_reactions.py +++ b/tests/channels/test_nbody_reactions.py @@ -165,7 +165,6 @@ def reduce_violated_rules( { "isospin_conservation", # Clebsch Gordan coefficient = 0 "c_parity_conservation", - "identical_particle_symmetrization", }, ), ( diff --git a/tests/channels/test_psi2s_to_eta_k_kstar.py b/tests/channels/test_psi2s_to_eta_k_kstar.py index fcc1133c..f4e8f8d5 100644 --- a/tests/channels/test_psi2s_to_eta_k_kstar.py +++ b/tests/channels/test_psi2s_to_eta_k_kstar.py @@ -18,7 +18,7 @@ ) def test_resonances(formalism: SpinFormalism, resonances, modified_pdg): reaction = qrules.generate_transitions( - initial_state=("psi(2S)", [+1, -1]), + initial_state="psi(2S)", final_state=["eta", "K-", "K*(892)+"], allowed_intermediate_particles=resonances, allowed_interaction_types=["em"], diff --git a/tests/channels/test_y_to_d0_d0bar_pi0_pi0.py b/tests/channels/test_y_to_d0_d0bar_pi0_pi0.py index 8b0049d2..7528bc18 100644 --- a/tests/channels/test_y_to_d0_d0bar_pi0_pi0.py +++ b/tests/channels/test_y_to_d0_d0bar_pi0_pi0.py @@ -8,13 +8,13 @@ def describe_reaction_generation(): @pytest.mark.parametrize( ("formalism", "n_solutions"), [ - ("helicity", 14), - ("canonical-helicity", 28), # two different LS couplings 2*14 = 28 + ("helicity", 1), + ("canonical-helicity", 3), # three different LS couplings ], ) def it_simple(formalism, n_solutions, particle_database): reaction = qrules.generate_transitions( - initial_state=[("Y(4260)", [-1, +1])], + initial_state=["Y(4260)"], final_state=["D*(2007)0", "D*(2007)~0"], particle_db=particle_database, formalism=formalism, @@ -27,13 +27,13 @@ def it_simple(formalism, n_solutions, particle_database): @pytest.mark.parametrize( ("formalism", "n_solutions"), [ - ("helicity", 14), - ("canonical-helicity", 28), # two different LS couplings 2*14 = 28 + ("helicity", 1), + ("canonical-helicity", 3), # three different LS couplings ], ) def it_full(formalism, n_solutions, particle_database): stm = StateTransitionManager( - initial_state=[("Y(4260)", [-1, +1])], + initial_state=["Y(4260)"], final_state=["D0", "D~0", "pi0", "pi0"], particle_db=particle_database, allowed_intermediate_particles=["D*"], @@ -49,7 +49,7 @@ def it_full(formalism, n_solutions, particle_database): def it_resonance_filter(particle_database): # https://github.com/ComPWA/qrules/issues/33 stm = StateTransitionManager( - initial_state=[("Y(4260)", [-1, +1])], + initial_state=["Y(4260)"], final_state=["D0", "D~0", "pi0", "pi0"], allowed_intermediate_particles=["D"], particle_db=particle_database, diff --git a/tests/unit/conftest.py b/tests/unit/conftest.py index a0d414e2..b570dd40 100644 --- a/tests/unit/conftest.py +++ b/tests/unit/conftest.py @@ -18,7 +18,7 @@ def reaction(request: SubRequest) -> ReactionInfo: formalism: SpinFormalism = request.param return qrules.generate_transitions( - initial_state=[("J/psi(1S)", [-1, 1])], + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"], allowed_interaction_types="strong", diff --git a/tests/unit/conservation_rules/test_c_parity.py b/tests/unit/conservation_rules/test_c_parity.py index 09f643db..df7ea7fb 100644 --- a/tests/unit/conservation_rules/test_c_parity.py +++ b/tests/unit/conservation_rules/test_c_parity.py @@ -4,6 +4,8 @@ import pytest from qrules.conservation_rules import ( + CParityCoupling, + CParityCouplingEdgeInput, CParityEdgeInput, CParityNodeInput, c_parity_conservation, @@ -100,3 +102,59 @@ def it_c_parity_multiparticle_boson(rule_input, expected): ) def it_c_parity_multiparticle_fermion(rule_input, expected): assert c_parity_conservation(*rule_input) is expected + + +def describe_CParityCoupling(): + pi_plus = CParityCouplingEdgeInput(spin_magnitude=0, parity=-1, pid=211) + pi_minus = CParityCouplingEdgeInput(spin_magnitude=0, parity=-1, pid=-211) + proton = CParityCouplingEdgeInput(spin_magnitude=0.5, parity=+1, pid=2212) + antiproton = CParityCouplingEdgeInput(spin_magnitude=0.5, parity=-1, pid=-2212) + + @pytest.mark.parametrize( + ("max_angular_momentum", "single", "expected"), + [ + (1, CParityCouplingEdgeInput(1, -1, 113, c_parity=-1), True), + (0, CParityCouplingEdgeInput(1, -1, 113, c_parity=-1), False), + (0, CParityCouplingEdgeInput(0, +1, 9010221, c_parity=+1), True), + (3, CParityCouplingEdgeInput(1, -1, 1, c_parity=+1), False), + (3, CParityCouplingEdgeInput(1, -1, 1, c_parity=None), True), + ], + ) + def it_couples_boson_pairs(max_angular_momentum, single, expected): + rule = CParityCoupling(max_angular_momentum) + assert rule([single], [pi_plus, pi_minus]) is expected + assert rule([pi_plus, pi_minus], [single]) is expected + + @pytest.mark.parametrize( + ("max_angular_momentum", "single", "expected"), + [ + # eta(c): 0-+ with L=0, S=0 + (0, CParityCouplingEdgeInput(0, -1, 441, c_parity=+1), True), + # J/psi: 1-- with L=0, S=1 + (0, CParityCouplingEdgeInput(1, -1, 443, c_parity=-1), True), + # chi(c1): 1++ needs L=1, S=1 + (0, CParityCouplingEdgeInput(1, +1, 20443, c_parity=+1), False), + (1, CParityCouplingEdgeInput(1, +1, 20443, c_parity=+1), True), + # exotic 1-+: parity allows L=0, but C-parity then requires S=0, so J=0 + (1, CParityCouplingEdgeInput(1, -1, 1, c_parity=+1), False), + ], + ) + def it_couples_fermion_pairs(max_angular_momentum, single, expected): + rule = CParityCoupling(max_angular_momentum) + assert rule([single], [proton, antiproton]) is expected + + def it_multiplies_defined_c_parities(): + rule = CParityCoupling(max_angular_momentum=0) + photon = CParityCouplingEdgeInput(1, -1, 22, c_parity=-1) + jpsi = CParityCouplingEdgeInput(1, -1, 443, c_parity=-1) + eta_c = CParityCouplingEdgeInput(0, -1, 441, c_parity=+1) + chi_c1 = CParityCouplingEdgeInput(1, +1, 20443, c_parity=+1) + assert rule([jpsi], [photon, eta_c]) is True + assert rule([jpsi], [photon, chi_c1]) is True + assert rule([eta_c], [photon, chi_c1]) is False + + def it_skips_pairs_that_are_not_particle_antiparticle(): + rule = CParityCoupling(max_angular_momentum=0) + exotic = CParityCouplingEdgeInput(1, -1, 1, c_parity=+1) + kaon = CParityCouplingEdgeInput(0, -1, pid=321) + assert rule([exotic], [pi_plus, kaon]) is True diff --git a/tests/unit/conservation_rules/test_duck_typing.py b/tests/unit/conservation_rules/test_duck_typing.py index 56856584..90634ebe 100644 --- a/tests/unit/conservation_rules/test_duck_typing.py +++ b/tests/unit/conservation_rules/test_duck_typing.py @@ -31,7 +31,9 @@ def test_protocol_compliance(): edge_input_classes = __get_duck_types(EdgeQuantumNumbers) assert edge_input_classes == { + conservation_rules.CParityCouplingEdgeInput, conservation_rules.CParityEdgeInput, + conservation_rules.GParityCouplingEdgeInput, conservation_rules.GParityEdgeInput, conservation_rules.GParityValidityInput, conservation_rules.HelicityParityEdgeInput, diff --git a/tests/unit/conservation_rules/test_g_parity.py b/tests/unit/conservation_rules/test_g_parity.py index 6bb690d9..31c4a999 100644 --- a/tests/unit/conservation_rules/test_g_parity.py +++ b/tests/unit/conservation_rules/test_g_parity.py @@ -4,6 +4,8 @@ import pytest from qrules.conservation_rules import ( + GParityCoupling, + GParityCouplingEdgeInput, GParityEdgeInput, GParityNodeInput, g_parity_conservation, @@ -96,3 +98,63 @@ def it_g_parity_all_defined(rule_input, expected): ) def it_g_parity_multiparticle_boson(rule_input, expected): assert g_parity_conservation(*rule_input) is expected + + +def describe_GParityCoupling(): + pi_plus = GParityCouplingEdgeInput( + 1, spin_magnitude=0, parity=-1, pid=211, g_parity=-1 + ) + pi_minus = GParityCouplingEdgeInput( + 1, spin_magnitude=0, parity=-1, pid=-211, g_parity=-1 + ) + k_plus = GParityCouplingEdgeInput(0.5, spin_magnitude=0, parity=-1, pid=321) + k_minus = GParityCouplingEdgeInput(0.5, spin_magnitude=0, parity=-1, pid=-321) + proton = GParityCouplingEdgeInput(0.5, spin_magnitude=0.5, parity=+1, pid=2212) + antiproton = GParityCouplingEdgeInput(0.5, spin_magnitude=0.5, parity=-1, pid=-2212) + + @pytest.mark.parametrize( + ("max_angular_momentum", "single", "expected"), + [ + # phi: I=0 and G=-1 with L=1 + (1, GParityCouplingEdgeInput(0, 1, -1, 333, g_parity=-1), True), + (0, GParityCouplingEdgeInput(0, 1, -1, 333, g_parity=-1), False), + # rho: I=1 and G=+1 with L=1 + (1, GParityCouplingEdgeInput(1, 1, -1, 113, g_parity=+1), True), + # f0: I=0 and G=+1 with L=0 + (0, GParityCouplingEdgeInput(0, 0, +1, 9010221, g_parity=+1), True), + # exotic 1-+ with I=1 and G=-1: parity requires odd L, but G-parity even L + (3, GParityCouplingEdgeInput(1, 1, -1, 1, g_parity=-1), False), + (3, GParityCouplingEdgeInput(1, 1, -1, 1, g_parity=None), True), + ], + ) + def it_couples_boson_pairs(max_angular_momentum, single, expected): + rule = GParityCoupling(max_angular_momentum) + assert rule([single], [k_plus, k_minus]) is expected + assert rule([k_plus, k_minus], [single]) is expected + + @pytest.mark.parametrize( + ("max_angular_momentum", "single", "expected"), + [ + # rho: 1-- with I=1 needs L=0, S=1, so C=-1 and G=+1 + (0, GParityCouplingEdgeInput(1, 1, -1, 113, g_parity=+1), True), + # exotic 1-+ with I=1: C=+1 requires S=0, so J=0 + (1, GParityCouplingEdgeInput(1, 1, -1, 1, g_parity=-1), False), + ], + ) + def it_couples_fermion_pairs(max_angular_momentum, single, expected): + rule = GParityCoupling(max_angular_momentum) + assert rule([single], [proton, antiproton]) is expected + + def it_multiplies_defined_g_parities(): + rule = GParityCoupling(max_angular_momentum=0) + rho = GParityCouplingEdgeInput(1, 1, -1, 113, g_parity=+1) + omega = GParityCouplingEdgeInput(0, 1, -1, 223, g_parity=-1) + assert rule([rho], [pi_plus, pi_minus]) is True + assert rule([omega], [pi_plus, pi_minus]) is False + assert rule([omega], [rho, pi_plus]) is True + assert rule([rho], [rho, pi_plus]) is False + + def it_skips_half_integer_isospin(): + rule = GParityCoupling(max_angular_momentum=0) + single = GParityCouplingEdgeInput(0.5, 1, -1, 1, g_parity=-1) + assert rule([single], [k_plus, k_minus]) is True diff --git a/tests/unit/io/conftest.py b/tests/unit/io/conftest.py index 647101fd..19eef62b 100644 --- a/tests/unit/io/conftest.py +++ b/tests/unit/io/conftest.py @@ -19,8 +19,8 @@ def particle_selection(particle_database: ParticleCollection): @pytest.fixture def stm() -> StateTransitionManager: stm = StateTransitionManager( - initial_state=[("J/psi(1S)", [+1])], - final_state=["K0", ("Sigma+", [+0.5]), ("p~", [+0.5])], + initial_state=["J/psi(1S)"], + final_state=["K0", "Sigma+", "p~"], allowed_intermediate_particles=["Sigma(1750)"], formalism="canonical-helicity", ) @@ -40,4 +40,4 @@ def qn_problem_and_result( ) -> tuple[QNProblemSet, QNResult]: qn_solutions = stm.find_quantum_number_transitions(problem_sets) strong_qn_solutions = qn_solutions[3600.0] - return strong_qn_solutions[1] + return next(pair for pair in strong_qn_solutions if pair[1].solutions) diff --git a/tests/unit/io/test_dot.py b/tests/unit/io/test_dot.py index 5b012536..da83dfd4 100644 --- a/tests/unit/io/test_dot.py +++ b/tests/unit/io/test_dot.py @@ -34,19 +34,19 @@ def it_exact_format(reaction: ReactionInfo): node [shape=none, width=0] edge [arrowhead=none] bgcolor=none - 0 [label="0: gamma[-1]"] - 1 [label="1: pi0[0]"] - 2 [label="2: pi0[0]"] - A [label="J/psi(1S)[-1]"] + 0 [label="0: gamma"] + 1 [label="1: pi0"] + 2 [label="2: pi0"] + A [label="J/psi(1S)"] { rank=same; A } { rank=same; 0 1 2 } A -> N0 - N0 -> N1 [label="f(0)(980)[0]"] + N0 -> N1 [label="f(0)(980)"] N0 -> 0 N1 -> 1 N1 -> 2 - N0 [label="P=+1"] - N1 [label="P=+1"] + N0 [label=""] + N1 [label=""] } """.replace("\n ", "\n") else: @@ -56,19 +56,19 @@ def it_exact_format(reaction: ReactionInfo): node [shape=none, width=0] edge [arrowhead=none] bgcolor=none - 0 [label="0: gamma[-1]"] - 1 [label="1: pi0[0]"] - 2 [label="2: pi0[0]"] - A [label="J/psi(1S)[-1]"] + 0 [label="0: gamma"] + 1 [label="1: pi0"] + 2 [label="2: pi0"] + A [label="J/psi(1S)"] { rank=same; A } { rank=same; 0 1 2 } A -> N0 - N0 -> N1 [label="f(0)(980)[0]"] + N0 -> N1 [label="f(0)(980)"] N0 -> 0 N1 -> 1 N1 -> 2 - N0 [label="L=|0,0⟩\nS=|1,-1⟩\nP=+1"] - N1 [label="L=|0,0⟩\nS=|0,0⟩\nP=+1"] + N0 [label="L=0\nS=1"] + N1 [label="L=0\nS=0"] } """.replace("\n ", "\n") assert src.strip() == expected_dot.strip() @@ -97,8 +97,8 @@ def it_with_styled_edges_and_nodes(reaction: ReactionInfo, output_dir): assert is_valid_dot(src) with open(output_dir + f"styled_{reaction.formalism}.gv", "w") as stream: stream.write(src) - assert '0 [fontcolor="blue", label="0: gamma[-1]"]' in src - assert 'N0 -> N1 [fontcolor="blue", label="f(0)(980)[0]"]' in src + assert '0 [fontcolor="blue", label="0: gamma"]' in src + assert 'N0 -> N1 [fontcolor="blue", label="f(0)(980)"]' in src assert 'N0 [fontcolor="darkgreen", shape="ellipse", label=""]' in src def it_no_label_overwriting(reaction: ReactionInfo): @@ -118,7 +118,7 @@ def it_no_label_overwriting(reaction: ReactionInfo): ) def it_problemset(formalism: SpinFormalism): stm = StateTransitionManager( - initial_state=[("J/psi(1S)", [+1])], + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], formalism=formalism, ) diff --git a/tests/unit/io/test_labels.py b/tests/unit/io/test_labels.py index 67dcf4e2..ec8f2600 100644 --- a/tests/unit/io/test_labels.py +++ b/tests/unit/io/test_labels.py @@ -18,7 +18,7 @@ from qrules.particle import Particle, ParticleCollection from qrules.quantum_numbers import InteractionProperties from qrules.solving import QNProblemSet, QNResult -from qrules.transition import ProblemSet, ReactionInfo, State +from qrules.transition import ProblemSet, ReactionInfo def describe_as_latex(): @@ -51,12 +51,9 @@ def __str__(self) -> str: "No LaTeX label renderer implemented type UnsupportedLabel" in caplog.text ) - def it_particle_and_state(particle_database: ParticleCollection): + def it_particle(particle_database: ParticleCollection): particle = particle_database["J/psi(1S)"] assert as_latex(particle) == R"J/\psi(1S)" - expected_state = R"J/\psi(1S)\left[\text{-}\frac{1}{2}\right]" - assert as_latex(State(particle, Fraction(-1, 2))) == expected_state - assert as_latex((particle, Fraction(-1, 2))) == expected_state particle_with_custom_latex = attrs.evolve( particle, @@ -166,9 +163,9 @@ def test_create_edge_label_accepts_renderer(reaction: ReactionInfo): render_label=as_latex, ) - assert plain_label.startswith(state.particle.name) - assert state.particle.latex is not None - assert latex_label.startswith(state.particle.latex) + assert plain_label.startswith(state.name) + assert state.latex is not None + assert latex_label.startswith(state.latex) multiline_label = create_edge_label( transition, @@ -198,7 +195,6 @@ def it_dict( print(src) expected_dot = dedent(""" RULES - spin_validity - 62 isospin_validity - 61 g_parity_validity - 60 gellmann_nishijima - 50 @@ -215,7 +211,6 @@ def it_dict( muon_lepton_number ∊ [0] parity ∊ [-1, +1] spin_magnitude ∊ [1/2] - spin_projection ∊ [-4, -7/2, -3, -5/2, -2, -3/2, -1, -1/2, 0, +1/2, +1, +3/2, +2, +5/2, +3, +7/2, +4] strangeness ∊ [-1, +1] tau_lepton_number ∊ [0] topness ∊ [0] @@ -229,7 +224,6 @@ def it_dict( RULES ChargeConservation - 100 BaryonNumberConservation - 90 - ls_spin_validity - 89 CharmConservation - 70 StrangenessConservation - 69 BottomnessConservation - 68 @@ -239,19 +233,12 @@ def it_dict( TauLNConservation - 43 MassConservation - 10 spin_magnitude_conservation - 8 - helicity_conservation - 7 parity_conservation - 6 c_parity_conservation - 5 - parity_conservation_helicity - 4 g_parity_conservation - 3 - identical_particle_symmetrization - 2 - clebsch_gordan_helicity_to_canonical - 1 DOMAINS l_magnitude ∊ [0, 1] - l_projection ∊ [0] - parity_prefactor ∊ [-1, +1] s_magnitude ∊ [0, 1/2, 1, 3/2, 2] - s_projection ∊ [-2, -3/2, -1, -1/2, 0, +1/2, +1, +3/2, +2] """).strip() assert src == expected_dot @@ -261,9 +248,6 @@ def it_dict( expected_dot = dedent(""" l_magnitude = 0 s_magnitude = 1/2 - l_projection = 0 - s_projection = -1/2 - parity_prefactor = +1 """).strip() assert src == expected_dot @@ -271,7 +255,6 @@ def it_dict( lines = set(src.splitlines()) expected_lines = { "spin_magnitude = 1/2", - "spin_projection = +1/2", "parity = +1", "isospin_magnitude = 1", "isospin_projection = -1", @@ -286,7 +269,7 @@ def it_dict( latex = as_latex(intermediate_setting) assert R"\text{RULES}" in latex - assert R"\text{spin\_validity - 62}" in latex + assert R"\text{isospin\_validity - 61}" in latex assert R"\text{DOMAINS}" in latex assert R"\text{spin\_magnitude} \in \left[\frac{1}{2}\right]" in latex @@ -294,20 +277,15 @@ def it_dict( assert R"\text{ChargeConservation - 100}" in latex assert R"\text{l\_magnitude} \in \left[0, 1\right]" in latex assert ( - R"\text{parity\_prefactor} \in \left[\text{-}1, \text{+}1\right]" in latex + R"\text{s\_magnitude} \in \left[0, \frac{1}{2}, 1, \frac{3}{2}, 2\right]" + in latex ) - assert ( - R"\text{s\_projection} \in \left[\text{-}2, " - R"\text{-}\frac{3}{2}, \text{-}1, \text{-}\frac{1}{2}, 0, " - R"\text{+}\frac{1}{2}, \text{+}1, \text{+}\frac{3}{2}, " - R"\text{+}2\right]" - ) in latex latex = as_latex(intermediate_state) assert R"\text{spin\_magnitude} = \frac{1}{2}" in latex assert R"\text{parity} = \text{+}1" in latex - def it_spin_tuple(particle_database: ParticleCollection): + def it_spin_tuple(): # non-spin src = as_string(("a", "b", "c")) assert src == "a\nb\nc" @@ -317,13 +295,8 @@ def it_spin_tuple(particle_database: ParticleCollection): # spin src = as_string((2, 1)) assert src == "|2,+1⟩" - - # particle with spin projection - pion = particle_database["J/psi(1S)"] - src = as_string((pion, 1)) - assert src == "J/psi(1S)[+1]" - src = as_string((pion, Fraction(-1))) - assert src == "J/psi(1S)[-1]" + src = as_string((Fraction(1, 2), Fraction(-1, 2))) + assert src == "|1/2,-1/2⟩" def test_collapse_graphs( @@ -400,26 +373,25 @@ def test_strip_projections(skh_particle_version: str): assert skh_particle_version is not None # skips test if particle version too low resonance = "Sigma(1670)~-" reaction = qrules.generate_transitions( - initial_state=[("J/psi(1S)", [+1])], - final_state=["K0", ("Sigma+", [+0.5]), ("p~", [+0.5])], + initial_state="J/psi(1S)", + final_state=["K0", "Sigma+", "p~"], allowed_intermediate_particles=[resonance], allowed_interaction_types="strong", ) - - assert len(reaction.transitions) == 5 transition = reaction.transitions[0] - - assert transition.intermediate_states[3].particle.name == resonance - assert transition.intermediate_states[3].spin_projection == -0.5 + assert transition.intermediate_states[3].name == resonance assert len(transition.interactions) == 2 - assert transition.interactions[0].s_projection == 1 - assert transition.interactions[0].l_projection == 0 - assert transition.interactions[1].s_projection == -0.5 - assert transition.interactions[1].l_projection == 0 - stripped_transition = strip_projections(transition) + # attach projections to the interactions, as a spin-projection extension would + transition_with_projections = transition.convert( + interaction_converter=lambda interaction: attrs.evolve( + interaction, l_projection=0, s_projection=interaction.s_magnitude + ) + ) + stripped_transition = strip_projections(transition_with_projections) assert stripped_transition.states[3].name == resonance - assert stripped_transition.interactions[0].s_projection is None - assert stripped_transition.interactions[0].l_projection is None - assert stripped_transition.interactions[1].s_projection is None - assert stripped_transition.interactions[1].l_projection is None + for interaction in stripped_transition.interactions.values(): + assert interaction.l_projection is None + assert interaction.s_projection is None + assert interaction.l_magnitude is not None + assert interaction.s_magnitude is not None diff --git a/tests/unit/io/test_mermaid.py b/tests/unit/io/test_mermaid.py index bee8182c..78c359d4 100644 --- a/tests/unit/io/test_mermaid.py +++ b/tests/unit/io/test_mermaid.py @@ -161,12 +161,11 @@ def it_latex_reaction(reaction: ReactionInfo): ) assert src.startswith("flowchart LR\n") assert not src.startswith("```mermaid") - assert R"J/\psi(1S)\left[" in src - assert R"f_{0}(980)\left[" in src - assert R"P = \text{+}1" in src + assert R"J/\psi(1S)$$" in src + assert R"f_{0}(980)$$" in src assert "
" not in src if reaction.formalism == "canonical-helicity": - assert R"$$\begin{gathered} L =" in src + assert R"$$\begin{gathered} L = 0" in src labeled_lines = [ line for line in src.splitlines() if '["' in line or '("' in line @@ -206,7 +205,7 @@ def it_reaction(reaction: ReactionInfo): initial_state_id = next(iter(transition.topology.incoming_edge_ids)) initial_node_id = transition.topology.edges[initial_state_id].ending_node_id initial_state = transition.states[initial_state_id] - assert f'N{initial_node_id}["{initial_state.particle.name}' in src + assert f'N{initial_node_id}["{initial_state.name}' in src assert f" A --- N{initial_node_id}" not in src src = io.asmermaid(reaction, latex=False) assert src.startswith("flowchart LR\n") @@ -221,9 +220,10 @@ def it_reaction(reaction: ReactionInfo): def it_reaction_with_node_labels(reaction: ReactionInfo): src = io.asmermaid(reaction.transitions[0], render_node=True, latex=False) assert src.startswith("flowchart LR\n") - assert "gamma[-1]" in src - assert "f(0)(980)[0]" in src - assert "P=+1" in src + assert "gamma" in src + assert "f(0)(980)" in src + if reaction.formalism == "canonical-helicity": + assert "L=0" in src assert " A --- N0" in src def it_latex_can_be_disabled(reaction: ReactionInfo): @@ -335,11 +335,11 @@ def it_qn_result(qn_problem_and_result: tuple[QNProblemSet, QNResult]): src = io.asmermaid(qn_result, render_node=True, latex=False) assert src.startswith("flowchart LR\n") assert " --- " in src - assert "parity_prefactor =" in src + assert "l_magnitude =" in src src = io.asmermaid(qn_result, render_node=True, latex=True) assert R"$$\begin{gathered}" in src - assert R"\text{parity\_prefactor} = \text{+}1" in src + assert R"\text{l\_magnitude} =" in src @pytest.mark.parametrize( "formalism", @@ -347,7 +347,7 @@ def it_qn_result(qn_problem_and_result: tuple[QNProblemSet, QNResult]): ) def it_problemset(formalism: SpinFormalism): stm = StateTransitionManager( - initial_state=[("J/psi(1S)", [+1])], + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], formalism=formalism, ) diff --git a/tests/unit/test_combinatorics.py b/tests/unit/test_combinatorics.py index bcbdd29c..84e2a960 100644 --- a/tests/unit/test_combinatorics.py +++ b/tests/unit/test_combinatorics.py @@ -22,38 +22,26 @@ def three_body_decay() -> Topology: def test_create_initial_facts(three_body_decay, particle_database): initial_facts = create_initial_facts( three_body_decay, - initial_state=[("J/psi(1S)", [-1, +1])], + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], particle_db=particle_database, ) - assert len(initial_facts) == 4 - - for fact in initial_facts: - edge_ids = sorted(fact.states) - assert edge_ids == [-1, 0, 1, 2] - particle_names = [fact.states[i][0].name for i in edge_ids] - assert particle_names == ["J/psi(1S)", "gamma", "pi0", "pi0"] - _, initial_polarization = fact.states[-1] - assert initial_polarization in {-1, +1} + edge_ids = sorted(initial_facts.states) + assert edge_ids == [-1, 0, 1, 2] + particle_names = [initial_facts.states[i].name for i in edge_ids] + assert particle_names == ["J/psi(1S)", "gamma", "pi0", "pi0"] -def test_create_initial_facts_without_spin_projections( +def test_create_initial_facts_wrong_number_of_states_raises( three_body_decay, particle_database ): - initial_facts = create_initial_facts( - three_body_decay, - initial_state=[("J/psi(1S)", [-1, +1])], - final_state=["gamma", "pi0", "pi0"], - particle_db=particle_database, - expand_spin_projections=False, - ) - assert len(initial_facts) == 1 - fact = initial_facts[0] - edge_ids = sorted(fact.states) - assert edge_ids == [-1, 0, 1, 2] - particle_names = [fact.states[i][0].name for i in edge_ids] - assert particle_names == ["J/psi(1S)", "gamma", "pi0", "pi0"] - assert all(projection is None for _, projection in fact.states.values()) + with pytest.raises(ValueError, match="Number of state definitions"): + create_initial_facts( + three_body_decay, + initial_state=["J/psi(1S)"], + final_state=["gamma", "pi0"], + particle_db=particle_database, + ) def describe_generate_kinematic_permutations(): @@ -150,7 +138,7 @@ def test_permutate_topology_kinematically( ): permutations = permutate_topology_kinematically( topology=three_body_decay, - initial_state=[("J/psi(1S)", [-1, +1])], + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], final_state_groupings=final_state_groupings, ) diff --git a/tests/unit/test_final_state_permutations.py b/tests/unit/test_final_state_permutations.py index 8dedb332..462bb608 100644 --- a/tests/unit/test_final_state_permutations.py +++ b/tests/unit/test_final_state_permutations.py @@ -24,7 +24,7 @@ def test_create_problem_sets(final_state_description: str): for problem_set in problem_sets.values(): for problem in problem_set: problem_final_state = [ - problem.initial_facts.states[i][0].name for i in range(3) + problem.initial_facts.states[i].name for i in range(3) ] assert problem_final_state == input_final_state @@ -36,7 +36,7 @@ def test_create_problem_sets(final_state_description: str): def test_generate_transitions(final_state_description: str): final_state = final_state_description.split(" ") reaction = qrules.generate_transitions( - initial_state=("J/psi(1S)", [-1, +1]), + initial_state="J/psi(1S)", final_state=final_state, allowed_intermediate_particles=["omega(782)"], allowed_interaction_types=["strong", "EM"], @@ -46,17 +46,16 @@ def test_generate_transitions(final_state_description: str): ] assert final_state == ordered_final_state - assert len(reaction.transitions) == 8 + assert len(reaction.transitions) == 1 for transition in reaction.transitions: ordered_final_state = [ - transition.final_states[i].particle.name - for i in sorted(transition.final_states) + transition.final_states[i].name for i in sorted(transition.final_states) ] assert final_state == ordered_final_state topology = transition.topology decay_products = { - transition.states[i].particle.name + transition.states[i].name for i in topology.get_edge_ids_outgoing_from_node(1) } assert decay_products == {"gamma", "pi0"} diff --git a/tests/unit/test_parity_prefactor.py b/tests/unit/test_parity_prefactor.py deleted file mode 100644 index 90b9aa48..00000000 --- a/tests/unit/test_parity_prefactor.py +++ /dev/null @@ -1,77 +0,0 @@ -from __future__ import annotations - -from typing import NamedTuple - -import pytest - -from qrules.settings import InteractionType -from qrules.transition import StateTransitionManager - - -class Input(NamedTuple): - initial_state: list - final_state: list - intermediate_states: list - final_state_grouping: list - - -@pytest.mark.parametrize( - ( - "test_input", - "ingoing_state", - "relative_parity_prefactor", - ), - [ - ( - Input( - [("J/psi(1S)", [1])], - [("gamma", [-1, 1]), ("pi0", [0]), ("pi0", [0])], - ["f(0)(980)"], - ["pi0", "pi0"], - ), - "J/psi(1S)", - 1.0, - ), - ( - Input( - [("J/psi(1S)", [1])], - [("pi0", [0]), ("pi+", [0]), ("pi-", [0])], - ["rho(770)"], - ["pi+", "pi-"], - ), - "J/psi(1S)", - -1.0, - ), - ], -) -def test_parity_prefactor( - test_input: Input, - ingoing_state: str, - relative_parity_prefactor: float, -) -> None: - stm = StateTransitionManager( - test_input.initial_state, - test_input.final_state, - allowed_intermediate_particles=test_input.intermediate_states, - ) - stm.add_final_state_grouping(test_input.final_state_grouping) - stm.set_allowed_interaction_types([InteractionType.EM]) - problem_sets = stm.create_problem_sets() - - reaction = stm.find_solutions(problem_sets) - - assert len(reaction.group_by_topology()) == 1 - for transition in reaction.transitions: - in_edges = [ - state_id - for state_id, state in transition.states.items() - if state.particle.name == ingoing_state - ] - assert len(in_edges) == 1 - - node_id = transition.topology.edges[in_edges[0]].ending_node_id - assert isinstance(node_id, int) - assert ( - relative_parity_prefactor - == transition.interactions[node_id].parity_prefactor - ) diff --git a/tests/unit/test_qrules.py b/tests/unit/test_qrules.py index c314921d..c3b87e5d 100644 --- a/tests/unit/test_qrules.py +++ b/tests/unit/test_qrules.py @@ -1,6 +1,7 @@ import pytest from qrules import generate_transitions +from qrules.particle import create_particle, load_pdg @pytest.mark.parametrize( @@ -27,6 +28,36 @@ def test_generate_transitions(resonance_names): final_state = dict(enumerate(final_state_names)) for transition in reaction.transitions: this_final_state = { - i: state.particle.name for i, state in transition.final_states.items() + i: state.name for i, state in transition.final_states.items() } assert final_state == this_final_state + + +@pytest.mark.filterwarnings( + "ignore:There are conservation rules that were not executed" +) +def test_ls_free_solving_conserves_pair_c_parity(): + """Existence rules prune particle-antiparticle pairs without LS couplings.""" + particle_db = load_pdg() + # 1-+ cannot decay to p pbar: parity requires L=0 and C-parity then requires S=0 + exotic = create_particle( + particle_db["rho(770)0"], + name="X(1-+)", + pid=9999999, + mass=2.5, + width=0.1, + c_parity=+1, + g_parity=-1, + ) + particle_db.add(exotic) + for ls_couplings in (True, False): + with pytest.raises(RuntimeError, match="No solutions were found"): + generate_transitions( + initial_state="J/psi(1S)", + final_state=["gamma", "p", "p~"], + particle_db=particle_db, + allowed_intermediate_particles=[exotic.name], + allowed_interaction_types=["strong", "em"], + max_angular_momentum=1, + ls_couplings=ls_couplings, + ) diff --git a/tests/unit/test_settings.py b/tests/unit/test_settings.py index 119c438b..a83f11ee 100644 --- a/tests/unit/test_settings.py +++ b/tests/unit/test_settings.py @@ -15,7 +15,6 @@ if TYPE_CHECKING: from qrules.particle import ParticleCollection - from qrules.transition import SpinFormalism def describe_InteractionType(): @@ -45,32 +44,26 @@ def test_create_domains(particle_database: ParticleCollection): pdg = particle_database pions = pdg.filter(lambda p: p.name.startswith("pi")) domains = _create_domains(pions) - assert len(domains) == 15 + assert len(domains) == 14 + assert EdgeQN.spin_projection not in domains assert domains[EdgeQN.baryon_number] == [0] assert domains[EdgeQN.strangeness] == [0] assert domains[EdgeQN.charmness] == [0] assert domains[EdgeQN.bottomness] == [0] assert domains[EdgeQN.charge] == [-1, 0, +1] assert domains[EdgeQN.spin_magnitude] == [0, 0.5, 1, 1.5, 2] - assert ( - domains[EdgeQN.spin_projection] - == [-2, -1.5, -1, -0.5] + domains[EdgeQN.spin_magnitude] - ) assert domains[EdgeQN.isospin_magnitude] == [0, 0.5, 1] assert domains[EdgeQN.isospin_projection] == [-1, -0.5, 0, 0.5, 1] @pytest.mark.parametrize("interaction_type", list(InteractionType)) @pytest.mark.parametrize("nbody_topology", [False, True]) -@pytest.mark.parametrize("formalism", ["canonical", "canonical-helicity", "helicity"]) def test_create_interaction_settings( particle_database: ParticleCollection, interaction_type: InteractionType, nbody_topology: bool, - formalism: SpinFormalism, ): settings = create_interaction_settings( - formalism, particle_db=particle_database, nbody_topology=nbody_topology, ) @@ -89,7 +82,6 @@ def test_create_interaction_settings( "c_parity": [-1, +1, None], "g_parity": [-1, +1, None], "spin_magnitude": _halves_domain(0, 4), - "spin_projection": _halves_domain(-4, 4), "charge": _int_domain(-2, 2), "isospin_magnitude": _halves_domain(0, 1.5), "isospin_projection": _halves_domain(-1.5, +1.5), @@ -102,19 +94,9 @@ def test_create_interaction_settings( "l_magnitude": _int_domain(0, 2), "s_magnitude": _halves_domain(0, 2), } - if "canonical" in formalism: - expected["l_projection"] = [-2, -1, 0, 1, 2] - expected["s_projection"] = _halves_domain(-2, 2) - if formalism == "canonical-helicity": - expected["l_projection"] = [0] - if "helicity" in formalism and interaction_type != InteractionType.WEAK: - expected["parity_prefactor"] = [-1, 1] if nbody_topology: expected["l_magnitude"] = [0] expected["s_magnitude"] = [0] - if nbody_topology and formalism != "helicity": - expected["l_projection"] = [0] - expected["s_projection"] = [0] node_qn_domains_str = { # strings are easier to compare with pytest qn_type.__name__: domain for qn_type, domain in node_settings.qn_domains.items() diff --git a/tests/unit/test_solving.py b/tests/unit/test_solving.py index 634ea17a..c1acf4ed 100644 --- a/tests/unit/test_solving.py +++ b/tests/unit/test_solving.py @@ -1,16 +1,18 @@ from __future__ import annotations +from fractions import Fraction from typing import TYPE_CHECKING +import attrs import pytest import qrules.particle import qrules.system_control from qrules.conservation_rules import ( c_parity_conservation, + helicity_conservation, parity_conservation, spin_magnitude_conservation, - spin_validity, ) from qrules.quantum_numbers import EdgeQuantumNumbers, NodeQuantumNumbers from qrules.solving import ( @@ -19,7 +21,9 @@ complete_intermediate_states, filter_quantum_number_problem_set, merge_qn_problem_sets, + remove_dominated_qn_problem_sets, ) +from qrules.topology import MutableTransition if TYPE_CHECKING: from qrules.argument_handling import GraphEdgePropertyMap @@ -35,35 +39,27 @@ def it_finds_solutions( result = complete_intermediate_states( qn_result, quantum_number_problem_set, all_particles ) - assert len(result.solutions) == 19 + assert len(result.solutions) == 24 - @pytest.mark.parametrize("with_spin_projection", [True, False]) def it_with_filtered_quantum_number_problem_set( all_particles: list[GraphEdgePropertyMap], quantum_number_problem_set: QNProblemSet, - with_spin_projection: bool, ) -> None: solver = CSPSolver() - parametrized_edge_properties_and_domains = { - EdgeQuantumNumbers.pid, # had to be added for c_parity_conservation to work - EdgeQuantumNumbers.spin_magnitude, - EdgeQuantumNumbers.parity, - EdgeQuantumNumbers.c_parity, - } - if with_spin_projection: - parametrized_edge_properties_and_domains.add( - EdgeQuantumNumbers.spin_projection - ) - new_quantum_number_problem_set = filter_quantum_number_problem_set( quantum_number_problem_set, - edge_rules={spin_validity}, + edge_rules=set(), node_rules={ spin_magnitude_conservation, parity_conservation, c_parity_conservation, }, - edge_properties=parametrized_edge_properties_and_domains, + edge_properties={ + EdgeQuantumNumbers.pid, # had to be added for c_parity_conservation to work + EdgeQuantumNumbers.spin_magnitude, + EdgeQuantumNumbers.parity, + EdgeQuantumNumbers.c_parity, + }, node_properties=( NodeQuantumNumbers.l_magnitude, NodeQuantumNumbers.s_magnitude, @@ -73,33 +69,72 @@ def it_with_filtered_quantum_number_problem_set( result = complete_intermediate_states( qn_result, new_quantum_number_problem_set, all_particles ) + assert len(result.solutions) == 127 + + +def test_merge_qn_problem_sets(quantum_number_problem_set: QNProblemSet) -> None: + """Merge problem sets that differ only in custom spin-projection facts.""" + + def with_initial_projection( + problem_set: QNProblemSet, projection: Fraction + ) -> QNProblemSet: + facts = problem_set.initial_facts + new_states = {i: dict(m) for i, m in facts.states.items()} + new_states[-1][EdgeQuantumNumbers.spin_projection] = projection + new_facts = MutableTransition( + facts.topology, + new_states, # type: ignore[arg-type] + dict(facts.interactions), # type: ignore[arg-type] + ) + return attrs.evolve(problem_set, initial_facts=new_facts) - if with_spin_projection: - assert len(result.solutions) == 319 - else: - assert len(result.solutions) == 127 - - -def test_merge_qn_problem_sets() -> None: - stm = qrules.StateTransitionManager( - initial_state=["psi(2S)"], - final_state=["gamma", "eta", "eta"], - formalism="helicity", - ) - problem_sets = stm.create_problem_sets() - qn_problem_sets = [ - p.to_qn_problem_set() - for strength in sorted(problem_sets) - for p in problem_sets[strength] + projections = [Fraction(-1), Fraction(+1)] + variants = [ + with_initial_projection(quantum_number_problem_set, p) for p in projections ] - merged_problem_sets = merge_qn_problem_sets(qn_problem_sets) - assert len(merged_problem_sets) < len(qn_problem_sets) - - def n_total_solutions(problem_sets: list[QNProblemSet]) -> int: - solver = CSPSolver() - return sum(len(solver.find_solutions(p).solutions) for p in problem_sets) - - assert n_total_solutions(merged_problem_sets) == n_total_solutions(qn_problem_sets) + merged = merge_qn_problem_sets( + variants, merge_qns={EdgeQuantumNumbers.spin_projection} + ) + assert len(merged) == 1 + merged_facts = merged[0].initial_facts.states[-1] + assert merged_facts[EdgeQuantumNumbers.spin_projection] == projections + + +def test_remove_dominated_qn_problem_sets( + quantum_number_problem_set: QNProblemSet, +) -> None: + def with_extra_node_rule(problem_set: QNProblemSet) -> QNProblemSet: + settings = problem_set.solving_settings + new_interactions = { + node_id: attrs.evolve( + node_settings, + conservation_rules={ + **node_settings.conservation_rules, + helicity_conservation: 1, + }, + ) + for node_id, node_settings in settings.interactions.items() + } + new_settings = MutableTransition( + settings.topology, + dict(settings.states), # type: ignore[arg-type] + new_interactions, # type: ignore[arg-type] + ) + return attrs.evolve(problem_set, solving_settings=new_settings) + + stricter = with_extra_node_rule(quantum_number_problem_set) + deduped = remove_dominated_qn_problem_sets({ + 1.0: [quantum_number_problem_set], + 60.0: [stricter], + }) + assert deduped == {1.0: [quantum_number_problem_set]} + + exact_duplicate = attrs.evolve(quantum_number_problem_set) + deduped = remove_dominated_qn_problem_sets({ + 1.0: [quantum_number_problem_set], + 60.0: [exact_duplicate], + }) + assert sum(len(group) for group in deduped.values()) == 1 @pytest.fixture(scope="session") diff --git a/tests/unit/test_system_control.py b/tests/unit/test_system_control.py index 5f8fcd04..2904c13a 100644 --- a/tests/unit/test_system_control.py +++ b/tests/unit/test_system_control.py @@ -1,7 +1,6 @@ from __future__ import annotations from copy import deepcopy -from fractions import Fraction from importlib.metadata import version import attrs @@ -13,7 +12,7 @@ match_external_edges, perform_external_edge_identical_particle_combinatorics, ) -from qrules.particle import Particle, ParticleWithSpin +from qrules.particle import Particle from qrules.quantum_numbers import ( EdgeQuantumNumbers, InteractionProperties, @@ -37,56 +36,44 @@ ), [ ( - [("Y(4260)", [-1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], + ["Y(4260)"], + ["D0", "D~0", "pi0", "pi0"], [[["D0", "pi0"], ["D~0", "pi0"]]], 1, ), ( - [("Y(4260)", [-1, 1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], - [[["D0", "pi0"], ["D~0", "pi0"]]], - 2, - ), - ( - [("Y(4260)", [1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], + ["Y(4260)"], + ["D0", "D~0", "pi0", "pi0"], [], 9, ), ( - [("Y(4260)", [-1, 1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], - [], - 18, - ), - ( - [("Y(4260)", [1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], + ["Y(4260)"], + ["D0", "D~0", "pi0", "pi0"], [[["D0", "pi0"], ["D~0", "pi0"]], ["D0", "pi0"]], 3, ), ( - [("J/psi(1S)", [-1, 1])], - [("gamma", [-1, 1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], [["pi0", "pi0"]], - 4, + 1, ), ( - [("J/psi(1S)", [-1, 1])], - [("gamma", [-1, 1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], [["pi0", "gamma"]], - 4, + 1, ), ( - [("J/psi(1S)", [-1, 1])], - [("gamma", [-1, 1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], [], - 8, + 2, ), ( - [("J/psi(1S)", [-1, 1])], - [("gamma", [-1, 1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], [["pi0", "pi-"]], 0, ), @@ -136,17 +123,15 @@ def __get_f2_1270_pos() -> tuple[float, float]: @pytest.mark.parametrize( - ("particle_name", "spin_projection", "expected_properties"), + ("particle_name", "expected_properties"), [ ( "pi0", - 0, { EdgeQuantumNumbers.pid: 111, EdgeQuantumNumbers.mass: 0.1349768, EdgeQuantumNumbers.width: get_pi0_width(), EdgeQuantumNumbers.spin_magnitude: 0.0, - EdgeQuantumNumbers.spin_projection: 0, EdgeQuantumNumbers.charge: 0, EdgeQuantumNumbers.isospin_magnitude: 1.0, EdgeQuantumNumbers.isospin_projection: 0.0, @@ -165,13 +150,11 @@ def __get_f2_1270_pos() -> tuple[float, float]: ), ( "D+", # no g and c parity - 0, { EdgeQuantumNumbers.pid: 411, EdgeQuantumNumbers.mass: __get_d_pos()[0], EdgeQuantumNumbers.width: __get_d_pos()[1], EdgeQuantumNumbers.spin_magnitude: 0.0, - EdgeQuantumNumbers.spin_projection: 0, EdgeQuantumNumbers.charge: 1, EdgeQuantumNumbers.isospin_magnitude: 0.5, EdgeQuantumNumbers.isospin_projection: 0.5, @@ -189,14 +172,12 @@ def __get_f2_1270_pos() -> tuple[float, float]: }, ), ( - "f(2)(1270)", # spin projection 1 - 1.0, + "f(2)(1270)", { EdgeQuantumNumbers.pid: 225, EdgeQuantumNumbers.mass: __get_f2_1270_pos()[0], EdgeQuantumNumbers.width: __get_f2_1270_pos()[1], EdgeQuantumNumbers.spin_magnitude: 2.0, - EdgeQuantumNumbers.spin_projection: 1.0, EdgeQuantumNumbers.charge: 0, EdgeQuantumNumbers.isospin_magnitude: 0.0, EdgeQuantumNumbers.isospin_projection: 0.0, @@ -217,13 +198,12 @@ def __get_f2_1270_pos() -> tuple[float, float]: ) def test_create_edge_properties( particle_name, - spin_projection, expected_properties, particle_database, skh_particle_version: str, ): particle = particle_database[particle_name] - assert create_edge_properties(particle, spin_projection) == expected_properties + assert create_edge_properties(particle) == expected_properties assert skh_particle_version is not None # dummy for skip tests @@ -240,7 +220,7 @@ def make_ls_test_graph( l_magnitude=angular_momentum_magnitude, ) } - states: dict[int, ParticleWithSpin] = {-1: (particle, Fraction(0))} + states: dict[int, Particle] = {-1: particle} return MutableTransition(topology, states, interactions) @@ -257,7 +237,7 @@ def make_ls_test_graph_scrambled( s_magnitude=coupled_spin_magnitude, ) } - states: dict[int, ParticleWithSpin] = {-1: (particle, Fraction(0))} + states: dict[int, Particle] = {-1: particle} return MutableTransition(topology, states, interactions) @@ -340,12 +320,7 @@ def it_filter_graphs_for_interaction_qns( tempgraph = make_ls_test_graph(value[1][0], value[1][1], pi0) tempgraph = attrs.evolve( tempgraph, - states={ - -1: ( - Particle(name=value[0], pid=0, mass=1.0, spin=1.0), - 0.0, - ) - }, + states={-1: Particle(name=value[0], pid=0, mass=1.0, spin=1.0)}, ) graphs.append(tempgraph) @@ -356,7 +331,7 @@ def it_filter_graphs_for_interaction_qns( def _create_graph( problem_set: ProblemSet, -) -> MutableTransition[ParticleWithSpin, InteractionProperties]: +) -> MutableTransition[Particle, InteractionProperties]: return MutableTransition( topology=problem_set.topology, interactions=problem_set.initial_facts.interactions, @@ -368,8 +343,8 @@ def _create_graph( ("initial_state", "final_state"), [ ( - [("Y(4260)", [-1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], + ["Y(4260)"], + ["D0", "D~0", "pi0", "pi0"], ), ], ) @@ -383,7 +358,7 @@ def test_edge_swap(particle_database, initial_state, final_state): stm.set_allowed_interaction_types([InteractionType.STRONG]) problem_sets = stm.create_problem_sets() - init_graphs: list[MutableTransition[ParticleWithSpin, InteractionProperties]] = [] + init_graphs: list[MutableTransition[Particle, InteractionProperties]] = [] for problem_set_list in problem_sets.values(): init_graphs.extend([_create_graph(x) for x in problem_set_list]) @@ -409,12 +384,12 @@ def test_edge_swap(particle_database, initial_state, final_state): ("initial_state", "final_state"), [ ( - [("Y(4260)", [-1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], + ["Y(4260)"], + ["D0", "D~0", "pi0", "pi0"], ), ( - [("J/psi(1S)", [-1, 1])], - [("gamma", [-1, 1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], ), ], ) @@ -429,7 +404,7 @@ def test_match_external_edges(particle_database, initial_state, final_state): stm.set_allowed_interaction_types([InteractionType.STRONG]) problem_sets = stm.create_problem_sets() - init_graphs: list[MutableTransition[ParticleWithSpin, InteractionProperties]] = [] + init_graphs: list[MutableTransition[Particle, InteractionProperties]] = [] for problem_set_list in problem_sets.values(): init_graphs.extend([_create_graph(x) for x in problem_set_list]) @@ -462,32 +437,32 @@ def test_match_external_edges(particle_database, initial_state, final_state): ), [ ( - [("Y(4260)", [1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], + ["Y(4260)"], + ["D0", "D~0", "pi0", "pi0"], [[["D0", "pi0"], ["D~0", "pi0"]]], 2, ), ( - [("Y(4260)", [1])], - [("D0", [0]), ("D~0", [0]), ("pi0", [0]), ("pi0", [0])], + ["Y(4260)"], + ["D0", "D~0", "pi0", "pi0"], [["D0", "pi0"]], 6, ), ( - [("J/psi(1S)", [1])], - [("gamma", [1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], [["pi0", "pi0"]], 1, ), ( - [("J/psi(1S)", [-1, 1])], - [("gamma", [-1, 1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], [], - 12, + 3, ), ( - [("J/psi(1S)", [1])], - [("gamma", [1]), ("pi0", [0]), ("pi0", [0])], + ["J/psi(1S)"], + ["gamma", "pi0", "pi0"], [["pi0", "gamma"]], 2, ), @@ -518,7 +493,7 @@ def test_external_edge_identical_particle_combinatorics( match_external_edges(init_graphs) - comb_graphs: list[MutableTransition[ParticleWithSpin, InteractionProperties]] = [] + comb_graphs: list[MutableTransition[Particle, InteractionProperties]] = [] for group in init_graphs: comb_graphs.extend( perform_external_edge_identical_particle_combinatorics(group) diff --git a/tests/unit/test_transition.py b/tests/unit/test_transition.py index 40789773..e066ee9e 100644 --- a/tests/unit/test_transition.py +++ b/tests/unit/test_transition.py @@ -17,7 +17,7 @@ MutableTransition, Topology, ) -from qrules.transition import ReactionInfo, SolvingMode, State, StateTransitionManager +from qrules.transition import ReactionInfo, SolvingMode, StateTransitionManager NAMESPACE_WITH_FRACTIONS = globals() NAMESPACE_WITH_FRACTIONS["Fraction"] = Fraction @@ -26,7 +26,7 @@ def describe_MutableTransition(): def it_intermediate_states(): stm = StateTransitionManager( - initial_state=[("J/psi(1S)", [-1, +1])], + initial_state=["J/psi(1S)"], final_state=["K0", "Sigma+", "p~"], allowed_intermediate_particles=["N(1700)", "Sigma(1750)"], formalism="helicity", @@ -48,9 +48,9 @@ def it_properties(reaction: ReactionInfo): assert reaction.final_state[2].name == "pi0" assert len(reaction.group_by_topology()) == 1 if reaction.formalism.startswith("cano"): - assert len(reaction.transitions) == 16 + assert len(reaction.transitions) == 4 else: - assert len(reaction.transitions) == 8 + assert len(reaction.transitions) == 2 for transition in reaction.transitions: assert isinstance(transition, FrozenTransition) @@ -66,49 +66,27 @@ def it_hash(reaction: ReactionInfo): def it_hash_value(reaction: ReactionInfo): if sys.version_info >= (3, 11) and not sys.version_info >= (3, 14): expected_hash = { - "canonical-helicity": "65106a44301f9340e633d09f66ad7d17", - "helicity": "9646d3ee5c5e8534deb8019435161f2e", + "canonical-helicity": "75f6d331aceefda11d14e61bea24b076", + "helicity": "9a5dad546caec7d2873ba7c745e8b321", }[reaction.formalism] elif sys.version_info >= (3, 14): expected_hash = { - "canonical-helicity": "762cc006a8c4c0a0a88fce934a32577d", - "helicity": "17fefe55a7da0810371e90bd762a176a", + "canonical-helicity": "9e6b4b5ea854785ac33e6049b6ab86d1", + "helicity": "a6880e15ca68c5d728574ffc3b5c59e4", }[reaction.formalism] else: expected_hash = { - "canonical-helicity": "0d8bc378677986e0dc2d3b02f5627e0b", - "helicity": "71404ad43550850a02109e8db044bd28", + "canonical-helicity": "995925fb32a00be871211bf9d3ac78bc", + "helicity": "1c1441c16cfbec426239b9c0e92d19c9", }[reaction.formalism] assert _compute_hash(reaction) == expected_hash -def describe_State(): - @pytest.mark.parametrize( - ("state_def_1", "state_def_2"), - [ - (("a", -1), ("a", +1)), - (("a", 0), ("a", 0)), - (("a", 0), ("b", 0)), - (("a", -1), ("b", +1)), - ], - ) - def it_ordering(state_def_1, state_def_2): - def create_state(state_def) -> State: - return State( - particle=Particle(name=state_def[0], pid=0, spin=0, mass=0), - spin_projection=state_def[1], - ) - - state1 = create_state(state_def_1) - state2 = create_state(state_def_2) - assert state2 >= state1 - - def describe_StateTransitionManager(): def it_allowed_intermediate_particles(): stm = StateTransitionManager( - initial_state=[("J/psi(1S)", list(map(Fraction, [-1, +1])))], + initial_state=["J/psi(1S)"], final_state=["p", "p~", "eta"], ) particle_name = "N(753)" @@ -142,8 +120,8 @@ def count_transitions(solving_mode: SolvingMode) -> int: reaction = stm.find_solutions(stm.create_problem_sets()) return len(reaction.transitions) - assert count_transitions(SolvingMode.FULL) == 294 - assert count_transitions(SolvingMode.FAST) == 90 + assert count_transitions(SolvingMode.FULL) == 32 + assert count_transitions(SolvingMode.FAST) == 10 def it_regex_pattern(): stm = StateTransitionManager( diff --git a/tests/unit/test_workflow.py b/tests/unit/test_workflow.py index 476c8fc8..4720f930 100644 --- a/tests/unit/test_workflow.py +++ b/tests/unit/test_workflow.py @@ -2,17 +2,23 @@ from fractions import Fraction from typing import Any +import attrs import pytest +from qrules.argument_handling import GraphEdgePropertyMap +from qrules.conservation_rules import helicity_conservation, spin_validity from qrules.io import asdict, asdot, asmermaid from qrules.particle import ParticleCollection, load_pdg from qrules.quantum_numbers import EdgeQuantumNumbers from qrules.settings import ( + CONSERVATION_LAW_PRIORITIES, DEFAULT_INTERACTION_TYPES, + EDGE_RULE_PRIORITIES, InteractionType, create_interaction_settings, ) -from qrules.solving import _create_merge_key +from qrules.solving import QNProblemSet +from qrules.topology import MutableTransition from qrules.transition import ReactionInfo, SolvingMode from qrules.workflow import ( InteractionConfig, @@ -23,7 +29,6 @@ find_qn_transitions, find_solutions, generate_qn_transitions, - strip_spin_projections, ) @@ -111,17 +116,15 @@ def count_transitions(solving_mode: SolvingMode) -> int: ) return len(reaction.transitions) - assert count_transitions(SolvingMode.FULL) == 294 - assert count_transitions(SolvingMode.FAST) == 90 + assert count_transitions(SolvingMode.FULL) == 32 + assert count_transitions(SolvingMode.FAST) == 10 def describe_InteractionConfig(): @pytest.fixture def config(particle_database: ParticleCollection) -> InteractionConfig: return InteractionConfig( - type_settings=create_interaction_settings( - "helicity", particle_db=particle_database - ) + type_settings=create_interaction_settings(particle_db=particle_database) ) def it_allows_all_types_by_default(config: InteractionConfig): @@ -167,26 +170,21 @@ def test_initial_state_restricts_interaction_types( assert sorted(qn_problem_sets.problem_sets) == expected_strengths -@pytest.mark.parametrize("merge_spin_projections", [False, True]) -def test_pipeline_reproduces_state_transition_manager( - reaction: ReactionInfo, merge_spin_projections: bool -): +def test_pipeline_reproduces_state_transition_manager(reaction: ReactionInfo): particle_db = load_pdg() qn_problem_sets = create_qn_problem_sets( - initial_state=[("J/psi(1S)", [-1, 1])], + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], particle_db=particle_db, allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"], interaction_config=InteractionConfig( type_settings=create_interaction_settings( - reaction.formalism, particle_db=particle_db, max_angular_momentum=2, ), allowed_types=[InteractionType.STRONG], ), formalism=reaction.formalism, - merge_spin_projections=merge_spin_projections, ) assert isinstance(qn_problem_sets, QNProblemSetCollection) assert qn_problem_sets.formalism == reaction.formalism @@ -200,27 +198,21 @@ def test_pipeline_reproduces_state_transition_manager( assert workflow_reaction == reaction -def test_projection_free_qn_transitions(): +def test_qn_transitions_are_projection_free(): particle_db = load_pdg() collection = create_qn_problem_sets( - initial_state=[("J/psi(1S)", [-1, 1])], + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], particle_db=particle_db, allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"], interaction_config=InteractionConfig( type_settings=create_interaction_settings( - "helicity", particle_db=particle_db, max_angular_momentum=2 + particle_db=particle_db, max_angular_momentum=2 ), allowed_types=[InteractionType.STRONG], ), ) - stripped = strip_spin_projections(collection) - assert isinstance(stripped, QNProblemSetCollection) - n_original = sum(map(len, collection.problem_sets.values())) - n_stripped = sum(map(len, stripped.problem_sets.values())) - assert n_stripped < n_original - - qn_transitions = find_qn_transitions(stripped) + qn_transitions = find_qn_transitions(collection) assert len(qn_transitions) > 0 qn_names = { qn_type.__name__ @@ -251,29 +243,138 @@ def test_projection_free_qn_transitions(): assert '"spin_projection"' not in serialized assert '"spin_magnitude"' in serialized - unexpanded = create_qn_problem_sets( - initial_state=[("J/psi(1S)", [-1, 1])], + +def test_spin_projections_reenabled_through_problem_sets(): + """Spin projections remain available as an extension of the QN problem sets. + + Spin projections are no longer part of the default workflow, but the `.CSPSolver` + handles any quantum number that the problem sets declare through facts and + domains. This test re-enables them for :math:`J/\\psi \\to \\gamma f_2(1270)` by + adding `~.EdgeQuantumNumbers.spin_projection` facts to the external edges (a + `list` fact is solved as a variable over that range), a projection domain plus + `.spin_validity` to the intermediate edges, and `.helicity_conservation` to the + interaction nodes, which prunes the helicity combinations to + :math:`|\\lambda_\\gamma - \\lambda_{f_2}| \\leq 1`. + """ + particle_db = load_pdg() + collection = create_qn_problem_sets( + initial_state=["J/psi(1S)"], final_state=["gamma", "pi0", "pi0"], particle_db=particle_db, - allowed_intermediate_particles=["f(0)(980)", "f(0)(1500)"], - interaction_config=InteractionConfig( - type_settings=create_interaction_settings( - "helicity", particle_db=particle_db, max_angular_momentum=2 - ), - allowed_types=[InteractionType.STRONG], - ), - spin_projections=False, + allowed_intermediate_particles=["f(2)(1270)"], + allowed_interaction_types=["strong", "EM"], + max_angular_momentum=2, + final_state_groupings=[[["pi0", "pi0"]]], ) - assert _to_merge_keys(unexpanded) == _to_merge_keys(stripped) - assert find_qn_transitions(unexpanded) == qn_transitions + spin_projections: dict[str, Fraction | list[Fraction]] = { + "J/psi(1S)": [Fraction(-1), Fraction(+1)], # from e+e- collision + "gamma": [Fraction(-1), Fraction(+1)], + "pi0": Fraction(0), + } + def get_particle_name(prop_map: GraphEdgePropertyMap) -> str: + pid = prop_map[EdgeQuantumNumbers.pid] + assert pid is not None + return particle_db.find(int(pid)).name + + def add_spin_projections(problem_set: QNProblemSet) -> QNProblemSet: + facts = problem_set.initial_facts + new_states = { + edge_id: { + **prop_map, + EdgeQuantumNumbers.spin_projection: spin_projections[ + get_particle_name(prop_map) + ], + } + for edge_id, prop_map in facts.states.items() + } + new_facts = MutableTransition( + facts.topology, + new_states, # type: ignore[arg-type] + dict(facts.interactions), # type: ignore[arg-type] + ) + settings = problem_set.solving_settings + new_edge_settings = {} + for edge_id, edge_settings in settings.states.items(): + if edge_id not in facts.topology.intermediate_edge_ids: + new_edge_settings[edge_id] = edge_settings + continue + max_spin = max(edge_settings.qn_domains[EdgeQuantumNumbers.spin_magnitude]) + projection_domain = [ + Fraction(x, 2) for x in range(-int(2 * max_spin), int(2 * max_spin) + 1) + ] + new_edge_settings[edge_id] = attrs.evolve( + edge_settings, + conservation_rules={ + **edge_settings.conservation_rules, + spin_validity: EDGE_RULE_PRIORITIES[spin_validity], + }, + qn_domains={ + **edge_settings.qn_domains, + EdgeQuantumNumbers.spin_projection: projection_domain, + }, + ) + new_node_settings = { + node_id: attrs.evolve( + node_settings, + conservation_rules={ + **node_settings.conservation_rules, + helicity_conservation: CONSERVATION_LAW_PRIORITIES[ + helicity_conservation + ], + }, + ) + for node_id, node_settings in settings.interactions.items() + } + new_settings = MutableTransition( + settings.topology, + new_edge_settings, # type: ignore[arg-type] + new_node_settings, # type: ignore[arg-type] + ) + return QNProblemSet(initial_facts=new_facts, solving_settings=new_settings) # type: ignore[arg-type] -def _to_merge_keys(collection: QNProblemSetCollection) -> set[tuple]: - return { - (strength, _create_merge_key(problem_set, set())) + collection.problem_sets = { + strength: [add_spin_projections(p) for p in problem_sets] for strength, problem_sets in collection.problem_sets.items() - for problem_set in problem_sets } + qn_transitions = find_qn_transitions(collection) + assert len(qn_transitions) == 36 + assert all( + EdgeQuantumNumbers.spin_projection in state + for transition in qn_transitions + for state in transition.states.values() + ) + + def get_projection(transition, edge_id: int) -> Fraction: + return transition.states[edge_id][EdgeQuantumNumbers.spin_projection] + + helicity_combinations = set() + for transition in qn_transitions: + topology = transition.topology + resonance_edge = next(iter(topology.intermediate_edge_ids)) + gamma_edge = next( + i + for i in topology.outgoing_edge_ids + if transition.states[i][EdgeQuantumNumbers.pid] == 22 + ) + helicity_combinations.add(( + get_projection(transition, gamma_edge), + get_projection(transition, resonance_edge), + )) + assert helicity_combinations == { + (Fraction(-1), Fraction(-2)), + (Fraction(-1), Fraction(-1)), + (Fraction(-1), Fraction(0)), + (Fraction(+1), Fraction(0)), + (Fraction(+1), Fraction(+1)), + (Fraction(+1), Fraction(+2)), + } + initial_state_projections = { + get_projection(transition, edge_id) + for transition in qn_transitions + for edge_id in transition.topology.incoming_edge_ids + } + assert initial_state_projections == {Fraction(-1), Fraction(+1)} def test_generate_qn_transitions(): @@ -345,6 +446,43 @@ def test_generate_qn_transitions_two_to_n(): assert len(reaction.group_by_topology()) > 1 +def test_generate_qn_transitions_without_ls_couplings(): + """LS-free solving must allow the same intermediate states (ComPWA/qrules#19).""" + particle_db = load_pdg() + reaction_kwargs: dict[str, Any] = dict( + initial_state=["gamma", "p"], + final_state=["p", "pi0"], + particle_db=particle_db, + allowed_intermediate_particles=[ + "Delta(1232)", + "N(1440)", + "rho(770)", + "omega(782)", + ], + allowed_interaction_types=["strong", "em"], + max_angular_momentum=2, + ) + signatures_by_mode = {} + n_transitions_by_mode = {} + for ls_couplings in (True, False): + reaction = generate_qn_transitions(ls_couplings=ls_couplings, **reaction_kwargs) + n_transitions_by_mode[ls_couplings] = len(reaction.transitions) + signatures_by_mode[ls_couplings] = { + ( + state[EdgeQuantumNumbers.spin_magnitude], + state[EdgeQuantumNumbers.parity], + ) + for transition in reaction.transitions + for state in transition.intermediate_states.values() + } + if not ls_couplings: + for transition in reaction.transitions: + for interactions in transition.interactions.values(): + assert not interactions + assert signatures_by_mode[False] == signatures_by_mode[True] + assert n_transitions_by_mode[False] < n_transitions_by_mode[True] + + def test_group_by_channel_and_channel_selection(): """Mandelstam channel encoding for 2-to-n reactions (ComPWA/qrules#29).""" particle_db = load_pdg() @@ -391,18 +529,7 @@ def test_qn_reaction_info_requires_particle_states(): final_state=["gamma", "pi0", "pi0"], particle_db=particle_db, allowed_intermediate_particles=["f(0)(980)"], - spin_projections=False, ) qn_transitions = find_qn_transitions(qn_problem_sets) with pytest.raises(TypeError, match="is of type FrozenDict, not Particle"): QNReactionInfo(qn_transitions) - - -def test_incompatible_spin_projection_flags_raise(): - with pytest.raises(ValueError, match="merge_spin_projections has no effect"): - create_qn_problem_sets( - initial_state=["J/psi(1S)"], - final_state=["gamma", "pi0", "pi0"], - merge_spin_projections=True, - spin_projections=False, - )