diff --git a/madgraph/core/base_objects.py b/madgraph/core/base_objects.py
index 94a0229c2..9e88ac149 100755
--- a/madgraph/core/base_objects.py
+++ b/madgraph/core/base_objects.py
@@ -34,6 +34,19 @@
 
 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'
+
+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
 
@@ -1956,8 +1969,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 \
@@ -1968,7 +1985,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
@@ -3389,6 +3406,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,
@@ -3738,6 +3759,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 +3855,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 +3949,55 @@ 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'] + \
+               [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
+        a version that did not know about interference processes)."""
+
+        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
+        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 ''
+        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) + ' '
+        return mystr
 
     def nice_string(self, indent=0, print_weighted=True, prefix=True, print_perturbated=True):
         """Returns a nicely formated string about current process
@@ -3970,7 +4054,7 @@ 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
                 value = int(self['orders'][key])
@@ -3997,8 +4081,10 @@ 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']:
@@ -4012,7 +4098,7 @@ 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 in self['constrained_orders']:
@@ -4047,6 +4133,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 +4192,7 @@ def input_string(self):
             prevleg = leg
 
         if self['orders']:
-            keys = list(self['orders'].keys())
+            keys = visible_orders(self['orders'])
             keys.sort(reverse=True)
             mystr = mystr + " ".join([key + '=' + repr(self['orders'][key]) \
                        for key in keys]) + ' '
@@ -4112,7 +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']]) + ' '
+                       for key in visible_orders(self['squared_orders'])]) + ' '
 
         # Add perturbation orders
         if self['perturbation_couplings']:
@@ -4923,15 +5011,17 @@ 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(visible_orders(self['orders']))]) + ' '
 
         if self['constrained_orders']:
             mystr = mystr + " ".join('%s%s%d' % (key, operator, value) for 
                                      (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']:
@@ -4945,9 +5035,11 @@ 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() \
+              for key in visible_orders(self['squared_orders']) \
                                     if print_weighted or key!='WEIGHTED']) + ' '
 
+        mystr = mystr + self.interference_string()
+
         # Remove last space
         mystr = mystr[:-1]
 
@@ -4989,7 +5081,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 fcd165322..e65aa07fc 100755
--- a/madgraph/core/diagram_generation.py
+++ b/madgraph/core/diagram_generation.py
@@ -41,6 +41,91 @@
 
 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): 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
+    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 +1352,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 +1395,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, diagram_filter)
+
         # Set diagrams to res if not asked to be returned
         if not returndiag:
            self['diagrams'] = res
@@ -1307,6 +1405,56 @@ def generate_diagrams(self, returndiag=False, diagram_filter=False):
         else:
            return not failed_crossing, res
 
+    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
+        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. 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')
+        register_interference_order(process.get('model'))
+        right_process = specialise_interference_process(process)
+
+        right_diagrams = base_objects.DiagramList()
+        if left_diagrams:
+            right_amplitude = Amplitude()
+            right_amplitude.set('process', right_process)
+            try:
+                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:
+            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 +2831,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
@@ -2842,7 +2993,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)
@@ -2850,7 +3005,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/core/helas_objects.py b/madgraph/core/helas_objects.py
index 3b3ae6a80..ef69162c8 100755
--- a/madgraph/core/helas_objects.py
+++ b/madgraph/core/helas_objects.py
@@ -4474,6 +4474,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/fks/fks_base.py b/madgraph/fks/fks_base.py
index 6041ca5bd..d290b054e 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 7f1eeac93..66295500f 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 42c481a73..7ae6c0a5b 100755
--- a/madgraph/interface/madgraph_interface.py
+++ b/madgraph/interface/madgraph_interface.py
@@ -144,6 +144,119 @@ 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. 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*[A-Za-z_]\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:
+        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 = rest
+        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 +939,27 @@ 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")
+        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' (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')
@@ -868,6 +1002,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.")
+        self._log_interference_syntax()
 
         logger.info("--  merge two model to create a new one", '$MG:color:BLUE')
         logger.info("syntax:",'$MG:BOLD')
@@ -1447,6 +1582,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 +3743,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 +4057,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 +4122,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:
@@ -4403,7 +4557,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.
@@ -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]:
@@ -6107,7 +6282,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' % 
@@ -6709,6 +6885,121 @@ 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=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
+        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.
+        """
+
+        # 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 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'):
+            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 +7034,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 +7079,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 +12392,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 0bfc1904c..d145a35a6 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 b84ff0c6b..e7e6472ea 100755
--- a/madgraph/iolibs/export_v4.py
+++ b/madgraph/iolibs/export_v4.py
@@ -8869,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
@@ -8895,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
@@ -8984,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
@@ -9047,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
@@ -9136,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
@@ -10562,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)
 
 
 
@@ -13826,6 +13857,82 @@ 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')]
+        # 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):
+            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(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]:
+                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:
+            return power // 2
+    return None
+
+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 11164c22f..65e53420a 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 17d365834..fc0534df0 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 """
         
@@ -1391,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 6fb8d51ea..91bbf6835 100755
--- a/madgraph/loop/loop_diagram_generation.py
+++ b/madgraph/loop/loop_diagram_generation.py
@@ -30,6 +30,13 @@
 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
+
 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 +765,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 +877,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
 
 
@@ -877,6 +892,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()
 
@@ -900,7 +921,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 +942,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",\
@@ -1029,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 = []
@@ -1069,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)
@@ -1204,6 +1236,76 @@ 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 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/madgraph/loop/loop_helas_objects.py b/madgraph/loop/loop_helas_objects.py
index f01eaeeba..1d599d4e4 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 3a7e9b9b0..dc6f056c3 100755
--- a/madgraph/various/banner.py
+++ b/madgraph/various/banner.py
@@ -5269,7 +5269,14 @@ 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 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/acceptance_tests/test_loop_induced_x_tree.py b/tests/acceptance_tests/test_loop_induced_x_tree.py
new file mode 100644
index 000000000..e306eafdf
--- /dev/null
+++ b/tests/acceptance_tests/test_loop_induced_x_tree.py
@@ -0,0 +1,182 @@
+################################################################################
+#
+# 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)
+
+    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 %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]
+        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,
+                                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)
+
+
+
+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_base_objects.py b/tests/unit_tests/core/test_base_objects.py
index d8c1d072d..c63393c59 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 000000000..972cf8f70
--- /dev/null
+++ b/tests/unit_tests/core/test_tree_interference.py
@@ -0,0 +1,293 @@
+################################################################################
+#
+# 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 contextlib
+import io
+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_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_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."""
+        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,
+                               '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 (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 ['matchbox', 'matrix', 'mg7_v5']:
+                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 000000000..689fc359a
--- /dev/null
+++ b/tests/unit_tests/interface/test_interference_syntax.py
@@ -0,0 +1,288 @@
+################################################################################
+#
+# 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_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',
+                     '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_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."""
+        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')
+        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 000000000..7eb9690fe
--- /dev/null
+++ b/tests/unit_tests/loop/test_loop_induced_x_tree.py
@@ -0,0 +1,234 @@
+################################################################################
+#
+# 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_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)
+
+    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):
+
+    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):
+
+    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]')
+
+    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()