Crossing symmetry

MadSpin/decay.py

@@ -66,6 +66,30 @@
 logger = logging.getLogger('decay.stdout') # -> stdout
 logger_stderr = logging.getLogger('decay.stderr') # ->stderr
 
+import contextlib
+
+@contextlib.contextmanager
+def _no_merge_crossing():
+    """Temporarily disable crossing-symmetry folding (merge_crossing='record')
+    for the enclosed MG5 generation.
+
+    MadSpin's legacy full_decay_chain (madspin_v1) and onshell_v1 paths evaluate
+    the production matrix element through a PDG-dispatch interface (the
+    standalone_msP fortran driver / smatrixhel) that cannot reach a crossing-
+    folded subprocess, so their generation must keep every subprocess on its own
+    (as before the crossing feature).  The density path is crossing-aware and
+    keeps crossing on.  MG_MERGE_CROSSING=off is the documented escape hatch
+    (see madgraph_interface.do_add)."""
+    saved = os.environ.get('MG_MERGE_CROSSING')
+    os.environ['MG_MERGE_CROSSING'] = 'off'
+    try:
+        yield
+    finally:
+        if saved is None:
+            os.environ.pop('MG_MERGE_CROSSING', None)
+        else:
+            os.environ['MG_MERGE_CROSSING'] = saved
+
 import random 
 import math
 from madgraph import MG5DIR, MadGraph5Error
@@ -3097,11 +3121,12 @@ def generate_all_matrix_element(self):
 
         commandline = commandline.replace('add process', 'generate',1)
         logger.info(commandline)
-
-        mgcmd.exec_cmd(commandline, precmd=True)
+
+        with _no_merge_crossing():
+            mgcmd.exec_cmd(commandline, precmd=True)
         commandline = 'output standalone_msP %s %s' % \
-        (pjoin(path_me,'production_me'), ' '.join(list(self.list_branches.keys())))        
-        mgcmd.exec_cmd(commandline, precmd=True)        
+        (pjoin(path_me,'production_me'), ' '.join(list(self.list_branches.keys())))
+        mgcmd.exec_cmd(commandline, precmd=True)
         logger.info('Done %.4g' % (time.time()-start))
 
         # 3. Create all_ME + topology objects ----------------------------------
@@ -3181,7 +3206,8 @@ def generate_all_matrix_element(self):
                         commandline += self.get_proc_with_decay(proc, one_decay, mgcmd._curr_model, self.options)
                 commandline = commandline.replace('add process', 'generate',1)
             logger.info(commandline)
-            mgcmd.exec_cmd(commandline, precmd=True)
+            with _no_merge_crossing():
+                mgcmd.exec_cmd(commandline, precmd=True)
             # remove decay with 0 branching ratio.
             mgcmd.remove_pointless_decay(self.banner.param_card)
             commandline = 'output standalone_msF %s %s' % (pjoin(path_me,'full_me'),
@@ -5123,7 +5149,17 @@ def generate_all_matrix_element(self):
 
         commandline = commandline.replace('add process', 'generate',1)
         mgcmd = self.mgcmd
-        mgcmd.exec_cmd(commandline, precmd=True)
+        # The legacy onshell_v1 path (mode=='onshell') evaluates the production
+        # ME through smatrixhel, which dispatches on concrete PDGs and returns 0
+        # for a crossing-folded subprocess (-> production_me==0 -> ZeroDivision).
+        # Keep every subprocess on its own for it, exactly like full_decay_chain.
+        # The density path (mode=='density') is crossing-aware (GET_DENSITY_IDX
+        # via _resolve_crossed) and keeps crossing on.
+        if self.mode == 'onshell':
+            with _no_merge_crossing():
+                mgcmd.exec_cmd(commandline, precmd=True)
+        else:
+            mgcmd.exec_cmd(commandline, precmd=True)
         # remove decay with 0 branching ratio.
         #mgcmd.remove_pointless_decay(self.banner.param_card)
         #
@@ -5135,16 +5171,52 @@ def generate_all_matrix_element(self):
         # store information about matrix element
         for matrix_element in mgcmd._curr_matrix_elements.get_matrix_elements():
             me_string = matrix_element.get('processes')[0].shell_string()
-            for me in matrix_element.get('processes'):
-                dirpath = pjoin(path_me, ms_me_subdir, 'SubProcesses', "P%s" % me_string)
-                # get the orignal order:
-                initial = []
-                final = [l.get('id') for l in me.get_legs_with_decays()\
-                          if l.get('state') or initial.append(l.get('id'))]
+
+            def register(leg_ids, _pdir="P%s" % me_string):
+                # ``leg_ids`` is the ordered list of external PDG ids (initial
+                # then final, in leg order).  Build the (sorted) lookup tag and
+                # keep the natural leg order for momentum extraction.  A base
+                # process always wins over a crossed one, so never overwrite.
+                initial = [i for i, l in zip(leg_ids, order_state) if not l]
+                final = [i for i, l in zip(leg_ids, order_state) if l]
                 order = (tuple(initial), tuple(final))
-                initial.sort(), final.sort()
-                tag = (tuple(initial), tuple(final))
-                self.all_me[tag] = {'pdir': "P%s" % me_string, 'order': order}
+                tag = (tuple(sorted(initial)), tuple(sorted(final)))
+                self.all_me.setdefault(tag, {'pdir': _pdir, 'order': order})
+
+            for me in matrix_element.get('processes'):
+                legs = me.get_legs_with_decays()
+                order_state = [l.get('state') for l in legs]
+                register([l.get('id') for l in legs])
+
+            # Crossed subprocesses folded into this matrix element with crossing
+            # symmetry on (merge_crossing='record') were NOT generated as their
+            # own directory: only the representative partner is on disk.  Its
+            # crossing-aware SMATRIX (smatrixhel dispatches on the signed PDGs,
+            # see interface_madspin.calculate_matrix_element) can still compute
+            # them, so register each crossed process against the SAME pdir.
+            # Without this the onshell/density lookup misses e.g. the recorded
+            # q q~ > t t~ g (partner of the kept q g > t t~ q) and raises
+            # KeyError.  The decays never cross (they ride on their production
+            # leg), so re-attach the base decay chains before expanding, exactly
+            # as the exporter does for the check_sa crossing demo.
+            try:
+                crossed_processes = matrix_element.get('crossed_processes')
+            except Exception:
+                crossed_processes = []
+            base_decays = matrix_element.get('processes')[0].get('decay_chains')
+            for cross_entry in crossed_processes:
+                proc = cross_entry[0]
+                if base_decays:
+                    proc = copy.copy(proc)
+                    proc.set('decay_chains', base_decays)
+                    # empty LegList (same class as proc's legs) forces
+                    # get_legs_with_decays to recompute with the re-attached decays
+                    proc.set('legs_with_decays', proc.get('legs').__class__())
+                    legs = proc.get_legs_with_decays()
+                else:
+                    legs = proc.get('legs')
+                order_state = [l.get('state') for l in legs]
+                register([l.get('id') for l in legs])
 
         return self.all_me