From cae3da848b41d3c160b3a01f54fc1fd0b17a598b Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Fri, 9 Oct 2026 06:24:02 +0200
Subject: [PATCH 1/5] Interference processes: [LIxtree=QCD] and [treextree]
 left/right syntax

Two processes on one line, of which only the interference
2 Re(A_left A_right^*) is computed:

  generate g g > h > t t~ [LIxtree=QCD] g g > t t~    (loop-induced x tree)
  generate g g > h > t t~ [LIxtree=QCD]               (right side implied)
  generate p p > z > e+ e- [treextree] p p > a > e+ e-  (tree x tree)

This automates the recipe of the wiki page LoopInducedTimesTree (UVtree
copies of the vertices in a special model, BKGQCD order, five source
patches) without exposing it.

Syntax (madgraph_interface, master_interface, base_objects)
- keyword case-insensitive; each side takes its own s-channels,
  exclusions and amplitude orders; squared-order constraints constrain
  the interference and go after ']' (refused on the right-hand side);
  both sides must have the same legs in the same order (right may be
  wider); LIxtree without right side = tree process of the same legs.
- Process keys interference_mode / interference_process, carried by
  get_process_with_legs, printed back in the input syntax.

Generation
- A hidden coupling order INTERF (hierarchy 0, added to the model in
  memory, never printed) tags the right-hand amplitudes; INTERF^2==1
  and INTERF as split order keep the left x right products only, also in
  the generated code.
- LIxtree: the right-hand trees are generated per subprocess and stored
  as loop_UVCT_diagrams (UVCT_orders {INTERF:1}, unit coupling); the
  ref_orders of a Born-less loop process now include the UVCT amplitudes
  (without it MadLoop stops: "Could not find squared orders").
  Subprocesses without surviving loops are dropped with their trees.
- treextree: the right-hand diagrams are appended after the left-hand
  generation and the squared-order filter is applied to the union; HELAS
  copies the tag onto the amplitude orders; no crossing reuse.
- treextree output limited to madevent and standalone_fortran (the C++
  outputs have no squared split orders and would return |L+R|^2).

Run card and systematics
- interference processes get the existing interference defaults
  (dynamical_scale_choice 3, no systematics by default);
- config_nqcd.inc: one alpha_s power for all channels, half the QCD order
  of the left x right products (QCD x EW u d > u d: 1, not 2 or 0);
- systematics summary: no scale line and wrong-signed percentages for a
  negative total cross-section (bounds started at 0/maxsize).

Validation (tests/acceptance_tests/test_loop_induced_x_tree.py)
- g g > h > t t~ [LIxtree=QCD] at a fixed point: interference equal to
  the wiki recipe, loop^2 to [sqrvirt=QCD], tree^2 (helicity filter off)
  to the tree-level ME, to 10 digits;
- u u~ Z x photon: equal to |Z+a|^2 - |Z|^2 - |a|^2 to 10 digits.
madevent: mH=500 GeV LIxtree -0.957(8) pb vs wiki recipe -0.965(2) pb;
p p > e+ e-, mll>200: treextree 0.0598(28) pb vs 0.0629(48) pb from the
three squared samples.
---
 madgraph/core/base_objects.py                 |  76 +++-
 madgraph/core/diagram_generation.py           | 147 +++++++-
 madgraph/core/helas_objects.py                |  16 +
 madgraph/interface/madgraph_interface.py      | 326 +++++++++++++++++-
 madgraph/interface/master_interface.py        |  16 +-
 madgraph/iolibs/export_v4.py                  |  79 +++++
 madgraph/loop/loop_base_objects.py            |  10 +
 madgraph/loop/loop_diagram_generation.py      |  64 +++-
 madgraph/loop/loop_helas_objects.py           |   8 +-
 madgraph/various/banner.py                    |   7 +-
 madgraph/various/systematics.py               |  55 +--
 .../test_loop_induced_x_tree.py               | 166 +++++++++
 tests/unit_tests/core/test_base_objects.py    |  20 +-
 .../unit_tests/core/test_tree_interference.py | 164 +++++++++
 .../interface/test_interference_syntax.py     | 247 +++++++++++++
 .../loop/test_loop_induced_x_tree.py          | 172 +++++++++
 16 files changed, 1527 insertions(+), 46 deletions(-)
 create mode 100644 tests/acceptance_tests/test_loop_induced_x_tree.py
 create mode 100644 tests/unit_tests/core/test_tree_interference.py
 create mode 100644 tests/unit_tests/interface/test_interference_syntax.py
 create mode 100644 tests/unit_tests/loop/test_loop_induced_x_tree.py

diff --git a/madgraph/core/base_objects.py b/madgraph/core/base_objects.py
index 94a0229c24..8f3b95edd0 100755
--- a/madgraph/core/base_objects.py
+++ b/madgraph/core/base_objects.py
@@ -34,6 +34,14 @@
 
 logger = logging.getLogger('madgraph.base_objects')
 pjoin = os.path.join
+
+# Coupling order carried only by the right-hand (tree-level) amplitudes of an
+# interference process ('left [LIxtree=QCD] right'): the squared-order
+# constraint INTERF^2==1 then keeps the products of one left-hand and one
+# right-hand amplitude, i.e. 2 Re(A_left A_right^*), and nothing else. It is
+# added to the model in memory with hierarchy 0, so that WEIGHTED is
+# unchanged, and is never printed in a process definition.
+INTERFERENCE_ORDER = 'INTERF'
 if madgraph.ordering:
     set = misc.OrderedSet
 
@@ -3738,6 +3746,13 @@ def default_setup(self):
         # contribution 'QCD=4 QED=0', the pure interference 'QCD=2 QED=2' and
         # the pure QED contribution of order 'QCD=0 QED=4'.
         self['split_orders'] = []
+        # Interference of two amplitude sets, 2Re(A_left A_right*), written
+        # as 'left [LIxtree=ORDERS] right' (loop-induced times tree) or
+        # 'left [treextree] right'. 'interference_mode' is '' for an ordinary
+        # process; 'interference_process' is then the right-hand process,
+        # carrying its own orders, s-channel requirements and exclusions.
+        self['interference_mode'] = ''
+        self['interference_process'] = None
 
     def filter(self, name, value):
         """Filter for valid process property values."""
@@ -3827,6 +3842,14 @@ def filter(self, name, value):
             import madgraph.interface.madgraph_interface as mg
             if value not in mg.MadGraphCmd._valid_nlo_modes:
                 raise self.PhysicsObjectError("%s is not a valid NLO_mode" % str(value))
+
+        if name == 'interference_mode':
+            if value not in ['', 'LIxtree', 'treextree']:
+                raise self.PhysicsObjectError("%s is not a valid interference_mode" % str(value))
+
+        if name == 'interference_process':
+            if value is not None and not isinstance(value, Process):
+                raise self.PhysicsObjectError("%s is not a valid Process" % str(value))
         return True
 
     def has_multiparticle_label(self):
@@ -3913,7 +3936,29 @@ def get_sorted_keys(self):
                 'forbidden_onsh_s_channels', 'forbidden_s_channels',
                 'forbidden_particles', 'is_decay_chain', 'decay_chains',
                 'legs_with_decays', 'perturbation_couplings', 'has_born', 
-                'NLO_mode', 'split_orders', 'born_sq_orders']
+                'NLO_mode', 'split_orders', 'born_sq_orders',
+                'interference_mode', 'interference_process']
+
+    def get_interference_mode(self):
+        """'LIxtree', 'treextree' or '' (also for a process unpickled from
+        a version that did not know about interference processes)."""
+
+        return dict.get(self, 'interference_mode', '') or ''
+
+    def interference_string(self):
+        """The tail of the process syntax for an interference process: the
+        '[ treextree ]' bracket (the LIxtree one is printed with the
+        perturbation orders) followed by the right-hand process, each piece
+        ending with a space. Empty for an ordinary process."""
+
+        mode = self.get_interference_mode()
+        if not mode:
+            return ''
+        mystr = '[ treextree ] ' if mode == 'treextree' else ''
+        right = self['interference_process']
+        if right is not None:
+            mystr += right.nice_string(prefix=False, print_weighted=False) + ' '
+        return mystr
 
     def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturbated=True):
         """Returns a nicely formated string about current process
@@ -3973,6 +4018,8 @@ def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturba
             for key in sorted(self['orders'].keys()):
                 if not print_weighted and key == 'WEIGHTED':
                     continue
+                if key == INTERFERENCE_ORDER:
+                    continue
                 value = int(self['orders'][key])
                 if key in self['squared_orders']:
                     if self.get_squared_order_type(key) in ['<=', '==', '='] and \
@@ -4003,6 +4050,8 @@ def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturba
             if self['NLO_mode']!='tree':
                 if self['NLO_mode']=='virt' and not self['has_born']:
                     mystr = mystr + 'sqrvirt = '
+                elif self.get_interference_mode() == 'LIxtree':
+                    mystr = mystr + 'LIxtree = '
                 else:
                     mystr = mystr + self['NLO_mode'] + ' = '
             for order in sorted(self['perturbation_couplings']):
@@ -4015,6 +4064,8 @@ def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturba
             for key in sorted(self['squared_orders'].keys()):
                 if not print_weighted and key == 'WEIGHTED':
                     continue
+                if key == INTERFERENCE_ORDER:
+                    continue
                 if key in self['constrained_orders']:
                     if self['constrained_orders'][key][0] == self['squared_orders'][key]/2 and \
                        self['constrained_orders'][key][1] == self.get_squared_order_type(key):
@@ -4047,6 +4098,8 @@ def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturba
                 forbpart = self['model'].get('particle_dict')[forb_id]
                 mystr = mystr + forbpart.get_name() + ' '
 
+        mystr = mystr + self.interference_string()
+
         # Remove last space
         mystr = mystr[:-1]
 
@@ -4104,7 +4157,7 @@ def input_string(self):
             prevleg = leg
 
         if self['orders']:
-            keys = list(self['orders'].keys())
+            keys = [k for k in self['orders'] if k != INTERFERENCE_ORDER]
             keys.sort(reverse=True)
             mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
                        for key in keys]) + ' '
@@ -4112,7 +4165,8 @@ def input_string(self):
         # Add squared orders
         if self['squared_orders']:
             mystr = mystr + " ".join([key + '^2=' + repr(self['squared_orders'][key]) \
-                       for key in self['squared_orders']]) + ' '
+                       for key in self['squared_orders']
+                       if key != INTERFERENCE_ORDER]) + ' '
 
         # Add perturbation orders
         if self['perturbation_couplings']:
@@ -4923,7 +4977,8 @@ def nice_string(self, indent=0, print_weighted=False, prefix=True):
 
         if self['orders']:
             mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
-                       for key in sorted(self['orders'])]) + ' '
+                       for key in sorted(self['orders'])
+                       if key != INTERFERENCE_ORDER]) + ' '
 
         if self['constrained_orders']:
             mystr = mystr + " ".join('%s%s%d' % (key, operator, value) for 
@@ -4936,6 +4991,8 @@ def nice_string(self, indent=0, print_weighted=False, prefix=True):
             if self['NLO_mode']!='tree':
                 if self['NLO_mode']=='virt' and not self['has_born']:
                     mystr = mystr + 'sqrvirt = '
+                elif self.get_interference_mode() == 'LIxtree':
+                    mystr = mystr + 'LIxtree = '
                 else:
                     mystr = mystr + self['NLO_mode'] + ' = '
             for order in self['perturbation_couplings']:
@@ -4946,7 +5003,10 @@ def nice_string(self, indent=0, print_weighted=False, prefix=True):
             mystr = mystr + " ".join([key + '^2%s%d'%\
                 (self.get_squared_order_type(key),self['squared_orders'][key]) \
               for key in self['squared_orders'].keys() \
-                                    if print_weighted or key!='WEIGHTED']) + ' '
+                    if (print_weighted or key!='WEIGHTED') and \
+                                            key != INTERFERENCE_ORDER]) + ' '
+
+        mystr = mystr + self.interference_string()
 
         # Remove last space
         mystr = mystr[:-1]
@@ -4989,7 +5049,11 @@ def get_process_with_legs(self, LegList):
             'overall_orders': self.get('overall_orders'),
             'split_orders': self.get('split_orders'),
             'born_sq_orders': self.get('born_sq_orders'),
-            'NLO_mode': self.get('NLO_mode')
+            'NLO_mode': self.get('NLO_mode'),
+            # the right-hand process stays the ProcessDefinition template:
+            # its legs are matched to these ones at generation time
+            'interference_mode': self.get_interference_mode(),
+            'interference_process': dict.get(self, 'interference_process')
             })
             
     def get_process(self, initial_state_ids, final_state_ids):
diff --git a/madgraph/core/diagram_generation.py b/madgraph/core/diagram_generation.py
index fcd1653229..7726faa968 100755
--- a/madgraph/core/diagram_generation.py
+++ b/madgraph/core/diagram_generation.py
@@ -41,6 +41,87 @@
 
 class NoDiagramException(InvalidCmd): pass
 
+#===============================================================================
+# Interference of two amplitude sets: 'left [LIxtree=QCD] right' (loop-induced
+# times tree, see loop_diagram_generation) and 'left [treextree] right'
+#===============================================================================
+INTERFERENCE_ORDER = base_objects.INTERFERENCE_ORDER
+# Output formats able to select squared split orders, hence to return the
+# interference alone (the others would give |A_left + A_right|^2)
+TREE_INTERFERENCE_FORMATS = ['madevent', 'standalone_fortran']
+
+def register_interference_order(model):
+    """Add the hidden coupling order tagging the right-hand amplitudes of an
+    interference process to the model, in memory and idempotently. Its
+    hierarchy is 0, so that the WEIGHTED order of every diagram is unchanged,
+    and no interaction carries it."""
+
+    if INTERFERENCE_ORDER not in model.get('order_hierarchy'):
+        model.get('order_hierarchy')[INTERFERENCE_ORDER] = 0
+    coupling_orders = model.get('coupling_orders')
+    if INTERFERENCE_ORDER not in coupling_orders:
+        coupling_orders.add(INTERFERENCE_ORDER)
+
+def prepare_interference_process(procdef):
+    """Complete the ProcessDefinition of an interference process
+    ('left [LIxtree=QCD] right' or 'left [treextree] right', see
+    MadGraphCmd.extract_interference_process) before its generation:
+
+    - the hidden order INTERFERENCE_ORDER is registered in the model;
+    - the left-hand process gets the squared-order constraint INTERF^2==1
+      and INTERF as split order, so that only the products of a left-hand
+      amplitude with a right-hand one are kept, also in the generated code;
+    - the right-hand process gets, when none is given, the minimal WEIGHTED
+      order exactly as 'generate <right-hand process>' would.
+
+    The right-hand diagrams themselves are generated per subprocess, with the
+    legs of that subprocess: Amplitude.add_interference_diagrams (treextree)
+    and LoopAmplitude.set_interference_trees (LIxtree)."""
+
+    mode = procdef.get_interference_mode()
+    if not mode:
+        return
+    if any(value < 0 for value in procdef.get('squared_orders').values()):
+        raise InvalidCmd("Negative squared-order constraints are not "
+                         "supported for %s interference processes." % mode)
+    model = procdef.get('model')
+    register_interference_order(model)
+
+    squared_orders = dict(procdef.get('squared_orders'))
+    squared_orders[INTERFERENCE_ORDER] = 1
+    procdef.set('squared_orders', squared_orders)
+    sqorders_types = dict(procdef.get('sqorders_types'))
+    sqorders_types[INTERFERENCE_ORDER] = '=='
+    procdef.set('sqorders_types', sqorders_types)
+    if INTERFERENCE_ORDER not in procdef.get('split_orders'):
+        procdef.set('split_orders',
+                    list(procdef.get('split_orders')) + [INTERFERENCE_ORDER])
+
+    right = procdef.get('interference_process')
+    if not right.get('orders') and not right.get('born_sq_orders'):
+        right.set('orders', MultiProcess.find_optimal_process_orders(right))
+    right.check_expansion_orders()
+
+def specialise_interference_process(process):
+    """The right-hand process of an interference subprocess, with the
+    external legs of that subprocess. It is also stored back in the
+    subprocess, which then shows the tree process it interferes with."""
+
+    right = process.get('interference_process')
+    legs = base_objects.LegList([copy.copy(leg) for leg in process.get('legs')])
+    if isinstance(right, base_objects.ProcessDefinition):
+        right_process = right.get_process_with_legs(legs)
+    else:
+        # already specialised to these legs (regeneration)
+        right_process = copy.copy(right)
+        right_process.set('legs', legs)
+    right_process.set('interference_mode', '')
+    right_process.set('interference_process', None)
+    # the process number is the one of the whole (left-hand) definition
+    right_process.set('id', 0)
+    process.set('interference_process', right_process)
+    return right_process
+
 #===============================================================================
 # DiagramTag mother class
 #===============================================================================
@@ -1267,7 +1348,17 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False):
         # in this way. We shall do this only if the diagrams are not asked to
         # be returned, as it is the case for NLO because it this case the
         # interference are not necessarily among the diagrams generated here only.
-        if not returndiag and len(res)>0:
+        # 'left [treextree] right': the squared-order constraints are about
+        # the products with the right-hand diagrams, added at the end (see
+        # add_interference_diagrams); only the amplitude-level constraints of
+        # this (left-hand) process apply here
+        treextree = not returndiag and \
+                         process.get_interference_mode() == 'treextree'
+        if treextree:
+            for name, (value, operator) in \
+                                  process.get('constrained_orders').items():
+                res.filter_constrained_orders(name, value, operator)
+        elif not returndiag and len(res)>0:
             res = self.apply_squared_order_constraints(res)
 
         if diagram_filter:
@@ -1300,6 +1391,9 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False):
             pertur = sorted(self.get('process')['perturbation_couplings'])[0]
         self.get('process').get('legs').sort(pert=pertur)
 
+        if treextree:
+            res = self.add_interference_diagrams(res)
+
         # Set diagrams to res if not asked to be returned
         if not returndiag:
            self['diagrams'] = res
@@ -1307,6 +1401,52 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False):
         else:
            return not failed_crossing, res
 
+    def add_interference_diagrams(self, left_diagrams):
+        """For an interference process 'left [treextree] right', add to the
+        diagrams of this (left-hand) process those of the right-hand process
+        with the same external legs, each tagged with one unit of the hidden
+        order INTERFERENCE_ORDER, and apply the squared-order constraints of
+        the process to the whole list: with INTERF^2==1 (see
+        prepare_interference_process) only the products of a left-hand and a
+        right-hand diagram, 2 Re(A_left A_right^*), contribute. All diagrams of
+        both sides are kept as integration channels. Returns the new diagram
+        list, empty if either side has no diagram."""
+
+        process = self.get('process')
+        register_interference_order(process.get('model'))
+        right_process = specialise_interference_process(process)
+
+        right_diagrams = base_objects.DiagramList()
+        if left_diagrams:
+            try:
+                right_diagrams = Amplitude(right_process).get('diagrams')
+            except InvalidCmd:
+                pass
+        if not left_diagrams or not right_diagrams:
+            logger.info("No diagram for %s: no interference with %s." % (
+                (right_process if left_diagrams else process).nice_string(
+                                                               prefix=False),
+                (process if left_diagrams else right_process).nice_string(
+                                                               prefix=False)))
+            return base_objects.DiagramList()
+
+        for diagram in right_diagrams:
+            orders = dict(diagram.get('orders'))
+            orders[INTERFERENCE_ORDER] = 1
+            diagram.set('orders', orders)
+
+        res = base_objects.DiagramList(list(left_diagrams) + list(right_diagrams))
+        while True:
+            new_res = res.apply_positive_sq_orders(res,
+                                          process.get('squared_orders'),
+                                          process.get('sqorders_types'))
+            if len(new_res) == len(res):
+                break
+            res = new_res
+        logger.info("Interference of %d with %d diagrams (%d contributing)" %
+                    (len(left_diagrams), len(right_diagrams), len(res)))
+        return res
+
     def apply_squared_order_constraints(self, diag_list):
         """Applies the user specified squared order constraints on the diagram
         list in argument."""
@@ -2683,10 +2823,13 @@ def get_fs_flavor(id, fsleg):
                         
                 # Check for successful crossings, unless we have specified
                 # properties that break crossing symmetry
+                # (an interference process carries a right-hand process with
+                # its own s-channel constraints: no reuse by crossing either)
                 if not process.get('required_s_channels') and \
                    not process.get('forbidden_onsh_s_channels') and \
                    not process.get('forbidden_s_channels') and \
-                   not process.get('is_decay_chain') and not diagram_filter:
+                   not process.get('is_decay_chain') and not diagram_filter \
+                   and not process.get_interference_mode():
                     try:
                         crossed_index = success_procs.index(sorted_legs)
                         # The relabeling of legs for loop amplitudes is cumbersome
diff --git a/madgraph/core/helas_objects.py b/madgraph/core/helas_objects.py
index 3d8bb1daaa..2cc22d46ed 100755
--- a/madgraph/core/helas_objects.py
+++ b/madgraph/core/helas_objects.py
@@ -4471,6 +4471,22 @@ def generate_helas_diagrams(self, amplitude, optimization=1,decay_ids=[]):
                     # Add amplitude to amplitdes in helas_diagram
                     helas_diagram.get('amplitudes').append(amp)
 
+            # A right-hand diagram of an interference process
+            # ('left [treextree] right') carries the hidden order tagging it,
+            # which no interaction has: pass it on to its amplitudes (on a
+            # copy, their orders being those of the interaction), so that the
+            # split orders of the generated code see it.
+            interference_tag = diagram.get('orders').get(
+                                    base_objects.INTERFERENCE_ORDER, 0) \
+                               if diagram.get('orders') else 0
+            if interference_tag:
+                for amp in helas_diagram.get('amplitudes'):
+                    orders = dict(amp.get('orders'))
+                    orders[base_objects.INTERFERENCE_ORDER] = \
+                        orders.get(base_objects.INTERFERENCE_ORDER, 0) + \
+                                                               interference_tag
+                    amp.set('orders', orders)
+
             # After generation of all wavefunctions and amplitudes,
             # first sort the wavefunctions according to number
             diagram_wavefunctions.sort(key=lambda wf:wf.get('number')) 
diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py
index 42c481a734..b6ad789886 100755
--- a/madgraph/interface/madgraph_interface.py
+++ b/madgraph/interface/madgraph_interface.py
@@ -144,6 +144,113 @@ class DuplicateParticle(madgraph.InvalidCmd):
     same misleading message.
     """
     
+#===============================================================================
+# Interference of two amplitude sets: 'left [LIxtree=QCD] right' and
+# 'left [treextree] right'
+#===============================================================================
+# The bracket of an interference process. Its keyword is matched without
+# regard to case, so that [LIxtree=QCD], [lixtree=QCD] and [TreeXTree] are all
+# accepted; INTERFERENCE_MODES maps the lower-case spelling on the one stored
+# in Process['interference_mode'].
+INTERFERENCE_MODES = {'lixtree': 'LIxtree', 'treextree': 'treextree'}
+_INTERFERENCE_BRACKET = re.compile(
+    r"\[\s*(?:(?P<option>\w+)\s*=)?\s*(?P<orders>[^\]\[]*?)\s*\]")
+# One coupling-order constraint at the start of a string: 'QED=0', 'QCD^2<=4',
+# 'QED == 2', ... The constraints written between the bracket and the
+# right-hand process belong to the left-hand one (the usual place of the
+# squared-order constraints of a loop process), those written after the
+# right-hand process belong to the right-hand one.
+_LEADING_ORDER_CONSTRAINT = re.compile(
+    r"^\s*(?:\w|\^2)+\s*(?:===|==|<=|>=|!=|=|<|>)\s*-?\d+(?=\s|$)")
+
+
+def interference_mode_of(option, orders):
+    """The interference mode named by a bracket, '' if it names none.
+    'option' is the word before '=' (None without '='), 'orders' the list of
+    words after it."""
+
+    if option and option.lower() in INTERFERENCE_MODES:
+        return INTERFERENCE_MODES[option.lower()]
+    # bare keyword: '[treextree]', or '[LIxtree]' (refused later: no orders)
+    if not option and len(orders) == 1 and orders[0].lower() in INTERFERENCE_MODES:
+        return INTERFERENCE_MODES[orders[0].lower()]
+    return ''
+
+
+def split_interference_line(line):
+    """Split a process definition of the form
+
+        left_process [LIxtree=ORDERS] left_constraints right_process
+        left_process [treextree] left_constraints right_process
+
+    into (left_line, mode, right_line).
+
+    'left_line' is the left-hand process written in the ordinary syntax: with
+    '[noborn=ORDERS]' for LIxtree (the loop-induced amplitude) and without any
+    bracket for treextree. 'mode' is 'LIxtree' or 'treextree', and
+    'right_line' the right-hand (tree-level) process, None when it is absent.
+    A trailing '@N' and what follows it stay with the left-hand line, which
+    carries the process number of the whole definition.
+
+    A line without an interference bracket is returned as (line, '', None).
+    Raises InvalidCmd for a malformed interference bracket.
+    """
+
+    match = _INTERFERENCE_BRACKET.search(line)
+    if not match:
+        return line, '', None
+    option = match.group('option')
+    orders = match.group('orders').split()
+    mode = interference_mode_of(option, orders)
+    if not mode:
+        return line, '', None
+    if not option:
+        # the only word in the bracket is the keyword itself
+        orders = []
+
+    head, tail = line[:match.start()], line[match.end():]
+    if '[' in tail or ']' in tail:
+        raise madgraph.InvalidCmd(
+            "Only one '[...]' is allowed in a %s process: the right-hand "
+            "process is a tree-level process." % mode)
+
+    if mode == 'LIxtree':
+        if not orders:
+            raise madgraph.InvalidCmd(
+                "The LIxtree mode needs the orders of the loop, as in "
+                "'g g > h > t t~ [LIxtree=QCD]'.")
+        bracket = ' [noborn= %s ] ' % ' '.join(orders)
+    else:
+        if option:
+            raise madgraph.InvalidCmd(
+                "The treextree mode takes no orders: write '[treextree]', "
+                "not '[%s=%s]'." % (option, ' '.join(orders)))
+        bracket = ' '
+
+    # '@N ...' closes the whole definition: keep it with the left-hand line
+    suffix = ''
+    at = re.search(r'@\s*\d+', tail)
+    if at:
+        tail, suffix = tail[:at.start()], ' ' + tail[at.start():]
+
+    left_constraints = []
+    tail = tail.strip()
+    constraint = _LEADING_ORDER_CONSTRAINT.match(tail)
+    while constraint:
+        left_constraints.append(constraint.group(0).strip())
+        tail = tail[constraint.end():].strip()
+        constraint = _LEADING_ORDER_CONSTRAINT.match(tail)
+
+    right = tail if tail else None
+    if right is None and mode == 'treextree':
+        raise madgraph.InvalidCmd(
+            "The treextree mode needs the right-hand process, as in "
+            "'p p > z > e+ e- [treextree] p p > a > e+ e-'.")
+
+    left = head.rstrip() + bracket + ' '.join(left_constraints) + suffix
+    return ' '.join(left.split()), mode, right
+
+
 # the same prompt without the colour escapes, for quoting commands inside
 # tutorial text and help messages
 MG7_PROMPT_TEXT = "MG7> "
@@ -826,6 +933,21 @@ def help_generate(self):
         logger.info(" > For processes without born amplitudes (i.e. loop-induced like g g > z), please use ")
         logger.info("   the 'virt=' NLO mode. aMC@NLO cannot integrate these processes, but standalone MadLoop5")
         logger.info("   can still handle these.")
+        logger.info("Interference syntax:",'$MG:BOLD')
+        logger.info(" o left process [LIxtree=LoopOrders] SQUAREDCOUPi=ORDERi right process")
+        logger.info(" o left process [treextree] SQUAREDCOUPi=ORDERi right process")
+        logger.info(" o Example: generate g g > h > t t~ [LIxtree=QCD] g g > t t~",'$MG:color:GREEN')
+        logger.info(" o Example: generate p p > z > e+ e- [treextree] p p > a > e+ e-",'$MG:color:GREEN')
+        logger.info(" > Only the interference 2 Re(A_left A_right^*) is computed. The left process is")
+        logger.info("   a loop-induced process for LIxtree, a tree-level one for treextree; the right")
+        logger.info("   process is always tree level. The keywords are case-insensitive.")
+        logger.info(" > Each side takes its own s-channel requirements, exclusions and amplitude orders.")
+        logger.info("   Squared-order constraints apply to the interference: give them right after ']'.")
+        logger.info(" > Both sides must have the same external legs, in the same order (the right-hand")
+        logger.info("   legs may be wider multiparticles). For LIxtree the right process is optional:")
+        logger.info("   by default it is the tree-level process with the same legs.")
+        logger.info(" > 'output madevent' integrates the interference (events with signed weights);")
+        logger.info("   for LIxtree, 'output standalone_fortran' evaluates it at a phase-space point.")
 
     def help_add(self):
         logger.info("-- generate diagrams for a process and add to existing processes",'$MG:color:BLUE')
@@ -868,6 +990,21 @@ def help_add(self):
         logger.info(" > For processes without born amplitudes (i.e. loop-induced like g g > z), please use ")
         logger.info("   the 'virt=' NLO mode. aMC@NLO cannot integrate these processes, but standalone MadLoop5")
         logger.info("   can still handle these.")
+        logger.info("Interference syntax:",'$MG:BOLD')
+        logger.info(" o left process [LIxtree=LoopOrders] SQUAREDCOUPi=ORDERi right process")
+        logger.info(" o left process [treextree] SQUAREDCOUPi=ORDERi right process")
+        logger.info(" o Example: generate g g > h > t t~ [LIxtree=QCD] g g > t t~",'$MG:color:GREEN')
+        logger.info(" o Example: generate p p > z > e+ e- [treextree] p p > a > e+ e-",'$MG:color:GREEN')
+        logger.info(" > Only the interference 2 Re(A_left A_right^*) is computed. The left process is")
+        logger.info("   a loop-induced process for LIxtree, a tree-level one for treextree; the right")
+        logger.info("   process is always tree level. The keywords are case-insensitive.")
+        logger.info(" > Each side takes its own s-channel requirements, exclusions and amplitude orders.")
+        logger.info("   Squared-order constraints apply to the interference: give them right after ']'.")
+        logger.info(" > Both sides must have the same external legs, in the same order (the right-hand")
+        logger.info("   legs may be wider multiparticles). For LIxtree the right process is optional:")
+        logger.info("   by default it is the tree-level process with the same legs.")
+        logger.info(" > 'output madevent' integrates the interference (events with signed weights);")
+        logger.info("   for LIxtree, 'output standalone_fortran' evaluates it at a phase-space point.")
 
         logger.info("--  merge two model to create a new one", '$MG:color:BLUE')
         logger.info("syntax:",'$MG:BOLD')
@@ -1447,6 +1584,13 @@ def check_generate(self, args):
     def check_process_format(self, process):
         """ check the validity of the string given to describe a format """
 
+        left, mode, right = split_interference_line(process)
+        if mode:
+            self.check_process_format(left)
+            if right is not None:
+                self.check_process_format(right)
+            return
+
         #check balance of paranthesis
         if process.count('(') != process.count(')'):
             raise self.InvalidCmd('Invalid Format, no balance between open and close parenthesis')
@@ -3601,7 +3745,8 @@ def _tutorial_opts(self):
                     'merge_quartic_vertices'
                     ]
     _valid_color_basis = ['auto', 'trace', 'ddm']
-    _valid_nlo_modes = ['all','real','virt','sqrvirt','tree','noborn','LOonly', 'only']
+    _valid_nlo_modes = ['all','real','virt','sqrvirt','tree','noborn','LOonly', 'only',
+                        'LIxtree']
     _valid_sqso_types = ['==','<=','=','>']
     _valid_amp_so_types = ['=','<=', '==', '>']
     _OLP_supported = ['MadLoop', 'GoSam']
@@ -3914,6 +4059,9 @@ def do_add(self, line, counter=0):
                 nb_proc = len([l for l in self.history if l.startswith(('generate','add process'))])
                 nb_proc += counter
                 myprocdef = self.extract_process(line, proc_number=nb_proc)
+                # 'left [treextree] right': hidden order selecting the
+                # interference, minimal orders of the right-hand process
+                diagram_generation.prepare_interference_process(myprocdef)
 
             
 
@@ -3976,14 +4124,22 @@ def do_add(self, line, counter=0):
                                "ignore_six_quark_processes" in self.options \
                                else []
     
-                myproc = diagram_generation.MultiProcess(myprocdef,
+                try:
+                    myproc = diagram_generation.MultiProcess(myprocdef,
                                          collect_mirror_procs = collect_mirror_procs,
                                          ignore_six_quark_processes = ignore_six_quark_processes,
                                          optimize=optimize, diagram_filter=diagram_filter,
                                          merge_crossing=merge_crossing)
+                    amplitudes = myproc.get('amplitudes')
+                except diagram_generation.NoDiagramException:
+                    if not myprocdef.get_interference_mode():
+                        raise
+                    raise self.InvalidCmd(
+                        "No subprocess of '%s' has diagrams on both sides of "
+                        "[treextree]: there is no interference to compute." %
+                        myprocdef.nice_string(prefix=False))
     
-    
-                for amp in myproc.get('amplitudes'):
+                for amp in amplitudes:
                     if amp not in self._curr_amps:
                         self._curr_amps.append(amp)
                     elif warning_duplicate:
@@ -5393,6 +5549,16 @@ def create_lambda_values_list(lower_bound, N):
         # Don't try to extract the process if just re-analyzing a saved run
         if not (args[0]=='cms' and options['analyze']!='None'):
             myprocdef = self.extract_process(proc_line)
+            if myprocdef and myprocdef.get_interference_mode():
+                raise self.InvalidCmd(
+                    "The check command does not support %s interference "
+                    "processes: generate the process and use 'output "
+                    "standalone_fortran', whose check executable prints the "
+                    "interference and, as cross-checks, the squared loop and "
+                    "tree amplitudes (complete only with HelicityFilterLevel "
+                    "0 in MadLoopParams.dat, the filter being set on the "
+                    "interference)."
+                    % myprocdef.get_interference_mode())
 
             # Check that we have something
             if not myprocdef:
@@ -5960,6 +6126,15 @@ def extract_process(self, line, proc_number = 0, overall_orders = {},
         """Extract a process definition from a string. Returns
         a ProcessDefinition."""
 
+        # 'left [LIxtree=QCD] right' / 'left [treextree] right': two processes
+        # on one line, see extract_interference_process
+        left, mode, right = split_interference_line(line)
+        if mode:
+            return self.extract_interference_process(left, mode, right,
+                                    proc_number=proc_number,
+                                    overall_orders=overall_orders,
+                                    avoid_squared_orders=avoid_squared_orders)
+
         orig_line = line
         # Check basic validity of the line
         if not len(re.findall(r'>\D', line)) in [1,2]:
@@ -6709,6 +6884,115 @@ def extract_process(self, line, proc_number = 0, overall_orders = {},
         #                       'is_decay_chain': decay_process\
 
 
+    def extract_interference_process(self, left, mode, right, proc_number=0,
+                                     overall_orders={},
+                                     avoid_squared_orders=False):
+        """Build the ProcessDefinition of an interference process, as split by
+        split_interference_line: 'left' is the left-hand process in the
+        ordinary syntax, 'mode' is 'LIxtree' or 'treextree' and 'right' the
+        right-hand (tree-level) process, None for the LIxtree default.
+
+        The left-hand process is returned, with 'interference_mode' set and
+        the right-hand ProcessDefinition in 'interference_process'. The
+        quantity computed is 2 Re(A_left A_right^*): the squared-order
+        constraints, which constrain that product, are given on the left-hand
+        side only. Without a right-hand process, LIxtree interferes with the
+        tree-level amplitudes of the same external legs, i.e. with what
+        'generate <legs>' alone would give: no s-channel requirement, no order
+        constraint and no exclusion is carried over from the loop process.
+        """
+
+        procdef = self.extract_process(left, proc_number=proc_number,
+                                       overall_orders=overall_orders,
+                                       avoid_squared_orders=avoid_squared_orders)
+        if not procdef:
+            raise self.InvalidCmd("Empty or wrong format process, please try again.")
+        if procdef.get('decay_chains'):
+            raise self.InvalidCmd("Decay chains are not supported for %s processes." % mode)
+
+        if right is None:
+            rightdef = base_objects.ProcessDefinition({
+                                        'legs': copy.deepcopy(procdef['legs']),
+                                        'model': procdef['model']})
+        else:
+            if ',' in right:
+                raise self.InvalidCmd("Decay chains are not supported for %s processes." % mode)
+            # the process number belongs to the whole definition, i.e. to
+            # the left-hand process: the right-hand one keeps id 0
+            rightdef = self.extract_process(right, avoid_squared_orders=True)
+            if not rightdef:
+                raise self.InvalidCmd("The right-hand process '%s' of the %s "
+                                      "process is not valid." % (right, mode))
+            if rightdef['squared_orders']:
+                raise self.InvalidCmd(
+                    "Squared-order constraints (%s) are not allowed on the "
+                    "right-hand process of a %s process: they constrain the "
+                    "interference, so they belong to the left-hand side, "
+                    "between ']' and the right-hand process, as in "
+                    "'g g > h > t t~ [LIxtree=QCD] QCD^2==6 g g > t t~'."
+                    % (', '.join('%s^2%s%d' % (k, rightdef.get_squared_order_type(k), v)
+                                 for k, v in rightdef['squared_orders'].items()),
+                       mode))
+
+        self.validate_interference_legs(procdef, rightdef, mode)
+
+        procdef.set('interference_mode', mode)
+        procdef.set('interference_process', rightdef)
+        return procdef
+
+    @staticmethod
+    def interference_leg_pdgs(leg, model):
+        """The set of signed pdg codes a MultiLeg stands for, with the merged
+        (flavour-grouped) particles expanded into their members, restricted
+        by the leg 'flavor' when it is set."""
+
+        merged = getattr(model, 'merged_particles', None) or {}
+        restriction = set(leg.get('flavor') or [])
+        pdgs = set()
+        for pid in leg.get('ids'):
+            if abs(pid) in merged:
+                sign = 1 if pid > 0 else -1
+                members = set(sign * m for m in merged[abs(pid)])
+                pdgs |= (members & restriction) or members
+            else:
+                pdgs.add(pid)
+        return pdgs
+
+    def validate_interference_legs(self, leftdef, rightdef, mode):
+        """The two sides of an interference process must describe the same
+        external legs, written in the same order: same number of initial and
+        final particles, the same polarisation, and at each position every
+        particle of the left-hand leg must be allowed by the right-hand one
+        (the right-hand legs may be wider, e.g. 'j' against 'g'). Raises
+        InvalidCmd otherwise."""
+
+        left_legs, right_legs = leftdef.get('legs'), rightdef.get('legs')
+        if leftdef.get_ninitial() != rightdef.get_ninitial() or \
+                                            len(left_legs) != len(right_legs):
+            raise self.InvalidCmd(
+                "The two sides of a %s process must have the same external "
+                "legs: '%s' and '%s' differ in the number of initial or final "
+                "particles." % (mode, leftdef.nice_string(prefix=False),
+                                rightdef.nice_string(prefix=False)))
+
+        model = leftdef.get('model')
+        names = lambda pdgs: ' '.join(sorted(model.get_particle(p).get_name()
+                                             for p in pdgs))
+        for i, (lleg, rleg) in enumerate(zip(left_legs, right_legs)):
+            missing = self.interference_leg_pdgs(lleg, model) - \
+                      self.interference_leg_pdgs(rleg, model)
+            if lleg.get('state') != rleg.get('state') or missing:
+                raise self.InvalidCmd(
+                    "The two sides of a %s process must have the same external "
+                    "legs, in the same order: leg %d of the right-hand process "
+                    "does not contain %s. Write the right-hand legs in the "
+                    "order of the left-hand ones." % (mode, i + 1,
+                    names(missing) if missing else 'the same state'))
+            if sorted(lleg.get('polarization')) != sorted(rleg.get('polarization')):
+                raise self.InvalidCmd(
+                    "The two sides of a %s process must have the same "
+                    "polarisation on each leg (leg %d differs)." % (mode, i + 1))
+
     def create_loop_induced(self, line, myprocdef=None):
         """ Routine to create the MultiProcess for the loop-induced case"""
         
@@ -6743,6 +7027,8 @@ def create_loop_induced(self, line, myprocdef=None):
         
         if not myprocdef:
             myprocdef = self.extract_process(' '.join(args))
+        # LIxtree: hidden order selecting the interference, right-hand orders
+        diagram_generation.prepare_interference_process(myprocdef)
         
         myprocdef.set('NLO_mode', 'noborn')
             
@@ -6786,13 +7072,23 @@ def create_loop_induced(self, line, myprocdef=None):
 
         # Decide here wether one needs a LoopMultiProcess or a MultiProcess
 
-        myproc = loop_diagram_generation.LoopInducedMultiProcess(myprocdef,
+        try:
+            myproc = loop_diagram_generation.LoopInducedMultiProcess(myprocdef,
                                  collect_mirror_procs = collect_mirror_procs,
                                  ignore_six_quark_processes = ignore_six_quark_processes,
                                  optimize=optimize,
                                  loop_filter=loop_filter)
+            amplitudes = myproc.get('amplitudes')
+        except diagram_generation.NoDiagramException:
+            if myprocdef.get_interference_mode() != 'LIxtree':
+                raise
+            raise self.InvalidCmd(
+                "No subprocess of '%s' has both loop-induced diagrams and "
+                "tree-level diagrams of '%s': there is no interference to "
+                "compute." % (myprocdef.nice_string(prefix=False).split('[')[0].strip(),
+                    myprocdef.get('interference_process').nice_string(prefix=False)))
 
-        for amp in myproc.get('amplitudes'):
+        for amp in amplitudes:
             if amp not in self._curr_amps:
                 self._curr_amps.append(amp)
                 if amp['has_born']:
@@ -12089,6 +12385,24 @@ def do_output(self, line):
                         export_v4.loop_induced_not_supported_msg(
                             format, self._curr_amps[0].get('process')))
 
+        # A tree x tree interference ('left [treextree] right') needs the
+        # selection of squared split orders in the generated code: other
+        # formats would silently return |A_left + A_right|^2.
+        # (a decay-chain amplitude has no 'process', and is never one)
+        interference = [amp for amp in self._curr_amps if 'process' in amp
+            and amp['process'].get_interference_mode() == 'treextree']
+        if interference:
+            for format in [self._export_format,
+                           options['me_exporter'].get('name')]:
+                if format and format not in \
+                           diagram_generation.TREE_INTERFERENCE_FORMATS:
+                    raise self.InvalidCmd(
+                        "The '%s' output format cannot select the "
+                        "interference of a treextree process (it would "
+                        "return |A_left + A_right|^2): use one of %s."
+                        % (format, ', '.join(
+                               diagram_generation.TREE_INTERFERENCE_FORMATS)))
+
         # check
         if os.path.realpath(self._export_dir) == os.getcwd():
             if len(args) == 0:
diff --git a/madgraph/interface/master_interface.py b/madgraph/interface/master_interface.py
index 0bfc1904cb..d145a35a6b 100755
--- a/madgraph/interface/master_interface.py
+++ b/madgraph/interface/master_interface.py
@@ -181,6 +181,14 @@ def extract_process_type(line):
         res=loopRE.search(re.split(r'\s\-\-', line,1)[0])
         if res:
             orders=res.group('orders').split() if res.group('orders') else []
+            # interference of two amplitude sets (keyword case-insensitive):
+            # 'left [LIxtree=QCD] right' is a loop-induced process, routed
+            # like [noborn=], 'left [treextree] right' a tree-level one.
+            interference = MGcmd.interference_mode_of(res.group('option'), orders)
+            if interference == 'LIxtree':
+                return ('NLO', 'LIxtree', orders if res.group('option') else [])
+            elif interference == 'treextree':
+                return ('tree', 'treextree', [])
             if res.group('option') and len(res.group('option').split())==1:
                 if res.group('option').split()[0]=='tree':
                     return ('tree',res.group('option').split()[0],orders)
@@ -219,11 +227,15 @@ def do_add(self, line, *args, **opts):
                     self.change_principal_cmd('aMC@NLO', allow_switch)
                 elif nlo_mode in ['virt', 'sqrvirt']:
                     self.change_principal_cmd('MadLoop', allow_switch)
-                elif nlo_mode == 'noborn': 
+                elif nlo_mode in ['noborn', 'LIxtree']:
+                    # LIxtree: loop-induced amplitude interfered with a tree.
+                    # Report a malformed line before the model is touched.
+                    if nlo_mode == 'LIxtree':
+                        MGcmd.split_interference_line(proc_line)
                     if self.current_interface == "MadGraph":
                         allow_switch = True
                     self.change_principal_cmd('MadLoop', allow_switch)
-                    self.cmd.validate_model(self, loop_type=nlo_mode,
+                    self.cmd.validate_model(self, loop_type='noborn',
                                                             coupling_type=orders)
                     self.change_principal_cmd('MadGraph', allow_switch)
                     return self.cmd.create_loop_induced(self, line, *args, **opts)
diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py
index b84ff0c6b3..2f75146242 100755
--- a/madgraph/iolibs/export_v4.py
+++ b/madgraph/iolibs/export_v4.py
@@ -7812,6 +7812,7 @@ def generate_subprocess_directory(self, matrix_element,
             matrix_element)
 
         filename = pjoin(Ppath, 'config_nqcd.inc')
+        nqcd_list = interference_nqcd_list(nqcd_list, [matrix_element])
         self.write_config_nqcd_file(writers.FortranWriter(filename),
                                nqcd_list)
 
@@ -10032,6 +10033,7 @@ def generate_subprocess_directory(self, subproc_group,
             subproc_diagrams_for_config)
 
         filename = 'config_nqcd.inc'
+        nqcd_list = interference_nqcd_list(nqcd_list, matrix_elements)
         self.write_config_nqcd_file(writers.FortranWriter(filename),
                                     nqcd_list)
 
@@ -13826,6 +13828,83 @@ def create_param_card(self, write_special=True):
 # standalone_msP / _msF / _rw. None of those has a loop backend.
 LOOP_INDUCED_FORMATS = ['madevent', 'plugin', 'standalone_fortran']
 
+def interference_alpha_s_power(matrix_elements):
+    """The power of alpha_s of the events of an interference process
+    ('left [LIxtree=QCD] right' or 'left [treextree] right').
+
+    madevent gives each event the QCD order of the diagram of its integration
+    channel (config_nqcd.inc), which is the power of alpha_s of |A|^2 for an
+    ordinary process. For an interference the quantity computed is
+    2 Re(A_left A_right^*): its power is half the QCD order of the products of
+    one left-hand and one right-hand amplitude passing the squared-order
+    constraints, whatever the channel. Returns:
+      False  for an ordinary process (keep the per-channel values),
+      n      when all those products share the power n,
+      None   when they do not, or when it is not an integer (the scale
+             variations of the systematics are then ill-defined).
+    """
+
+    import madgraph.loop.loop_helas_objects as loop_helas_objects
+    powers = set()
+    for me in matrix_elements:
+        process = me.get('processes')[0]
+        if not process.get_interference_mode():
+            return False
+        hierarchy = process.get('model').get('order_hierarchy')
+        if isinstance(me, loop_helas_objects.LoopHelasMatrixElement):
+            amplitude = me.get('base_amplitude')
+            orders = [d.get('orders') for d in amplitude.get('loop_diagrams') +
+                      amplitude.get('loop_UVCT_diagrams')]
+        else:
+            orders = [d.calculate_orders() for d in me.get('diagrams')]
+        sides = {0: [], 1: []}
+        for signature in set(tuple(sorted((k, v) for k, v in o.items()
+                                          if k != 'WEIGHTED')) for o in orders):
+            side = dict(signature).get(base_objects.INTERFERENCE_ORDER, 0)
+            if side in sides:
+                sides[side].append(dict(signature))
+        squared_orders = process.get('squared_orders')
+        for left in sides[0]:
+            for right in sides[1]:
+                product = dict((k, left.get(k, 0) + right.get(k, 0))
+                               for k in set(left) | set(right))
+                product['WEIGHTED'] = sum(hierarchy.get(k, 0) * v
+                                          for k, v in product.items())
+                passed = True
+                for order, value in squared_orders.items():
+                    if value < 0:
+                        continue
+                    kind = process.get_squared_order_type(order)
+                    combined = product.get(order, 0)
+                    if (kind == '==' and combined != value) or \
+                       (kind in ['=', '<='] and combined > value) or \
+                       (kind == '>' and combined <= value):
+                        passed = False
+                        break
+                if passed:
+                    powers.add(product.get('QCD', 0))
+    if len(powers) == 1:
+        power = powers.pop()
+        if power % 2 == 0:
+            return power // 2
+    return None
+
+def interference_nqcd_list(nqcd_list, matrix_elements):
+    """config_nqcd.inc entries of an interference process: the same power of
+    alpha_s for every channel (see interference_alpha_s_power)."""
+
+    power = interference_alpha_s_power(matrix_elements)
+    if power is False:
+        return nqcd_list
+    if power is None:
+        logger.warning("The power of alpha_s of this interference is not "
+                       "unique: the scale variations of the systematics "
+                       "(off by default for interference processes) would be "
+                       "wrong. Constrain the QCD order of the interference "
+                       "(e.g. QCD^2==n after ']') to make it unique.")
+        return nqcd_list
+    return [power] * len(nqcd_list)
+
 def loop_induced_not_supported_msg(format, process=None):
     """Refusal text for a format that cannot write a LoopHelasMatrixElement."""
 
diff --git a/madgraph/loop/loop_base_objects.py b/madgraph/loop/loop_base_objects.py
index 17d3658341..e0bd17c1d0 100755
--- a/madgraph/loop/loop_base_objects.py
+++ b/madgraph/loop/loop_base_objects.py
@@ -1366,6 +1366,16 @@ def get_sorted_keys(self):
         
         return ['vertices', 'UVCT_couplings', 'UVCT_orders', 'type', 'orders']
 
+    def get_contracted_loop_diagram(self, model, struct_rep=None):
+        """A UVCT diagram has no loop to contract: return a copy of it, as
+        for the tree-level diagrams interfered with a loop-induced amplitude
+        (LIxtree), which are stored as UVCT diagrams."""
+        return copy.copy(self)
+
+    def get_contracted_loop_diagram_without_tag(self, struct_rep=None):
+        """See get_contracted_loop_diagram."""
+        return copy.copy(self)
+
     def get_UVCTinteraction(self, model):
         """ Finds the UV counter-term interaction present in this UVCTDiagram """
         
diff --git a/madgraph/loop/loop_diagram_generation.py b/madgraph/loop/loop_diagram_generation.py
index 6fb8d51eae..e7f0cbf3a8 100755
--- a/madgraph/loop/loop_diagram_generation.py
+++ b/madgraph/loop/loop_diagram_generation.py
@@ -30,6 +30,14 @@
 from madgraph import InvalidCmd
 logger = logging.getLogger('madgraph.loop_diagram_generation')
 
+#===============================================================================
+# Interference of a loop-induced amplitude with tree-level ones (LIxtree): see
+# the helpers in diagram_generation, shared with the tree x tree case
+#===============================================================================
+INTERFERENCE_ORDER = base_objects.INTERFERENCE_ORDER
+register_interference_order = diagram_generation.register_interference_order
+prepare_interference_process = diagram_generation.prepare_interference_process
+
 def ldg_debug_info(msg,val, force=False):
     # This subroutine has typically quite large DEBUG info.
     # So even in debug mode, they are turned off by default.
@@ -758,6 +766,11 @@ def generate_diagrams(self, loop_filter=None, diagram_filter=None):
             self['process']['orders'].update(user_orders)
             return False
 
+        # 'left [LIxtree=QCD] right': the tree-level amplitudes of the
+        # right-hand process, interfered with the loop-induced ones
+        if self['process'].get_interference_mode() == 'LIxtree':
+            self.set_interference_trees()
+
         # We add here the UV renormalization contribution built in
         # LoopUVCTDiagram. It is done before the squared order selection because
         # it is possible that some UV-renorm. diagrams are removed as well.
@@ -865,6 +878,9 @@ def generate_diagrams(self, loop_filter=None, diagram_filter=None):
         if not self['process']['has_born'] and not self['loop_diagrams']:
             self['process']['squared_orders'].clear()
             self['process']['squared_orders'].update(user_squared_orders)
+            # nothing left to interfere with the trees of a LIxtree process
+            # (they would otherwise keep this amplitude alive)
+            self['loop_UVCT_diagrams'] = base_objects.DiagramList()
             return False
 
 
@@ -900,7 +916,11 @@ def generate_diagrams(self, loop_filter=None, diagram_filter=None):
         # Give some info about the run
         nLoopDiag = 0
         nCT={'UV':0,'R2':0}
+        nInterferenceTrees = 0
         for ldiag in self['loop_UVCT_diagrams']:
+            if ldiag['UVCT_orders'].get(INTERFERENCE_ORDER):
+                nInterferenceTrees += 1
+                continue
             nCT[ldiag['type'][:2]]+=len(ldiag['UVCT_couplings'])
         for ldiag in self['loop_diagrams']:
             nLoopDiag+=1
@@ -917,7 +937,9 @@ def generate_diagrams(self, loop_filter=None, diagram_filter=None):
         logger.info("Contributing diagrams generated: "+\
           "%d Born, %d%s loops, %d R2, %d UV"%(len(self['born_diagrams']),
                     len(self['loop_diagrams']),'(+%d)'%nLoopsIdentified \
-                            if nLoopsIdentified>0 else '' ,nCT['R2'],nCT['UV']))
+                            if nLoopsIdentified>0 else '' ,nCT['R2'],nCT['UV'])+\
+          (", %d interfering tree%s" % (nInterferenceTrees,
+            's' if nInterferenceTrees > 1 else '') if nInterferenceTrees else ''))
         
         ldg_debug_info("#Diags after filtering",len(self['loop_diagrams']))
         ldg_debug_info("# of different structures identified",\
@@ -1204,6 +1226,46 @@ def generate_loop_diagrams(self):
         return totloopsuccessful
 
 
+    def set_interference_trees(self):
+        """For an interference process 'left [LIxtree=QCD] right', generate
+        the tree-level diagrams of the right-hand process with the external
+        legs of this (loop-induced) process, and store them among the
+        'loop_UVCT_diagrams': like a UV counterterm, a tree diagram is a
+        loop-less amplitude that the loop amplitudes are squared against.
+        Each one carries one unit of the hidden order INTERFERENCE_ORDER in
+        its 'UVCT_orders' and a unit 'UVCT_couplings' factor, so that the
+        squared-order constraint INTERF^2==1 (see
+        prepare_interference_process) keeps 2 Re(A_loop A_tree^*) only.
+        Returns the number of tree diagrams added (0 if the right-hand
+        process has none for these legs, in which case the squared-order
+        selection then discards this subprocess)."""
+
+        process = self['process']
+        model = process.get('model')
+        register_interference_order(model)
+        tree_process = diagram_generation.specialise_interference_process(process)
+
+        try:
+            tree_amplitude = diagram_generation.Amplitude(tree_process)
+            tree_diagrams = tree_amplitude.get('diagrams')
+        except (InvalidCmd, diagram_generation.NoDiagramException):
+            tree_diagrams = []
+
+        for diagram in tree_diagrams:
+            tree = loop_base_objects.LoopUVCTDiagram({
+                            'vertices': copy.deepcopy(diagram.get('vertices'))})
+            tree.set('type', 'UVtree')
+            tree.set('UVCT_couplings', [1])
+            tree.set('UVCT_orders', {INTERFERENCE_ORDER: 1})
+            tree.calculate_orders(model)
+            self['loop_UVCT_diagrams'].append(tree)
+
+        if not tree_diagrams:
+            logger.info("No tree-level diagram of %s: no interference."
+                        % tree_process.nice_string(prefix=False))
+        ldg_debug_info("#Interfering tree diagrams", len(tree_diagrams))
+        return len(tree_diagrams)
+
     def set_Born_CT(self):
         """ Scan all born diagrams and add for each all the corresponding UV 
         counterterms. It creates one LoopUVCTDiagram per born diagram and set
diff --git a/madgraph/loop/loop_helas_objects.py b/madgraph/loop/loop_helas_objects.py
index f01eaeebae..1d599d4e4e 100755
--- a/madgraph/loop/loop_helas_objects.py
+++ b/madgraph/loop/loop_helas_objects.py
@@ -1918,7 +1918,13 @@ def get_split_orders_mapping(self):
         if process.get('has_born'):
             ref_orders = [bao[0] for bao in born_amp_orders]
         else:
-            ref_orders = [lao[0] for lao in loop_orders+ct_amp_orders]
+            # Without Born, the loop amplitudes are squared against all the
+            # loop-less ones too: the tree-level amplitudes of an interference
+            # process (LIxtree) are stored as UVCT amplitudes. Their squared
+            # orders must be known even if only the interference is selected,
+            # since the generated code evaluates every product.
+            ref_orders = [lao[0] for lao in loop_orders+ct_amp_orders+
+                                                               uvct_amp_orders]
         
         # Temporarily we set squared_orders to be a dictionary with keys being
         # the actual contributing squared_orders and the values are the list 
diff --git a/madgraph/various/banner.py b/madgraph/various/banner.py
index 39e0f02520..878f235a55 100755
--- a/madgraph/various/banner.py
+++ b/madgraph/various/banner.py
@@ -5269,7 +5269,12 @@ def create_default_for_process(self, proc_characteristic, history, proc_def):
         interference = False
         for proc in proc_def:
             for oneproc in proc:
-                if '^2' in oneproc.nice_string():
+                # a squared-order constraint, or an interference process
+                # ('left [LIxtree=QCD] right' / '[treextree]'), whose hidden
+                # squared-order constraint is not printed
+                if '^2' in oneproc.nice_string() or \
+                        (hasattr(oneproc, 'get_interference_mode') and
+                                            oneproc.get_interference_mode()):
                     interference = True
                     break
             else:
diff --git a/madgraph/various/systematics.py b/madgraph/various/systematics.py
index 71c9b317b3..99f1845eb4 100644
--- a/madgraph/various/systematics.py
+++ b/madgraph/various/systematics.py
@@ -490,9 +490,14 @@ def print_cross_sections(self, all_cross, nb_event, stdout):
         #print "normalisation is ", norm
         #print "nb_event is ", nb_event
     
-        max_scale, min_scale = 0,sys.maxsize
-        max_alps, min_alps = 0, sys.maxsize
-        max_dyn, min_dyn = 0,sys.maxsize
+        # Bounds of each envelope, None until a variation is met: a
+        # cross-section can be negative (an interference, for instance), so
+        # neither 0 nor sys.maxsize is a neutral starting value.
+        max_scale, min_scale = None, None
+        max_alps, min_alps = None, None
+        max_dyn, min_dyn = None, None
+        upper = lambda bound, value: value if bound is None else max(bound, value)
+        lower = lambda bound, value: value if bound is None else min(bound, value)
         pdfs = {}
         dyns = {} # dyn : {'max': , 'min':}
 
@@ -520,28 +525,28 @@ def print_cross_sections(self, all_cross, nb_event, stdout):
             mur, muf, alps, dyn, pdf = arg[:5]
             if pdf == self.orig_pdf and (dyn==self.orig_dyn or dyn==-1)\
                                               and (mur!=1 or muf!=1 or alps!=1):
-                max_scale = max(max_scale,all_cross[i])
-                min_scale = min(min_scale,all_cross[i])
+                max_scale = upper(max_scale,all_cross[i])
+                min_scale = lower(min_scale,all_cross[i])
             if pdf == self.orig_pdf and mur==1 and muf==1 and \
                                     (dyn==self.orig_dyn or dyn==-1) and alps!=1:
-                max_alps = max(max_alps,all_cross[i])
-                min_alps = min(min_alps,all_cross[i]) 
+                max_alps = upper(max_alps,all_cross[i])
+                min_alps = lower(min_alps,all_cross[i])
             
             if pdf == self.orig_pdf and mur==1 and muf==1 and alps==1:
-                max_dyn = max(max_dyn,all_cross[i])
-                min_dyn = min(min_dyn,all_cross[i])
+                max_dyn = upper(max_dyn,all_cross[i])
+                min_dyn = lower(min_dyn,all_cross[i])
                                             
             if pdf == self.orig_pdf and (alps!=1 or mur!=1 or muf!=1) and \
                                                 (dyn!=self.orig_dyn or dyn!=-1):
                 if dyn not in dyns:
-                    dyns[dyn] = {'max':0, 'min':sys.maxsize,'central':0}
+                    dyns[dyn] = {'max':None, 'min':None,'central':0}
                 curr = dyns[dyn]
-                curr['max'] = max(curr['max'],all_cross[i])
-                curr['min'] = min(curr['min'],all_cross[i])
+                curr['max'] = upper(curr['max'],all_cross[i])
+                curr['min'] = lower(curr['min'],all_cross[i])
             if pdf == self.orig_pdf and (alps==1 and mur==1 and muf==1) and \
                                                 (dyn!=self.orig_dyn or dyn!=-1):
                 if dyn not in dyns:
-                    dyns[dyn] = {'max':0, 'min':sys.maxsize,'central':all_cross[i]}
+                    dyns[dyn] = {'max':None, 'min':None,'central':all_cross[i]}
                 else:
                     dyns[dyn]['central'] = all_cross[i]          
                 
@@ -561,12 +566,16 @@ def print_cross_sections(self, all_cross, nb_event, stdout):
         resume.write( '#***************************************************************************\n')
         resume.write( "#\n")
         resume.write( '# original cross-section: %s\n' % all_cross[0])
-        if max_scale:
-            resume.write( '#     scale variation: +%2.3g%% -%2.3g%%\n' % ((max_scale-all_cross[0])/all_cross[0]*100,(all_cross[0]-min_scale)/all_cross[0]*100))
-        if max_alps:
-            resume.write( '#     emission scale variation: +%2.3g%% -%2.3g%%\n' % ((max_alps-all_cross[0])/all_cross[0]*100,(all_cross[0]-min_alps)/all_cross[0]*100))
-        if max_dyn and (max_dyn!= all_cross[0] or min_dyn != all_cross[0]):
-            resume.write( '#     central scheme variation: +%2.3g%% -%2.3g%%\n' % ((max_dyn-all_cross[0])/all_cross[0]*100,(all_cross[0]-min_dyn)/all_cross[0]*100))
+        # variations in percent of the size of the central value: '+' is the
+        # increase of the cross-section, also when it is negative
+        central_size = abs(all_cross[0]) if all_cross[0] else 1.
+        percent = lambda value: value / central_size * 100
+        if max_scale is not None:
+            resume.write( '#     scale variation: +%2.3g%% -%2.3g%%\n' % (percent(max_scale-all_cross[0]),percent(all_cross[0]-min_scale)))
+        if max_alps is not None:
+            resume.write( '#     emission scale variation: +%2.3g%% -%2.3g%%\n' % (percent(max_alps-all_cross[0]),percent(all_cross[0]-min_alps)))
+        if max_dyn is not None and (max_dyn!= all_cross[0] or min_dyn != all_cross[0]):
+            resume.write( '#     central scheme variation: +%2.3g%% -%2.3g%%\n' % (percent(max_dyn-all_cross[0]),percent(all_cross[0]-min_dyn)))
         if self.banner.run_card['pdlabel'] in ['eva']:
             resume.write( '# PDF variation not available for EVA.\n')
         elif self.orig_pdf.lhapdfID in pdfs:
@@ -578,7 +587,7 @@ def print_cross_sections(self, all_cross, nb_event, stdout):
             except RuntimeError:
                 resume.write( '# PDF variation: missing combination\n')
             else:
-                resume.write( '# PDF variation: +%2.3g%% -%2.3g%%\n' % (pdferr.errplus*100/all_cross[0], pdferr.errminus*100/all_cross[0]))       
+                resume.write( '# PDF variation: +%2.3g%% -%2.3g%%\n' % (percent(pdferr.errplus), percent(pdferr.errminus)))       
         # report error/central not directly linked to the central
         resume.write( "#\n")        
         for lhapdfid,values in pdfs.items():
@@ -599,7 +608,7 @@ def print_cross_sections(self, all_cross, nb_event, stdout):
                 # the same error can happend to some other type of error like custom.
                 pass
             else:
-                resume.write( '#PDF %s: %g +%2.3g%% -%2.3g%%\n' % (pdfset.name, pdferr.central,pdferr.errplus*100/all_cross[0], pdferr.errminus*100/all_cross[0]))
+                resume.write( '#PDF %s: %g +%2.3g%% -%2.3g%%\n' % (pdfset.name, pdferr.central,percent(pdferr.errplus), percent(pdferr.errminus)))
 
         dyn_name = {1: r'\sum ET', 2:r'\sum\sqrt{m^2+pt^2}', 3:r'0.5 \sum\sqrt{m^2+pt^2}',4:r'\sqrt{\hat s}' }
         for key, curr in dyns.items():
@@ -608,10 +617,10 @@ def print_cross_sections(self, all_cross, nb_event, stdout):
             central, maxvalue, minvalue = curr['central'], curr['max'], curr['min']
             if central == 0:
                 continue
-            if maxvalue == 0:
+            if maxvalue is None:
                 resume.write("# dynamical scheme # %s : %g # %s\n" %(key, central, dyn_name[key]))
             else:
-                resume.write("# dynamical scheme # %s : %g +%2.3g%% -%2.3g%% # %s\n" %(key, central, (maxvalue-central)/central*100,(central-minvalue)/central*100, dyn_name[key]))
+                resume.write("# dynamical scheme # %s : %g +%2.3g%% -%2.3g%% # %s\n" %(key, central, (maxvalue-central)/abs(central)*100,(central-minvalue)/abs(central)*100, dyn_name[key]))
       
         resume.write('\n'.join(to_report))
 
diff --git a/tests/acceptance_tests/test_loop_induced_x_tree.py b/tests/acceptance_tests/test_loop_induced_x_tree.py
new file mode 100644
index 0000000000..c7ca8112c4
--- /dev/null
+++ b/tests/acceptance_tests/test_loop_induced_x_tree.py
@@ -0,0 +1,166 @@
+################################################################################
+#
+# Copyright (c) 2026 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+"""Numerical checks of the interference processes. First the loop-induced
+x tree interference
+
+    generate g g > h > t t~ [LIxtree=QCD]
+    output standalone_fortran
+
+at one fixed phase-space point (sqrt(s) = 1 TeV, default loop_sm parameters).
+The standalone MadLoop check prints the three products it can form:
+
+  INTERF=1  2 Re(A_loop A_tree^*)   the selected result
+  INTERF=0  |A_loop|^2              must equal 'g g > h > t t~ [sqrvirt=QCD]'
+  INTERF=2  |A_tree|^2              must equal the tree-level matrix element
+                                    of 'g g > t t~' (with the helicity filter
+                                    off: it is set on the interference)
+
+The references were obtained independently: the tree-level matrix element
+from 'output standalone_fortran' of 'g g > t t~', the loop squared from
+'[sqrvirt=QCD]', and the interference from the original recipe of the wiki
+page LoopInducedTimesTree (tree vertices copied as UVtree counterterms with a
+'BKGQCD' order), all at this same point.
+"""
+
+from __future__ import absolute_import
+
+import os
+import re
+import shutil
+import subprocess
+import tempfile
+
+import tests.unit_tests as unittest
+
+import madgraph.interface.master_interface as MGCmd
+
+pjoin = os.path.join
+
+PS_POINT = """5.0000000000000000e+02  0.0000000000000000e+00  0.0000000000000000e+00  5.0000000000000000e+02
+5.0000000000000000e+02  0.0000000000000000e+00  0.0000000000000000e+00 -5.0000000000000000e+02
+4.9999999999999989e+02  1.0407299191319960e+02  4.1735559881462342e+02 -1.8722744588420289e+02
+4.9999999999999989e+02 -1.0407299191319960e+02 -4.1735559881462342e+02  1.8722744588420281e+02
+"""
+
+REFERENCES = {1: -2.3256785052943376e-04,   # wiki recipe
+              0: 8.1044796923476000e-05,    # [sqrvirt=QCD]
+              2: 5.9262610090352252e-01}    # tree-level g g > t t~
+
+
+class LoopInducedXTreeStandaloneTest(unittest.TestCase):
+
+    def setUp(self):
+        self.tmpdir = tempfile.mkdtemp(prefix='lixtree_')
+
+    def tearDown(self):
+        shutil.rmtree(self.tmpdir, ignore_errors=True)
+
+    def test_loop_induced_x_tree_standalone(self):
+        out_dir = pjoin(self.tmpdir, 'gg_h_ttx')
+        interface = MGCmd.MasterCmd()
+        for command in ['import model loop_sm',
+                        'generate g g > h > t t~ [LIxtree=QCD]',
+                        'output standalone_fortran %s -f' % out_dir]:
+            interface.exec_cmd(command, printcmd=False, precmd=True,
+                               errorhandling=False)
+
+        proc_dir = pjoin(out_dir, 'SubProcesses', 'P0_gg_h_ttx')
+        # read the point above instead of the RAMBO one
+        check_sa = pjoin(proc_dir, 'check_sa.f')
+        text = open(check_sa).read()
+        self.assertIn('PARAMETER (READPS = .FALSE.)', text)
+        open(check_sa, 'w').write(text.replace('PARAMETER (READPS = .FALSE.)',
+                                               'PARAMETER (READPS = .TRUE.)'))
+        open(pjoin(proc_dir, 'PS.input'), 'w').write(PS_POINT)
+        # no helicity filter, so that |A_tree|^2 is complete as well
+        card = pjoin(proc_dir, 'MadLoop5_resources', 'MadLoopParams.dat')
+        lines = open(card).read().split('\n')
+        index = lines.index('#HelicityFilterLevel')
+        lines[index + 1] = '0'
+        open(card, 'w').write('\n'.join(lines))
+
+        subprocess.check_call(['make', 'check'], cwd=proc_dir,
+                              stdout=subprocess.DEVNULL,
+                              stderr=subprocess.DEVNULL)
+        output = subprocess.run(['./check'], cwd=proc_dir, timeout=600,
+                                capture_output=True, text=True).stdout
+
+        results = {}
+        for match in re.finditer(r'Loop ME for orders \(INTERF=(\d) QCD=4\) :\s*'
+                                 r'> accuracy\s*=\s*\S+\s*> finite\s*=\s*(\S+)',
+                                 output):
+            results[int(match.group(1))] = float(match.group(2).replace('D', 'E'))
+        self.assertEqual(sorted(results), [0, 1, 2], output[-2000:])
+        for interf, reference in REFERENCES.items():
+            self.assertAlmostEqual(results[interf] / reference, 1., places=10,
+                                   msg='INTERF=%d: %r instead of %r' %
+                                       (interf, results[interf], reference))
+
+        # the result returned is the interference alone
+        total = re.search(r'Matrix element finite\s*=\s*(\S+)', output)
+        self.assertAlmostEqual(float(total.group(1)) / REFERENCES[1], 1.,
+                               places=10)
+
+
+
+class TreeXTreeStandaloneTest(unittest.TestCase):
+    """'u u~ > z > e+ e- [treextree] u u~ > a > e+ e-' at the phase-space
+    point of the standalone check, against |A_Z + A_a|^2 - |A_Z|^2 - |A_a|^2
+    from three ordinary standalone outputs evaluated at the same point."""
+
+    def setUp(self):
+        self.tmpdir = tempfile.mkdtemp(prefix='treextree_')
+
+    def tearDown(self):
+        shutil.rmtree(self.tmpdir, ignore_errors=True)
+
+    def evaluate(self, name, process):
+        """(momenta, matrix element) printed by the standalone check."""
+        out_dir = pjoin(self.tmpdir, name)
+        interface = MGCmd.MasterCmd()
+        for command in ['import model sm', 'generate %s' % process,
+                        'output standalone_fortran %s -f' % out_dir]:
+            interface.exec_cmd(command, printcmd=False, precmd=True,
+                               errorhandling=False)
+        sub_dir = pjoin(out_dir, 'SubProcesses')
+        proc_dir = [pjoin(sub_dir, d) for d in os.listdir(sub_dir)
+                    if d.startswith('P')][0]
+        subprocess.check_call(['make', 'check'], cwd=proc_dir,
+                              stdout=subprocess.DEVNULL,
+                              stderr=subprocess.DEVNULL)
+        output = subprocess.run(['./check'], cwd=proc_dir, timeout=300,
+                                capture_output=True, text=True).stdout
+        momenta = re.findall(r'^\s*[1-4]\s+(\S+\s+\S+\s+\S+\s+\S+)', output,
+                             re.MULTILINE)
+        value = re.search(r'Matrix element =\s*(\S+)', output).group(1)
+        return momenta, float(value.replace('D', 'E'))
+
+    def test_tree_x_tree_standalone(self):
+        point, interference = self.evaluate('interf',
+                        'u u~ > z > e+ e- [treextree] u u~ > a > e+ e-')
+        self.assertEqual(len(point), 4)
+        references = {}
+        for name, process in [('za', 'u u~ > e+ e- / h'),
+                              ('z', 'u u~ > z > e+ e-'),
+                              ('a', 'u u~ > a > e+ e-')]:
+            momenta, references[name] = self.evaluate(name, process)
+            self.assertEqual(momenta, point)
+        expected = references['za'] - references['z'] - references['a']
+        self.assertTrue(abs(expected) > 1e-3 * references['za'])
+        self.assertAlmostEqual(interference / expected, 1., places=10)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/unit_tests/core/test_base_objects.py b/tests/unit_tests/core/test_base_objects.py
index d8c1d072d2..c63393c59c 100755
--- a/tests/unit_tests/core/test_base_objects.py
+++ b/tests/unit_tests/core/test_base_objects.py
@@ -1739,7 +1739,9 @@ def setUp(self):
                        'overall_orders': {},
                        'born_sq_orders': {},
                        'NLO_mode':'tree',
-                       'split_orders':[]}
+                       'split_orders':[],
+                       'interference_mode': '',
+                       'interference_process': None}
 
         self.myprocess = base_objects.Process(self.mydict)
 
@@ -1827,7 +1829,9 @@ def test_representation(self):
         goal = goal + "    \'has_born\': True,\n"
         goal = goal + "    \'NLO_mode\': 'tree',\n"
         goal = goal + "    \'split_orders\': [],\n"
-        goal = goal + "    \'born_sq_orders\': {}\n}"
+        goal = goal + "    \'born_sq_orders\': {},\n"
+        goal = goal + "    \'interference_mode\': '',\n"
+        goal = goal + "    \'interference_process\': None\n}"
 
         for a, b in zip(goal.split('\n'), str(self.myprocess).split('\n')):
             self.assertEqual(a,b)
@@ -2160,7 +2164,9 @@ def setUp(self):
                        'overall_orders':{},
                        'sqorders_types':{},
                        'NLO_mode':'tree',
-                       'split_orders':[]}
+                       'split_orders':[],
+                       'interference_mode': '',
+                       'interference_process': None}
 
         self.my_process_definition = base_objects.ProcessDefinition(self.mydict)
 
@@ -2226,7 +2232,11 @@ def test_get_process_with_legs(self):
         myleglist = base_objects.LegList(mylist)
         my_new_process_definition = copy.copy(self.my_process_definition)
         my_new_process_definition['born_sq_orders'] = {'QCD':99, 'QED':99}
+        my_new_process_definition['interference_mode'] = 'treextree'
+        my_new_process_definition['interference_process'] = \
+                                       base_objects.ProcessDefinition({'id': 0})
         testproc = my_new_process_definition.get_process_with_legs(myleglist)
+        self.assertEqual(testproc['interference_mode'], 'treextree')
 
         for (k, v) in testproc.items():
             if k not in list(self.my_process_definition.keys()): continue
@@ -2403,7 +2413,9 @@ def test_representation(self):
         goal = goal + "    \'has_born\': True,\n"
         goal = goal + "    \'NLO_mode\': 'tree',\n"
         goal = goal + "    \'split_orders\': [],\n"                
-        goal = goal + "    \'born_sq_orders\': {}\n}"                
+        goal = goal + "    \'born_sq_orders\': {},\n"
+        goal = goal + "    \'interference_mode\': '',\n"
+        goal = goal + "    \'interference_process\': None\n}"                
         self.assertEqual(goal, str(self.my_process_definition))
 
 #===============================================================================
diff --git a/tests/unit_tests/core/test_tree_interference.py b/tests/unit_tests/core/test_tree_interference.py
new file mode 100644
index 0000000000..0ddaec4ee1
--- /dev/null
+++ b/tests/unit_tests/core/test_tree_interference.py
@@ -0,0 +1,164 @@
+################################################################################
+#
+# Copyright (c) 2026 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+"""Generation of the interference of two tree-level amplitude sets:
+
+    generate u u~ > z > e+ e- [treextree] u u~ > a > e+ e-
+
+The diagrams of the right-hand process, generated with the legs of each
+subprocess, are added to those of the left-hand one with one unit of the
+hidden order base_objects.INTERFERENCE_ORDER, and the squared-order constraint
+INTERF^2==1 keeps 2 Re(A_left A_right^*) only. The numerical check (against
+|A_Z + A_a|^2 - |A_Z|^2 - |A_a|^2) is the acceptance test
+test_tree_x_tree_standalone in tests/acceptance_tests/test_loop_induced_x_tree.py.
+"""
+
+from __future__ import absolute_import
+
+import os
+import re
+import shutil
+import sys
+import tempfile
+
+root_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0]
+sys.path.append(os.path.join(root_path, os.path.pardir, os.path.pardir))
+
+import tests.unit_tests as unittest
+
+import madgraph.core.base_objects as base_objects
+import madgraph.core.helas_objects as helas_objects
+import madgraph.interface.master_interface as MGCmd
+
+from madgraph import InvalidCmd
+
+INTERF = base_objects.INTERFERENCE_ORDER
+
+
+def generate(line, model='sm'):
+    interface = MGCmd.MasterCmd()
+    interface.exec_cmd('import model %s' % model, printcmd=False, precmd=True)
+    interface.exec_cmd('generate %s' % line, printcmd=False, precmd=True,
+                       errorhandling=False)
+    return interface
+
+
+class TreeInterferenceGenerationTest(unittest.TestCase):
+
+    def test_diagrams_and_orders(self):
+        interface = generate('u u~ > z > e+ e- [treextree] u u~ > a > e+ e-')
+        self.assertEqual(len(interface._curr_amps), 1)
+        amp = interface._curr_amps[0]
+        diagrams = amp.get('diagrams')
+        self.assertEqual(len(diagrams), 2)
+        # left-hand (Z) diagram first, then the tagged right-hand (photon) one
+        self.assertEqual([d.get_order(INTERF) for d in diagrams], [0, 1])
+        self.assertEqual([d.get_order('WEIGHTED') for d in diagrams], [4, 4])
+
+        process = amp.get('process')
+        self.assertEqual(process['squared_orders'], {INTERF: 1})
+        self.assertEqual(process['split_orders'], [INTERF])
+        self.assertNotIn(INTERF, process.nice_string())
+        right = process['interference_process']
+        self.assertIsInstance(right, base_objects.Process)
+        self.assertEqual(right['required_s_channels'], [[22]])
+        self.assertEqual(right['id'], 0)
+
+        # the tag reaches the HELAS amplitudes (on a copy: the photon vertex
+        # of the model is untouched)
+        me = helas_objects.HelasMatrixElement(amp)
+        squared, amp_orders = me.get_split_orders_mapping()
+        self.assertEqual(sorted(so[0][0] for so in amp_orders), [0, 1])
+        self.assertEqual(sorted(so[0] for so in squared), [0, 1, 2])
+        model = process['model']
+        for inter in model['interactions']:
+            self.assertNotIn(INTERF, inter['orders'])
+
+    def test_right_hand_constraints_apply_to_the_right_only(self):
+        """'/ z' on the right does not remove the left-hand Z diagram; the
+        left-hand squared-order constraint applies to the product."""
+        interface = generate('u u~ > e+ e- / h [treextree] QED^2==4 '
+                             'u u~ > e+ e- / h z')
+        diagrams = interface._curr_amps[0].get('diagrams')
+        # left: photon and Z; right: photon only
+        self.assertEqual(sorted(d.get_order(INTERF) for d in diagrams),
+                         [0, 0, 1])
+
+    def test_multiprocess_without_crossing_reuse(self):
+        """Every subprocess gets the right-hand diagrams of its own legs."""
+        interface = generate('p p > z > e+ e- [treextree] p p > a > e+ e-')
+        self.assertTrue(interface._curr_amps)
+        for amp in interface._curr_amps:
+            right = amp.get('process')['interference_process']
+            self.assertEqual([l['id'] for l in right['legs']],
+                             [l['id'] for l in amp.get('process')['legs']])
+            self.assertEqual(sorted(d.get_order(INTERF)
+                                    for d in amp.get('diagrams')), [0, 1])
+
+    def test_no_interference(self):
+        self.assertRaisesRegex(InvalidCmd, 'no interference',
+                               generate,
+                               'u u~ > z > e+ e- [treextree] u u~ > h > e+ e-')
+
+
+class TreeInterferenceOutputTest(unittest.TestCase):
+    """Only the formats selecting squared split orders return the
+    interference alone: the others are refused before anything is written."""
+
+    def test_unsupported_formats_are_refused(self):
+        interface = generate('u u~ > z > e+ e- [treextree] u u~ > a > e+ e-')
+        tmpdir = tempfile.mkdtemp(prefix='treextree_')
+        try:
+            for fmt in ['mg7', 'standalone', 'matchbox']:
+                out = os.path.join(tmpdir, fmt)
+                self.assertRaisesRegex(InvalidCmd, 'cannot select the interference',
+                                       interface.exec_cmd,
+                                       'output %s %s -f' % (fmt, out),
+                                       printcmd=False, precmd=True,
+                                       errorhandling=False)
+                self.assertFalse(os.path.exists(out))
+        finally:
+            shutil.rmtree(tmpdir, ignore_errors=True)
+
+
+    def test_madevent_output(self):
+        """madevent: the interference defaults of the run card (HT/2 scale,
+        no systematics by default) and one power of alpha_s for all the
+        channels, the one of the product of a left- and a right-hand
+        amplitude: QCD (QCD=2) x electroweak (QCD=0) is alpha_s^1, whereas
+        the channels of the two sides have QCD orders 2 and 0."""
+        interface = generate('u d > u d QED=0 [treextree] u d > u d QCD=0')
+        tmpdir = tempfile.mkdtemp(prefix='treextree_')
+        try:
+            out = os.path.join(tmpdir, 'me')
+            interface.exec_cmd('output madevent %s -f' % out, printcmd=False,
+                               precmd=True, errorhandling=False)
+            run_card = open(os.path.join(out, 'Cards', 'run_card.dat')).read()
+            self.assertRegex(run_card, r'\n\s*3\s*=\s*dynamical_scale_choice')
+            self.assertRegex(run_card, r'\n\s*False\s*=\s*use_syst')
+            self.assertRegex(run_card, r'\n\s*none\s*=\s*systematics_program')
+            sub = os.path.join(out, 'SubProcesses')
+            powers = set()
+            for proc_dir in os.listdir(sub):
+                if not proc_dir.startswith('P'):
+                    continue
+                text = open(os.path.join(sub, proc_dir, 'config_nqcd.inc')).read()
+                powers.update(re.findall(r'NQCD\(\d+\)/(\d+)/', text))
+            self.assertEqual(powers, set(['1']))
+        finally:
+            shutil.rmtree(tmpdir, ignore_errors=True)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/unit_tests/interface/test_interference_syntax.py b/tests/unit_tests/interface/test_interference_syntax.py
new file mode 100644
index 0000000000..a9db26338e
--- /dev/null
+++ b/tests/unit_tests/interface/test_interference_syntax.py
@@ -0,0 +1,247 @@
+################################################################################
+#
+# Copyright (c) 2026 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+"""Syntax of the interference processes
+
+    generate g g > h > t t~ [LIxtree=QCD] g g > t t~     (loop-induced x tree)
+    generate p p > z > e+ e- [treextree] p p > a > e+ e-  (tree x tree)
+
+Only the parsing is tested here (the generation is tested in
+tests/unit_tests/core/test_tree_interference.py and
+tests/unit_tests/loop/test_loop_induced_x_tree.py): the line is split in a
+left-hand and a right-hand process, each with its own constraints, the two sides are checked
+for consistency, and the command is routed to the right interface. The
+keywords are case-insensitive and squared-order constraints are refused on
+the right-hand process.
+"""
+
+from __future__ import absolute_import
+
+import os
+import sys
+
+root_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0]
+sys.path.append(os.path.join(root_path, os.path.pardir, os.path.pardir))
+
+import tests.unit_tests as unittest
+
+import madgraph.core.base_objects as base_objects
+import madgraph.interface.madgraph_interface as mg_interface
+import madgraph.interface.master_interface as MGCmd
+
+from madgraph import InvalidCmd
+
+split = mg_interface.split_interference_line
+
+
+class InterferenceLineSplitTest(unittest.TestCase):
+    """split_interference_line, without any model."""
+
+    def test_lixtree_with_right_process(self):
+        self.assertEqual(split('g g > h > t t~ [LIxtree=QCD] g g > t t~'),
+                         ('g g > h > t t~ [noborn= QCD ]', 'LIxtree', 'g g > t t~'))
+
+    def test_lixtree_default_right_process(self):
+        self.assertEqual(split('g g > h > t t~ [LIxtree=QCD]'),
+                         ('g g > h > t t~ [noborn= QCD ]', 'LIxtree', None))
+
+    def test_keywords_are_case_insensitive(self):
+        for keyword in ['LIxtree', 'lixtree', 'LIXTREE', 'LiXTree']:
+            self.assertEqual(
+                split('g g > h > t t~ [%s=QCD] g g > t t~' % keyword)[1:],
+                ('LIxtree', 'g g > t t~'))
+        for keyword in ['treextree', 'TreeXTree', 'TREEXTREE']:
+            self.assertEqual(
+                split('p p > z > e+ e- [%s] p p > a > e+ e-' % keyword),
+                ('p p > z > e+ e-', 'treextree', 'p p > a > e+ e-'))
+
+    def test_constraints_go_to_their_side(self):
+        """Constraints between ']' and the right-hand process belong to the
+        left-hand one, those after it to the right-hand one; '@N' closes the
+        whole definition and stays with the left-hand line."""
+        self.assertEqual(
+            split('g g > h > t t~ [lixtree = QCD] QCD^2==6 QED = 2 '
+                  'g g > t t~ QED=0 / b @3'),
+            ('g g > h > t t~ [noborn= QCD ] QCD^2==6 QED = 2 @3', 'LIxtree',
+             'g g > t t~ QED=0 / b'))
+
+    def test_ordinary_lines_are_untouched(self):
+        for line in ['g g > z z [noborn=QCD]', 'p p > t t~ [QCD]',
+                     'p p > t t~ [real=QCD] QCD^2<=4', 'p p > t t~ QED=0',
+                     'p p > e+ e- [LOonly]']:
+            self.assertEqual(split(line), (line, '', None))
+
+    def test_malformed_brackets_are_refused(self):
+        for line in ['g g > h [LIxtree=] g g > h',
+                     'g g > h [LIxtree] g g > h',
+                     'p p > z > e+ e- [treextree=QED] p p > a > e+ e-',
+                     'p p > z > e+ e- [treextree]',
+                     'g g > h > t t~ [LIxtree=QCD] g g > t t~ [QCD]']:
+            self.assertRaises(InvalidCmd, split, line)
+
+    def test_routing(self):
+        """The master interface sends LIxtree where [noborn=] goes (MadLoop
+        model validation, then the loop-induced generation) and treextree to
+        the tree-level MadGraph interface."""
+        process_type = MGCmd.MasterCmd.extract_process_type
+        self.assertEqual(process_type('g g > h > t t~ [LIxtree] g g > t t~'),
+                         ('NLO', 'LIxtree', []))
+        self.assertEqual(process_type('g g > h > t t~ [LIxtree=QCD] g g > t t~'),
+                         ('NLO', 'LIxtree', ['QCD']))
+        self.assertEqual(process_type('g g > h > t t~ [lixtree= QCD] QCD^2==6 g g > t t~'),
+                         ('NLO', 'LIxtree', ['QCD']))
+        self.assertEqual(process_type('p p > z > e+ e- [TreeXTree] p p > a > e+ e-'),
+                         ('tree', 'treextree', []))
+        # unchanged
+        self.assertEqual(process_type('g g > z z [noborn=QCD]'), ('NLO', 'noborn', ['QCD']))
+        self.assertEqual(process_type('p p > t t~ [QCD]'), ('NLO', 'all', ['QCD']))
+        self.assertEqual(process_type('p p > t t~'), ('tree', None, []))
+
+
+class _StubModel(object):
+    merged_particles = {81: [1, 3, 5]}
+
+
+class InterferenceLegPdgsTest(unittest.TestCase):
+    """Merged (flavour-grouped) legs are compared through their members."""
+
+    def leg(self, ids, flavor=[]):
+        return base_objects.MultiLeg({'ids': ids, 'flavor': list(flavor)})
+
+    def test_expansion(self):
+        pdgs = mg_interface.MadGraphCmd.interference_leg_pdgs
+        self.assertEqual(pdgs(self.leg([81]), _StubModel), set([1, 3, 5]))
+        self.assertEqual(pdgs(self.leg([-81]), _StubModel), set([-1, -3, -5]))
+        self.assertEqual(pdgs(self.leg([81], [3]), _StubModel), set([3]))
+        self.assertEqual(pdgs(self.leg([21, 81, 2]), _StubModel), set([21, 1, 3, 5, 2]))
+
+
+class InterferenceProcessTest(unittest.TestCase):
+    """extract_process on interference lines, with loop_sm."""
+
+    interface = None
+
+    def setUp(self):
+        if InterferenceProcessTest.interface is None:
+            interface = MGCmd.MasterCmd()
+            interface.exec_cmd('import model loop_sm', printcmd=False,
+                               precmd=True)
+            InterferenceProcessTest.interface = interface
+        self.cmd = InterferenceProcessTest.interface
+
+    def test_process_defaults(self):
+        """An ordinary process carries no interference."""
+        procdef = self.cmd.extract_process('g g > t t~')
+        self.assertEqual(procdef['interference_mode'], '')
+        self.assertIsNone(procdef['interference_process'])
+        self.assertEqual(base_objects.Process()['interference_mode'], '')
+
+    def test_lixtree(self):
+        procdef = self.cmd.extract_process(
+            'g g > h > t t~ [LIxtree=QCD] QCD^2==6 g g > t t~ QED=0 / b @2')
+        self.assertEqual(procdef['interference_mode'], 'LIxtree')
+        self.assertEqual(procdef['NLO_mode'], 'noborn')
+        self.assertFalse(procdef['has_born'])
+        self.assertEqual(procdef['perturbation_couplings'], ['QCD'])
+        self.assertEqual(procdef['squared_orders'], {'QCD': 6})
+        self.assertEqual(procdef['id'], 2)
+        self.assertEqual(procdef['required_s_channels'], [[25]])
+
+        right = procdef['interference_process']
+        self.assertIsInstance(right, base_objects.ProcessDefinition)
+        self.assertEqual(right['NLO_mode'], 'tree')
+        self.assertTrue(right['has_born'])
+        self.assertEqual(right['interference_mode'], '')
+        self.assertEqual(right['orders'], {'QED': 0})
+        self.assertEqual(right['squared_orders'], {})
+        self.assertEqual(right['required_s_channels'], [])
+        self.assertEqual(right['forbidden_particles'], [5])
+        self.assertEqual(right['id'], 0)
+        self.assertEqual([l['ids'] for l in right['legs']],
+                         [[21], [21], [6], [-6]])
+
+    def test_lixtree_default_right_process(self):
+        """Without a right-hand process, the trees of the same external legs:
+        the '> h >' requirement and the loop orders are not carried over."""
+        procdef = self.cmd.extract_process('g g > h > t t~ QED=2 [lixtree=QCD]')
+        right = procdef['interference_process']
+        self.assertEqual([l['ids'] for l in right['legs']],
+                         [l['ids'] for l in procdef['legs']])
+        self.assertEqual(right['required_s_channels'], [])
+        self.assertEqual(right['orders'], {})
+        self.assertEqual(right['perturbation_couplings'], [])
+        self.assertEqual(right['NLO_mode'], 'tree')
+
+    def test_treextree(self):
+        procdef = self.cmd.extract_process(
+            'u u~ > z > e+ e- [TreeXTree] QED^2==4 p p > a > e+ e-')
+        self.assertEqual(procdef['interference_mode'], 'treextree')
+        self.assertEqual(procdef['NLO_mode'], 'tree')
+        self.assertEqual(procdef['perturbation_couplings'], [])
+        self.assertEqual(procdef['squared_orders'], {'QED': 4})
+        self.assertEqual(procdef['required_s_channels'], [[23]])
+        right = procdef['interference_process']
+        self.assertEqual(right['required_s_channels'], [[22]])
+
+    def test_display_round_trip(self):
+        """The printed definition is the input syntax, and parses back to the
+        same interference."""
+        # (multiparticles are printed as 'u/d/...', which is not input
+        # syntax: hence particles only)
+        for line in ['g g > h > t t~ [LIxtree=QCD] QCD^2==6 g g > t t~ QED=0',
+                     'u u~ > z > e+ e- [treextree] u u~ > a > e+ e- / h']:
+            procdef = self.cmd.extract_process(line)
+            printed = procdef.nice_string(prefix=False)
+            self.assertIn('[ %s' % procdef['interference_mode'], printed)
+            again = self.cmd.extract_process(printed)
+            self.assertEqual(again['interference_mode'], procdef['interference_mode'])
+            self.assertEqual(again['squared_orders'], procdef['squared_orders'])
+            self.assertEqual(again['interference_process']['legs'],
+                             procdef['interference_process']['legs'])
+            self.assertEqual(again['interference_process']['required_s_channels'],
+                             procdef['interference_process']['required_s_channels'])
+
+    def test_right_hand_squared_orders_are_refused(self):
+        for line in ['g g > h > t t~ [LIxtree=QCD] g g > t t~ QCD^2==4',
+                     'u u~ > z > e+ e- [treextree] p p > a > e+ e- QED^2<=4']:
+            self.assertRaisesRegex(InvalidCmd, 'right-hand process',
+                                   self.cmd.extract_process, line)
+
+    def test_inconsistent_legs_are_refused(self):
+        for line in [
+                # one more final-state particle
+                'g g > h > t t~ [LIxtree=QCD] g g > t t~ g',
+                # same legs, other order
+                'g g > h > t t~ [LIxtree=QCD] g g > t~ t',
+                # other initial state
+                'g g > h > t t~ [LIxtree=QCD] u u~ > t t~',
+                # a right-hand leg narrower than the left-hand one
+                'p p > z > e+ e- [treextree] u u~ > a > e+ e-',
+                # other polarisation
+                'u u~ > z > e+{L} e- [treextree] u u~ > a > e+ e-']:
+            self.assertRaisesRegex(InvalidCmd, 'two sides',
+                                   self.cmd.extract_process, line)
+
+    def test_wider_right_hand_legs_are_accepted(self):
+        procdef = self.cmd.extract_process(
+            'u u~ > z > e+ e- [treextree] p p > a > e+ e-')
+        self.assertEqual(procdef['interference_mode'], 'treextree')
+
+    def test_decay_chains_are_refused(self):
+        self.assertRaises(InvalidCmd, self.cmd.extract_process,
+                    'u u~ > z > e+ e- [treextree] u u~ > a > e+ e-, a > e+ e-')
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/tests/unit_tests/loop/test_loop_induced_x_tree.py b/tests/unit_tests/loop/test_loop_induced_x_tree.py
new file mode 100644
index 0000000000..ccb185a2ab
--- /dev/null
+++ b/tests/unit_tests/loop/test_loop_induced_x_tree.py
@@ -0,0 +1,172 @@
+################################################################################
+#
+# Copyright (c) 2026 The MadGraph5_aMC@NLO Development team and Contributors
+#
+# This file is a part of the MadGraph5_aMC@NLO project, an application which
+# automatically generates Feynman diagrams and matrix elements for arbitrary
+# high-energy processes in the Standard Model and beyond.
+#
+# It is subject to the MadGraph5_aMC@NLO license which should accompany this
+# distribution.
+#
+# For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
+#
+################################################################################
+"""Generation of the interference of a loop-induced amplitude with tree-level
+ones:
+
+    generate g g > h > t t~ [LIxtree=QCD]                 (right side implied)
+    generate g g > h > t t~ [LIxtree=QCD] g g > t t~ QED=0
+
+The tree-level diagrams of the right-hand process are stored among the
+'loop_UVCT_diagrams' of the loop-induced amplitude, each with one unit of the
+hidden order base_objects.INTERFERENCE_ORDER, and the squared-order constraint
+INTERF^2==1 keeps only the loop x tree products. The numerical check of the
+result (against the tree-level matrix element, [sqrvirt=] and the original
+wiki recipe) is the acceptance test test_loop_induced_x_tree_standalone in
+tests/acceptance_tests/test_cmd_madloop.py.
+"""
+
+from __future__ import absolute_import
+
+import os
+import sys
+
+root_path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0]
+sys.path.append(os.path.join(root_path, os.path.pardir, os.path.pardir))
+
+import tests.unit_tests as unittest
+
+import madgraph.core.base_objects as base_objects
+import madgraph.interface.master_interface as MGCmd
+import madgraph.loop.loop_base_objects as loop_base_objects
+import madgraph.loop.loop_helas_objects as loop_helas_objects
+
+from madgraph import InvalidCmd
+
+INTERF = base_objects.INTERFERENCE_ORDER
+
+
+def generate(line):
+    """A fresh interface (the hidden order is added to its model) with 'line'
+    generated."""
+    interface = MGCmd.MasterCmd()
+    interface.exec_cmd('import model loop_sm', printcmd=False, precmd=True)
+    interface.exec_cmd('generate %s' % line, printcmd=False, precmd=True,
+                       errorhandling=False)
+    return interface
+
+
+class LoopInducedXTreeGenerationTest(unittest.TestCase):
+
+    gg_h_ttx = None
+
+    def setUp(self):
+        if LoopInducedXTreeGenerationTest.gg_h_ttx is None:
+            LoopInducedXTreeGenerationTest.gg_h_ttx = \
+                                generate('g g > h > t t~ [LIxtree=QCD]')
+        self.interface = LoopInducedXTreeGenerationTest.gg_h_ttx
+
+    def test_diagrams(self):
+        """4 loops (the top triangles of g g > h, h > t t~) and the 3 QCD
+        trees of g g > t t~, which carry one unit of the hidden order."""
+        self.assertEqual(len(self.interface._curr_amps), 1)
+        amp = self.interface._curr_amps[0]
+        self.assertFalse(amp['has_born'])
+        self.assertEqual(len(amp['loop_diagrams']), 4)
+        self.assertEqual(len(amp['loop_UVCT_diagrams']), 3)
+        for loop in amp['loop_diagrams']:
+            self.assertEqual(loop.get_order(INTERF), 0)
+        for tree in amp['loop_UVCT_diagrams']:
+            self.assertIsInstance(tree, loop_base_objects.LoopUVCTDiagram)
+            self.assertEqual(tree['UVCT_orders'], {INTERF: 1})
+            self.assertEqual(tree['UVCT_couplings'], [1])
+            self.assertEqual(tree.get_order(INTERF), 1)
+            self.assertEqual(tree.get_order('QCD'), 2)
+            self.assertEqual(tree.get_order('QED'), 0)
+            # the hidden order does not change the weight
+            self.assertEqual(tree.get_order('WEIGHTED'), 2)
+
+    def test_process(self):
+        """The left-hand process selects the interference with the hidden
+        order, which is not shown; the subprocess keeps the right-hand tree
+        process of its own legs, without the '> h >' requirement."""
+        process = self.interface._curr_amps[0]['process']
+        self.assertEqual(process['squared_orders'][INTERF], 1)
+        self.assertEqual(process.get_squared_order_type(INTERF), '==')
+        self.assertIn(INTERF, process['split_orders'])
+        self.assertIn('QCD', process['split_orders'])
+        self.assertNotIn(INTERF, process.nice_string())
+        self.assertIn('[ LIxtree = QCD ] g g > t t~', process.nice_string())
+        model = process['model']
+        self.assertEqual(model['order_hierarchy'][INTERF], 0)
+        self.assertIn(INTERF, model['coupling_orders'])
+
+        right = process['interference_process']
+        self.assertIsInstance(right, base_objects.Process)
+        self.assertEqual([l['id'] for l in right['legs']], [21, 21, 6, -6])
+        self.assertEqual(right['required_s_channels'], [])
+        # minimal order, as 'generate g g > t t~' would choose
+        self.assertEqual(right['orders'], {'WEIGHTED': 2})
+
+    def test_helas_split_orders(self):
+        """In the generated code each tree amplitude has INTERF=1 (not 2: the
+        UVCT orders are added once), and the squared orders of the three
+        products (loop x tree, loop x loop, tree x tree) are all defined, the
+        first one being the selected one."""
+        me = loop_helas_objects.LoopHelasMatrixElement(
+                     self.interface._curr_amps[0], optimized_output=True)
+        squared, amp_orders = me.get_split_orders_mapping()
+        split_orders = me.get('processes')[0]['split_orders']
+        interf = split_orders.index(INTERF)
+        self.assertEqual(sorted(so[interf] for so, _ in squared), [0, 1, 2])
+        loop_amp_orders = sorted(so[interf] for so, _ in
+                                 amp_orders['loop_amp_orders'])
+        self.assertEqual(loop_amp_orders, [0, 1])
+        for diagram in me.get_loop_UVCT_diagrams():
+            for amp in diagram.get_loop_UVCTamplitudes():
+                self.assertEqual(amp.get('orders')[INTERF], 1)
+
+    def test_alpha_s_power(self):
+        """loop (QCD=2) x tree (QCD=2): alpha_s^2 for every event, which
+        madevent writes in config_nqcd.inc for the systematics."""
+        import madgraph.iolibs.export_v4 as export_v4
+        me = loop_helas_objects.LoopHelasMatrixElement(
+                     self.interface._curr_amps[0], optimized_output=True)
+        self.assertEqual(export_v4.interference_alpha_s_power([me]), 2)
+        self.assertEqual(export_v4.interference_nqcd_list([5, 0, 3], [me]),
+                         [2, 2, 2])
+
+
+class LoopInducedXTreeMultiProcessTest(unittest.TestCase):
+
+    def test_subprocesses_without_loop_are_dropped(self):
+        """p p > h > t t~: the light-quark channels have trees but no loop,
+        so no interference: only g g survives (the trees alone would have
+        kept the subprocess alive)."""
+        interface = generate('p p > h > t t~ [LIxtree=QCD]')
+        self.assertEqual([[l['id'] for l in amp['process']['legs']]
+                          for amp in interface._curr_amps],
+                         [[21, 21, 6, -6]])
+
+    def test_right_hand_constraints(self):
+        """The constraints of the right-hand process apply to the trees only,
+        those between ']' and the right-hand process to the interference."""
+        interface = generate('g g > h > t t~ [LIxtree=QCD] QCD^2==4 '
+                             'g g > t t~ QED=0 / b')
+        amp = interface._curr_amps[0]
+        self.assertEqual(len(amp['loop_UVCT_diagrams']), 3)
+        self.assertEqual(amp['process']['squared_orders'],
+                         {'QCD': 4, INTERF: 1})
+        right = amp['process']['interference_process']
+        self.assertEqual(right['orders'], {'QED': 0})
+        self.assertEqual(right['forbidden_particles'], [5])
+
+    def test_no_tree_level_interference(self):
+        """g g > z z has no tree-level diagram: nothing to interfere with."""
+        self.assertRaisesRegex(InvalidCmd, 'no interference',
+                               generate, 'g g > h > z z [LIxtree=QCD]')
+
+
+if __name__ == '__main__':
+    unittest.main()

From 08080763238572e48c901c9adb89af5dd0dd3934 Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Fri, 9 Oct 2026 09:05:38 +0200
Subject: [PATCH 2/5] Interference processes: fixes from review

- Left-hand amplitude constraints ('==', '>') no longer become squared-order
  constraints on the left x right product: 'u d > u d QED==0 [treextree]
  u d > u d QCD==0' got QED^2==0 and no interference at all.
- Order constraints after the bracket must be identifiers followed by a
  process: '25 > 6 -6' or 'h > 5 -5' as right-hand process were eaten as
  'ORDER > n'.
- The LIxtree bracket is printed with the right-hand process, after the
  left-hand exclusions, so that a subprocess string parses back.
- alpha_s power: the configs.inc writers take one power per subprocess and
  use that of the subprocess owning each channel; an interference process
  never shares a subprocess group with an ordinary one (the interference
  mode is part of the process class; ordinary grouping unchanged); warning
  when the subprocesses of one group have different powers.
- The hidden order no longer leaks into later processes of the session: no
  zero INTERF in diagram orders, no extra level in get_particles_hierarchy.
- --diagram_filter also filters the right-hand diagrams.
- A Process pickled before the interference keys prints again.
- LIxtree: trees whose loop partners are removed by the loop filters run
  after the squared-order selection are dropped.

Tests for each; full unit suite 2464 OK, acceptance interference tests OK.
---
 madgraph/core/base_objects.py                 |  49 ++++++---
 madgraph/core/diagram_generation.py           |  12 +-
 madgraph/interface/madgraph_interface.py      |  20 +++-
 madgraph/iolibs/export_v4.py                  | 104 ++++++++++++------
 madgraph/iolibs/group_subprocs.py             |  10 +-
 madgraph/loop/loop_base_objects.py            |   3 +
 madgraph/loop/loop_diagram_generation.py      |  36 ++++++
 .../unit_tests/core/test_tree_interference.py |  77 +++++++++++++
 .../interface/test_interference_syntax.py     |  22 ++++
 .../loop/test_loop_induced_x_tree.py          |  36 +++++-
 10 files changed, 306 insertions(+), 63 deletions(-)

diff --git a/madgraph/core/base_objects.py b/madgraph/core/base_objects.py
index 8f3b95edd0..621195cb4e 100755
--- a/madgraph/core/base_objects.py
+++ b/madgraph/core/base_objects.py
@@ -1964,8 +1964,12 @@ def get_particles_hierarchy(self):
         with lowest value) where it has an interaction.
         """
         
-        # Find coupling orders in model
-        coupling_orders = self.get('coupling_orders')
+        # Find coupling orders in model (not the hidden order of interference
+        # processes, which no interaction carries: with its hierarchy 0 it
+        # would otherwise add a level and move the particles of order-less
+        # interactions there)
+        coupling_orders = [k for k in self.get('coupling_orders')
+                           if k != INTERFERENCE_ORDER]
         # Loop through the different coupling hierarchy values, so we
         # start with the most dominant and proceed to the least dominant
         hierarchy = sorted(list(set([self.get('order_hierarchy')[k] for \
@@ -1976,7 +1980,7 @@ def get_particles_hierarchy(self):
         for value in hierarchy:
             orders.append([ k for (k, v) in \
                             self.get('order_hierarchy').items() if \
-                            v == value ])
+                            v == value and k != INTERFERENCE_ORDER ])
 
         # Extract the interaction that correspond to the different
         # coupling hierarchies, and the corresponding particles
@@ -3397,6 +3401,10 @@ def calculate_orders(self, model):
             weight += sum([model.get('order_hierarchy')[c]*n for \
                               (c,n) in couplings.items()])
         coupling_orders['WEIGHTED'] = weight
+        # the hidden order of interference processes, once in the model, is
+        # only listed by the diagrams which carry it
+        if not coupling_orders.get(INTERFERENCE_ORDER):
+            coupling_orders.pop(INTERFERENCE_ORDER, None)
         self.set('orders', coupling_orders)
 
     def pass_squared_order_constraints(self, diag_multiplier, squared_orders,
@@ -3936,8 +3944,9 @@ def get_sorted_keys(self):
                 'forbidden_onsh_s_channels', 'forbidden_s_channels',
                 'forbidden_particles', 'is_decay_chain', 'decay_chains',
                 'legs_with_decays', 'perturbation_couplings', 'has_born', 
-                'NLO_mode', 'split_orders', 'born_sq_orders',
-                'interference_mode', 'interference_process']
+                'NLO_mode', 'split_orders', 'born_sq_orders'] + \
+               [key for key in ('interference_mode', 'interference_process')
+                if key in self]   # absent from a process pickled before them
 
     def get_interference_mode(self):
         """'LIxtree', 'treextree' or '' (also for a process unpickled from
@@ -3947,14 +3956,20 @@ def get_interference_mode(self):
 
     def interference_string(self):
         """The tail of the process syntax for an interference process: the
-        '[ treextree ]' bracket (the LIxtree one is printed with the
-        perturbation orders) followed by the right-hand process, each piece
-        ending with a space. Empty for an ordinary process."""
+        '[ treextree ]' or '[ LIxtree = ORDERS ]' bracket followed by the
+        right-hand process, each piece ending with a space. It comes after
+        everything that belongs to the left-hand process (orders, squared
+        orders, exclusions), since what follows the bracket is read back as
+        the right-hand process. Empty for an ordinary process."""
 
         mode = self.get_interference_mode()
         if not mode:
             return ''
-        mystr = '[ treextree ] ' if mode == 'treextree' else ''
+        if mode == 'treextree':
+            mystr = '[ treextree ] '
+        else:
+            mystr = '[ LIxtree = %s ] ' % ' '.join(sorted(
+                                              self['perturbation_couplings']))
         right = self['interference_process']
         if right is not None:
             mystr += right.nice_string(prefix=False, print_weighted=False) + ' '
@@ -4044,14 +4059,14 @@ def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturba
               self['constrained_orders'][key][1], self['constrained_orders'][key][0]) 
                     for key in sorted(self['constrained_orders'].keys()))  + ' '
 
-        # Add perturbation_couplings
-        if print_perturbated and self['perturbation_couplings']:
+        # Add perturbation_couplings (for LIxtree: with the right-hand process,
+        # see interference_string)
+        if print_perturbated and self['perturbation_couplings'] and \
+                                   self.get_interference_mode() != 'LIxtree':
             mystr = mystr + '[ '
             if self['NLO_mode']!='tree':
                 if self['NLO_mode']=='virt' and not self['has_born']:
                     mystr = mystr + 'sqrvirt = '
-                elif self.get_interference_mode() == 'LIxtree':
-                    mystr = mystr + 'LIxtree = '
                 else:
                     mystr = mystr + self['NLO_mode'] + ' = '
             for order in sorted(self['perturbation_couplings']):
@@ -4985,14 +5000,14 @@ def nice_string(self, indent=0, print_weighted=False, prefix=True):
                                      (key,(value, operator)) 
                                    in self['constrained_orders'].items()) + ' '
 
-        # Add perturbation_couplings
-        if self['perturbation_couplings']:
+        # Add perturbation_couplings (for LIxtree: with the right-hand process,
+        # see interference_string)
+        if self['perturbation_couplings'] and \
+                                   self.get_interference_mode() != 'LIxtree':
             mystr = mystr + '[ '
             if self['NLO_mode']!='tree':
                 if self['NLO_mode']=='virt' and not self['has_born']:
                     mystr = mystr + 'sqrvirt = '
-                elif self.get_interference_mode() == 'LIxtree':
-                    mystr = mystr + 'LIxtree = '
                 else:
                     mystr = mystr + self['NLO_mode'] + ' = '
             for order in self['perturbation_couplings']:
diff --git a/madgraph/core/diagram_generation.py b/madgraph/core/diagram_generation.py
index 7726faa968..97a4f066e2 100755
--- a/madgraph/core/diagram_generation.py
+++ b/madgraph/core/diagram_generation.py
@@ -1392,7 +1392,7 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False):
         self.get('process').get('legs').sort(pert=pertur)
 
         if treextree:
-            res = self.add_interference_diagrams(res)
+            res = self.add_interference_diagrams(res, diagram_filter)
 
         # Set diagrams to res if not asked to be returned
         if not returndiag:
@@ -1401,7 +1401,7 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False):
         else:
            return not failed_crossing, res
 
-    def add_interference_diagrams(self, left_diagrams):
+    def add_interference_diagrams(self, left_diagrams, diagram_filter=False):
         """For an interference process 'left [treextree] right', add to the
         diagrams of this (left-hand) process those of the right-hand process
         with the same external legs, each tagged with one unit of the hidden
@@ -1409,7 +1409,8 @@ def add_interference_diagrams(self, left_diagrams):
         the process to the whole list: with INTERF^2==1 (see
         prepare_interference_process) only the products of a left-hand and a
         right-hand diagram, 2 Re(A_left A_right^*), contribute. All diagrams of
-        both sides are kept as integration channels. Returns the new diagram
+        both sides are kept as integration channels. A user diagram filter
+        ('--diagram_filter') applies to both sides. Returns the new diagram
         list, empty if either side has no diagram."""
 
         process = self.get('process')
@@ -1418,8 +1419,11 @@ def add_interference_diagrams(self, left_diagrams):
 
         right_diagrams = base_objects.DiagramList()
         if left_diagrams:
+            right_amplitude = Amplitude()
+            right_amplitude.set('process', right_process)
             try:
-                right_diagrams = Amplitude(right_process).get('diagrams')
+                right_amplitude.generate_diagrams(diagram_filter=diagram_filter)
+                right_diagrams = right_amplitude.get('diagrams')
             except InvalidCmd:
                 pass
         if not left_diagrams or not right_diagrams:
diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py
index b6ad789886..5781be38c1 100755
--- a/madgraph/interface/madgraph_interface.py
+++ b/madgraph/interface/madgraph_interface.py
@@ -159,9 +159,10 @@ class DuplicateParticle(madgraph.InvalidCmd):
 # 'QED == 2', ... The constraints written between the bracket and the
 # right-hand process belong to the left-hand one (the usual place of the
 # squared-order constraints of a loop process), those written after the
-# right-hand process belong to the right-hand one.
+# right-hand process belong to the right-hand one. An order name is an
+# identifier, never a pdg code: '25 > 6 -6' is a process, not 'ORDER > 6'.
 _LEADING_ORDER_CONSTRAINT = re.compile(
-    r"^\s*(?:\w|\^2)+\s*(?:===|==|<=|>=|!=|=|<|>)\s*-?\d+(?=\s|$)")
+    r"^\s*[A-Za-z_]\w*(?:\^2)?\s*(?:===|==|<=|>=|!=|=|<|>)\s*-?\d+(?=\s|$)")
 
 
 def interference_mode_of(option, orders):
@@ -237,8 +238,13 @@ def split_interference_line(line):
     tail = tail.strip()
     constraint = _LEADING_ORDER_CONSTRAINT.match(tail)
     while constraint:
+        rest = tail[constraint.end():].strip()
+        # a constraint is followed by more constraints or by the right-hand
+        # process, which has an arrow: in 'h > 5 -5', 'h > 5' is no constraint
+        if rest and not re.search(r'>\D', rest):
+            break
         left_constraints.append(constraint.group(0).strip())
-        tail = tail[constraint.end():].strip()
+        tail = rest
         constraint = _LEADING_ORDER_CONSTRAINT.match(tail)
 
     right = tail if tail else None
@@ -6902,9 +6908,15 @@ def extract_interference_process(self, left, mode, right, proc_number=0,
         constraint and no exclusion is carried over from the loop process.
         """
 
+        # An amplitude-level '==' or '>' constraint on the left-hand process
+        # constrains the left-hand amplitudes only: unlike for an ordinary
+        # process, it must not turn into a squared-order constraint, which
+        # here would act on the left x right products (u d > u d QED==0
+        # [treextree] u d > u d QCD==0 would get QED^2==0 and no interference
+        # at all). Explicit '^2' constraints are kept.
         procdef = self.extract_process(left, proc_number=proc_number,
                                        overall_orders=overall_orders,
-                                       avoid_squared_orders=avoid_squared_orders)
+                                       avoid_squared_orders=True)
         if not procdef:
             raise self.InvalidCmd("Empty or wrong format process, please try again.")
         if procdef.get('decay_chains'):
diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py
index 2f75146242..e7cfa11a86 100755
--- a/madgraph/iolibs/export_v4.py
+++ b/madgraph/iolibs/export_v4.py
@@ -7812,7 +7812,6 @@ def generate_subprocess_directory(self, matrix_element,
             matrix_element)
 
         filename = pjoin(Ppath, 'config_nqcd.inc')
-        nqcd_list = interference_nqcd_list(nqcd_list, [matrix_element])
         self.write_config_nqcd_file(writers.FortranWriter(filename),
                                nqcd_list)
 
@@ -8870,7 +8869,8 @@ def write_configs_file(self, writer, matrix_element):
                                                             [[c[1]] for c in configs],
                                                             mapconfigs,
                                                             nexternal, ninitial,
-                                                            model)
+                                                            model,
+                    alpha_s_powers=interference_alpha_s_powers([matrix_element]))
 
     #===========================================================================
     # write_run_configs_file
@@ -8896,7 +8896,8 @@ def write_run_config_file(self, writer):
     # write_configs_file_from_diagrams
     #===========================================================================
     def write_configs_file_from_diagrams(self, writer, configs, mapconfigs,
-                                         nexternal, ninitial, model):
+                                         nexternal, ninitial, model,
+                                         alpha_s_powers=None):
         """Write the actual configs.inc file.
         
         configs is the diagrams corresponding to configs (each
@@ -8985,14 +8986,22 @@ def write_configs_file_from_diagrams(self, writer, configs, mapconfigs,
             lines.append("data mapconfig(%d)/%d/" % (nconfigs,
                                                      mapconfigs[iconfig]))
             lines.append("data tstrategy(%d)/%d/" % (nconfigs, tchannels_strategy))
-            # Number of QCD couplings in this diagram
+            # Number of QCD couplings in this diagram, i.e. the power of
+            # alpha_s of the events of this channel; for an interference
+            # process the power of the left x right product instead, given
+            # per subprocess in alpha_s_powers (see
+            # interference_alpha_s_power)
             nqcd = 0
-            for h in helas_diags:
+            for isub, h in enumerate(helas_diags):
                 if h:
-                    try:
-                        nqcd = h.calculate_orders()['QCD']
-                    except KeyError:
-                        pass
+                    power = alpha_s_powers[isub] if alpha_s_powers else None
+                    if power is not None:
+                        nqcd = power
+                    else:
+                        try:
+                            nqcd = h.calculate_orders()['QCD']
+                        except KeyError:
+                            pass
                     break
                 else:
                     continue
@@ -9048,7 +9057,7 @@ def write_configs_file_from_diagrams(self, writer, configs, mapconfigs,
     #===========================================================================
     def write_configs_file_from_onia_diagrams(self, writer, configs, mapconfigs,
                                          nexternal, ninitial, model,
-                                         onia_pairs=None):
+                                         onia_pairs=None, alpha_s_powers=None):
         """Write the actual configs.inc file.
 
         configs is the diagrams corresponding to configs (each
@@ -9137,14 +9146,22 @@ def write_configs_file_from_onia_diagrams(self, writer, configs, mapconfigs,
             lines.append("data mapconfig(%d)/%d/" % (nconfigs,
                                                      mapconfigs[iconfig]))
             lines.append("data tstrategy(%d)/%d/" % (nconfigs, tchannels_strategy))
-            # Number of QCD couplings in this diagram
+            # Number of QCD couplings in this diagram, i.e. the power of
+            # alpha_s of the events of this channel; for an interference
+            # process the power of the left x right product instead, given
+            # per subprocess in alpha_s_powers (see
+            # interference_alpha_s_power)
             nqcd = 0
-            for h in helas_diags:
+            for isub, h in enumerate(helas_diags):
                 if h:
-                    try:
-                        nqcd = h.calculate_orders()['QCD']
-                    except KeyError:
-                        pass
+                    power = alpha_s_powers[isub] if alpha_s_powers else None
+                    if power is not None:
+                        nqcd = power
+                    else:
+                        try:
+                            nqcd = h.calculate_orders()['QCD']
+                        except KeyError:
+                            pass
                     break
                 else:
                     continue
@@ -10033,7 +10050,6 @@ def generate_subprocess_directory(self, subproc_group,
             subproc_diagrams_for_config)
 
         filename = 'config_nqcd.inc'
-        nqcd_list = interference_nqcd_list(nqcd_list, matrix_elements)
         self.write_config_nqcd_file(writers.FortranWriter(filename),
                                     nqcd_list)
 
@@ -10564,18 +10580,31 @@ def write_configs_file(self, writer, subproc_group, diagrams_for_config):
         (nexternal, ninitial) = subproc_group.get_nexternal_ninitial()
         onia_pairs = matrix_elements[0].get_onia_pairs()
 
+        # power of alpha_s per subprocess for the interference processes
+        alpha_s_powers = interference_alpha_s_powers(matrix_elements)
+        if len(set(p for p in alpha_s_powers if p is not None)) > 1:
+            logger.warning("The interference subprocesses of %s have different "
+                           "powers of alpha_s, but share the integration "
+                           "channels: each channel takes the power of its first "
+                           "subprocess, so the scale variations of the "
+                           "systematics (off by default) would be approximate. "
+                           "Generate these subprocesses with separate 'add "
+                           "process' lines to avoid it."
+                           % subproc_group.get('name'))
         if onia_pairs:
             return len(diagrams), \
                self.write_configs_file_from_onia_diagrams(writer, diagrams,
                                                 config_numbers,
                                                 nexternal, ninitial,
-                                                     model,onia_pairs)
+                                                     model,onia_pairs,
+                                                alpha_s_powers=alpha_s_powers)
         else:
             return len(diagrams), \
                self.write_configs_file_from_diagrams(writer, diagrams,
                                                 config_numbers,
                                                 nexternal, ninitial,
-                                                     model)
+                                                     model,
+                                                alpha_s_powers=alpha_s_powers)
 
 
 
@@ -13889,21 +13918,28 @@ def interference_alpha_s_power(matrix_elements):
             return power // 2
     return None
 
-def interference_nqcd_list(nqcd_list, matrix_elements):
-    """config_nqcd.inc entries of an interference process: the same power of
-    alpha_s for every channel (see interference_alpha_s_power)."""
-
-    power = interference_alpha_s_power(matrix_elements)
-    if power is False:
-        return nqcd_list
-    if power is None:
-        logger.warning("The power of alpha_s of this interference is not "
-                       "unique: the scale variations of the systematics "
-                       "(off by default for interference processes) would be "
-                       "wrong. Constrain the QCD order of the interference "
-                       "(e.g. QCD^2==n after ']') to make it unique.")
-        return nqcd_list
-    return [power] * len(nqcd_list)
+def interference_alpha_s_powers(matrix_elements):
+    """The power of alpha_s to give the events of each matrix element in
+    config_nqcd.inc (see interference_alpha_s_power): the power of the left x
+    right products for an interference process, None for an ordinary one,
+    whose events keep the QCD order of the diagram of their channel. Kept
+    per matrix element, so that a subprocess group mixing interference and
+    ordinary processes gets the right value for each channel."""
+
+    powers = []
+    for me in matrix_elements:
+        power = interference_alpha_s_power([me])
+        if power is None:
+            logger.warning("The power of alpha_s of the interference %s is not "
+                           "unique: the scale variations of the systematics "
+                           "(off by default for interference processes) would "
+                           "be wrong. Constrain the QCD order of the "
+                           "interference (e.g. QCD^2==n after ']') to make it "
+                           "unique." % me.get('processes')[0].nice_string(
+                                                                 prefix=False))
+        # 'is False': a power of 0 is a valid one
+        powers.append(None if power is False else power)
+    return powers
 
 def loop_induced_not_supported_msg(format, process=None):
     """Refusal text for a format that cannot write a LoopHelasMatrixElement."""
diff --git a/madgraph/iolibs/group_subprocs.py b/madgraph/iolibs/group_subprocs.py
index 11164c22f6..65e53420aa 100755
--- a/madgraph/iolibs/group_subprocs.py
+++ b/madgraph/iolibs/group_subprocs.py
@@ -517,6 +517,12 @@ def find_process_classes(amplitudes, criteria):
             pols = tuple(base_objects.canonical_polarization(
                              l.get('polarization'))
                          for l in process.get('legs'))
+            # Likewise an interference process ('[LIxtree=QCD]', '[treextree]')
+            # never shares a group with an ordinary one, even with the same
+            # '@N': the events of a group share the power of alpha_s of each
+            # channel (config_nqcd.inc), which differs between the two.
+            # Ordinary processes all have '' here: their grouping is unchanged.
+            interference = process.get_interference_mode()
 
             # This is where the requirements for which particles to
             # combine are defined. Include p.get('is_part') in
@@ -531,7 +537,7 @@ def find_process_classes(amplitudes, criteria):
                              abs(p.get('color')),l.get('onshell'),l.get('onium').get('id')) for (p, l) \
                              in zip(is_parts + fs_parts, process.get('legs'))],
                            amplitude.get('process').get('id'),
-                           process.get('id'), pols]
+                           process.get('id'), pols, interference]
             if (criteria=="madweight"):
               proc_class = [ [(abs(p.get('pdg_code'))==5, abs(p.get('pdg_code'))==11, 
                            abs(p.get('pdg_code'))==13, abs(p.get('pdg_code'))==15) for p in \
@@ -546,7 +552,7 @@ def find_process_classes(amplitudes, criteria):
                              abs(p.get('color')),l.get('onshell'),l.get('onium').get('id')) for (p, l) \
                              in zip(is_parts + fs_parts, process.get('legs'))],
                            amplitude.get('process').get('id'),
-                           process.get('id'), pols]
+                           process.get('id'), pols, interference]
 
             try:
                 amplitude_classes[iamp] = proc_classes.index(proc_class)
diff --git a/madgraph/loop/loop_base_objects.py b/madgraph/loop/loop_base_objects.py
index e0bd17c1d0..fc0534df00 100755
--- a/madgraph/loop/loop_base_objects.py
+++ b/madgraph/loop/loop_base_objects.py
@@ -1401,6 +1401,9 @@ def calculate_orders(self, model):
             weight += sum([model.get('order_hierarchy')[c]*n for \
                               (c,n) in couplings.items()])
         coupling_orders['WEIGHTED'] = weight
+        # only the interfering trees (LIxtree) list the hidden order
+        if not coupling_orders.get(base_objects.INTERFERENCE_ORDER):
+            coupling_orders.pop(base_objects.INTERFERENCE_ORDER, None)
         self.set('orders', coupling_orders)
 
     def nice_string(self):
diff --git a/madgraph/loop/loop_diagram_generation.py b/madgraph/loop/loop_diagram_generation.py
index e7f0cbf3a8..a27f7bcd0d 100755
--- a/madgraph/loop/loop_diagram_generation.py
+++ b/madgraph/loop/loop_diagram_generation.py
@@ -893,6 +893,12 @@ def generate_diagrams(self, loop_filter=None, diagram_filter=None):
         # will find examples of common filters.
         self.user_filter(model,self['structure_repository'], filter=loop_filter)
 
+        # LIxtree: the filters above may have removed every loop that some
+        # interfering trees were paired with; such trees no longer contribute
+        # to the interference and are dropped
+        if self['process'].get_interference_mode() == 'LIxtree':
+            self.drop_orphan_interference_trees()
+
         # Set the necessary UV/R2 CounterTerms for each loop diagram generated
         self.set_LoopCT_vertices()
 
@@ -1266,6 +1272,36 @@ def set_interference_trees(self):
         ldg_debug_info("#Interfering tree diagrams", len(tree_diagrams))
         return len(tree_diagrams)
 
+    def drop_orphan_interference_trees(self):
+        """Remove the interfering trees of a LIxtree process (see
+        set_interference_trees) of which no remaining loop diagram forms a
+        product passing the squared-order constraints. The squared-order
+        selection keeps a tree as long as one loop partner exists, but the loop
+        filters applied after it (Furry, wave-function corrections, vanishing
+        tadpoles, perturbative orders, user loop filter) can remove that
+        partner. Returns the number of trees removed."""
+
+        squared_orders = dict((order, value) for order, value in
+                              self['process']['squared_orders'].items()
+                              if value >= 0)
+        sqorders_types = copy.copy(self['process'].get('sqorders_types'))
+        if 'WEIGHTED' not in sqorders_types:
+            sqorders_types['WEIGHTED'] = '<='
+        kept = base_objects.DiagramList()
+        removed = 0
+        for diagram in self['loop_UVCT_diagrams']:
+            if not diagram['UVCT_orders'].get(INTERFERENCE_ORDER) or \
+                   any(diagram.pass_squared_order_constraints(loop,
+                                          squared_orders, sqorders_types)
+                       for loop in self['loop_diagrams']):
+                kept.append(diagram)
+            else:
+                removed += 1
+        self['loop_UVCT_diagrams'] = kept
+        if removed:
+            ldg_debug_info("#Orphan interfering trees removed", removed)
+        return removed
+
     def set_Born_CT(self):
         """ Scan all born diagrams and add for each all the corresponding UV 
         counterterms. It creates one LoopUVCTDiagram per born diagram and set
diff --git a/tests/unit_tests/core/test_tree_interference.py b/tests/unit_tests/core/test_tree_interference.py
index 0ddaec4ee1..73c9d285cb 100644
--- a/tests/unit_tests/core/test_tree_interference.py
+++ b/tests/unit_tests/core/test_tree_interference.py
@@ -106,6 +106,83 @@ def test_multiprocess_without_crossing_reuse(self):
             self.assertEqual(sorted(d.get_order(INTERF)
                                     for d in amp.get('diagrams')), [0, 1])
 
+    def test_left_hand_amplitude_constraints(self):
+        """'==' on the left constrains the left-hand amplitudes only: it does
+        not become a squared-order constraint on the product (QED^2==0 would
+        remove every QCD x EW product here)."""
+        interface = generate('u d > u d QED==0 [treextree] u d > u d QCD==0')
+        amp = interface._curr_amps[0]
+        process = amp.get('process')
+        self.assertEqual(process['squared_orders'], {INTERF: 1})
+        self.assertEqual(process['constrained_orders'], {'QED': (0, '==')})
+        sides = [(d.get_order(INTERF), d.get_order('QED')) for d in amp.get('diagrams')]
+        self.assertTrue(sides)
+        self.assertTrue(all(qed == 0 for side, qed in sides if side == 0))
+        self.assertTrue(all(qed == 2 for side, qed in sides if side == 1))
+        self.assertEqual(set(side for side, _ in sides), set([0, 1]))
+
+    def test_user_diagram_filter_applies_to_both_sides(self):
+        """'--diagram_filter' runs the user filter on the right-hand diagrams
+        too."""
+        import madgraph.core.diagram_generation as diagram_generation
+        filtered = []
+        original = diagram_generation.Amplitude.apply_user_filter
+        def record(amplitude, diag_list):
+            filtered.append(amplitude.get('process').get('required_s_channels'))
+            return diag_list
+        diagram_generation.Amplitude.apply_user_filter = record
+        try:
+            generate('u u~ > z > e+ e- [treextree] u u~ > a > e+ e- '
+                     '--diagram_filter')
+        finally:
+            diagram_generation.Amplitude.apply_user_filter = original
+        self.assertIn([[23]], filtered)
+        self.assertIn([[22]], filtered)
+
+    def test_hidden_order_does_not_leak(self):
+        """Once registered by an interference process, the hidden order
+        neither appears in the diagrams of later ordinary processes nor adds a
+        level to the coupling hierarchy used for the default orders."""
+        interface = MGCmd.MasterCmd()
+        interface.exec_cmd('import model sm', printcmd=False, precmd=True)
+        model = interface._curr_model
+        before = model.get_particles_hierarchy()
+        interface.exec_cmd('generate u u~ > z > e+ e- [treextree] u u~ > a > e+ e-',
+                           printcmd=False, precmd=True, errorhandling=False)
+        self.assertIn(INTERF, model.get('coupling_orders'))
+        self.assertEqual(model.get_particles_hierarchy(), before)
+        interface.exec_cmd('generate u u~ > e+ e-', printcmd=False,
+                           precmd=True, errorhandling=False)
+        for diagram in interface._curr_amps[0].get('diagrams'):
+            self.assertNotIn(INTERF, diagram.get('orders'))
+
+    def test_never_grouped_with_an_ordinary_process(self):
+        """Same '@N', same legs: still two subprocess groups, since the
+        events of a group share the power of alpha_s of each channel."""
+        import madgraph.iolibs.group_subprocs as group_subprocs
+        import madgraph.core.diagram_generation as diagram_generation
+        interface = generate('u u~ > z > e+ e- [treextree] u u~ > a > e+ e- @1')
+        interface.exec_cmd('add process u u~ > z > e+ e- @1', printcmd=False,
+                           precmd=True, errorhandling=False)
+        amps = diagram_generation.AmplitudeList(interface._curr_amps)
+        self.assertEqual(len(amps), 2)
+        classes = group_subprocs.SubProcessGroup.find_process_classes(amps,
+                                                                 'madevent')
+        self.assertNotEqual(classes[0], classes[1])
+
+    def test_alpha_s_power_per_subprocess(self):
+        """One power per matrix element: None for an ordinary process (its
+        channels keep the QCD order of their diagram), the power of the
+        left x right products for an interference."""
+        import madgraph.iolibs.export_v4 as export_v4
+        interface = generate('u d > u d QED=0 [treextree] u d > u d QCD=0')
+        interference = helas_objects.HelasMatrixElement(interface._curr_amps[0])
+        interface.exec_cmd('add process u d > u d', printcmd=False,
+                           precmd=True, errorhandling=False)
+        ordinary = helas_objects.HelasMatrixElement(interface._curr_amps[-1])
+        self.assertEqual(export_v4.interference_alpha_s_powers(
+                                     [ordinary, interference]), [None, 1])
+
     def test_no_interference(self):
         self.assertRaisesRegex(InvalidCmd, 'no interference',
                                generate,
diff --git a/tests/unit_tests/interface/test_interference_syntax.py b/tests/unit_tests/interface/test_interference_syntax.py
index a9db26338e..f4b0fb01c1 100644
--- a/tests/unit_tests/interface/test_interference_syntax.py
+++ b/tests/unit_tests/interface/test_interference_syntax.py
@@ -76,6 +76,19 @@ def test_constraints_go_to_their_side(self):
             ('g g > h > t t~ [noborn= QCD ] QCD^2==6 QED = 2 @3', 'LIxtree',
              'g g > t t~ QED=0 / b'))
 
+    def test_pdg_code_processes_are_not_constraints(self):
+        """An order name is an identifier, and a constraint is followed by a
+        process: '25 > 6' or 'h > 5' starting a right-hand decay is no
+        'ORDER > n' constraint."""
+        self.assertEqual(split('25 > 6 -6 [treextree] 25 > 6 -6'),
+                         ('25 > 6 -6', 'treextree', '25 > 6 -6'))
+        self.assertEqual(split('h > b b~ [treextree] h > 5 -5'),
+                         ('h > b b~', 'treextree', 'h > 5 -5'))
+        self.assertEqual(split('u u~ > e+ e- [treextree] QED>2 u u~ > e+ e-'),
+                         ('u u~ > e+ e- QED>2', 'treextree', 'u u~ > e+ e-'))
+        self.assertEqual(split('g g > h > t t~ [LIxtree=QCD] QCD^2==6'),
+                         ('g g > h > t t~ [noborn= QCD ] QCD^2==6', 'LIxtree', None))
+
     def test_ordinary_lines_are_untouched(self):
         for line in ['g g > z z [noborn=QCD]', 'p p > t t~ [QCD]',
                      'p p > t t~ [real=QCD] QCD^2<=4', 'p p > t t~ QED=0',
@@ -147,6 +160,15 @@ def test_process_defaults(self):
         self.assertIsNone(procdef['interference_process'])
         self.assertEqual(base_objects.Process()['interference_mode'], '')
 
+    def test_process_from_an_older_pickle(self):
+        """A process unpickled from a version without the interference keys
+        still prints."""
+        process = base_objects.Process()
+        for key in ('interference_mode', 'interference_process'):
+            dict.__delitem__(process, key)
+        self.assertNotIn('interference_mode', str(process))
+        self.assertEqual(process.get_interference_mode(), '')
+
     def test_lixtree(self):
         procdef = self.cmd.extract_process(
             'g g > h > t t~ [LIxtree=QCD] QCD^2==6 g g > t t~ QED=0 / b @2')
diff --git a/tests/unit_tests/loop/test_loop_induced_x_tree.py b/tests/unit_tests/loop/test_loop_induced_x_tree.py
index ccb185a2ab..838586557a 100644
--- a/tests/unit_tests/loop/test_loop_induced_x_tree.py
+++ b/tests/unit_tests/loop/test_loop_induced_x_tree.py
@@ -134,8 +134,40 @@ def test_alpha_s_power(self):
         me = loop_helas_objects.LoopHelasMatrixElement(
                      self.interface._curr_amps[0], optimized_output=True)
         self.assertEqual(export_v4.interference_alpha_s_power([me]), 2)
-        self.assertEqual(export_v4.interference_nqcd_list([5, 0, 3], [me]),
-                         [2, 2, 2])
+        self.assertEqual(export_v4.interference_alpha_s_powers([me]), [2])
+
+    def test_orphan_trees_are_dropped(self):
+        """Trees whose loop partners are all removed after the squared-order
+        selection no longer contribute and are dropped."""
+        import copy
+        import madgraph.core.base_objects as base_objects
+        amp = copy.copy(self.interface._curr_amps[0])
+        amp['loop_UVCT_diagrams'] = base_objects.DiagramList(
+                                            amp['loop_UVCT_diagrams'])
+        self.assertEqual(amp.drop_orphan_interference_trees(), 0)
+        self.assertEqual(len(amp['loop_UVCT_diagrams']), 3)
+        amp['loop_diagrams'] = base_objects.DiagramList()
+        self.assertEqual(amp.drop_orphan_interference_trees(), 3)
+        self.assertEqual(len(amp['loop_UVCT_diagrams']), 0)
+
+
+class LoopInducedXTreeDisplayTest(unittest.TestCase):
+
+    def test_subprocess_string_parses_back(self):
+        """The string of a subprocess puts the left-hand exclusions before
+        the bracket, so that they stay on the left-hand side when the string
+        is read back."""
+        interface = generate('g g > h > t t~ / b [LIxtree=QCD] g g > t t~')
+        process = interface._curr_amps[0]['process']
+        printed = process.nice_string(prefix=False)
+        self.assertLess(printed.index('/ b'), printed.index('[ LIxtree'))
+        again = interface.extract_process(printed)
+        self.assertEqual(again['interference_mode'], 'LIxtree')
+        self.assertEqual(again['forbidden_particles'], [5])
+        right = again['interference_process']
+        self.assertEqual(right['forbidden_particles'], [])
+        self.assertEqual([l['ids'] for l in right['legs']],
+                         [[21], [21], [6], [-6]])
 
 
 class LoopInducedXTreeMultiProcessTest(unittest.TestCase):

From dfbd39e870a68b5a5ab569ab88bf66a068c69464 Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Fri, 9 Oct 2026 09:16:33 +0200
Subject: [PATCH 3/5] Interference processes: cleanup from review

- The run card asks the process whether it keeps only some squared orders
  (Process.has_squared_order_selection: an interference, or a squared-order
  constraint not implied by an amplitude '=='/'>'), instead of searching
  '^2' in its printed form; same answer for ordinary processes (tested).
- The squared-order test of interference_alpha_s_power reuses
  Diagram.pass_squared_order_constraints instead of a copy of its rules.
- One place for hiding the internal order in process strings
  (base_objects.visible_orders).
- The interference help is written once, for 'help generate' and 'help add'.
- No mutable default argument in extract_interference_process; unused alias
  removed from loop_diagram_generation.

Full unit suite 2465 OK, acceptance interference tests OK.
---
 madgraph/core/base_objects.py                 | 45 +++++++++++++------
 madgraph/interface/madgraph_interface.py      | 26 ++++-------
 madgraph/iolibs/export_v4.py                  | 30 +++++--------
 madgraph/loop/loop_diagram_generation.py      |  1 -
 madgraph/various/banner.py                    | 14 +++---
 .../interface/test_interference_syntax.py     | 19 ++++++++
 6 files changed, 77 insertions(+), 58 deletions(-)

diff --git a/madgraph/core/base_objects.py b/madgraph/core/base_objects.py
index 621195cb4e..9e88ac1490 100755
--- a/madgraph/core/base_objects.py
+++ b/madgraph/core/base_objects.py
@@ -42,6 +42,11 @@
 # added to the model in memory with hierarchy 0, so that WEIGHTED is
 # unchanged, and is never printed in a process definition.
 INTERFERENCE_ORDER = 'INTERF'
+
+def visible_orders(orders):
+    """The keys of a dictionary of (squared) coupling orders that a process
+    description shows: all but the hidden INTERFERENCE_ORDER."""
+    return [key for key in orders if key != INTERFERENCE_ORDER]
 if madgraph.ordering:
     set = misc.OrderedSet
 
@@ -3954,6 +3959,25 @@ def get_interference_mode(self):
 
         return dict.get(self, 'interference_mode', '') or ''
 
+    def has_squared_order_selection(self):
+        """Whether only some squared orders of this process are kept: an
+        interference process ('[LIxtree=QCD]', '[treextree]'), or a
+        squared-order constraint other than the one implied by an amplitude
+        '==' / '>' constraint, here or in a decay chain. This is exactly when
+        nice_string shows a '^2' (or would, but for the hidden order), and it
+        is what the run card treats as an interference."""
+
+        if self.get_interference_mode():
+            return True
+        for key, value in self['squared_orders'].items():
+            constrained = self['constrained_orders'].get(key)
+            if constrained and constrained[0] == value / 2 and \
+                     constrained[1] == self.get_squared_order_type(key):
+                continue
+            return True
+        return any(decay.has_squared_order_selection()
+                   for decay in self['decay_chains'])
+
     def interference_string(self):
         """The tail of the process syntax for an interference process: the
         '[ treextree ]' or '[ LIxtree = ORDERS ]' bracket followed by the
@@ -4030,11 +4054,9 @@ def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturba
         # Add orders
         if self['orders']:
             to_add = []
-            for key in sorted(self['orders'].keys()):
+            for key in sorted(visible_orders(self['orders'])):
                 if not print_weighted and key == 'WEIGHTED':
                     continue
-                if key == INTERFERENCE_ORDER:
-                    continue
                 value = int(self['orders'][key])
                 if key in self['squared_orders']:
                     if self.get_squared_order_type(key) in ['<=', '==', '='] and \
@@ -4076,11 +4098,9 @@ def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturba
         # Add squared orders
         if self['squared_orders']:
             to_add = []
-            for key in sorted(self['squared_orders'].keys()):
+            for key in sorted(visible_orders(self['squared_orders'])):
                 if not print_weighted and key == 'WEIGHTED':
                     continue
-                if key == INTERFERENCE_ORDER:
-                    continue
                 if key in self['constrained_orders']:
                     if self['constrained_orders'][key][0] == self['squared_orders'][key]/2 and \
                        self['constrained_orders'][key][1] == self.get_squared_order_type(key):
@@ -4172,7 +4192,7 @@ def input_string(self):
             prevleg = leg
 
         if self['orders']:
-            keys = [k for k in self['orders'] if k != INTERFERENCE_ORDER]
+            keys = visible_orders(self['orders'])
             keys.sort(reverse=True)
             mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
                        for key in keys]) + ' '
@@ -4180,8 +4200,7 @@ def input_string(self):
         # Add squared orders
         if self['squared_orders']:
             mystr = mystr + " ".join([key + '^2=' + repr(self['squared_orders'][key]) \
-                       for key in self['squared_orders']
-                       if key != INTERFERENCE_ORDER]) + ' '
+                       for key in visible_orders(self['squared_orders'])]) + ' '
 
         # Add perturbation orders
         if self['perturbation_couplings']:
@@ -4992,8 +5011,7 @@ def nice_string(self, indent=0, print_weighted=False, prefix=True):
 
         if self['orders']:
             mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
-                       for key in sorted(self['orders'])
-                       if key != INTERFERENCE_ORDER]) + ' '
+                       for key in sorted(visible_orders(self['orders']))]) + ' '
 
         if self['constrained_orders']:
             mystr = mystr + " ".join('%s%s%d' % (key, operator, value) for 
@@ -5017,9 +5035,8 @@ def nice_string(self, indent=0, print_weighted=False, prefix=True):
         if self['squared_orders']:
             mystr = mystr + " ".join([key + '^2%s%d'%\
                 (self.get_squared_order_type(key),self['squared_orders'][key]) \
-              for key in self['squared_orders'].keys() \
-                    if (print_weighted or key!='WEIGHTED') and \
-                                            key != INTERFERENCE_ORDER]) + ' '
+              for key in visible_orders(self['squared_orders']) \
+                                    if print_weighted or key!='WEIGHTED']) + ' '
 
         mystr = mystr + self.interference_string()
 
diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py
index 5781be38c1..e6cb32b0db 100755
--- a/madgraph/interface/madgraph_interface.py
+++ b/madgraph/interface/madgraph_interface.py
@@ -939,6 +939,10 @@ def help_generate(self):
         logger.info(" > For processes without born amplitudes (i.e. loop-induced like g g > z), please use ")
         logger.info("   the 'virt=' NLO mode. aMC@NLO cannot integrate these processes, but standalone MadLoop5")
         logger.info("   can still handle these.")
+        self._log_interference_syntax()
+
+    def _log_interference_syntax(self):
+        """The interference syntax, for 'help generate' and 'help add'."""
         logger.info("Interference syntax:",'$MG:BOLD')
         logger.info(" o left process [LIxtree=LoopOrders] SQUAREDCOUPi=ORDERi right process")
         logger.info(" o left process [treextree] SQUAREDCOUPi=ORDERi right process")
@@ -996,21 +1000,7 @@ def help_add(self):
         logger.info(" > For processes without born amplitudes (i.e. loop-induced like g g > z), please use ")
         logger.info("   the 'virt=' NLO mode. aMC@NLO cannot integrate these processes, but standalone MadLoop5")
         logger.info("   can still handle these.")
-        logger.info("Interference syntax:",'$MG:BOLD')
-        logger.info(" o left process [LIxtree=LoopOrders] SQUAREDCOUPi=ORDERi right process")
-        logger.info(" o left process [treextree] SQUAREDCOUPi=ORDERi right process")
-        logger.info(" o Example: generate g g > h > t t~ [LIxtree=QCD] g g > t t~",'$MG:color:GREEN')
-        logger.info(" o Example: generate p p > z > e+ e- [treextree] p p > a > e+ e-",'$MG:color:GREEN')
-        logger.info(" > Only the interference 2 Re(A_left A_right^*) is computed. The left process is")
-        logger.info("   a loop-induced process for LIxtree, a tree-level one for treextree; the right")
-        logger.info("   process is always tree level. The keywords are case-insensitive.")
-        logger.info(" > Each side takes its own s-channel requirements, exclusions and amplitude orders.")
-        logger.info("   Squared-order constraints apply to the interference: give them right after ']'.")
-        logger.info(" > Both sides must have the same external legs, in the same order (the right-hand")
-        logger.info("   legs may be wider multiparticles). For LIxtree the right process is optional:")
-        logger.info("   by default it is the tree-level process with the same legs.")
-        logger.info(" > 'output madevent' integrates the interference (events with signed weights);")
-        logger.info("   for LIxtree, 'output standalone_fortran' evaluates it at a phase-space point.")
+        self._log_interference_syntax()
 
         logger.info("--  merge two model to create a new one", '$MG:color:BLUE')
         logger.info("syntax:",'$MG:BOLD')
@@ -6891,7 +6881,7 @@ def extract_process(self, line, proc_number = 0, overall_orders = {},
 
 
     def extract_interference_process(self, left, mode, right, proc_number=0,
-                                     overall_orders={},
+                                     overall_orders=None,
                                      avoid_squared_orders=False):
         """Build the ProcessDefinition of an interference process, as split by
         split_interference_line: 'left' is the left-hand process in the
@@ -6915,8 +6905,8 @@ def extract_interference_process(self, left, mode, right, proc_number=0,
         # [treextree] u d > u d QCD==0 would get QED^2==0 and no interference
         # at all). Explicit '^2' constraints are kept.
         procdef = self.extract_process(left, proc_number=proc_number,
-                                       overall_orders=overall_orders,
-                                       avoid_squared_orders=True)
+                    overall_orders=overall_orders if overall_orders else {},
+                    avoid_squared_orders=True)
         if not procdef:
             raise self.InvalidCmd("Empty or wrong format process, please try again.")
         if procdef.get('decay_chains'):
diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py
index e7cfa11a86..e7e6472ea8 100755
--- a/madgraph/iolibs/export_v4.py
+++ b/madgraph/iolibs/export_v4.py
@@ -13886,32 +13886,24 @@ def interference_alpha_s_power(matrix_elements):
                       amplitude.get('loop_UVCT_diagrams')]
         else:
             orders = [d.calculate_orders() for d in me.get('diagrams')]
+        # one diagram per distinct set of orders, on each side
         sides = {0: [], 1: []}
         for signature in set(tuple(sorted((k, v) for k, v in o.items()
                                           if k != 'WEIGHTED')) for o in orders):
-            side = dict(signature).get(base_objects.INTERFERENCE_ORDER, 0)
+            diagram_orders = dict(signature)
+            diagram_orders['WEIGHTED'] = sum(hierarchy.get(k, 0) * v
+                                             for k, v in signature)
+            side = diagram_orders.get(base_objects.INTERFERENCE_ORDER, 0)
             if side in sides:
-                sides[side].append(dict(signature))
+                sides[side].append(base_objects.Diagram({'orders': diagram_orders}))
         squared_orders = process.get('squared_orders')
+        types = dict((order, process.get_squared_order_type(order))
+                     for order in squared_orders)
         for left in sides[0]:
             for right in sides[1]:
-                product = dict((k, left.get(k, 0) + right.get(k, 0))
-                               for k in set(left) | set(right))
-                product['WEIGHTED'] = sum(hierarchy.get(k, 0) * v
-                                          for k, v in product.items())
-                passed = True
-                for order, value in squared_orders.items():
-                    if value < 0:
-                        continue
-                    kind = process.get_squared_order_type(order)
-                    combined = product.get(order, 0)
-                    if (kind == '==' and combined != value) or \
-                       (kind in ['=', '<='] and combined > value) or \
-                       (kind == '>' and combined <= value):
-                        passed = False
-                        break
-                if passed:
-                    powers.add(product.get('QCD', 0))
+                if left.pass_squared_order_constraints(right, squared_orders,
+                                                       types):
+                    powers.add(left.get_order('QCD') + right.get_order('QCD'))
     if len(powers) == 1:
         power = powers.pop()
         if power % 2 == 0:
diff --git a/madgraph/loop/loop_diagram_generation.py b/madgraph/loop/loop_diagram_generation.py
index a27f7bcd0d..b54e8b391a 100755
--- a/madgraph/loop/loop_diagram_generation.py
+++ b/madgraph/loop/loop_diagram_generation.py
@@ -36,7 +36,6 @@
 #===============================================================================
 INTERFERENCE_ORDER = base_objects.INTERFERENCE_ORDER
 register_interference_order = diagram_generation.register_interference_order
-prepare_interference_process = diagram_generation.prepare_interference_process
 
 def ldg_debug_info(msg,val, force=False):
     # This subroutine has typically quite large DEBUG info.
diff --git a/madgraph/various/banner.py b/madgraph/various/banner.py
index 878f235a55..2057159201 100755
--- a/madgraph/various/banner.py
+++ b/madgraph/various/banner.py
@@ -5269,12 +5269,14 @@ def create_default_for_process(self, proc_characteristic, history, proc_def):
         interference = False
         for proc in proc_def:
             for oneproc in proc:
-                # a squared-order constraint, or an interference process
-                # ('left [LIxtree=QCD] right' / '[treextree]'), whose hidden
-                # squared-order constraint is not printed
-                if '^2' in oneproc.nice_string() or \
-                        (hasattr(oneproc, 'get_interference_mode') and
-                                            oneproc.get_interference_mode()):
+                # a squared-order selection: an explicit squared-order
+                # constraint or an interference process ('[LIxtree=QCD]',
+                # '[treextree]', whose constraint is hidden)
+                if hasattr(oneproc, 'has_squared_order_selection'):
+                    selection = oneproc.has_squared_order_selection()
+                else:
+                    selection = '^2' in oneproc.nice_string()
+                if selection:
                     interference = True
                     break
             else:
diff --git a/tests/unit_tests/interface/test_interference_syntax.py b/tests/unit_tests/interface/test_interference_syntax.py
index f4b0fb01c1..689fc359ae 100644
--- a/tests/unit_tests/interface/test_interference_syntax.py
+++ b/tests/unit_tests/interface/test_interference_syntax.py
@@ -160,6 +160,25 @@ def test_process_defaults(self):
         self.assertIsNone(procdef['interference_process'])
         self.assertEqual(base_objects.Process()['interference_mode'], '')
 
+    def test_squared_order_selection(self):
+        """What the run card treats as an interference: for ordinary
+        processes exactly when the process string shows a '^2' (the squared
+        order implied by an amplitude '==' is not one), and every interference
+        process, whose constraint is hidden."""
+        for line, expected in [('u u~ > e+ e-', False),
+                               ('p p > j j QED==2', False),
+                               ('p p > j j QCD^2==2', True),
+                               ('p p > j j QED^2<=2', True),
+                               ('u u~ > z > e+ e- [treextree] u u~ > a > e+ e-', True)]:
+            procdef = self.cmd.extract_process(line)
+            if procdef.get_interference_mode():
+                import madgraph.core.diagram_generation as diagram_generation
+                diagram_generation.prepare_interference_process(procdef)
+            process = next(iter(procdef))
+            self.assertEqual(process.has_squared_order_selection(), expected, line)
+            if not procdef.get_interference_mode():
+                self.assertEqual('^2' in process.nice_string(), expected, line)
+
     def test_process_from_an_older_pickle(self):
         """A process unpickled from a version without the interference keys
         still prints."""

From 95d265d2515fa39f64198bf5b4cd59dec16dd2ff Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Fri, 9 Oct 2026 10:24:06 +0200
Subject: [PATCH 4/5] Interference processes: never print the hidden order
 INTERF

The hidden order INTERFERENCE_ORDER (hierarchy 0), registered in the model
by an interference process, still leaked into:
- the WEIGHTED definition of the minimal-order search ("WEIGTHED IS
  QCD+2*QED+INTERF", which also reads as weight 1): orders of hierarchy 0
  are skipped, and the WEIGTHED typo is fixed;
- 'display coupling_order' (" INTERF : weight = 0");
- the orders set by default_unset_couplings (INTERF<=N was imposed);
- the debug summary of the loop squared-order selection (an INTERF column);
- a later NLO process of the same session: squared-order warning, split
  orders and perturbation couplings ("[ all = INTERF QCD QED ]").
---
 madgraph/core/diagram_generation.py           |  8 +++--
 madgraph/fks/fks_base.py                      |  2 +-
 madgraph/interface/amcatnlo_interface.py      | 10 ++++--
 madgraph/interface/madgraph_interface.py      |  7 ++--
 madgraph/loop/loop_diagram_generation.py      | 23 +++++++-----
 .../unit_tests/core/test_tree_interference.py | 35 +++++++++++++++++++
 .../loop/test_loop_induced_x_tree.py          | 30 ++++++++++++++++
 7 files changed, 98 insertions(+), 17 deletions(-)

diff --git a/madgraph/core/diagram_generation.py b/madgraph/core/diagram_generation.py
index 97a4f066e2..a08cac5455 100755
--- a/madgraph/core/diagram_generation.py
+++ b/madgraph/core/diagram_generation.py
@@ -2989,7 +2989,11 @@ def find_optimal_process_orders(process_definition, diagram_filter=False):
         hierarchy = list(hierarchydef.items())
         hierarchy.sort()
         for key, value in hierarchydef.items():
-            if value>1:
+            # an order of hierarchy 0 (like the hidden INTERFERENCE_ORDER)
+            # does not enter WEIGHTED
+            if value == 0:
+                continue
+            elif value>1:
                 tmp.append('%s*%s' % (value,key))
             else:
                 tmp.append('%s' % key)
@@ -2997,7 +3001,7 @@ def find_optimal_process_orders(process_definition, diagram_filter=False):
         # Run diagram generation with increasing max_order_now until
         # we manage to get diagrams
         while max_order_now < max_WEIGHTED_order:
-            logger.info("Trying coupling order WEIGHTED<=%d: WEIGTHED IS %s" % (max_order_now, wgtdef))
+            logger.info("Trying coupling order WEIGHTED<=%d: WEIGHTED IS %s" % (max_order_now, wgtdef))
 
             oldloglevel = logger.level
             logger.setLevel(logging.WARNING)
diff --git a/madgraph/fks/fks_base.py b/madgraph/fks/fks_base.py
index 6041ca5bd1..d290b054e9 100755
--- a/madgraph/fks/fks_base.py
+++ b/madgraph/fks/fks_base.py
@@ -920,7 +920,7 @@ def find_reals(self, pert_orders = []):
         # are imposed)
         if not pert_orders:
             if not self.born_amp['process']['orders']:
-                pert_orders = model['coupling_orders']
+                pert_orders = MG.visible_orders(model['coupling_orders'])
             else:
                 pert_orders = self.born_amp['process']['perturbation_couplings']
 
diff --git a/madgraph/interface/amcatnlo_interface.py b/madgraph/interface/amcatnlo_interface.py
index 7f1eeac939..66295500f8 100755
--- a/madgraph/interface/amcatnlo_interface.py
+++ b/madgraph/interface/amcatnlo_interface.py
@@ -52,6 +52,7 @@
 import madgraph.iolibs.export_v4 as export_v4
 import madgraph.iolibs.helas_call_writers as helas_call_writers
 import madgraph.loop.loop_base_objects as loop_base_objects
+import madgraph.core.base_objects as base_objects
 import madgraph.core.diagram_generation as diagram_generation
 import madgraph.core.helas_objects as helas_objects
 
@@ -752,7 +753,10 @@ def do_add(self, line, *args,**opt):
 
         # check that all the couplings of the model have been constrained
         # in the squared orders, otherwise set the others to zero
-        for o in myprocdef['model'].get('coupling_orders'):
+        # (not the hidden order an interference process adds to the model)
+        model_orders = base_objects.visible_orders(
+                                     myprocdef['model'].get('coupling_orders'))
+        for o in model_orders:
             if o not in myprocdef['squared_orders'].keys():
                 logger.warning('No squared order constraint for order %s. Setting to 0' % o)
                 myprocdef['squared_orders'][o] = 0 
@@ -761,7 +765,7 @@ def do_add(self, line, *args,**opt):
         # split all orders in the model, for the moment it's the simplest solution
         # mz02/2014
         #if proc_type[1] != 'only':
-        myprocdef['split_orders'] += [o for o in myprocdef['model'].get('coupling_orders') \
+        myprocdef['split_orders'] += [o for o in model_orders \
                 if o not in myprocdef['split_orders']]
 
         # now set the squared orders
@@ -809,7 +813,7 @@ def do_add(self, line, *args,**opt):
         # squared-orders constraints. In that case, all kind of splittings/loop-particles
         # must be included
         if not myprocdef['orders'] and self.options['nlo_mixed_expansion']:
-            myprocdef['perturbation_couplings'] = list(myprocdef['model']['coupling_orders'])
+            myprocdef['perturbation_couplings'] = list(model_orders)
 
         self._curr_proc_defs.append(myprocdef)
 
diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py
index e6cb32b0db..f287ed0599 100755
--- a/madgraph/interface/madgraph_interface.py
+++ b/madgraph/interface/madgraph_interface.py
@@ -4555,7 +4555,9 @@ def key_sort(x):
                 print(self.boundstate_string(key))
 
         elif args[0] == 'coupling_order':
-            hierarchy = list(self._curr_model['order_hierarchy'].items())
+            hierarchy = [(order, weight) for order, weight in
+                         self._curr_model['order_hierarchy'].items()
+                         if order != base_objects.INTERFERENCE_ORDER]
             hierarchy.sort(key=operator.itemgetter(1))
             # an order declared by the model can have no interaction left
             # carrying it -- a restriction card typically removes all of them.
@@ -6278,7 +6280,8 @@ def extract_process(self, line, proc_number = 0, overall_orders = {},
         if self.options['default_unset_couplings'] != 99 and \
                                                      (orders or squared_orders): 
                            
-                to_set = [name for name in self._curr_model.get('coupling_orders')
+                to_set = [name for name in base_objects.visible_orders(
+                                       self._curr_model.get('coupling_orders'))
                           if name not in orders and name not in squared_orders]
                 if to_set:
                     logger.info('the following coupling will be allowed up to the maximal value of %s: %s' % 
diff --git a/madgraph/loop/loop_diagram_generation.py b/madgraph/loop/loop_diagram_generation.py
index b54e8b391a..91bbf68359 100755
--- a/madgraph/loop/loop_diagram_generation.py
+++ b/madgraph/loop/loop_diagram_generation.py
@@ -1056,9 +1056,16 @@ def print_split_order_infos(self):
                                                            key=lambda el: el[1])
         
         
+        # the hidden INTERFERENCE_ORDER is not printed
+        shown = [i for i, order in enumerate(sorted_hierarchy[:-1])
+                                               if order != INTERFERENCE_ORDER]
+        def SO_string(SO):
+            return '(%s,W%d)' % (','.join('%d' % SO[i] for i in shown), SO[-1])
+
         logger.debug("Coupling order combinations considered:"+\
-                                            " (%s)"%','.join(sorted_hierarchy))
-        
+                    " (%s)"%','.join([sorted_hierarchy[i] for i in shown] +
+                                                                 ['WEIGHTED']))
+
         # Now check what is left
         born_considered = []
         loop_considered = []
@@ -1096,14 +1103,12 @@ def print_split_order_infos(self):
             if len(considered)==0:
                 logger.debug(" > %s : None"%name)
             else:
-                logger.debug(" > %s : %s"%(name,' '.join(['(%s,W%d)'%(
-                            ','.join(list('%d'%s for s in c[:-1])),c[-1]) 
-                                                         for c in considered])))
-            
+                logger.debug(" > %s : %s"%(name,' '.join(
+                                        [SO_string(c) for c in considered])))
+
             if len(extra)!=0:
-                logger.debug(" > %s (not selected but available): %s"%(name,' '.
-                    join(['(%s,W%d)'%(','.join(list('%d'%s for s in e[:-1])),
-                                                       e[-1]) for e in extra])))
+                logger.debug(" > %s (not selected but available): %s"%(name,
+                                    ' '.join([SO_string(e) for e in extra])))
                 
         # In case it is needed, the considered orders are returned 
         # (it is used by some of the unit tests)
diff --git a/tests/unit_tests/core/test_tree_interference.py b/tests/unit_tests/core/test_tree_interference.py
index 73c9d285cb..90df5e30b8 100644
--- a/tests/unit_tests/core/test_tree_interference.py
+++ b/tests/unit_tests/core/test_tree_interference.py
@@ -26,6 +26,8 @@
 
 from __future__ import absolute_import
 
+import contextlib
+import io
 import os
 import re
 import shutil
@@ -156,6 +158,39 @@ def test_hidden_order_does_not_leak(self):
         for diagram in interface._curr_amps[0].get('diagrams'):
             self.assertNotIn(INTERF, diagram.get('orders'))
 
+    def test_hidden_order_is_never_printed(self):
+        """Registered in the model by an interference process, the hidden
+        order (hierarchy 0) is neither part of the WEIGHTED definition printed
+        by the minimal-order search, nor listed by 'display coupling_order',
+        nor one of the orders that default_unset_couplings sets."""
+        interface = MGCmd.MasterCmd()
+        interface.exec_cmd('import model sm', printcmd=False, precmd=True)
+        # (the left-hand process goes through WEIGHTED<=2 and 3 before 4)
+        with self.assertLogs('madgraph.diagram_generation', level='INFO') as log:
+            interface.exec_cmd('generate u d > u d [treextree] u d > u d QCD=0',
+                               printcmd=False, precmd=True, errorhandling=False)
+        self.assertEqual(interface._curr_model['order_hierarchy'][INTERF], 0)
+        searches = [line for line in log.output
+                    if 'Trying coupling order WEIGHTED' in line]
+        self.assertTrue(searches)
+        for line in searches:
+            self.assertEqual(sorted(line.rsplit(' IS ', 1)[1].split('+')),
+                             ['2*QED', 'QCD'])
+            self.assertIn('WEIGHTED IS ', line)
+
+        stdout = io.StringIO()
+        with contextlib.redirect_stdout(stdout):
+            interface.exec_cmd('display coupling_order', printcmd=False,
+                               precmd=True, errorhandling=False)
+        self.assertIn('QED : weight = 2', stdout.getvalue())
+        self.assertNotIn(INTERF, stdout.getvalue())
+
+        interface.options['default_unset_couplings'] = 0
+        interface.exec_cmd('generate u u~ > e+ e- QED=2', printcmd=False,
+                           precmd=True, errorhandling=False)
+        process = interface._curr_amps[0].get('process')
+        self.assertEqual(process.get('orders'), {'QED': 2, 'QCD': 0})
+
     def test_never_grouped_with_an_ordinary_process(self):
         """Same '@N', same legs: still two subprocess groups, since the
         events of a group share the power of alpha_s of each channel."""
diff --git a/tests/unit_tests/loop/test_loop_induced_x_tree.py b/tests/unit_tests/loop/test_loop_induced_x_tree.py
index 838586557a..7eb9690fef 100644
--- a/tests/unit_tests/loop/test_loop_induced_x_tree.py
+++ b/tests/unit_tests/loop/test_loop_induced_x_tree.py
@@ -150,6 +150,22 @@ def test_orphan_trees_are_dropped(self):
         self.assertEqual(amp.drop_orphan_interference_trees(), 3)
         self.assertEqual(len(amp['loop_UVCT_diagrams']), 0)
 
+    def test_split_order_infos_do_not_print_the_hidden_order(self):
+        """The debug summary of the squared-order selection leaves the hidden
+        order out: loop x tree (INTERF=1) is selected, loop x loop (0) and
+        tree x tree (2) are not. The returned combinations keep it."""
+        amp = self.interface._curr_amps[0]
+        with self.assertLogs('madgraph.loop_diagram_generation',
+                             level='DEBUG') as log:
+            infos = amp.print_split_order_infos()
+        text = '\n'.join(log.output)
+        self.assertNotIn(INTERF, text)
+        self.assertIn('considered: (QCD,QED,WEIGHTED)', text)
+        self.assertIn(' > loop : (4,2,W8)', text)
+        self.assertIn('(4,4,W12)', text)
+        self.assertIn('(4,0,W4)', text)
+        self.assertEqual(infos[2], [(1, 4, 2, 8)])
+
 
 class LoopInducedXTreeDisplayTest(unittest.TestCase):
 
@@ -199,6 +215,20 @@ def test_no_tree_level_interference(self):
         self.assertRaisesRegex(InvalidCmd, 'no interference',
                                generate, 'g g > h > z z [LIxtree=QCD]')
 
+    def test_later_nlo_process_ignores_the_hidden_order(self):
+        """The hidden order the interference process added to the model is
+        not one of the orders of a later NLO process: neither constrained,
+        split nor perturbed."""
+        interface = generate('g g > h > t t~ [LIxtree=QCD]')
+        interface.exec_cmd('generate u u~ > e+ e- [QCD]', printcmd=False,
+                           precmd=True, errorhandling=False)
+        procdef = interface._curr_proc_defs[-1]
+        self.assertEqual(sorted(procdef['perturbation_couplings']),
+                         ['QCD', 'QED'])
+        self.assertNotIn(INTERF, procdef['split_orders'])
+        self.assertNotIn(INTERF, procdef['squared_orders'])
+        self.assertNotIn(INTERF, procdef.nice_string())
+
 
 if __name__ == '__main__':
     unittest.main()

From 68ca0109187601243ee847ba54f76e9b18429212 Mon Sep 17 00:00:00 2001
From: Olivier Mattelaer <olivier.mattelaer@uclouvain.be>
Date: Sat, 10 Oct 2026 23:40:22 +0200
Subject: [PATCH 5/5] Interference processes: [treextree] with output mg7 and
 the C++ standalone

Since PR #244 (now in main) the madmatrix color sum pairs the amplitude split
orders on every backend, CPU and GPU, and madspace integrates and unweights
the signed weights of an interference. Allow 'output' (mg7) and 'output
standalone' for treextree; madevent and standalone_fortran stay, and LIxtree
still needs a loop backend (madevent, standalone_fortran). The madmatrix
message announcing CPU-only builds comes corrected with the merge of main.

Checked:
- u u~ > z > e+ e- [treextree] u u~ > a > e+ e-, C++ standalone built in
  double precision: equal to |Z+a|^2 - |Z|^2 - |a|^2 from the same code to
  2e-16, and to the Fortran value to 1e-15 (new acceptance test; the default
  mixed-precision build agrees to 7e-8 only);
- mg7, p p > e+ e-, m(ee) > 200 GeV: interference 0.0650(9) pb over four
  seeds, |Z+a|^2 - |Z|^2 - |a|^2 0.0649(9) pb; the Z+a, Z and a cross
  sections agree with madevent. madevent's interference comes out 7.5%
  lower (0.0601(17) pb, 2.6 sigma), its known iteration-combination bias on
  sign-changing integrands; same matrix element (agrees to 1e-15).

Full unit suite 2499 OK, interference acceptance tests OK.
---
 madgraph/core/diagram_generation.py           |  8 +++++--
 madgraph/interface/madgraph_interface.py      |  6 +++--
 .../test_loop_induced_x_tree.py               | 22 ++++++++++++++++---
 .../unit_tests/core/test_tree_interference.py | 21 ++++++++++++++++--
 4 files changed, 48 insertions(+), 9 deletions(-)

diff --git a/madgraph/core/diagram_generation.py b/madgraph/core/diagram_generation.py
index a08cac5455..e65aa07fcc 100755
--- a/madgraph/core/diagram_generation.py
+++ b/madgraph/core/diagram_generation.py
@@ -47,8 +47,12 @@ class NoDiagramException(InvalidCmd): pass
 #===============================================================================
 INTERFERENCE_ORDER = base_objects.INTERFERENCE_ORDER
 # Output formats able to select squared split orders, hence to return the
-# interference alone (the others would give |A_left + A_right|^2)
-TREE_INTERFERENCE_FORMATS = ['madevent', 'standalone_fortran']
+# interference alone (the others would give |A_left + A_right|^2): the Fortran
+# ones (GET_MATRIX) and the madmatrix ones (mg7, and standalone for the C++
+# matrix element), whose color sum pairs the amplitude split orders on every
+# backend, CPU and GPU (PR #244)
+TREE_INTERFERENCE_FORMATS = ['madevent', 'standalone_fortran', 'mg7',
+                             'standalone']
 
 def register_interference_order(model):
     """Add the hidden coupling order tagging the right-hand amplitudes of an
diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py
index f287ed0599..7ae6c0a5b7 100755
--- a/madgraph/interface/madgraph_interface.py
+++ b/madgraph/interface/madgraph_interface.py
@@ -956,8 +956,10 @@ def _log_interference_syntax(self):
         logger.info(" > Both sides must have the same external legs, in the same order (the right-hand")
         logger.info("   legs may be wider multiparticles). For LIxtree the right process is optional:")
         logger.info("   by default it is the tree-level process with the same legs.")
-        logger.info(" > 'output madevent' integrates the interference (events with signed weights);")
-        logger.info("   for LIxtree, 'output standalone_fortran' evaluates it at a phase-space point.")
+        logger.info(" > 'output' (mg7) and 'output madevent' integrate the interference (events with")
+        logger.info("   signed weights); 'output standalone' and 'output standalone_fortran' evaluate")
+        logger.info("   it at a phase-space point. LIxtree needs a loop backend: 'output madevent' or")
+        logger.info("   'output standalone_fortran' only.")
 
     def help_add(self):
         logger.info("-- generate diagrams for a process and add to existing processes",'$MG:color:BLUE')
diff --git a/tests/acceptance_tests/test_loop_induced_x_tree.py b/tests/acceptance_tests/test_loop_induced_x_tree.py
index c7ca8112c4..e306eafdff 100644
--- a/tests/acceptance_tests/test_loop_induced_x_tree.py
+++ b/tests/acceptance_tests/test_loop_induced_x_tree.py
@@ -126,21 +126,29 @@ def setUp(self):
     def tearDown(self):
         shutil.rmtree(self.tmpdir, ignore_errors=True)
 
+    output_format = 'standalone_fortran'
+
     def evaluate(self, name, process):
         """(momenta, matrix element) printed by the standalone check."""
         out_dir = pjoin(self.tmpdir, name)
         interface = MGCmd.MasterCmd()
         for command in ['import model sm', 'generate %s' % process,
-                        'output standalone_fortran %s -f' % out_dir]:
+                        'output %s %s -f' % (self.output_format, out_dir)]:
             interface.exec_cmd(command, printcmd=False, precmd=True,
                                errorhandling=False)
         sub_dir = pjoin(out_dir, 'SubProcesses')
         proc_dir = [pjoin(sub_dir, d) for d in os.listdir(sub_dir)
                     if d.startswith('P')][0]
-        subprocess.check_call(['make', 'check'], cwd=proc_dir,
+        if self.output_format == 'standalone_fortran':
+            build, check = ['make', 'check'], ['./check']
+        else:
+            # madmatrix: the default build mixes precisions (float color
+            # sum), which caps the comparison at ~1e-7; build it in double
+            build, check = ['make', 'FPTYPE=d', '-j4'], ['./check_sa.exe']
+        subprocess.check_call(build, cwd=proc_dir,
                               stdout=subprocess.DEVNULL,
                               stderr=subprocess.DEVNULL)
-        output = subprocess.run(['./check'], cwd=proc_dir, timeout=300,
+        output = subprocess.run(check, cwd=proc_dir, timeout=300,
                                 capture_output=True, text=True).stdout
         momenta = re.findall(r'^\s*[1-4]\s+(\S+\s+\S+\s+\S+\s+\S+)', output,
                              re.MULTILINE)
@@ -162,5 +170,13 @@ def test_tree_x_tree_standalone(self):
         self.assertAlmostEqual(interference / expected, 1., places=10)
 
 
+
+class TreeXTreeMadmatrixStandaloneTest(TreeXTreeStandaloneTest):
+    """The same comparison with the C++ matrix element of madmatrix ('output
+    standalone', the code mg7 runs), built in double precision."""
+
+    output_format = 'standalone'
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/unit_tests/core/test_tree_interference.py b/tests/unit_tests/core/test_tree_interference.py
index 90df5e30b8..972cf8f70f 100644
--- a/tests/unit_tests/core/test_tree_interference.py
+++ b/tests/unit_tests/core/test_tree_interference.py
@@ -226,13 +226,30 @@ def test_no_interference(self):
 
 class TreeInterferenceOutputTest(unittest.TestCase):
     """Only the formats selecting squared split orders return the
-    interference alone: the others are refused before anything is written."""
+    interference alone (madevent, standalone_fortran, mg7, standalone): the
+    others are refused before anything is written."""
+
+    def test_madmatrix_formats_are_accepted(self):
+        """mg7 and the C++ standalone: the madmatrix color sum pairs the
+        amplitude split orders (CPU and GPU backends)."""
+        interface = generate('u u~ > z > e+ e- [treextree] u u~ > a > e+ e-')
+        tmpdir = tempfile.mkdtemp(prefix='treextree_')
+        try:
+            for fmt in ['mg7', 'standalone']:
+                out = os.path.join(tmpdir, fmt)
+                interface.exec_cmd('output %s %s -f' % (fmt, out),
+                                   printcmd=False, precmd=True,
+                                   errorhandling=False)
+                self.assertTrue(os.path.isdir(os.path.join(out, 'SubProcesses')),
+                                fmt)
+        finally:
+            shutil.rmtree(tmpdir, ignore_errors=True)
 
     def test_unsupported_formats_are_refused(self):
         interface = generate('u u~ > z > e+ e- [treextree] u u~ > a > e+ e-')
         tmpdir = tempfile.mkdtemp(prefix='treextree_')
         try:
-            for fmt in ['mg7', 'standalone', 'matchbox']:
+            for fmt in ['matchbox', 'matrix', 'mg7_v5']:
                 out = os.path.join(tmpdir, fmt)
                 self.assertRaisesRegex(InvalidCmd, 'cannot select the interference',
                                        interface.exec_cmd,