Crossing symmetry

Template/LO/SubProcesses/makefile

@@ -38,8 +38,16 @@ endif
 
 MATRIX_HEL = $(patsubst %.f,%.o,$(wildcard matrix*_orig.f))
 MATRIX = $(patsubst %.f,%.o,$(wildcard matrix*_optim.f))
+# Crossing-symmetry routers: matrix*_router.f share a base subprocess's matrix
+# element and are never recycled, so they are compiled into both the forhel and
+# the optimized binaries alongside the recycled bases. When recycling is off the
+# matrix*.f glob below already covers them.
+ROUTER = $(patsubst %.f,%.o,$(wildcard matrix*_router.f))
 ifeq ($(strip $(MATRIX_HEL)),)
         MATRIX = $(patsubst %.f,%.o,$(wildcard matrix*.f))
+else
+        MATRIX += $(ROUTER)
+        MATRIX_HEL += $(ROUTER)
 endif
 
 
@@ -102,3 +110,10 @@ initcluster.o: message.inc
 
 clean:
 	$(RM) *.o gensym madevent madevent_forhel
+
+# Track B cross-group crossing: present only in a dependent P directory, it makes
+# the base group's matrix<b> object(s) be symlinked from the base directory
+# instead of recompiled here (see write_crossgroup_mk). Included LAST so its
+# specific rules override the makefile's %.o pattern / $(MATRIX) static-pattern
+# rules. Absent (a no-op) elsewhere.
+-include crossgroup.mk

aloha/__init__.py

@@ -1,4 +1,12 @@
 complex_mass = False  # Tag for activating the complex mass scheme
+t_channel_width = False  # Whether to keep the width i*M*Gamma in the propagator
+                         # denominator for spacelike (t-channel, P^2<0) momenta.
+                         # False (default): drop it there -- the correct tree-level
+                         # treatment outside the complex-mass scheme (a t-channel
+                         # propagator has no pole to regulate, and the spurious
+                         # width breaks gauge cancellations). True: keep the width
+                         # in every propagator (legacy behaviour). Ignored when
+                         # complex_mass is True (the width lives in the mass then).
 unitary_gauge = True  # Tag choosing between Feynman Gauge or unitary gauge
                       # 0/False: Feynman
                       # 1/True: unitary

aloha/aloha_writers.py

@@ -956,8 +956,25 @@ def sort_fct(a, b):
                             out.write('    denom = %(COUP)s/(%(denom)s)\n' % {'COUP': coup_name,\
                                 'denom':self.write_obj(self.routine.denominator)}) 
                     else:
-                        out.write('    denom = %(COUP)s/(P%(i)s(0)**2-P%(i)s(1)**2-P%(i)s(2)**2-P%(i)s(3)**2 - M%(i)s * (M%(i)s -CI* W%(i)s))\n' % \
-                                  {'i': self.outgoing, 'COUP': coup_name})
+                        p2 = 'P%(i)s(0)**2-P%(i)s(1)**2-P%(i)s(2)**2-P%(i)s(3)**2' \
+                              % {'i': self.outgoing}
+                        wdenom = '%(p2)s - M%(i)s * (M%(i)s -CI* W%(i)s)' \
+                              % {'i': self.outgoing, 'p2': p2}
+                        if aloha.t_channel_width:
+                            out.write('    denom = %(COUP)s/(%(wd)s)\n' %
+                                      {'COUP': coup_name, 'wd': wdenom})
+                        else:
+                            # spacelike (t-channel) propagator: no pole to regulate,
+                            # so drop the width unless in the complex-mass scheme.
+                            out.write('    if (dble(%(p2)s).gt.0d0) then\n'
+                                      % {'p2': p2})
+                            out.write('       denom = %(COUP)s/(%(wd)s)\n' %
+                                      {'COUP': coup_name, 'wd': wdenom})
+                            out.write('    else\n')
+                            out.write('       denom = %(COUP)s/(%(p2)s - M%(i)s**2)\n'
+                                      % {'i': self.outgoing, 'COUP': coup_name,
+                                         'p2': p2})
+                            out.write('    endif\n')
                 else:
                     if self.routine.denominator:
                         if 'P1N' not in self.tag:
@@ -2050,8 +2067,19 @@ def sort_fct(a, b):
                             out.write('    denom = %(pre_coup)s%(coup)s%(post_coup)s/(%(denom)s);\n' % \
                                   mydict) 
                     else:
-                        out.write('    denom = %(pre_coup)s%(coup)s%(post_coup)s/((P%(i)s[0]*P%(i)s[0])-(P%(i)s[1]*P%(i)s[1])-(P%(i)s[2]*P%(i)s[2])-(P%(i)s[3]*P%(i)s[3]) - M%(i)s * (M%(i)s -cI* W%(i)s));\n' % \
-                                  mydict)
+                        p2 = '(P%(i)s[0]*P%(i)s[0])-(P%(i)s[1]*P%(i)s[1])-(P%(i)s[2]*P%(i)s[2])-(P%(i)s[3]*P%(i)s[3])' % mydict
+                        wd = '%(p2)s - M%(i)s * (M%(i)s -cI* W%(i)s)' % dict(mydict, p2=p2)
+                        if aloha.t_channel_width:
+                            out.write('    denom = %(pre_coup)s%(coup)s%(post_coup)s/(%(wd)s);\n' % dict(mydict, wd=wd))
+                        else:
+                            # spacelike (t-channel) propagator: no pole to regulate,
+                            # so drop the width unless in the complex-mass scheme.
+                            p2sign = '(%s)%s' % (p2, self.realoperator) if aloha.loop_mode else p2
+                            out.write('    if (%s > 0.){\n' % p2sign)
+                            out.write('       denom = %(pre_coup)s%(coup)s%(post_coup)s/(%(wd)s);\n' % dict(mydict, wd=wd))
+                            out.write('    } else {\n')
+                            out.write('       denom = %(pre_coup)s%(coup)s%(post_coup)s/(%(p2)s - M%(i)s*M%(i)s);\n' % dict(mydict, p2=p2))
+                            out.write('    }\n')
                 else:
                     if self.routine.denominator:
                         raise Exception('modify denominator are not compatible with complex mass scheme')                
@@ -2583,8 +2611,17 @@ def sort_fct(a, b):
                         out.write('    denom = %(COUP)s/(%(denom)s)\n' % {'COUP': coup_name,\
                                 'denom':self.write_obj(self.routine.denominator)}) 
                     else:
-                        out.write('    denom = %(coup)s/(P%(i)s[0]**2-P%(i)s[1]**2-P%(i)s[2]**2-P%(i)s[3]**2 - M%(i)s * (M%(i)s -1j* W%(i)s))\n' % 
-                          {'i': self.outgoing,'coup':coup_name})
+                        p2 = 'P%(i)s[0]**2-P%(i)s[1]**2-P%(i)s[2]**2-P%(i)s[3]**2' % {'i': self.outgoing}
+                        wd = '%(p2)s - M%(i)s * (M%(i)s -1j* W%(i)s)' % {'i': self.outgoing, 'p2': p2}
+                        if aloha.t_channel_width:
+                            out.write('    denom = %(coup)s/(%(wd)s)\n' % {'coup': coup_name, 'wd': wd})
+                        else:
+                            # spacelike (t-channel) propagator: no pole to regulate,
+                            # so drop the width unless in the complex-mass scheme.
+                            out.write('    if (%s).real > 0:\n' % p2)
+                            out.write('        denom = %(coup)s/(%(wd)s)\n' % {'coup': coup_name, 'wd': wd})
+                            out.write('    else:\n')
+                            out.write('        denom = %(coup)s/(%(p2)s - M%(i)s**2)\n' % {'i': self.outgoing, 'coup': coup_name, 'p2': p2})
                 else:
                     if self.routine.denominator:
                         raise Exception('modify denominator are not compatible with complex mass scheme')                

madgraph/core/diagram_generation.py

@@ -444,6 +444,13 @@ def default_setup(self):
         # has_mirror_process is True if the same process but with the
         # two incoming particles interchanged has been generated
         self['has_mirror_process'] = False
+        # Crossed subprocesses folded into this amplitude and NOT generated on
+        # their own (merge_crossing='record'): each entry is
+        # (crossed Process, base_permutation, crossed_permutation), enough for
+        # the exporter to reach the crossed process through this amplitude's
+        # crossing-aware SMATRIX (see MultiProcess.cross_amplitude for the same
+        # permutation pair). Empty in the historical modes.
+        self['crossed_processes'] = []
 
     def __init__(self, argument=None):
         """Allow initialization with Process"""
@@ -470,6 +477,9 @@ def filter(self, name, value):
         if name == 'has_mirror_process':
             if not isinstance(value, bool):
                 raise self.PhysicsObjectError("%s is not a valid boolean" % str(value))
+        if name == 'crossed_processes':
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError("%s is not a valid list" % str(value))
         return True
 
     def get(self, name):
@@ -487,7 +497,8 @@ def get(self, name):
     def get_sorted_keys(self):
         """Return diagram property names as a nicely sorted list."""
 
-        return ['process', 'diagrams', 'has_mirror_process']
+        return ['process', 'diagrams', 'has_mirror_process',
+                'crossed_processes']
 
     def get_number_of_diagrams(self):
         """Returns number of diagrams for this amplitude"""
@@ -1944,6 +1955,20 @@ def get_flavor(id, fsleg):
                             non_permuted_procs.append(fast_proc)
                             logger.info("Crossed process found for %s, reuse diagrams." % \
                                         process.base_string())
+                        elif merge_crossing == 'record':
+                            # Found crossing - do NOT generate a separate
+                            # amplitude, but record the crossed process on the
+                            # base so the exporter can still reach it through the
+                            # base's crossing-aware SMATRIX (its partonic
+                            # contribution is not lost, unlike merge_crossing=True).
+                            amplitudes[crossed_index].get('crossed_processes')\
+                                .append((process, permutations[crossed_index],
+                                         permutation))
+                            logger.info("Crossed process %s recorded on %s "
+                                        "(not generated)." %
+                                        (process.base_string(),
+                                         amplitudes[crossed_index].get('process')
+                                         .base_string()))
                         else:
                             logger.info("Crossed process found for %s, do not generate diagrams." % \
                                         process.base_string())

madgraph/core/helas_objects.py

@@ -3984,6 +3984,11 @@ def default_setup(self):
         # has_mirror_process is True if the same process but with the
         # two incoming particles interchanged has been generated
         self['has_mirror_process'] = False
+        # Crossed subprocesses folded into this ME and NOT generated on their
+        # own (merge_crossing='record'); carried over from the base amplitude so
+        # the exporter can reach them through this ME's crossing-aware SMATRIX.
+        # Each entry is (crossed Process, base_permutation, crossed_permutation).
+        self['crossed_processes'] = []
         self['allowed_flavors'] = [] # list of all allowed flavors for the process
         self['allowed_flavors_with_iden'] = [] # list of all allowed flavors for the process but grouped by identical matrix-element
         self['allowed_flavors_with_iden_sign'] = [] # list of all allowed flavors for the process but grouped by identical matrix-element
@@ -4023,14 +4028,18 @@ def filter(self, name, value):
         if name == 'has_mirror_process':
             if not isinstance(value, bool):
                 raise self.PhysicsObjectError("%s is not a valid boolean" % str(value))
+        if name == 'crossed_processes':
+            if not isinstance(value, list):
+                raise self.PhysicsObjectError("%s is not a valid list" % str(value))
         return True
 
     def get_sorted_keys(self):
         """Return particle property names as a nicely sorted list."""
 
         return ['processes', 'identical_particle_factor',
                 'diagrams', 'color_basis', 'color_matrix',
-                'base_amplitude', 'has_mirror_process']
+                'base_amplitude', 'has_mirror_process',
+                'crossed_processes']
 
     # Enhanced get function
     def get(self, name):
@@ -4055,6 +4064,9 @@ def __init__(self, amplitude=None, optimization=1,
                 self.get('processes').append(amplitude.get('process'))
                 self.set('has_mirror_process',
                          amplitude.get('has_mirror_process'))
+                if amplitude.get('crossed_processes'):
+                    self.set('crossed_processes',
+                             list(amplitude.get('crossed_processes')))
                 self.generate_helas_diagrams(amplitude, optimization, decay_ids)
                 self.calculate_fermionfactors()
                 self.calculate_identical_particle_factor()

madgraph/interface/common_run_interface.py

@@ -218,6 +218,14 @@ def check_set(self, args):
                 self.help_set()
                 raise self.InvalidCmd('set needs an option and an argument')
 
+        if args[0] == 'zerowidth_tchannel':
+            raise self.InvalidCmd(
+                "'zerowidth_tchannel' is a generation-time option: the T-channel "
+                "width treatment is now baked into the matrix element (ALOHA) at "
+                "'output' time and cannot be changed at run time. Choose it in MG5 "
+                "before output ('set zerowidth_tchannel True|False') and regenerate "
+                "the process.")
+
         if args[0] not in self._set_options + list(self.options.keys()):
             self.help_set()
             raise self.InvalidCmd('Possible options for set are %s' % \
@@ -6860,7 +6868,28 @@ def check_card_consistency(self):
             if 'dressed_ee' in  proc_charac['limitations']:
                 if self.run_card['lpp1'] not in [0,1,-1] or self.run_card['lpp1'] not in [0,1,-1]:
                     raise InvalidCmd("dressed lepton mode is not available for this process (see warning associated to the code generation to understand why)")
-            # 
+
+            if 'crossing' in proc_charac['limitations']:
+                # Crossing reuses one matrix element across physically distinct
+                # (crossed) initial states, so a per-beam property is ambiguous.
+                if self.run_card['polbeam1'] or self.run_card['polbeam2']:
+                    raise InvalidCmd(
+                        "Beam polarisation is not compatible with crossing symmetry:\n"
+                        "this process reuses a matrix element across crossed initial\n"
+                        "states, for which a per-beam polarisation is ill-defined.\n"
+                        "Regenerate the process with crossing disabled, e.g.\n"
+                        "  generate <process> --use_crossing=False\n"
+                        "and 'output' again, to run polarised beams.")
+                if 'eva' in (self.run_card['pdlabel'],
+                             self.run_card['pdlabel1'], self.run_card['pdlabel2']):
+                    raise InvalidCmd(
+                        "The EVA luminosity is not compatible with crossing symmetry:\n"
+                        "this process reuses a matrix element across crossed initial\n"
+                        "states, for which the per-beam EVA density is ill-defined.\n"
+                        "Regenerate the process with crossing disabled, e.g.\n"
+                        "  generate <process> --use_crossing=False\n"
+                        "and 'output' again, to use EVA.")
+            #
             if 'fix_scale' in proc_charac['limitations']:
                 if not self.run_card['fixed_fac_scale'] or not self.run_card['fixed_ren_scale']:
                     raise InvalidCmd("Your model is identified as having not SM running of the strong coupling.\n"+\