From 0a941e4d0c125852b9dcaba339c6003b603587f4 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 04:13:42 +0200
Subject: [PATCH 1/9] FEAT: build topologies for 2-to-n reactions

---
 src/qrules/topology.py      | 211 +++++++++++++++++++++++++++++++-----
 src/qrules/transition.py    |   2 +-
 src/qrules/workflow.py      |   7 +-
 tests/unit/test_topology.py |  45 +++++++-
 tests/unit/test_workflow.py |  32 ++++++
 5 files changed, 260 insertions(+), 37 deletions(-)

diff --git a/src/qrules/topology.py b/src/qrules/topology.py
index 846eb5da..bdf6db31 100644
--- a/src/qrules/topology.py
+++ b/src/qrules/topology.py
@@ -20,7 +20,7 @@
 import itertools
 import logging
 from abc import ABC, abstractmethod
-from collections import abc
+from collections import abc, defaultdict
 from functools import total_ordering
 from typing import TYPE_CHECKING, Any, Generic, TypeVar, overload
 
@@ -343,7 +343,11 @@ def get_originating_node_list(
     def __get_originating_node(edge_id: int) -> int | None:
         return topology.edges[edge_id].originating_node_id
 
-    return [node_id for node_id in map(__get_originating_node, edge_ids) if node_id]
+    return [
+        node_id
+        for node_id in map(__get_originating_node, edge_ids)
+        if node_id is not None
+    ]
 
 
 def _to_mutable_topology_nodes(inst: Iterable[int]) -> set[int]:
@@ -551,15 +555,21 @@ def build(
             active_graph_list = extendable_graph_list
             extendable_graph_list = []
             for active_graph in active_graph_list:
-                # check if finished
+                topology, open_end_edges = active_graph
+                connected_edge_groups = _group_connected_edge_ids(topology)
                 if (
-                    len(active_graph[1]) == number_of_final_edges
-                    and len(active_graph[0].nodes) > 0
+                    len(open_end_edges) == number_of_final_edges
+                    and len(topology.nodes) > 0
+                    and len(connected_edge_groups) == 1
                 ):
                     graph_tuple_list.append(active_graph)
                     continue
-
-                extendable_graph_list.extend(self._extend_graph(active_graph))
+                if self.__may_still_reach_final_state(
+                    open_end_edges, len(connected_edge_groups), number_of_final_edges
+                ):
+                    extendable_graph_list.extend(
+                        self._extend_graph(active_graph, connected_edge_groups)
+                    )
 
         _LOGGER.info("finished building topology graphs...")
         # strip the current open end edges list from the result graph tuples
@@ -570,50 +580,177 @@ def build(
             topologies.append(topology.freeze())
         return tuple(topologies)
 
+    def __may_still_reach_final_state(
+        self,
+        open_end_edges: Sequence[int],
+        number_of_connected_groups: int,
+        number_of_final_edges: int,
+    ) -> bool:
+        r"""Check whether extending a graph can still lead to a valid final state.
+
+        Only interaction nodes with more than one ingoing edge can reduce the number
+        of open edges, and — since attaching such a node to already-connected edges is
+        forbidden (see `._extend_graph`) — each one merges at least two of the
+        graph's connected groups. Graphs whose open edges exceed the number of final
+        edges by more than these merges can compensate are dead ends. Pruning them
+        keeps the recursive search finite when the interaction node set contains
+        :math:`2 \to 1` building blocks.
+        """
+        merging_nodes = [
+            node
+            for node in self.interaction_node_set
+            if node.number_of_ingoing_edges > 1
+        ]
+        if number_of_connected_groups > 1 and not merging_nodes:
+            return False
+        max_reduction_per_merge = max(
+            (
+                node.number_of_ingoing_edges - node.number_of_outgoing_edges
+                for node in merging_nodes
+            ),
+            default=0,
+        )
+        max_reduction = max(max_reduction_per_merge, 0) * (
+            number_of_connected_groups - 1
+        )
+        return len(open_end_edges) - max_reduction <= number_of_final_edges
+
     def _extend_graph(
-        self, pair: tuple[MutableTopology, Sequence[int]]
+        self,
+        pair: tuple[MutableTopology, Sequence[int]],
+        connected_edge_groups: list[set[int]],
     ) -> list[tuple[MutableTopology, list[int]]]:
         extended_graph_list: list[tuple[MutableTopology, list[int]]] = []
 
         topology, current_open_end_edges = pair
+        group_index = {
+            edge_id: i
+            for i, group in enumerate(connected_edge_groups)
+            for edge_id in group
+        }
 
         # Try to extend the graph with interaction nodes
         # that have equal or less ingoing lines than active lines
         for interaction_node in self.interaction_node_set:
-            if interaction_node.number_of_ingoing_edges <= len(current_open_end_edges):
-                # make all combinations
-                combis = list(
-                    itertools.combinations(
-                        current_open_end_edges,
-                        interaction_node.number_of_ingoing_edges,
-                    )
+            number_of_ingoing_edges = interaction_node.number_of_ingoing_edges
+            if number_of_ingoing_edges > len(current_open_end_edges):
+                continue
+            seen_originating_nodes: list[list[int]] = []
+            for combination in itertools.combinations(
+                current_open_end_edges, number_of_ingoing_edges
+            ):
+                groups = {group_index[edge_id] for edge_id in combination}
+                if len(groups) < number_of_ingoing_edges:
+                    # attaching already-connected edges to one node creates a cycle
+                    continue
+                originating_nodes = get_originating_node_list(
+                    topology, edge_ids=combination
+                )
+                # skip combinations that are symmetric to an earlier one
+                if originating_nodes in seen_originating_nodes:
+                    continue
+                seen_originating_nodes.append(originating_nodes)
+                extended_graph_list.append(
+                    _attach_node_to_edges(pair, interaction_node, combination)
                 )
-                # remove all combinations that originate from the same nodes
-                for comb1, comb2 in itertools.combinations(combis, 2):
-                    if get_originating_node_list(
-                        topology, edge_ids=comb1
-                    ) == get_originating_node_list(topology, edge_ids=comb2):
-                        combis.remove(comb2)
-
-                for combi in combis:
-                    new_graph = _attach_node_to_edges(pair, interaction_node, combi)
-                    extended_graph_list.append(new_graph)
 
         return extended_graph_list
 
 
+def _group_connected_edge_ids(topology: MutableTopology) -> list[set[int]]:
+    """Group the edge IDs of a topology into connected components.
+
+    Two edges are connected when they are attached to a common node. Edges that are
+    not attached to any node form a component of their own.
+
+    >>> topology = MutableTopology()
+    >>> topology.add_edges([0, 1, 2])
+    >>> topology.add_node(0)
+    >>> topology.attach_edges_to_node_ingoing([0, 1], node_id=0)
+    >>> _group_connected_edge_ids(topology)
+    [{0, 1}, {2}]
+    """
+    edge_ids_per_node: defaultdict[int, set[int]] = defaultdict(set)
+    for edge_id, edge in topology.edges.items():
+        for node_id in (edge.originating_node_id, edge.ending_node_id):
+            if node_id is not None:
+                edge_ids_per_node[node_id].add(edge_id)
+    components: list[set[int]] = []
+    remaining_edge_ids = set(topology.edges)
+    while remaining_edge_ids:
+        component = {remaining_edge_ids.pop()}
+        queue = list(component)
+        while queue:
+            edge = topology.edges[queue.pop()]
+            for node_id in (edge.originating_node_id, edge.ending_node_id):
+                if node_id is None:
+                    continue
+                new_neighbors = edge_ids_per_node[node_id] - component
+                component |= new_neighbors
+                queue += new_neighbors
+        remaining_edge_ids -= component
+        components.append(component)
+    return sorted(components, key=min)
+
+
+def _create_isomorphism_invariant(topology: Topology) -> tuple[tuple[int, int], ...]:
+    """Compute a canonical form of a topology under node and edge relabeling.
+
+    Two topologies have the same invariant if and only if they are isomorphic as
+    directed graphs, i.e. equal up to a relabeling of their nodes and edges. The
+    invariant is the lexicographically smallest sorted tuple of
+    ``(originating_node_id, ending_node_id)`` pairs over all node relabelings, with
+    ``-1`` for open edge ends.
+    """
+    nodes = sorted(topology.nodes)
+    return min(
+        tuple(
+            sorted(
+                (
+                    -1
+                    if e.originating_node_id is None
+                    else relabel[e.originating_node_id],
+                    -1 if e.ending_node_id is None else relabel[e.ending_node_id],
+                )
+                for e in topology.edges.values()
+            )
+        )
+        for permutation in itertools.permutations(nodes)
+        for relabel in [dict(zip(nodes, permutation, strict=True))]
+    )
+
+
+def _remove_isomorphic_topologies(
+    topologies: Iterable[Topology],
+) -> tuple[Topology, ...]:
+    unique_topologies: dict[tuple, Topology] = {}
+    for topology in sorted(topologies):
+        unique_topologies.setdefault(_create_isomorphism_invariant(topology), topology)
+    return tuple(unique_topologies.values())
+
+
 def create_isobar_topologies(
     number_of_final_states: int,
+    number_of_initial_states: int = 1,
 ) -> tuple[Topology, ...]:
-    """Builder function to create a set of unique isobar decay topologies.
+    r"""Builder function to create a set of unique isobar decay topologies.
+
+    All topologies are built from two-body interaction nodes only: :math:`1 \to 2`
+    decay nodes and — when there is more than one initial state — :math:`2 \to 1`
+    production nodes. For two initial states, this covers both :math:`s`-channel
+    topologies (the initial states annihilate into a single intermediate edge) and
+    :math:`t`-channel-like topologies (the initial states are connected through an
+    exchange edge). See `ComPWA/qrules#29 <https://github.com/ComPWA/qrules/issues/29>`_.
 
     Args:
         number_of_final_states: The number of `~Topology.outgoing_edge_ids`
             (`~.Transition.final_states`).
+        number_of_initial_states: The number of `~Topology.incoming_edge_ids`
+            (`~.Transition.initial_states`).
 
     Returns:
         A sorted `tuple` of non-isomorphic `Topology` instances, all with the same
-        number of final states.
+        number of initial and final states.
 
     Example:
         >>> topologies = create_isobar_topologies(number_of_final_states=4)
@@ -625,16 +762,30 @@ def create_isobar_topologies(
         2
         >>> list(topologies) == sorted(topologies)  # ordered
         True
+        >>> topologies = create_isobar_topologies(
+        ...     number_of_final_states=2,
+        ...     number_of_initial_states=2,
+        ... )
+        >>> len(topologies)
+        2
+        >>> {len(t.incoming_edge_ids) for t in topologies}
+        {2}
     """
+    if number_of_initial_states < 1:
+        msg = "At least one initial state required"
+        raise ValueError(msg)
     if number_of_final_states < 2:
         msg = "At least two final states required for an isobar decay"
         raise ValueError(msg)
-    builder = SimpleStateTransitionTopologyBuilder([InteractionNode(1, 2)])
+    interaction_nodes = [InteractionNode(1, 2)]
+    if number_of_initial_states > 1:
+        interaction_nodes.append(InteractionNode(2, 1))
+    builder = SimpleStateTransitionTopologyBuilder(interaction_nodes)
     topologies = builder.build(
-        number_of_initial_edges=1,
+        number_of_initial_edges=number_of_initial_states,
         number_of_final_edges=number_of_final_states,
     )
-    return tuple(sorted(topologies))
+    return _remove_isomorphic_topologies(topologies)
 
 
 def create_n_body_topology(
diff --git a/src/qrules/transition.py b/src/qrules/transition.py
index da0f8d94..035f0b7d 100644
--- a/src/qrules/transition.py
+++ b/src/qrules/transition.py
@@ -252,7 +252,7 @@ def __init__(  # ruff: ignore[too-many-positional-arguments]
         """`.Topology` instances over which the STM propagates quantum numbers."""
         # turn off mass conservation, in case more than one initial state
         # particle is present
-        if use_nbody_topology and len(initial_state) > 1:
+        if len(initial_state) > 1:
             mass_conservation_factor = None
 
         if reload_pdg or len(self.__particles) == 0:
diff --git a/src/qrules/workflow.py b/src/qrules/workflow.py
index f9387ffd..80bfa5d2 100644
--- a/src/qrules/workflow.py
+++ b/src/qrules/workflow.py
@@ -763,7 +763,7 @@ def create_qn_problem_sets(  # ruff: ignore[too-many-positional-arguments]
     )
     # turn off mass conservation, in case more than one initial state
     # particle is present
-    if use_nbody_topology and len(initial_state) > 1:
+    if len(initial_state) > 1:
         mass_conservation_factor = None
     if interaction_config is None:
         interaction_config = InteractionConfig(
@@ -1172,7 +1172,10 @@ def _create_topologies(
 ) -> tuple[tuple[Topology, ...], bool]:
     topology_building = topology_building.lower()
     if topology_building == "isobar":
-        return create_isobar_topologies(number_of_final_states), False
+        return (
+            create_isobar_topologies(number_of_final_states, number_of_initial_states),
+            False,
+        )
     if "n-body" in topology_building or "nbody" in topology_building:
         return (
             create_n_body_topology(number_of_initial_states, number_of_final_states),
diff --git a/tests/unit/test_topology.py b/tests/unit/test_topology.py
index a09e7eb9..e4d439c0 100644
--- a/tests/unit/test_topology.py
+++ b/tests/unit/test_topology.py
@@ -285,10 +285,12 @@ def it_unique_ordering(n_final_states):
         (2, 1, None),
         (3, 1, None),
         (4, 2, None),
-        (5, 5, None),
-        (6, 16, None),
-        (7, 61, None),
-        (8, 272, None),
+        # the number of topologies follows the Wedderburn-Etherington numbers
+        # (number of unordered rooted binary trees, https://oeis.org/A001190)
+        (5, 3, None),
+        (6, 6, None),
+        (7, 11, None),
+        (8, 23, None),
     ],
 )
 def test_create_isobar_topologies(
@@ -310,6 +312,41 @@ def test_create_isobar_topologies(
             assert len(topology.nodes) == n_expected_nodes
 
 
+@pytest.mark.parametrize(
+    ("n_final", "n_topologies"),
+    [
+        (2, 2),
+        (3, 5),
+        (4, 12),
+        (5, 30),
+    ],
+)
+def test_create_isobar_topologies_two_initial_states(n_final: int, n_topologies: int):
+    topologies = create_isobar_topologies(n_final, number_of_initial_states=2)
+    assert len(topologies) == n_topologies
+    for topology in topologies:
+        assert len(topology.incoming_edge_ids) == 2
+        assert len(topology.outgoing_edge_ids) == n_final
+        assert len(topology.intermediate_edge_ids) == n_final - 1
+        assert len(topology.nodes) == n_final
+
+
+def test_topology_builder_two_to_n():
+    """The builder terminates on 2-to-1 building blocks (ComPWA/qrules#29)."""
+    builder = SimpleStateTransitionTopologyBuilder([
+        InteractionNode(1, 2),
+        InteractionNode(2, 1),
+    ])
+    topologies = builder.build(
+        number_of_initial_edges=2,
+        number_of_final_edges=3,
+    )
+    assert len(topologies) == 6
+    for topology in topologies:
+        assert len(topology.incoming_edge_ids) == 2
+        assert len(topology.outgoing_edge_ids) == 3
+
+
 @pytest.mark.parametrize(
     ("n_initial", "n_final", "exception"),
     [
diff --git a/tests/unit/test_workflow.py b/tests/unit/test_workflow.py
index 1a993dc9..9a85373a 100644
--- a/tests/unit/test_workflow.py
+++ b/tests/unit/test_workflow.py
@@ -1,4 +1,5 @@
 import json
+from fractions import Fraction
 
 import pytest
 
@@ -312,6 +313,37 @@ def test_generate_qn_transitions():
     assert "0^{+}(0^{++})" in collapsed_mermaid
 
 
+def test_generate_qn_transitions_two_to_n():
+    """Production reactions with two initial states (ComPWA/qrules#29)."""
+    particle_db = load_pdg()
+    reaction = generate_qn_transitions(
+        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"],
+    )
+    assert len(reaction.transitions) > 0
+    assert {p.name for p in reaction.initial_state.values()} == {"gamma", "p"}
+    assert {p.name for p in reaction.final_state.values()} == {"p", "pi0"}
+    intermediate_signatures = {
+        (
+            state[EdgeQuantumNumbers.baryon_number],
+            state[EdgeQuantumNumbers.spin_magnitude],
+        )
+        for transition in reaction.transitions
+        for state in transition.intermediate_states.values()
+    }
+    assert (1, Fraction(3, 2)) in intermediate_signatures  # s-channel Delta(1232)
+    assert (0, Fraction(1)) in intermediate_signatures  # t-channel vector exchange
+    assert len(reaction.group_by_topology()) > 1
+
+
 def test_qn_reaction_info_requires_particle_states():
     particle_db = load_pdg()
     qn_problem_sets = create_qn_problem_sets(

From 1ce8132dd70d03cd7e986ac98dffa0693e7627dc Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 04:14:02 +0200
Subject: [PATCH 2/9] FIX: generalize conservation rules to production nodes

---
 src/qrules/conservation_rules.py              | 126 +++++++++++-------
 tests/channels/test_gammap_to_p_pi0.py        |  32 +++++
 .../test_helicity_conservation.py             |  41 ++++--
 .../test_parity_conservation.py               |  24 ++++
 4 files changed, 170 insertions(+), 53 deletions(-)
 create mode 100644 tests/channels/test_gammap_to_p_pi0.py

diff --git a/src/qrules/conservation_rules.py b/src/qrules/conservation_rules.py
index c77c1ebd..f528f256 100644
--- a/src/qrules/conservation_rules.py
+++ b/src/qrules/conservation_rules.py
@@ -229,10 +229,15 @@ def parity_conservation(
     outgoing_edge_qns: list[EdgeQN.parity],
     l_magnitude: NodeQN.l_magnitude,
 ) -> bool:
-    r"""Implement :math:`P_{in} = P_{out} \cdot (-1)^L`."""
+    r"""Implement :math:`P_{in} = P_{out} \cdot (-1)^L`.
+
+    The :math:`(-1)^L` factor originates from the two-body side of the interaction
+    node, so the same equation applies to both decay (:math:`1 \to 2\,3`) and
+    production (:math:`1\,2 \to 3`) nodes.
+    """
     if any(p is None for p in [*ingoing_edge_qns, *outgoing_edge_qns]):
         return False
-    if len(ingoing_edge_qns) == 1 and len(outgoing_edge_qns) == 2:
+    if (len(ingoing_edge_qns), len(outgoing_edge_qns)) in {(1, 2), (2, 1)}:
         parity_in = reduce(lambda x, y: x * y.value, ingoing_edge_qns, 1)
         parity_out = reduce(lambda x, y: x * y.value, outgoing_edge_qns, 1)
         return parity_in == (parity_out * (-1) ** l_magnitude)
@@ -262,24 +267,34 @@ def parity_conservation_helicity(
 
     .. note:: Only the special case :math:`\lambda_1=\lambda_2=0` may return `False`
         independent on the parity prefactor.
+
+    For a production node :math:`1\,2 \to 3`, the same check applies with the roles of
+    the single state and the two-body pair mirrored.
     """
-    if len(ingoing_edge_qns) == 1 and len(outgoing_edge_qns) == 2:
-        out_spins = [x.spin_magnitude for x in outgoing_edge_qns]
-        parity_product = reduce(
-            lambda x, y: x * y.parity.value if y.parity else x,
-            ingoing_edge_qns + outgoing_edge_qns,
-            1,
-        )
+    arity = (len(ingoing_edge_qns), len(outgoing_edge_qns))
+    if arity == (1, 2):
+        single_state = ingoing_edge_qns[0]
+        two_body_states = outgoing_edge_qns
+    elif arity == (2, 1):
+        single_state = outgoing_edge_qns[0]
+        two_body_states = ingoing_edge_qns
+    else:
+        return True
+    two_body_spins = [x.spin_magnitude for x in two_body_states]
+    parity_product = reduce(
+        lambda x, y: x * y.parity.value if y.parity else x,
+        ingoing_edge_qns + outgoing_edge_qns,
+        1,
+    )
 
-        prefactor = parity_product * (-1.0) ** (
-            sum(out_spins) - ingoing_edge_qns[0].spin_magnitude
-        )
+    prefactor = parity_product * (-1.0) ** (
+        sum(two_body_spins) - single_state.spin_magnitude
+    )
 
-        if all(x.spin_projection == 0.0 for x in outgoing_edge_qns) and prefactor == -1:  # ruff: ignore[float-equality-comparison]
-            return False
+    if all(x.spin_projection == 0.0 for x in two_body_states) and prefactor == -1:  # ruff: ignore[float-equality-comparison]
+        return False
 
-        return prefactor == parity_prefactor
-    return True
+    return prefactor == parity_prefactor
 
 
 @frozen
@@ -523,6 +538,7 @@ class SpinMagnitudeFacts(TypedDict):
 class HelicityFacts(TypedDict):
     """Facts required by `helicity_conservation`; a subset of `.EdgeFacts`."""
 
+    spin_magnitude: Fraction
     spin_projection: Fraction
 
 
@@ -783,7 +799,7 @@ def clebsch_gordan_helicity_to_canonical(
     outgoing_spins: list[SpinEdgeInput],
     interaction_qns: SpinNodeInput,
 ) -> bool:
-    """Implement Clebsch-Gordan checks.
+    r"""Implement Clebsch-Gordan checks.
 
     For :math:`S_1, S_2` to :math:`S` and the :math:`L,S` to :math:`J` coupling based on
     the conversion of helicity to canonical amplitude sums.
@@ -791,52 +807,72 @@ def clebsch_gordan_helicity_to_canonical(
     .. note:: This rule does not check that the spin magnitudes couple correctly to
         :math:`L` and :math:`S`, as this is already performed by
         `.spin_magnitude_conservation`.
+
+    For a production node :math:`1\,2 \to 3`, the same couplings are checked with the
+    roles of the single state and the two-body pair mirrored.
     """
     if len(ingoing_spins) == 1 and len(outgoing_spins) == 2:
-        out_spin1 = _Spin(
-            outgoing_spins[0].spin_magnitude,
-            outgoing_spins[0].spin_projection,
+        return __check_two_body_helicity_to_canonical_coupling(
+            ingoing_spins[0].spin_magnitude, outgoing_spins, interaction_qns
         )
-        out_spin2 = _Spin(
-            outgoing_spins[1].spin_magnitude,
-            -outgoing_spins[1].spin_projection,
+    if len(ingoing_spins) == 2 and len(outgoing_spins) == 1:
+        return __check_two_body_helicity_to_canonical_coupling(
+            outgoing_spins[0].spin_magnitude, ingoing_spins, interaction_qns
         )
+    return False
 
-        helicity_diff = out_spin1.projection + out_spin2.projection
-        if helicity_diff != interaction_qns.s_projection:
-            return False
 
-        ang_mom = _Spin(interaction_qns.l_magnitude, interaction_qns.l_projection)
-        coupled_spin = _Spin(interaction_qns.s_magnitude, interaction_qns.s_projection)
-        parent_spin = ingoing_spins[0].spin_magnitude
+def __check_two_body_helicity_to_canonical_coupling(
+    single_spin_magnitude: Fraction,
+    two_body_spins: list[SpinEdgeInput],
+    interaction_qns: SpinNodeInput,
+) -> bool:
+    spin1 = _Spin(
+        two_body_spins[0].spin_magnitude,
+        two_body_spins[0].spin_projection,
+    )
+    spin2 = _Spin(
+        two_body_spins[1].spin_magnitude,
+        -two_body_spins[1].spin_projection,
+    )
 
-        coupled_spin = _Spin(coupled_spin.magnitude, helicity_diff)
-        if not _check_spin_valid(coupled_spin.magnitude, coupled_spin.projection):
-            return False
-        in_spin = _Spin(parent_spin, helicity_diff)
-        if not _check_spin_valid(in_spin.magnitude, in_spin.projection):
-            return False
+    helicity_diff = spin1.projection + spin2.projection
+    if helicity_diff != interaction_qns.s_projection:
+        return False
 
-        if _is_clebsch_gordan_coefficient_zero(out_spin1, out_spin2, coupled_spin):
-            return False
+    ang_mom = _Spin(interaction_qns.l_magnitude, interaction_qns.l_projection)
 
-        return not _is_clebsch_gordan_coefficient_zero(ang_mom, coupled_spin, in_spin)
-    return False
+    coupled_spin = _Spin(interaction_qns.s_magnitude, helicity_diff)
+    if not _check_spin_valid(coupled_spin.magnitude, coupled_spin.projection):
+        return False
+    single_spin = _Spin(single_spin_magnitude, helicity_diff)
+    if not _check_spin_valid(single_spin.magnitude, single_spin.projection):
+        return False
+
+    if _is_clebsch_gordan_coefficient_zero(spin1, spin2, coupled_spin):
+        return False
+
+    return not _is_clebsch_gordan_coefficient_zero(ang_mom, coupled_spin, single_spin)
 
 
 def helicity_conservation(
-    ingoing_spin_mags: list[SpinMagnitudeFacts],
+    ingoing_helicities: list[HelicityFacts],
     outgoing_helicities: list[HelicityFacts],
 ) -> bool:
     r"""Implementation of helicity conservation.
 
-    Check for :math:`|\lambda_2-\lambda_3| \leq S_1`.
+    For a decay node :math:`1 \to 2\,3`, check :math:`|\lambda_2-\lambda_3| \leq S_1`.
+    For a production node :math:`1\,2 \to 3`, check the mirrored condition
+    :math:`|\lambda_1-\lambda_2| \leq S_3`.
     """
-    if len(ingoing_spin_mags) == 1 and len(outgoing_helicities) == 2:
-        mother_spin = ingoing_spin_mags[0]["spin_magnitude"]
+    if len(ingoing_helicities) == 1 and len(outgoing_helicities) == 2:
+        parent_spin = ingoing_helicities[0]["spin_magnitude"]
         helicities = [x["spin_projection"] for x in outgoing_helicities]
-        if mother_spin >= abs(helicities[0] - helicities[1]):
-            return True
+        return abs(helicities[0] - helicities[1]) <= parent_spin
+    if len(ingoing_helicities) == 2 and len(outgoing_helicities) == 1:
+        daughter_spin = outgoing_helicities[0]["spin_magnitude"]
+        helicities = [x["spin_projection"] for x in ingoing_helicities]
+        return abs(helicities[0] - helicities[1]) <= daughter_spin
     return False
 
 
diff --git a/tests/channels/test_gammap_to_p_pi0.py b/tests/channels/test_gammap_to_p_pi0.py
new file mode 100644
index 00000000..a9c4c7cc
--- /dev/null
+++ b/tests/channels/test_gammap_to_p_pi0.py
@@ -0,0 +1,32 @@
+"""Test 2-to-n production reactions, https://github.com/ComPWA/qrules/issues/29."""
+
+import pytest
+
+import qrules
+from qrules.transition import SpinFormalism
+
+
+@pytest.mark.parametrize("formalism", ["helicity", "canonical-helicity"])
+def test_pi0_photoproduction(formalism: SpinFormalism, particle_database):
+    reaction = qrules.generate_transitions(
+        initial_state=["gamma", "p"],
+        final_state=["p", "pi0"],
+        allowed_intermediate_particles=[
+            "Delta(1232)",
+            "N(1440)",
+            "rho(770)",
+            "omega(782)",
+        ],
+        allowed_interaction_types=["strong", "em"],
+        formalism=formalism,
+        particle_db=particle_database,
+    )
+    assert reaction.get_intermediate_particles().names == [
+        "Delta(1232)~-",  # u-channel baryon exchange
+        "Delta(1232)+",  # s-channel resonance
+        "N(1440)~-",
+        "N(1440)+",
+        "omega(782)",  # t-channel meson exchange
+        "rho(770)0",
+    ]
+    assert {p.name for p in reaction.initial_state.values()} == {"gamma", "p"}
diff --git a/tests/unit/conservation_rules/test_helicity_conservation.py b/tests/unit/conservation_rules/test_helicity_conservation.py
index 4ba63398..cdaa1a61 100644
--- a/tests/unit/conservation_rules/test_helicity_conservation.py
+++ b/tests/unit/conservation_rules/test_helicity_conservation.py
@@ -3,22 +3,47 @@
 
 import pytest
 
-from qrules.conservation_rules import (
-    HelicityFacts,
-    SpinMagnitudeFacts,
-    helicity_conservation,
+from qrules.conservation_rules import HelicityFacts, helicity_conservation
+
+
+def _helicity_facts(magnitude: float, projection: float) -> HelicityFacts:
+    return HelicityFacts(
+        spin_magnitude=Fraction(magnitude),
+        spin_projection=Fraction(projection),
+    )
+
+
+@pytest.mark.parametrize(
+    ("in_edge_qns", "out_edge_qns", "expected"),
+    [
+        (
+            [_helicity_facts(s_magnitude, min(s_magnitude, 1))],
+            [
+                _helicity_facts(abs(lambda1), lambda1),
+                _helicity_facts(abs(lambda2), lambda2),
+            ],
+            abs(lambda1 - lambda2) <= s_magnitude,
+        )
+        for s_magnitude, lambda1, lambda2 in product(
+            [0, 0.5, 1, 1.5, 2],
+            [-2, -1.5, -1.0, -0.5, 0, 0.5, 1, 1.5, 2],
+            [-1, 0, 1],
+        )
+    ],
 )
+def test_helicity_conservation_decay(in_edge_qns, out_edge_qns, expected):
+    assert helicity_conservation(in_edge_qns, out_edge_qns) is expected
 
 
 @pytest.mark.parametrize(
     ("in_edge_qns", "out_edge_qns", "expected"),
     [
         (
-            [SpinMagnitudeFacts(spin_magnitude=Fraction(s_magnitude))],
             [
-                HelicityFacts(spin_projection=Fraction(lambda1)),
-                HelicityFacts(spin_projection=Fraction(lambda2)),
+                _helicity_facts(abs(lambda1), lambda1),
+                _helicity_facts(abs(lambda2), lambda2),
             ],
+            [_helicity_facts(s_magnitude, min(s_magnitude, 1))],
             abs(lambda1 - lambda2) <= s_magnitude,
         )
         for s_magnitude, lambda1, lambda2 in product(
@@ -28,5 +53,5 @@
         )
     ],
 )
-def test_helicity_conservation(in_edge_qns, out_edge_qns, expected):
+def test_helicity_conservation_production(in_edge_qns, out_edge_qns, expected):
     assert helicity_conservation(in_edge_qns, out_edge_qns) is expected
diff --git a/tests/unit/conservation_rules/test_parity_conservation.py b/tests/unit/conservation_rules/test_parity_conservation.py
index 566b05da..33d1f7e1 100644
--- a/tests/unit/conservation_rules/test_parity_conservation.py
+++ b/tests/unit/conservation_rules/test_parity_conservation.py
@@ -34,6 +34,30 @@ def describe_parity_conservation():
     def it_conserves_parity(in_parities, out_parities, l_magnitude, expected):
         assert parity_conservation(in_parities, out_parities, l_magnitude) is expected
 
+    @pytest.mark.parametrize(
+        ("in_parities", "out_parities", "l_magnitude", "expected"),
+        [
+            (
+                [
+                    Parity(parity_in1),
+                    Parity(1),
+                ],
+                [
+                    Parity(parity_out),
+                ],
+                NodeQuantumNumbers.l_magnitude(Fraction(l_magnitude)),
+                parity_in1 == parity_out * (-1) ** (l_magnitude),
+            )
+            for parity_in1, parity_out, l_magnitude in product(
+                [-1, 1], [-1, 1], range(5)
+            )
+        ],
+    )
+    def it_conserves_parity_in_production(
+        in_parities, out_parities, l_magnitude, expected
+    ):
+        assert parity_conservation(in_parities, out_parities, l_magnitude) is expected
+
     @pytest.mark.parametrize(
         ("in_parities", "out_parities", "l_magnitude", "expected"),
         [

From b516d954b18210a07629ec507bc63ed1311cbed2 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 04:14:02 +0200
Subject: [PATCH 3/9] DOC: demonstrate production reactions in example notebook

---
 .cspell.json                    |   1 +
 docs/usage/qn-transitions.ipynb | 100 ++++++++++++++++++++++++++++++++
 2 files changed, 101 insertions(+)

diff --git a/.cspell.json b/.cspell.json
index d81641b9..b3689f2b 100644
--- a/.cspell.json
+++ b/.cspell.json
@@ -155,6 +155,7 @@
     "bottomness",
     "breit",
     "charmness",
+    "charmonium",
     "clebsch",
     "combi",
     "conda",
diff --git a/docs/usage/qn-transitions.ipynb b/docs/usage/qn-transitions.ipynb
index 9c2ffb31..9c85c2ac 100644
--- a/docs/usage/qn-transitions.ipynb
+++ b/docs/usage/qn-transitions.ipynb
@@ -34,6 +34,7 @@
     "\n",
     "import qrules.io\n",
     "from qrules.quantum_numbers import EdgeQuantumNumbers\n",
+    "from qrules.topology import create_isobar_topologies\n",
     "from qrules.workflow import (\n",
     "    create_qn_problem_sets,\n",
     "    find_qn_transitions,\n",
@@ -376,6 +377,105 @@
     "f\"{len(reaction_5body.transitions):,} transitions over {len(reaction_5body.group_by_topology())} topologies\""
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Production reactions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Reactions with **two initial states** are supported as well (see [ComPWA/qrules#29](https://github.com/ComPWA/qrules/issues/29)). For more than one initial state, {func}`.create_isobar_topologies` builds the topologies from $2 \\to 1$ production nodes in addition to the usual $1 \\to 2$ decay nodes, which produces both $s$-channel topologies (the initial states annihilate into a single edge) and $t$-channel-like topologies (the initial states are connected through an exchange edge). For a $2 \\to 2$ reaction:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "topologies = create_isobar_topologies(\n",
+    "    number_of_final_states=2,\n",
+    "    number_of_initial_states=2,\n",
+    ")\n",
+    "Markdown(qrules.io.asmermaid(topologies, markdown=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This makes it possible to describe production processes like $\\pi^0$ photoproduction, $\\gamma p \\to p\\pi^0$ (GlueX). The allowed transitions contain the $s$-channel $\\Delta$ and $N^*$ resonances as well as $t$-channel $\\rho/\\omega$ exchange and $u$-channel-like baryon exchange:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "reaction_2to2 = generate_qn_transitions(\n",
+    "    initial_state=[\"gamma\", \"p\"],\n",
+    "    final_state=[\"p\", \"pi0\"],\n",
+    "    particle_db=PDG,\n",
+    "    allowed_intermediate_particles=[\"Delta(1232)\", \"N(1440)\", \"rho(770)\", \"omega(782)\"],\n",
+    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
+    ")\n",
+    "Markdown(qrules.io.asmermaid(reaction_2to2, collapse_graphs=True, markdown=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When there is more than one initial state, the {class}`.MassConservation` rule is switched off, so the collision energy does not enter the problem. Model a specific energy by restricting :code:`allowed_intermediate_particles` to the resonances that are accessible at that energy. For instance, $e^+e^- \\to p\\bar p\\eta$ at BESIII charmonium energies, where the annihilation proceeds through a $J^{PC}=1^{--}$ state and the $p\\eta$ subsystem can form $N^*$ resonances:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "reaction_ee = generate_qn_transitions(\n",
+    "    initial_state=[\"e+\", \"e-\"],\n",
+    "    final_state=[\"p\", \"p~\", \"eta\"],\n",
+    "    particle_db=PDG,\n",
+    "    allowed_intermediate_particles=[\"J/psi(1S)\", \"psi(2S)\", \"N(1535)\"],\n",
+    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
+    ")\n",
+    "len(reaction_ee.transitions)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Only $s$-channel topologies survive here: a $t$-channel exchange would require an intermediate state that carries lepton number. The spins, parities, and baryon numbers of the intermediate edges confirm the photon-like annihilation state and the nucleon resonances:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "{\n",
+    "    (\n",
+    "        state[EdgeQuantumNumbers.spin_magnitude],\n",
+    "        state[EdgeQuantumNumbers.parity],\n",
+    "        state[EdgeQuantumNumbers.baryon_number],\n",
+    "    )\n",
+    "    for transition in reaction_ee.transitions\n",
+    "    for state in transition.intermediate_states.values()\n",
+    "}"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

From cdbd308aa8de7194ad48a0a16395701067728477 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 09:02:11 +0200
Subject: [PATCH 4/9] FIX: assign initial facts to sorted external edge IDs

---
 src/qrules/combinatorics.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/qrules/combinatorics.py b/src/qrules/combinatorics.py
index 32194c77..24f02680 100644
--- a/src/qrules/combinatorics.py
+++ b/src/qrules/combinatorics.py
@@ -236,7 +236,7 @@ def create_initial_facts(
     :math:`J^{P(C)}` level (see `.strip_spin_projections`).
     """
     states = __create_states_with_spin_projections(
-        list(topology.incoming_edge_ids) + list(topology.outgoing_edge_ids),
+        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,

From 50c746035b3c285d540a96e978cd4623afadeb88 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 09:02:29 +0200
Subject: [PATCH 5/9] BEHAVIOR: validate G-parity against isospin and C-parity

---
 .cspell.json                                  |  2 +
 src/qrules/conservation_rules.py              | 56 +++++++++++++++++++
 src/qrules/settings.py                        |  3 +
 .../conservation_rules/test_duck_typing.py    |  1 +
 tests/unit/io/test_labels.py                  |  1 +
 5 files changed, 63 insertions(+)

diff --git a/.cspell.json b/.cspell.json
index b3689f2b..cd1d3817 100644
--- a/.cspell.json
+++ b/.cspell.json
@@ -188,7 +188,9 @@
     "lambdify",
     "lineshape",
     "lineshapes",
+    "mandelstam",
     "mathbb",
+    "multiplet",
     "nishijima",
     "pydot",
     "PyPA",
diff --git a/src/qrules/conservation_rules.py b/src/qrules/conservation_rules.py
index f528f256..8d11cdf3 100644
--- a/src/qrules/conservation_rules.py
+++ b/src/qrules/conservation_rules.py
@@ -366,6 +366,62 @@ class GParityNodeInput:
     s_magnitude: NodeSMagnitude = field(converter=to_fraction)
 
 
+@frozen
+class GParityValidityInput:
+    isospin_magnitude: EdgeIsospinMagnitude = field(converter=to_fraction)
+    charge: EdgeQN.charge = field(converter=EdgeQN.charge)
+    c_parity: EdgeCParity | None = field(converter=optional(to_parity), default=None)
+    g_parity: EdgeGParity | None = field(converter=optional(to_parity), default=None)
+
+
+def g_parity_validity(state: GParityValidityInput) -> bool:
+    r"""Check that the :math:`G`-parity is consistent with isospin and :math:`C`-parity.
+
+    :math:`G`-parity is only defined for states with integer isospin and satisfies
+    :math:`G = C \cdot (-1)^I`, with :math:`C` the :math:`C`-parity of the neutral
+    member of the isospin multiplet. A neutral state that carries a :math:`G`-parity
+    is such a neutral member, so its :math:`C`-parity must be defined and satisfy this
+    relation. For charged states, the relation cannot be checked, since their
+    :math:`C`-parity is undefined.
+
+    >>> from fractions import Fraction
+    >>> from qrules.quantum_numbers import Parity
+    >>> rho_meson = GParityValidityInput(
+    ...     isospin_magnitude=Fraction(1),
+    ...     charge=0,
+    ...     c_parity=Parity(-1),
+    ...     g_parity=Parity(+1),
+    ... )
+    >>> g_parity_validity(rho_meson)
+    True
+    >>> undefined_c_parity = GParityValidityInput(
+    ...     isospin_magnitude=Fraction(1),
+    ...     charge=0,
+    ...     g_parity=Parity(-1),
+    ... )
+    >>> g_parity_validity(undefined_c_parity)
+    False
+    >>> charged_rho_meson = GParityValidityInput(
+    ...     isospin_magnitude=Fraction(1),
+    ...     charge=+1,
+    ...     g_parity=Parity(+1),
+    ... )
+    >>> g_parity_validity(charged_rho_meson)
+    True
+    """
+    if state.g_parity is None:
+        return True
+    if state.isospin_magnitude.denominator != 1:
+        return False
+    if state.charge != 0:
+        return True
+    if state.c_parity is None:
+        return False
+    return state.g_parity.value == state.c_parity.value * (-1) ** int(
+        state.isospin_magnitude
+    )
+
+
 def g_parity_conservation(  # ruff: ignore[complex-structure]
     ingoing_edge_qns: list[GParityEdgeInput],
     outgoing_edge_qns: list[GParityEdgeInput],
diff --git a/src/qrules/settings.py b/src/qrules/settings.py
index 61c18124..be252a30 100644
--- a/src/qrules/settings.py
+++ b/src/qrules/settings.py
@@ -30,6 +30,7 @@
     c_parity_conservation,
     clebsch_gordan_helicity_to_canonical,
     g_parity_conservation,
+    g_parity_validity,
     gellmann_nishijima,
     helicity_conservation,
     identical_particle_symmetrization,
@@ -85,6 +86,7 @@
 
 EDGE_RULE_PRIORITIES: dict[RuleKey, int] = {
     gellmann_nishijima: 50,
+    g_parity_validity: 60,
     isospin_validity: 61,
     spin_validity: 62,
 }
@@ -149,6 +151,7 @@ def create_interaction_settings(  # ruff: ignore[too-many-positional-arguments]
         conservation_rules=_with_priorities(
             {
                 isospin_validity,
+                g_parity_validity,
                 gellmann_nishijima,
                 spin_validity,
             },
diff --git a/tests/unit/conservation_rules/test_duck_typing.py b/tests/unit/conservation_rules/test_duck_typing.py
index 79ca597a..56856584 100644
--- a/tests/unit/conservation_rules/test_duck_typing.py
+++ b/tests/unit/conservation_rules/test_duck_typing.py
@@ -33,6 +33,7 @@ def test_protocol_compliance():
     assert edge_input_classes == {
         conservation_rules.CParityEdgeInput,
         conservation_rules.GParityEdgeInput,
+        conservation_rules.GParityValidityInput,
         conservation_rules.HelicityParityEdgeInput,
         conservation_rules.IdenticalParticleSymmetryOutEdgeInput,
         conservation_rules.IsoSpinEdgeInput,
diff --git a/tests/unit/io/test_labels.py b/tests/unit/io/test_labels.py
index 09bcab8e..67dcf4e2 100644
--- a/tests/unit/io/test_labels.py
+++ b/tests/unit/io/test_labels.py
@@ -200,6 +200,7 @@ def it_dict(
             RULES
             spin_validity - 62
             isospin_validity - 61
+            g_parity_validity - 60
             gellmann_nishijima - 50
             DOMAINS
             baryon_number ∊ [-1, +1]

From 587a834414d828bee5601855335f3a119ad98be4 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 09:02:44 +0200
Subject: [PATCH 6/9] ENH: render minimal quantum-number signatures when
 collapsing

---
 src/qrules/io/_labels.py | 62 ++++++++++++++++++++++++++++++++++++----
 1 file changed, 56 insertions(+), 6 deletions(-)

diff --git a/src/qrules/io/_labels.py b/src/qrules/io/_labels.py
index afb2c3aa..2e306d66 100644
--- a/src/qrules/io/_labels.py
+++ b/src/qrules/io/_labels.py
@@ -19,7 +19,13 @@
     NodeSettings,
     QNProblemSet,
 )
-from qrules.topology import FrozenTransition, MutableTransition, Topology, Transition
+from qrules.topology import (
+    FrozenDict,
+    FrozenTransition,
+    MutableTransition,
+    Topology,
+    Transition,
+)
 from qrules.transition import ProblemSet, State
 
 if TYPE_CHECKING:
@@ -562,6 +568,9 @@ class QuantumNumberSignature:
 
     `collapse_graphs` converts states that are quantum-number property maps to this
     compact form, so that a collapsed edge lists signatures instead of complete maps.
+    The :math:`C`-parity is only taken over for self-conjugate states (zero charge,
+    baryon number, and strangeness) and the :math:`G`-parity only for nonstrange
+    non-baryonic states, since the quantum numbers are undefined otherwise.
 
     >>> from qrules.quantum_numbers import EdgeQuantumNumbers as EQN
     >>> signature = QuantumNumberSignature.from_property_map({
@@ -577,6 +586,16 @@ class QuantumNumberSignature:
     '1^{+}(1^{--})'
     >>> as_string(QuantumNumberSignature.from_property_map({EQN.spin_magnitude: 0.5}))
     '1/2'
+    >>> baryon = QuantumNumberSignature.from_property_map({
+    ...     EQN.spin_magnitude: 1.5,
+    ...     EQN.parity: +1,
+    ...     EQN.c_parity: -1,
+    ...     EQN.isospin_magnitude: 1.5,
+    ...     EQN.g_parity: +1,
+    ...     EQN.baryon_number: 1,
+    ... })
+    >>> as_string(baryon)
+    '3/2(3/2⁺)'
     """
 
     spin_magnitude: Fraction | None = None
@@ -587,16 +606,25 @@ class QuantumNumberSignature:
 
     @classmethod
     def from_property_map(cls, qn_map: Mapping[Any, Any]) -> QuantumNumberSignature:
+        baryon_number = qn_map.get(EdgeQuantumNumbers.baryon_number) or 0
+        charge = qn_map.get(EdgeQuantumNumbers.charge) or 0
+        strangeness = qn_map.get(EdgeQuantumNumbers.strangeness) or 0
+        c_parity = None
+        if baryon_number == 0 and charge == 0 and strangeness == 0:
+            c_parity = _to_optional_int(qn_map.get(EdgeQuantumNumbers.c_parity))
+        g_parity = None
+        if baryon_number == 0 and strangeness == 0:
+            g_parity = _to_optional_int(qn_map.get(EdgeQuantumNumbers.g_parity))
         return cls(
             spin_magnitude=_to_optional_fraction(
                 qn_map.get(EdgeQuantumNumbers.spin_magnitude)
             ),
             parity=_to_optional_int(qn_map.get(EdgeQuantumNumbers.parity)),
-            c_parity=_to_optional_int(qn_map.get(EdgeQuantumNumbers.c_parity)),
+            c_parity=c_parity,
             isospin_magnitude=_to_optional_fraction(
                 qn_map.get(EdgeQuantumNumbers.isospin_magnitude)
             ),
-            g_parity=_to_optional_int(qn_map.get(EdgeQuantumNumbers.g_parity)),
+            g_parity=g_parity,
         )
 
 
@@ -757,8 +785,8 @@ def collapse_graphs(
             FrozenTransition(
                 topology,
                 states={
-                    i: tuple(sorted(particles, key=_sorting_key))
-                    for i, particles in group.states.items()
+                    i: tuple(sorted(_summarize_property_maps(states), key=_sorting_key))
+                    for i, states in group.states.items()
                 },
                 interactions=group.interactions,
             )
@@ -770,10 +798,32 @@ def _strip_properties(state: Any) -> Any:
     if isinstance(state, State):
         return state.particle
     if isinstance(state, abc.Mapping):
-        return QuantumNumberSignature.from_property_map(state)
+        return FrozenDict(state)
     return state
 
 
+def _summarize_property_maps(states: Iterable[Any]) -> set[Any]:
+    """Replace collapsed quantum-number property maps by unique signatures.
+
+    Property maps in which some quantum numbers are merely unassigned (`None`) are
+    dropped when a more determined map with the same assigned values is present, so
+    that the collapsed edge label stays minimal.
+    """
+    property_maps = [state for state in states if isinstance(state, abc.Mapping)]
+    summarized: set[Any] = {
+        state for state in states if not isinstance(state, abc.Mapping)
+    }
+    assigned_values = [
+        {key: value for key, value in qn_map.items() if value is not None}
+        for qn_map in property_maps
+    ]
+    for qn_map, assigned in zip(property_maps, assigned_values, strict=True):
+        is_subsumed = any(assigned.items() < other.items() for other in assigned_values)
+        if not is_subsumed:
+            summarized.add(QuantumNumberSignature.from_property_map(qn_map))
+    return summarized
+
+
 def _sorting_key(obj: Any) -> Any:
     if isinstance(obj, State):
         return obj.particle.name

From 0485aed33e284b674a879053037d14edbfcb2014 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 09:02:59 +0200
Subject: [PATCH 7/9] FEAT: classify transitions by Mandelstam channel

---
 src/qrules/topology.py                 | 97 +++++++++++++++++++++++++-
 src/qrules/transition.py               | 17 ++++-
 src/qrules/workflow.py                 | 60 +++++++++++++++-
 tests/channels/test_gammap_to_p_pi0.py | 29 ++++++++
 tests/unit/test_workflow.py            | 40 +++++++++++
 5 files changed, 239 insertions(+), 4 deletions(-)

diff --git a/src/qrules/topology.py b/src/qrules/topology.py
index bdf6db31..ae195ec6 100644
--- a/src/qrules/topology.py
+++ b/src/qrules/topology.py
@@ -22,7 +22,7 @@
 from abc import ABC, abstractmethod
 from collections import abc, defaultdict
 from functools import total_ordering
-from typing import TYPE_CHECKING, Any, Generic, TypeVar, overload
+from typing import TYPE_CHECKING, Any, Generic, Literal, TypeVar, overload
 
 import attrs
 from attrs import define, field, frozen
@@ -788,6 +788,101 @@ def create_isobar_topologies(
     return _remove_isomorphic_topologies(topologies)
 
 
+def determine_mandelstam_channel(
+    topology: Topology, intermediate_edge_id: int
+) -> Literal["s", "t", "u"]:
+    r"""Determine which Mandelstam channel an intermediate edge represents.
+
+    Cutting the intermediate edge splits a tree topology into two sides. If one side
+    contains all initial states, the edge carries the invariant mass of the final
+    states on the other side, that is, an :math:`s`-type channel. Otherwise, the edge
+    is an exchange between the initial states. Following the Mandelstam convention
+    :math:`t = (p_1 - p_3)^2` and :math:`u = (p_1 - p_4)^2` for :math:`1\,2 \to 3\,4`,
+    the exchange is labeled :math:`t` if the first initial state and the first final
+    state lie on the same side of the cut and :math:`u` otherwise. Reorder the final
+    states to switch between the two.
+
+    >>> topologies = create_isobar_topologies(
+    ...     number_of_final_states=2,
+    ...     number_of_initial_states=2,
+    ... )
+    >>> [determine_mandelstam_channel(t, intermediate_edge_id=2) for t in topologies]
+    ['s', 't']
+    >>> u_channel_topology = topologies[1].relabel_edges({0: 1, 1: 0})
+    >>> determine_mandelstam_channel(u_channel_topology, intermediate_edge_id=2)
+    'u'
+    >>> determine_mandelstam_channel(topologies[0], intermediate_edge_id=0)
+    Traceback (most recent call last):
+        ...
+    ValueError: Edge 0 is not an intermediate edge
+    """
+    if intermediate_edge_id not in topology.intermediate_edge_ids:
+        msg = f"Edge {intermediate_edge_id} is not an intermediate edge"
+        raise ValueError(msg)
+    cut_edge = topology.edges[intermediate_edge_id]
+    originating_side_node_ids = __collect_node_ids_on_side(
+        topology,
+        start_node_id=min(cut_edge.get_connected_nodes()),
+        cut_edge_id=intermediate_edge_id,
+    )
+
+    def is_on_originating_side(external_edge_id: int) -> bool:
+        (node_id,) = topology.edges[external_edge_id].get_connected_nodes()
+        return node_id in originating_side_node_ids
+
+    initial_state_sides = {
+        is_on_originating_side(i) for i in topology.incoming_edge_ids
+    }
+    if len(initial_state_sides) == 1:
+        return "s"
+    beam_side = is_on_originating_side(min(topology.incoming_edge_ids))
+    first_final_state_side = is_on_originating_side(min(topology.outgoing_edge_ids))
+    if beam_side == first_final_state_side:
+        return "t"
+    return "u"
+
+
+def __collect_node_ids_on_side(
+    topology: Topology, start_node_id: int, cut_edge_id: int
+) -> set[int]:
+    connected_node_ids = {start_node_id}
+    queue = [start_node_id]
+    while queue:
+        node_id = queue.pop()
+        for edge_id, edge in topology.edges.items():
+            if edge_id == cut_edge_id or node_id not in edge.get_connected_nodes():
+                continue
+            new_node_ids = edge.get_connected_nodes() - connected_node_ids
+            connected_node_ids |= new_node_ids
+            queue += new_node_ids
+    return connected_node_ids
+
+
+def determine_reaction_channel(topology: Topology) -> str:
+    """Summarize the Mandelstam channels of a topology as a single label.
+
+    The label is :code:`"s"` if none of the intermediate edges is an exchange between
+    the initial states, and otherwise the alphabetically sorted concatenation of the
+    `.determine_mandelstam_channel` letters of the exchange edges — :code:`"t"` or
+    :code:`"u"` for a single exchange, :code:`"tt"` etc. for double exchange.
+
+    >>> topologies = create_isobar_topologies(
+    ...     number_of_final_states=2,
+    ...     number_of_initial_states=2,
+    ... )
+    >>> [determine_reaction_channel(t) for t in topologies]
+    ['s', 't']
+    """
+    channels = [
+        determine_mandelstam_channel(topology, edge_id)
+        for edge_id in topology.intermediate_edge_ids
+    ]
+    exchange_channels = [channel for channel in channels if channel != "s"]
+    if not exchange_channels:
+        return "s"
+    return "".join(sorted(exchange_channels))
+
+
 def create_n_body_topology(
     number_of_initial_states: int, number_of_final_states: int
 ) -> Topology:
diff --git a/src/qrules/transition.py b/src/qrules/transition.py
index 035f0b7d..6a592c12 100644
--- a/src/qrules/transition.py
+++ b/src/qrules/transition.py
@@ -41,7 +41,13 @@
     create_edge_properties,
     create_node_properties,
 )
-from qrules.topology import FrozenDict, FrozenTransition, MutableTransition, Topology
+from qrules.topology import (
+    FrozenDict,
+    FrozenTransition,
+    MutableTransition,
+    Topology,
+    determine_reaction_channel,
+)
 
 if TYPE_CHECKING:
     from collections.abc import Iterable, Sequence
@@ -430,3 +436,12 @@ def group_by_topology(self) -> dict[Topology, list[StateTransition]]:
         for transition in self.transitions:
             groupings[transition.topology].append(transition)
         return dict(groupings)
+
+    def group_by_channel(self) -> dict[str, list[StateTransition]]:
+        """Group transitions by Mandelstam channel (`.determine_reaction_channel`)."""
+        groupings = defaultdict(list)
+        for transition in self.transitions:
+            groupings[determine_reaction_channel(transition.topology)].append(
+                transition
+            )
+        return dict(sorted(groupings.items()))
diff --git a/src/qrules/workflow.py b/src/qrules/workflow.py
index 80bfa5d2..b39efbb7 100644
--- a/src/qrules/workflow.py
+++ b/src/qrules/workflow.py
@@ -74,6 +74,7 @@
     MutableTransition,
     create_isobar_topologies,
     create_n_body_topology,
+    determine_reaction_channel,
 )
 from qrules.transition import (
     ExecutionInfo,
@@ -438,8 +439,15 @@ def create_problem_sets(  # ruff: ignore[too-many-positional-arguments]
     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."""
+    """Create a `.ProblemSet` collection over all topologies, grouped by strength.
+
+    With :code:`allowed_channels`, only the kinematic permutations whose
+    `.determine_reaction_channel` label is in the given selection (e.g.
+    :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 = [
@@ -451,6 +459,8 @@ def create_problem_sets(  # ruff: ignore[too-many-positional-arguments]
             final_state,
             final_state_groupings,
         )
+        if allowed_channels is None
+        or determine_reaction_channel(permutation) in allowed_channels
         for initial_facts in create_initial_facts(
             permutation,
             initial_state,
@@ -465,6 +475,36 @@ def create_problem_sets(  # ruff: ignore[too-many-positional-arguments]
     return _group_by_strength(problem_sets)
 
 
+def _validate_channels(allowed_channels: Iterable[str] | None) -> set[str] | None:
+    """Normalize a Mandelstam channel selection to `.determine_reaction_channel` labels.
+
+    >>> _validate_channels("t") == {"t"}
+    True
+    >>> _validate_channels(["s", "ut"]) == {"s", "tu"}
+    True
+    >>> _validate_channels(None) is None
+    True
+    >>> _validate_channels(["x"])
+    Traceback (most recent call last):
+        ...
+    ValueError: Invalid Mandelstam channel 'x'. Valid examples: 's', 't', 'u', 'tt'
+    """
+    if allowed_channels is None:
+        return None
+    if isinstance(allowed_channels, str):
+        allowed_channels = [allowed_channels]
+    channels = set()
+    for channel in allowed_channels:
+        if channel != "s" and (not channel or set(channel) - {"t", "u"}):
+            msg = (
+                f"Invalid Mandelstam channel {channel!r}."
+                " Valid examples: 's', 't', 'u', 'tt'"
+            )
+            raise ValueError(msg)
+        channels.add("".join(sorted(channel)))
+    return channels
+
+
 def _group_by_strength(
     problem_sets: list[ProblemSet],
 ) -> dict[float, list[ProblemSet]]:
@@ -726,6 +766,7 @@ def create_qn_problem_sets(  # ruff: ignore[too-many-positional-arguments]
     max_angular_momentum: int = 1,
     max_spin_magnitude: float = 2,
     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:
@@ -750,7 +791,10 @@ def create_qn_problem_sets(  # ruff: ignore[too-many-positional-arguments]
 
     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`.
+    :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`).
     """
     if not spin_projections and merge_spin_projections:
         msg = "merge_spin_projections has no effect when spin_projections=False"
@@ -793,6 +837,7 @@ def create_qn_problem_sets(  # ruff: ignore[too-many-positional-arguments]
         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:
@@ -948,6 +993,15 @@ def group_by_topology(self) -> dict[Topology, list[QNTransition]]:
             groupings[transition.topology].append(transition)
         return dict(groupings)
 
+    def group_by_channel(self) -> dict[str, list[QNTransition]]:
+        """Group transitions by Mandelstam channel (`.determine_reaction_channel`)."""
+        groupings = defaultdict(list)
+        for transition in self.transitions:
+            groupings[determine_reaction_channel(transition.topology)].append(
+                transition
+            )
+        return dict(sorted(groupings.items()))
+
 
 def generate_qn_transitions(  # ruff: ignore[too-many-positional-arguments]
     initial_state: StateDefinitionInput | Sequence[StateDefinitionInput],
@@ -960,6 +1014,7 @@ def generate_qn_transitions(  # ruff: ignore[too-many-positional-arguments]
     max_angular_momentum: int = 1,
     max_spin_magnitude: float = 2,
     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.
@@ -986,6 +1041,7 @@ def generate_qn_transitions(  # ruff: ignore[too-many-positional-arguments]
         max_angular_momentum=max_angular_momentum,
         max_spin_magnitude=max_spin_magnitude,
         final_state_groupings=final_state_groupings,
+        allowed_channels=allowed_channels,
         spin_projections=False,
     )
     transitions = find_qn_transitions(qn_problem_sets, particle_db)
diff --git a/tests/channels/test_gammap_to_p_pi0.py b/tests/channels/test_gammap_to_p_pi0.py
index a9c4c7cc..d195119a 100644
--- a/tests/channels/test_gammap_to_p_pi0.py
+++ b/tests/channels/test_gammap_to_p_pi0.py
@@ -30,3 +30,32 @@ def test_pi0_photoproduction(formalism: SpinFormalism, particle_database):
         "rho(770)0",
     ]
     assert {p.name for p in reaction.initial_state.values()} == {"gamma", "p"}
+
+
+def test_group_by_channel(particle_database):
+    reaction = qrules.generate_transitions(
+        initial_state=["gamma", "p"],
+        final_state=["pi0", "p"],
+        allowed_intermediate_particles=[
+            "Delta(1232)",
+            "N(1440)",
+            "rho(770)",
+            "omega(782)",
+        ],
+        allowed_interaction_types=["strong", "em"],
+        formalism="helicity",
+        particle_db=particle_database,
+    )
+    channels = {
+        channel: {
+            state.particle.name
+            for transition in transitions
+            for state in transition.intermediate_states.values()
+        }
+        for channel, transitions in reaction.group_by_channel().items()
+    }
+    assert channels == {
+        "s": {"Delta(1232)+", "N(1440)+"},
+        "t": {"omega(782)", "rho(770)0"},
+        "u": {"Delta(1232)+", "Delta(1232)~-", "N(1440)+", "N(1440)~-"},
+    }
diff --git a/tests/unit/test_workflow.py b/tests/unit/test_workflow.py
index 9a85373a..476c8fc8 100644
--- a/tests/unit/test_workflow.py
+++ b/tests/unit/test_workflow.py
@@ -1,5 +1,6 @@
 import json
 from fractions import Fraction
+from typing import Any
 
 import pytest
 
@@ -344,6 +345,45 @@ def test_generate_qn_transitions_two_to_n():
     assert len(reaction.group_by_topology()) > 1
 
 
+def test_group_by_channel_and_channel_selection():
+    """Mandelstam channel encoding for 2-to-n reactions (ComPWA/qrules#29)."""
+    particle_db = load_pdg()
+    reaction_kwargs: dict[str, Any] = dict(
+        initial_state=["gamma", "p"],
+        final_state=["pi0", "p"],
+        particle_db=particle_db,
+        allowed_intermediate_particles=[
+            "Delta(1232)",
+            "N(1440)",
+            "rho(770)",
+            "omega(782)",
+        ],
+        allowed_interaction_types=["strong", "em"],
+    )
+    reaction = generate_qn_transitions(**reaction_kwargs)
+    channels = reaction.group_by_channel()
+    assert sorted(channels) == ["s", "t", "u"]
+    assert sum(map(len, channels.values())) == len(reaction.transitions)
+
+    def get_baryon_numbers(channel: str) -> set:
+        return {
+            state[EdgeQuantumNumbers.baryon_number]
+            for transition in channels[channel]
+            for state in transition.intermediate_states.values()
+        }
+
+    assert get_baryon_numbers("s") == {+1}  # Delta/N* resonances
+    assert get_baryon_numbers("t") == {0}  # meson exchange
+    assert get_baryon_numbers("u") == {-1, +1}  # baryon exchange
+
+    t_channel_only = generate_qn_transitions(**reaction_kwargs, allowed_channels="t")
+    assert sorted(t_channel_only.group_by_channel()) == ["t"]
+    assert len(t_channel_only.transitions) == len(channels["t"])
+
+    with pytest.raises(ValueError, match="Invalid Mandelstam channel 'x'"):
+        generate_qn_transitions(**reaction_kwargs, allowed_channels=["x"])
+
+
 def test_qn_reaction_info_requires_particle_states():
     particle_db = load_pdg()
     qn_problem_sets = create_qn_problem_sets(

From faf8fc352b9038f70d4f1f95d44f22950085fd40 Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 09:03:20 +0200
Subject: [PATCH 8/9] DOC: add notebook on production reactions and Mandelstam
 channels

---
 docs/usage.ipynb                |   1 +
 docs/usage/production.ipynb     | 389 ++++++++++++++++++++++++++++++++
 docs/usage/qn-transitions.ipynb | 102 +--------
 3 files changed, 391 insertions(+), 101 deletions(-)
 create mode 100644 docs/usage/production.ipynb

diff --git a/docs/usage.ipynb b/docs/usage.ipynb
index 1a6f9dc3..cee4d0d2 100644
--- a/docs/usage.ipynb
+++ b/docs/usage.ipynb
@@ -204,6 +204,7 @@
     "---\n",
     "usage/reaction\n",
     "usage/qn-transitions\n",
+    "usage/production\n",
     "usage/particle\n",
     "usage/visualize\n",
     "usage/conservation\n",
diff --git a/docs/usage/production.ipynb b/docs/usage/production.ipynb
new file mode 100644
index 00000000..bdfb28a1
--- /dev/null
+++ b/docs/usage/production.ipynb
@@ -0,0 +1,389 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Production reactions\n",
+    "\n",
+    ":::{autolink-concat}\n",
+    ":::"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In addition to decay processes with a single initial state, qrules can generate transitions for **production reactions with two initial states** (see [ComPWA/qrules#29](https://github.com/ComPWA/qrules/issues/29)) — think of photoproduction at GlueX ($\\gamma p$), fixed-target proton scattering at HADES ($pp$), antiproton annihilation at PANDA ($p\\bar p$), or $e^+e^-$ annihilation at BESIII. This page shows how such reactions are built up from topologies, how the resulting transitions are classified into **Mandelstam channels** ($s$, $t$, and $u$), and how to visualize each channel as a single collapsed diagram."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    },
+    "tags": [
+     "hide-cell"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from IPython.display import Markdown\n",
+    "\n",
+    "import qrules\n",
+    "import qrules.io\n",
+    "from qrules.quantum_numbers import EdgeQuantumNumbers\n",
+    "from qrules.topology import create_isobar_topologies, determine_reaction_channel\n",
+    "from qrules.workflow import generate_qn_transitions\n",
+    "\n",
+    "PDG = qrules.load_pdg()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Topologies with two initial states"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For more than one initial state, {func}`.create_isobar_topologies` builds the topologies from $2 \\to 1$ *production* nodes in addition to the usual $1 \\to 2$ *decay* nodes. For a $2 \\to 2$ reaction, this results in two shapes: one where the initial states annihilate into a single intermediate edge and one where they are connected through an exchange edge."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "topologies = create_isobar_topologies(\n",
+    "    number_of_final_states=2,\n",
+    "    number_of_initial_states=2,\n",
+    ")\n",
+    "Markdown(qrules.io.asmermaid(topologies, markdown=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The number of topologies grows quickly with the number of final states, because the exchange edge can attach at several places in the decay chain:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "{\n",
+    "    n_final: len(create_isobar_topologies(n_final, number_of_initial_states=2))\n",
+    "    for n_final in (2, 3, 4, 5)\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Mandelstam channels"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Which Mandelstam variable an intermediate edge carries follows from cutting the topology at that edge. If one side of the cut contains *both* initial states, the edge carries the invariant mass of an $s$-type channel. Otherwise the edge is an *exchange* between the initial states: following the convention $t = (p_1 - p_3)^2$ and $u = (p_1 - p_4)^2$ for $1\\,2 \\to 3\\,4$, it is a $t$-channel if the first initial state and the first final state lie on the same side of the cut and a $u$-channel otherwise. {func}`.determine_mandelstam_channel` implements this classification per edge and {func}`.determine_reaction_channel` summarizes a whole topology:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "[determine_reaction_channel(topology) for topology in topologies]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "As an example, take $\\pi^0$ photoproduction, $\\gamma p \\to \\pi^0 p$ (GlueX). Note how the ordering of the final state matches the Mandelstam convention: the photon pairs with the $\\pi^0$ in the $t$-channel and with the recoil proton in the $u$-channel. The reaction is generated just like a decay — here at the quantum-number level with {func}`.generate_qn_transitions` — and {meth}`.QNReactionInfo.group_by_channel` classifies the transitions:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "reaction = generate_qn_transitions(\n",
+    "    initial_state=[\"gamma\", \"p\"],\n",
+    "    final_state=[\"pi0\", \"p\"],\n",
+    "    particle_db=PDG,\n",
+    "    allowed_intermediate_particles=[\n",
+    "        \"Delta(1232)\",\n",
+    "        \"N(1440)\",\n",
+    "        \"rho(770)\",\n",
+    "        \"omega(782)\",\n",
+    "    ],\n",
+    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
+    ")\n",
+    "{\n",
+    "    channel: len(transitions)\n",
+    "    for channel, transitions in reaction.group_by_channel().items()\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Combined with {func}`.asdot`'s :code:`collapse=\"topology\"` option, each channel reduces to a few minimal diagrams that together still represent the full reaction object. The intermediate edges list the allowed exchanged or resonant states in $I^G(J^{PC})$ notation. The $s$-channel contains the baryonic $\\Delta$ and $N^*$ resonances:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "channels = reaction.group_by_channel()\n",
+    "Markdown(qrules.io.asmermaid(channels[\"s\"], collapse=\"topology\", markdown=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "the $t$-channel proceeds through vector and axial-vector meson exchange ($\\omega/\\rho$-like and $h_1/b_1$-like quantum numbers):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Markdown(qrules.io.asmermaid(channels[\"t\"], collapse=\"topology\", markdown=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "and the $u$-channel exchanges a baryon between the two vertices — the two diagrams are the two orientations of the exchange edge, corresponding to baryon and antibaryon quantum numbers:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Markdown(qrules.io.asmermaid(channels[\"u\"], collapse=\"topology\", markdown=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Selecting channels"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If only certain channels are of interest, pass :code:`allowed_channels` to restrict the problem sets *before* they are solved. This speeds up the generation, since fewer constraint problems have to be solved:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "t_channel_reaction = generate_qn_transitions(\n",
+    "    initial_state=[\"gamma\", \"p\"],\n",
+    "    final_state=[\"pi0\", \"p\"],\n",
+    "    particle_db=PDG,\n",
+    "    allowed_intermediate_particles=[\n",
+    "        \"Delta(1232)\",\n",
+    "        \"N(1440)\",\n",
+    "        \"rho(770)\",\n",
+    "        \"omega(782)\",\n",
+    "    ],\n",
+    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
+    "    allowed_channels=\"t\",\n",
+    ")\n",
+    "len(t_channel_reaction.transitions)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Particle-level transitions"
+   ]
+  },
+  {
+   "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:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "particle_reaction = qrules.generate_transitions(\n",
+    "    initial_state=[\"gamma\", \"p\"],\n",
+    "    final_state=[\"pi0\", \"p\"],\n",
+    "    allowed_intermediate_particles=[\n",
+    "        \"Delta(1232)\",\n",
+    "        \"N(1440)\",\n",
+    "        \"rho(770)\",\n",
+    "        \"omega(782)\",\n",
+    "    ],\n",
+    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
+    "    formalism=\"helicity\",\n",
+    "    particle_db=PDG,\n",
+    ")\n",
+    "{\n",
+    "    channel: sorted({\n",
+    "        state.particle.name\n",
+    "        for transition in transitions\n",
+    "        for state in transition.intermediate_states.values()\n",
+    "    })\n",
+    "    for channel, transitions in particle_reaction.group_by_channel().items()\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Collapsing the full particle-level reaction now labels the intermediate edges by particle name:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Markdown(qrules.io.asmermaid(particle_reaction, collapse=\"topology\", markdown=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Reactions at different energies"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "When there is more than one initial state, the {class}`.MassConservation` rule is switched off, so the collision energy does not enter the problem. Model a specific energy by restricting :code:`allowed_intermediate_particles` to the resonances that are accessible at that energy. For instance, $e^+e^- \\to p\\bar p\\eta$ at BESIII charmonium energies, where the annihilation proceeds through a $J^{PC}=1^{--}$ state and the $p\\eta$ subsystem can form $N^*$ resonances:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "reaction_ee = generate_qn_transitions(\n",
+    "    initial_state=[\"e+\", \"e-\"],\n",
+    "    final_state=[\"p\", \"p~\", \"eta\"],\n",
+    "    particle_db=PDG,\n",
+    "    allowed_intermediate_particles=[\"J/psi(1S)\", \"psi(2S)\", \"N(1535)\"],\n",
+    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
+    ")\n",
+    "len(reaction_ee.transitions)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Only the $s$-channel survives here: a $t$- or $u$-channel exchange would require an intermediate state that carries lepton number. The spins, parities, and baryon numbers of the intermediate edges confirm the photon-like annihilation state and the nucleon resonances:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sorted(reaction_ee.group_by_channel())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "{\n",
+    "    (\n",
+    "        state[EdgeQuantumNumbers.spin_magnitude],\n",
+    "        state[EdgeQuantumNumbers.parity],\n",
+    "        state[EdgeQuantumNumbers.baryon_number],\n",
+    "    )\n",
+    "    for transition in reaction_ee.transitions\n",
+    "    for state in transition.intermediate_states.values()\n",
+    "}"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "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",
+    ":::"
+   ]
+  }
+ ],
+ "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.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/docs/usage/qn-transitions.ipynb b/docs/usage/qn-transitions.ipynb
index 9c85c2ac..950ab4e4 100644
--- a/docs/usage/qn-transitions.ipynb
+++ b/docs/usage/qn-transitions.ipynb
@@ -34,7 +34,6 @@
     "\n",
     "import qrules.io\n",
     "from qrules.quantum_numbers import EdgeQuantumNumbers\n",
-    "from qrules.topology import create_isobar_topologies\n",
     "from qrules.workflow import (\n",
     "    create_qn_problem_sets,\n",
     "    find_qn_transitions,\n",
@@ -377,111 +376,12 @@
     "f\"{len(reaction_5body.transitions):,} transitions over {len(reaction_5body.group_by_topology())} topologies\""
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Production reactions"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Reactions with **two initial states** are supported as well (see [ComPWA/qrules#29](https://github.com/ComPWA/qrules/issues/29)). For more than one initial state, {func}`.create_isobar_topologies` builds the topologies from $2 \\to 1$ production nodes in addition to the usual $1 \\to 2$ decay nodes, which produces both $s$-channel topologies (the initial states annihilate into a single edge) and $t$-channel-like topologies (the initial states are connected through an exchange edge). For a $2 \\to 2$ reaction:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "topologies = create_isobar_topologies(\n",
-    "    number_of_final_states=2,\n",
-    "    number_of_initial_states=2,\n",
-    ")\n",
-    "Markdown(qrules.io.asmermaid(topologies, markdown=True))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This makes it possible to describe production processes like $\\pi^0$ photoproduction, $\\gamma p \\to p\\pi^0$ (GlueX). The allowed transitions contain the $s$-channel $\\Delta$ and $N^*$ resonances as well as $t$-channel $\\rho/\\omega$ exchange and $u$-channel-like baryon exchange:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "%%time\n",
-    "reaction_2to2 = generate_qn_transitions(\n",
-    "    initial_state=[\"gamma\", \"p\"],\n",
-    "    final_state=[\"p\", \"pi0\"],\n",
-    "    particle_db=PDG,\n",
-    "    allowed_intermediate_particles=[\"Delta(1232)\", \"N(1440)\", \"rho(770)\", \"omega(782)\"],\n",
-    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
-    ")\n",
-    "Markdown(qrules.io.asmermaid(reaction_2to2, collapse_graphs=True, markdown=True))"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "When there is more than one initial state, the {class}`.MassConservation` rule is switched off, so the collision energy does not enter the problem. Model a specific energy by restricting :code:`allowed_intermediate_particles` to the resonances that are accessible at that energy. For instance, $e^+e^- \\to p\\bar p\\eta$ at BESIII charmonium energies, where the annihilation proceeds through a $J^{PC}=1^{--}$ state and the $p\\eta$ subsystem can form $N^*$ resonances:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "%%time\n",
-    "reaction_ee = generate_qn_transitions(\n",
-    "    initial_state=[\"e+\", \"e-\"],\n",
-    "    final_state=[\"p\", \"p~\", \"eta\"],\n",
-    "    particle_db=PDG,\n",
-    "    allowed_intermediate_particles=[\"J/psi(1S)\", \"psi(2S)\", \"N(1535)\"],\n",
-    "    allowed_interaction_types=[\"strong\", \"em\"],\n",
-    ")\n",
-    "len(reaction_ee.transitions)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Only $s$-channel topologies survive here: a $t$-channel exchange would require an intermediate state that carries lepton number. The spins, parities, and baryon numbers of the intermediate edges confirm the photon-like annihilation state and the nucleon resonances:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "{\n",
-    "    (\n",
-    "        state[EdgeQuantumNumbers.spin_magnitude],\n",
-    "        state[EdgeQuantumNumbers.parity],\n",
-    "        state[EdgeQuantumNumbers.baryon_number],\n",
-    "    )\n",
-    "    for transition in reaction_ee.transitions\n",
-    "    for state in transition.intermediate_states.values()\n",
-    "}"
-   ]
-  },
   {
    "cell_type": "markdown",
    "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`.\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",
     ":::"
    ]
   }

From 662eb3fff52007ee7b60d8eac3bce2fd185c512c Mon Sep 17 00:00:00 2001
From: Remco de Boer <29308176+redeboer@users.noreply.github.com>
Date: Mon, 13 Jul 2026 09:12:09 +0200
Subject: [PATCH 9/9] ENH: align exchange vertices at the same rank in DOT
 rendering

---
 docs/usage/production.ipynb | 11 ++++++++++-
 src/qrules/io/_dot.py       | 19 ++++++++++++++++++-
 2 files changed, 28 insertions(+), 2 deletions(-)

diff --git a/docs/usage/production.ipynb b/docs/usage/production.ipynb
index bdfb28a1..af0a8e4d 100644
--- a/docs/usage/production.ipynb
+++ b/docs/usage/production.ipynb
@@ -319,6 +319,15 @@
     "len(reaction_ee.transitions)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "Markdown(qrules.io.asmermaid(reaction_ee, collapse=\"topology\", markdown=True))"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -381,7 +390,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.12"
+   "version": "3.13.14"
   }
  },
  "nbformat": 4,
diff --git a/src/qrules/io/_dot.py b/src/qrules/io/_dot.py
index 412a05e4..7d2f3249 100644
--- a/src/qrules/io/_dot.py
+++ b/src/qrules/io/_dot.py
@@ -15,7 +15,7 @@
 
 from qrules.io import _labels
 from qrules.solving import QNProblemSet, QNResult
-from qrules.topology import Topology, Transition
+from qrules.topology import Topology, Transition, determine_mandelstam_channel
 from qrules.transition import ProblemSet, ReactionInfo
 
 if TYPE_CHECKING:
@@ -112,6 +112,7 @@ def _render_transition(  # ruff: ignore[complex-structure, too-many-branches]
             lines += [self._create_graphviz_node(graphviz_node, label, self.edge_style)]
         lines += [_create_same_rank_line(topology.incoming_edge_ids, prefix)]
         lines += [_create_same_rank_line(topology.outgoing_edge_ids, prefix)]
+        lines += _create_exchange_same_rank_lines(topology, prefix)
         for i, edge in topology.edges.items():
             j, k = edge.ending_node_id, edge.originating_node_id
             from_node = prefix + _get_graphviz_node(i, k)
@@ -222,3 +223,19 @@ def _create_same_rank_line(node_edge_ids: Iterable[int], prefix: str = "") -> st
     name_list = [f"{prefix}{_get_graphviz_node(i)}" for i in node_edge_ids]
     name_string = " ".join(name_list)
     return f"{{ rank=same; {name_string} }}"
+
+
+def _create_exchange_same_rank_lines(topology: Topology, prefix: str = "") -> list[str]:
+    """Align the two nodes of each exchange edge at the same rank.
+
+    Without this constraint, Graphviz places the vertices of :math:`t`- and
+    :math:`u`-channel exchanges at subsequent ranks, which makes them look
+    time-ordered like a decay chain.
+    """
+    return [
+        f"{{ rank=same; {prefix}N{edge.originating_node_id}"
+        f" {prefix}N{edge.ending_node_id} }}"
+        for edge_id in sorted(topology.intermediate_edge_ids)
+        for edge in [topology.edges[edge_id]]
+        if determine_mandelstam_channel(topology, edge_id) != "s"
+    ]