From 86f17bb64b95f4a289734884bdcb8a312293dab6 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Fri, 21 Aug 2026 22:56:26 +0200 Subject: [PATCH 1/6] get_external_flavors: L-cut wavefunctions are not external legs `output standalone` on a loop-induced process died with an IndexError in export_v4.write_check_sa, on `all_pdgs[j-1]`, because the flavor tuples from `get_external_flavors()` were longer than the process' leg list. The reported symptom pointed at the flavor-grouping block, but that block is innocent: `legs_with_decays` was correct all along (4 entries for g g > z z), and the flavor tuples were wrong (6 entries). Root cause is in `HelasMatrixElement._flavor_enumeration_context`, the single place that decides what the external legs of a matrix element are. It collects every motherless wavefunction and walks `number_external` upwards. In a LoopHelasMatrixElement the L-cut wavefunctions are motherless too and carry number_external = nexternal+1 / nexternal+2, so two fictitious "external legs" (the cut quark and antiquark of the first loop) were appended to `pdgs`. Every consumer of `allowed_flavors` then saw nexternal+2 slots. `LoopHelasMatrixElement.get_nexternal_ninitial` already overrides the base method with exactly the `not wf.get('is_loop')` filter for this reason; the flavor enumeration simply never got the same treatment. Apply it there too. No effect on tree-level matrix elements: no tree wavefunction has is_loop set, so `pdgs` and hence the emitted FLAVOR / PDG_FOR_FLAVOR tables are unchanged (verified byte-identical for a plain and for a merged-particle process). For a loop ME the only shipped consumer today is the madevent loop-induced exporter, where the fix shows up as the `C FLAVOR = [...]` comment in auto_dsig1.f finally having nexternal entries. Co-Authored-By: Claude Opus 5 --- madgraph/core/helas_objects.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/madgraph/core/helas_objects.py b/madgraph/core/helas_objects.py index b93648282..8159fd2a3 100755 --- a/madgraph/core/helas_objects.py +++ b/madgraph/core/helas_objects.py @@ -5479,8 +5479,11 @@ def _flavor_enumeration_context(self, model): """ pdgs = [] pdg_signs = [] + # L-cut wavefunctions of a loop ME are also motherless, but they are not + # external legs of the process (same filter as get_nexternal_ninitial). external_wfs = sorted([wf for wf in self.get_all_wavefunctions() - if len(wf.get('mothers')) == 0], + if len(wf.get('mothers')) == 0 + and not wf.get('is_loop')], key=lambda w: w['number_external']) external_number = 1 id_to_wf = collections.defaultdict(list) From 78ce3f34ed541124ae39d53e2d53b1794c7c3bfa Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Fri, 21 Aug 2026 22:56:40 +0200 Subject: [PATCH 2/6] output standalone: use the MadLoop exporter for a loop-induced process Fixing the flavor tuples above lets write_check_sa run, and `output standalone` then dies one step further on, in the tree-level write_matrix_element_v4, with "wavefunction_rank has not been computed". Same underlying cause: the exporter itself is wrong. For `[noborn=]`, master_interface borrows the MadLoop interface only to validate the model, then switches back to 'MadGraph' and calls create_loop_induced. So the export goes through ExportV4Factory with output_type='default', whose `format.startswith('standalone')` branch returns the tree-level ProcessExporterFortranSA unconditionally -- unlike the two madevent branches right below it, which do check for a LoopAmplitude. ProcessExporterFortranSA cannot write a LoopHelasMatrixElement. Route that branch to the MadLoop standalone exporter when the amplitude is a LoopAmplitude, mirroring what the madevent branches do. The one obstacle was LoopProcessExporterFortranSA.generate_subprocess_directory naming its third positional parameter `second_exporter` and asserting it is None: the base-class convention for that slot is the subprocess number, which is what madgraph_interface.export passes (loop_interface.ML5export omits it). Renamed; `second_exporter` stays available as a keyword and keeps its assert. The result is the ordinary MadLoop standalone layout. For g g > z z it is identical, modulo the process comment line, to what `[sqrvirt=QCD]` -- which reaches the MadLoop interface and so was never broken -- produces on main, and `./check` in the generated P0 dir returns bit-identical values from both. Other formats with no MadLoop backend (matrix, standalone_cpp, mg7) still fail on a loop-induced process; they are out of scope here. Co-Authored-By: Claude Opus 5 --- madgraph/iolibs/export_v4.py | 22 +++++++++++++++++++++- madgraph/loop/loop_exporters.py | 6 ++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py index 459e43a1d..b7de6c95c 100755 --- a/madgraph/iolibs/export_v4.py +++ b/madgraph/iolibs/export_v4.py @@ -11652,8 +11652,28 @@ def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True opt['madanalysis5'] = cmd.options['madanalysis5_path'] if format == 'matrix' or format.startswith('standalone'): + # A loop-induced ([noborn=]) process reaches this factory through the + # MadGraph interface (master_interface switches back to it after + # generation), but ProcessExporterFortranSA cannot write a + # LoopHelasMatrixElement. Use the MadLoop standalone exporter, as + # the madevent branches below already do for their own format. + if format == 'standalone' and cmd._curr_amps and isinstance( + cmd._curr_amps[0], loop_diagram_generation.LoopAmplitude): + import madgraph.loop.loop_exporters as loop_exporters + if not os.path.isdir(os.path.join(cmd._mgme_dir, + 'Template/loop_material')): + raise MadGraph5Error( + 'MG5_aMC cannot find the \'loop_material\' directory' + ' in %s' % str(cmd._mgme_dir)) + if cmd.options['loop_optimized_output']: + MadLoop_SA_options['export_format'] = 'madloop_optimized' + ExporterClass = \ + loop_exporters.LoopProcessOptimizedExporterFortranSA + else: + ExporterClass = loop_exporters.LoopProcessExporterFortranSA + return ExporterClass(cmd._export_dir, MadLoop_SA_options) return ProcessExporterFortranSA(cmd._export_dir, opt, format=format) - + elif format in ['madevent'] and group_subprocesses: if isinstance(cmd._curr_amps[0], loop_diagram_generation.LoopAmplitude): diff --git a/madgraph/loop/loop_exporters.py b/madgraph/loop/loop_exporters.py index c856b53f8..4e2aa9a84 100755 --- a/madgraph/loop/loop_exporters.py +++ b/madgraph/loop/loop_exporters.py @@ -1200,10 +1200,12 @@ def write_process_info_file(self, writer, matrix_element): writer.writelines(proc_include) - def generate_subprocess_directory(self, matrix_element, fortran_model, second_exporter=None): + def generate_subprocess_directory(self, matrix_element, fortran_model, + me_number=None, second_exporter=None): """ To overload the default name for this function such that the correct function is used when called from the command interface """ - + # 3rd positional slot is the subprocess number (base-class convention); + # loop_interface.ML5export omits it, madgraph_interface.export passes it. assert second_exporter is None self.unique_id +=1 return self.generate_loop_subprocess(matrix_element,fortran_model, From 017217d64ff3de080e6176531a27426e25cee98a Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Fri, 21 Aug 2026 23:16:41 +0200 Subject: [PATCH 3/6] refuse loop-induced output for the formats with no MadLoop backend 'standalone' now has a MadLoop backend to route to, but the other formats a loop-induced ([noborn=]) process can reach through the tree-level do_output still have none, and still die deep inside a tree-level exporter: output matrix -> "wavefunction_rank has not been computed" output mg7 / mg7_v5 -> KeyError on the first loop leg, which the mg7 exporter's edge-name map does not contain output standalone_cpp -> same shape; no C++ exporter has a loop backend All pre-existing. Refuse them up front and point at [sqrvirt=], which stays in the MadLoop interface and yields the very same matrix element: g g > h h [sqrvirt=QCD] + output standalone gives a directory whose ./check returns exactly the 3.2829343688358318E-005 of the [noborn=] output. The allow-list is LOOP_INDUCED_FORMATS = ['madevent', 'plugin', 'standalone'], and membership is tested EXACTLY, never with startswith: 'standalone' is a prefix of standalone_cpp / _mg7 / _msP / _msF / _rw, none of which has a loop backend, and a startswith test would silently re-open all five. Checked in three places: MadGraphCmd.do_output, ahead of the rmtree that cleans an existing output directory so a guaranteed refusal never deletes one first; plus ExportV4Factory and ExportCPPFactory as backstops for direct callers. [virt=]/[sqrvirt=] are untouched -- they go through loop_interface.do_output (output_type='madloop') and reach none of the three sites -- and 'output aloha' returns before the check. Co-Authored-By: Claude Opus 5 --- madgraph/interface/madgraph_interface.py | 12 +++++++ madgraph/iolibs/export_cpp.py | 14 +++++++- madgraph/iolibs/export_v4.py | 46 +++++++++++++++++++++--- 3 files changed, 67 insertions(+), 5 deletions(-) diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py index e7cd08fcf..5e5ff0c5a 100755 --- a/madgraph/interface/madgraph_interface.py +++ b/madgraph/interface/madgraph_interface.py @@ -9806,6 +9806,18 @@ def do_output(self, line): # merge_quartic_vertices can be resolved -- before anything is built self.apply_quartic_diagram_order(options) + # A loop-induced process is exported by this tree-level do_output (see + # create_loop_induced), but only the formats in LOOP_INDUCED_FORMATS + # have a loop backend to route it to. Refuse the others here, ahead of + # the directory cleaning just below, so that a guaranteed refusal never + # deletes an existing output directory first. The exporter factories + # carry the same check as a backstop. + if self._export_format not in export_v4.LOOP_INDUCED_FORMATS and \ + self._curr_amps and isinstance(self._curr_amps[0], + loop_diagram_generation.LoopAmplitude): + raise self.InvalidCmd(export_v4.loop_induced_not_supported_msg( + self._export_format, self._curr_amps[0].get('process'))) + # check if os.path.realpath(self._export_dir) == os.getcwd(): if len(args) == 0: diff --git a/madgraph/iolibs/export_cpp.py b/madgraph/iolibs/export_cpp.py index 1b83b6f67..7561e961c 100755 --- a/madgraph/iolibs/export_cpp.py +++ b/madgraph/iolibs/export_cpp.py @@ -40,10 +40,12 @@ import madgraph.iolibs.file_writers as writers import madgraph.iolibs.template_files as template_files import madgraph.iolibs.ufo_expression_parsers as parsers +import madgraph.loop.loop_diagram_generation as loop_diagram_generation import madgraph.various.banner as banner_mod from madgraph import MadGraph5Error, InvalidCmd, MG5DIR from madgraph.iolibs.files import cp, ln, mv +import madgraph.iolibs.export_v4 as export_v4 from madgraph.iolibs.export_v4 import VirtualExporter, ProcessExporterFortran import madgraph.various.misc as misc @@ -3559,7 +3561,17 @@ def ExportCPPFactory(cmd, group_subprocesses=False, cmd_options={}): opt = dict(cmd.options) opt['output_options'] = cmd_options cformat = cmd._export_format - + + # None of the C++ exporters below has a MadLoop backend, so a loop-induced + # process would reach them as a LoopHelasMatrixElement whose loop legs they + # cannot even index (the mg7 exporter builds its edge names from the + # external legs alone). Refuse it here instead. Plugins are left alone: + # they are free to implement their own loop support. + if cformat not in export_v4.LOOP_INDUCED_FORMATS and cmd._curr_amps and \ + isinstance(cmd._curr_amps[0], loop_diagram_generation.LoopAmplitude): + raise InvalidCmd(export_v4.loop_induced_not_supported_msg( + cformat, cmd._curr_amps[0].get('process'))) + if cformat == 'pythia8': return ProcessExporterPythia8(cmd._export_dir, opt) elif cformat == 'standalone_cpp': diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py index b7de6c95c..b86ff5e1d 100755 --- a/madgraph/iolibs/export_v4.py +++ b/madgraph/iolibs/export_v4.py @@ -11502,7 +11502,41 @@ def create_param_card(self, write_special=True): rule_card_path=rule_card, mssm_convert=True, write_special=write_special) - + +# The output formats that can serve a loop-induced ([noborn=]) process coming +# through the tree-level do_output: 'madevent' has the LoopInducedExporterME* +# exporters, 'standalone' is routed to the MadLoop standalone exporter, and a +# plugin is free to bring its own. Every other format sends the +# LoopHelasMatrixElement to a tree-level exporter that cannot write it. +# Membership must be tested exactly: 'standalone' is a prefix of +# standalone_cpp / _mg7 / _msP / _msF / _rw, none of which has a loop backend. +LOOP_INDUCED_FORMATS = ['madevent', 'plugin', 'standalone'] + +def loop_induced_not_supported_msg(format, process=None): + """Error text for an output format that has no MadLoop backend. + + A loop-induced ([noborn=]) process is exported by the *tree-level* output + machinery: master_interface only borrows the MadLoop interface to validate + the model, then switches back to 'MadGraph' and calls create_loop_induced. + So a format whose exporter cannot write a LoopHelasMatrixElement has to say + so here rather than let the tree-level exporter fail deep inside. + + The same matrix element is available through [sqrvirt=], which does stay in + the MadLoop interface and therefore reaches the MadLoop exporters. + """ + + orders = 'QCD' + if process: + try: + orders = ' '.join(process.get('perturbation_couplings')) or orders + except Exception: + pass + + return """The '%(format)s' output format does not support loop-induced processes. +Generate the process with [sqrvirt=%(orders)s] rather than [noborn=%(orders)s] to obtain the +same matrix element as a standalone MadLoop output, or use 'output madevent' +to integrate it.""" % {'format': format, 'orders': orders} + def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True, cmd_options={}): """ Determine which Export_v4 class is required. cmd is the command interface containing all potential usefull information. @@ -11655,10 +11689,14 @@ def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True # A loop-induced ([noborn=]) process reaches this factory through the # MadGraph interface (master_interface switches back to it after # generation), but ProcessExporterFortranSA cannot write a - # LoopHelasMatrixElement. Use the MadLoop standalone exporter, as - # the madevent branches below already do for their own format. - if format == 'standalone' and cmd._curr_amps and isinstance( + # LoopHelasMatrixElement. + if cmd._curr_amps and isinstance( cmd._curr_amps[0], loop_diagram_generation.LoopAmplitude): + # Only plain 'standalone' has a MadLoop backend to route to, as + # the madevent branches below have for their own format. + if format not in LOOP_INDUCED_FORMATS: + raise InvalidCmd( + loop_induced_not_supported_msg(format, curr_proc)) import madgraph.loop.loop_exporters as loop_exporters if not os.path.isdir(os.path.join(cmd._mgme_dir, 'Template/loop_material')): From 79ab2288b8b962bbb592b5cc5c427f914fafa983 Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Fri, 21 Aug 2026 23:16:41 +0200 Subject: [PATCH 4/6] validate_model: bootstrap the model for [noborn=] `generate g g > h [noborn=QCD]` as the first command of a session died with AttributeError: 'NoneType' object has no attribute 'merged_particles'. master_interface calls validate_model directly for noborn, before create_loop_induced runs check_add. The [virt=]/[real=] paths go through check_add -> check_generate first, which imports the Standard Model when none is active; noborn never gets that far. A guard on the merged_particles line alone would not do: the rest of validate_model dereferences _curr_model too (perturbation_couplings, get('name'), get('gauge')). So import the default model at the top instead, exactly as check_generate does -- the existing sm -> loop_sm upgrade below then runs unchanged. Co-Authored-By: Claude Opus 5 --- madgraph/interface/loop_interface.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/madgraph/interface/loop_interface.py b/madgraph/interface/loop_interface.py index f3e53c7a5..e67159093 100755 --- a/madgraph/interface/loop_interface.py +++ b/madgraph/interface/loop_interface.py @@ -298,11 +298,18 @@ def proc_validity(self, proc, mode): def validate_model(self, loop_type='virtual',coupling_type=['QCD'], stop=True): """ Upgrade the model sm to loop_sm if needed """ - # Allow to call this function with a string instead of a list of + # Allow to call this function with a string instead of a list of # perturbation orders. if isinstance(coupling_type,str): coupling_type = [coupling_type,] + # Everything below assumes a model. [virt=]/[real=] get one from + # check_generate before we are called, but master_interface calls us + # directly for [noborn=], ahead of create_loop_induced's check_add. + if not self._curr_model: + logger.info("No model currently active, so we import the Standard Model") + self.do_import('model sm') + active_interface = getattr(self, 'current_interface', None) ## if coupling_type!= ['QCD'] and loop_type not in ['virtual','noborn']: From c1a8eb5ae5bf78e14b8754f546996e7865f3bffa Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Fri, 21 Aug 2026 23:17:10 +0200 Subject: [PATCH 5/6] tests: cover loop-induced output, allowed and refused Nothing covered this before: the only test of non-optimized loop-induced output was the IOTest short_ML_SMQCD_LoopInduced/gg_hh, which drives the exporter directly and never reaches ExportV4Factory. tests/unit_tests/loop/test_loop_induced_output.py, 15 tests, ~22 s: - the flavor tuples of a loop ME have nexternal entries with the right PDGs (3 for g g > h, 4 for g g > z z); - 'output standalone' goes through and writes a MadLoop directory, and the factory hands it a MadLoop exporter while a tree process keeps the tree-level one; - matrix / mg7 / standalone_cpp / standalone_msP refuse it and name [sqrvirt=], the refusal does not wipe an existing output directory, and LOOP_INDUCED_FORMATS is asserted not to contain the standalone_* formats so nobody can turn the membership test into a startswith; - [noborn=] works with no model imported; - madevent on [noborn=], and standalone on [sqrvirt=] and [virt=], still work -- these are the controls: the refusal keys off the amplitude being a LoopAmplitude, which those are too, and they are spared only because they never reach the two factories. 11 of the 15 fail on a clean tree; the 4 that pass are exactly those controls. testIO_loop_induced_standalone_output joins the existing MadLoop_output_from_the_interface group: runs the real `output standalone` on g g > h [noborn=QCD] through MasterCmd -- with no `import model`, so it covers the validate_model bootstrap too -- and compares check_sa.f and loop_matrix.f. On a clean tree it fails with the original IndexError. ~4 s. One trap worth knowing: building a loop ME with compute_loop_nc=False leaves colour data that a later loop-induced madevent output reuses, and it then dies in get_icolamp_lines on a None Nc power. That is pre-existing and reproducible on a clean tree; the tests build their matrix elements the way group_subprocs does (compute_loop_nc=True) to stay clear of it. Co-Authored-By: Claude Opus 5 --- tests/acceptance_tests/test_cmd_madloop.py | 22 + ...hLI_IOTest%SubProcesses%P0_gg_h%check_sa.f | 729 ++++ ..._IOTest%SubProcesses%P0_gg_h%loop_matrix.f | 3180 +++++++++++++++++ .../loop/test_loop_induced_output.py | 335 ++ 4 files changed, 4266 insertions(+) create mode 100644 tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%check_sa.f create mode 100644 tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%loop_matrix.f create mode 100644 tests/unit_tests/loop/test_loop_induced_output.py diff --git a/tests/acceptance_tests/test_cmd_madloop.py b/tests/acceptance_tests/test_cmd_madloop.py index ee46b977e..48188e4ee 100755 --- a/tests/acceptance_tests/test_cmd_madloop.py +++ b/tests/acceptance_tests/test_cmd_madloop.py @@ -1055,6 +1055,28 @@ def run_cmd(cmd): IOTests.IOTest.remove_f77_function_from_file( pjoin(self.IOpath,'ggttx_IOTest', 'SubProcesses','MadLoopCommons.f'), 'PRINT_MADLOOP_BANNER') + + @IOTests.createIOTest(groupName='MadLoop_output_from_the_interface') + def testIO_loop_induced_standalone_output(self): + r""" target: gghLI_IOTest/SubProcesses/P0_gg_h/[(check_sa|loop_matrix)\.f] + """ + # A loop-induced ([noborn=]) process is exported from the MadGraph + # interface, which used to hand the loop matrix element to the + # tree-level standalone exporter and crash. + interface = MGCmd.MasterCmd() + interface.no_notification() + + # Select the Tensor Integral to include in the test + misc.deactivate_dependence('pjfry', cmd = interface, log='stdout') + misc.deactivate_dependence('samurai', cmd = interface, log='stdout') + misc.deactivate_dependence('golem', cmd = interface, log='stdout') + misc.activate_dependence('ninja', cmd = interface, log='stdout',MG5dir=MG5DIR) + + # no 'import model': validate_model must bootstrap sm -> loop_sm itself + interface.exec_cmd('generate g g > h [noborn=QCD]', errorhandling=False, + printcmd=False, precmd=True, postcmd=True) + interface.onecmd('output standalone %s -f' % + str(pjoin(self.IOpath,'gghLI_IOTest'))) diff --git a/tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%check_sa.f b/tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%check_sa.f new file mode 100644 index 000000000..0df98dd0d --- /dev/null +++ b/tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%check_sa.f @@ -0,0 +1,729 @@ + PROGRAM DRIVER +C ***************************************************************** +C ******** +C THIS IS THE DRIVER FOR CHECKING THE STANDALONE MATRIX ELEMENT. +C IT USES A SIMPLE PHASE SPACE GENERATOR +C ***************************************************************** +C ******** + IMPLICIT NONE +C +C CONSTANTS +C + REAL*8 ZERO + PARAMETER (ZERO=0D0) + + LOGICAL READPS + PARAMETER (READPS = .FALSE.) + + INTEGER NPSPOINTS + PARAMETER (NPSPOINTS = 10) + +C integer nexternal C number particles (incoming+outgoing) in the +C me + INTEGER NEXTERNAL, NINCOMING + PARAMETER (NEXTERNAL=3,NINCOMING=2) + + CHARACTER(512) MADLOOPRESOURCEPATH + +C +C INCLUDE FILES +C +C the include file with the values of the parameters and masses +C + + INCLUDE 'coupl.inc' +C particle masses + REAL*8 PMASS(NEXTERNAL) +C integer n_max_cg + INCLUDE 'ngraphs.inc' + INCLUDE 'nsquaredSO.inc' + + INCLUDE 'MadLoopParams.inc' + +C +C LOCAL +C + INTEGER I,J,K,L +C four momenta. Energy is the zeroth component. + REAL*8 P(0:3,NEXTERNAL) + INTEGER MATELEM_ARRAY_DIM + REAL*8 , ALLOCATABLE :: MATELEM(:,:) + REAL*8 SQRTS,AO2PI,TOTMASS +C sqrt(s)= center of mass energy + REAL*8 PIN(0:3), POUT(0:3) + CHARACTER*120 BUFF(NEXTERNAL) + INTEGER RETURNCODE, UNITS, TENS, HUNDREDS + INTEGER NSQUAREDSO_LOOP + REAL*8 , ALLOCATABLE :: PREC_FOUND(:) + INTEGER NB_INTER + LOGICAL INIT + DATA INIT/.TRUE./ + COMMON/INITCHECKSA/INIT + + INTEGER NLOOPFLOWS + PARAMETER (NLOOPFLOWS=1) + INTEGER NCOMB + PARAMETER (NCOMB=4) + INTEGER N_CHANGING + PARAMETER (N_CHANGING=2) + INTEGER N_COMB_RHO + PARAMETER (N_COMB_RHO=2) + INTEGER H + INTEGER HELC(NEXTERNAL,NCOMB) + DOUBLE PRECISION ALPHAS, MU_R2 + + DOUBLE COMPLEX INTER_DENS(0:3,0:NSQUAREDSO) + DOUBLE COMPLEX RMATRIX((N_COMB_RHO*(N_COMB_RHO+1))/2, + $ 0:NSQUAREDSO) + DOUBLE PRECISION RES(0:3,0:NSQUAREDSO), SUM_INT + INTEGER HEL_MULT + PARAMETER (HEL_MULT = 1) + INTEGER H1, H2 + +C +C GLOBAL VARIABLES +C +C This is from ML code for the list of split orders selected by +C the process definition +C + INTEGER NLOOPCHOSEN + CHARACTER*20 CHOSEN_LOOP_SO_INDICES(NSQUAREDSO) + LOGICAL CHOSEN_LOOP_SO_CONFIGS(NSQUAREDSO) + COMMON/ML5_0_CHOSEN_LOOP_SQSO/CHOSEN_LOOP_SO_CONFIGS + +C +C EXTERNAL +C + REAL*8 DOT + EXTERNAL DOT + + INTEGER NHEL(NEXTERNAL) + INTEGER, ALLOCATABLE :: POS(:) + INTEGER, ALLOCATABLE :: ALLOW_HEL(:) + DOUBLE COMPLEX, ALLOCATABLE :: INTER(:) + INTEGER SOL + +C +C BEGIN CODE +C + IF (INIT) THEN + INIT=.FALSE. + CALL ML5_0_GET_ANSWER_DIMENSION(MATELEM_ARRAY_DIM) + ALLOCATE(MATELEM(0:3,0:MATELEM_ARRAY_DIM)) + CALL ML5_0_GET_NSQSO_LOOP(NSQUAREDSO_LOOP) + ALLOCATE(PREC_FOUND(0:NSQUAREDSO_LOOP)) +C +C INITIALIZATION CALLS +C +C Call to initialize the values of the couplings, masses and +C widths +C used in the evaluation of the matrix element. The primary +C parameters of the +C models are read from Cards/param_card.dat. The secondary +C parameters are calculated +C in Source/MODEL/couplings.f. The values are stored in common +C blocks that are listed +C in coupl.inc . +C first call to setup the paramaters + CALL SETPARA('param_card.dat') +C set up masses + INCLUDE 'pmass.inc' + + ENDIF + +C Start by initializing what is the squared split orders indices +C chosen + NLOOPCHOSEN=0 + DO I=1,NSQUAREDSO + IF (CHOSEN_LOOP_SO_CONFIGS(I)) THEN + NLOOPCHOSEN=NLOOPCHOSEN+1 + WRITE(CHOSEN_LOOP_SO_INDICES(NLOOPCHOSEN),'(I3,A2)') I,'L)' + ENDIF + ENDDO + + AO2PI=G**2/(8.D0*(3.14159265358979323846D0**2)) + + WRITE(*,*) 'AO2PI=',AO2PI +C Now use a simple multipurpose PS generator (RAMBO) just to get a +C RANDOM set of four momenta of given masses pmass(i) to be used +C to evaluate +C the madgraph matrix-element. +C Alternatevely, here the user can call or set the four momenta at +C his will, see below. +C + IF(NINCOMING.EQ.1) THEN + SQRTS=PMASS(1) + ELSE + TOTMASS = 0.0D0 + DO I=1,NEXTERNAL + TOTMASS = TOTMASS + PMASS(I) + ENDDO +C CMS energy in GEV + SQRTS=MAX(1000D0,2.0D0*TOTMASS) + ENDIF + + CALL PRINTOUT() + + DO K=1,NPSPOINTS + + IF(READPS) THEN + OPEN(967, FILE='PS.input', ERR=976, STATUS='OLD', + $ ACTION='READ') + DO I=1,NEXTERNAL + READ(967,*,END=978) P(0,I),P(1,I),P(2,I),P(3,I) + ENDDO + GOTO 978 + 976 CONTINUE + STOP 'Could not read the PS.input phase-space point.' + 978 CONTINUE + CLOSE(967) + ELSE + IF ((NINCOMING.EQ.2).AND.((NEXTERNAL - NINCOMING .EQ.1))) + $ THEN + IF (PMASS(3).EQ.0.0D0) THEN + STOP 'Cannot generate 2>1 kin. config. with m3=0.0d0' + ELSE +C deal with the case of only one particle in the final +C state + P(0,1) = PMASS(3)/2D0 + P(1,1) = 0D0 + P(2,1) = 0D0 + P(3,1) = PMASS(3)/2D0 + IF (PMASS(1).GT.0D0) THEN + P(3,1) = DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2) + ENDIF + P(0,2) = PMASS(3)/2D0 + P(1,2) = 0D0 + P(2,2) = 0D0 + P(3,2) = -PMASS(3)/2D0 + IF (PMASS(2) > 0D0) THEN + P(3,2) = -DSQRT(PMASS(3)**2/4D0 - PMASS(1)**2) + ENDIF + P(0,3) = PMASS(3) + P(1,3) = 0D0 + P(2,3) = 0D0 + P(3,3) = 0D0 + ENDIF + ELSE + CALL GET_MOMENTA(SQRTS,PMASS,P) + ENDIF + ENDIF + + DO I=0,3 + PIN(I)=0.0D0 + DO J=1,NINCOMING + PIN(I)=PIN(I)+P(I,J) + ENDDO + ENDDO + +C In standalone mode, always use sqrt_s as the renormalization +C scale. + SQRTS=DSQRT(DABS(DOT(PIN(0),PIN(0)))) + MU_R=SQRTS + +C Update the couplings with the new MU_R + CALL UPDATE_AS_PARAM(1) + +C Optionally the user can set where to find the +C MadLoop5_resources folder. +C Otherwise it will look for it automatically and find it if it +C has not +C been moved +C MadLoopResourcePath = '' +C CALL SETMADLOOPPATH(MadLoopResourcePath) +C To force the stabiliy check to also be performed in the +C initialization phase +C CALL ML5_0_FORCE_STABILITY_CHECK(.TRUE.) +C To chose a particular tartget split order, SOTARGET is an +C integer labeling +C the possible squared order couplings contributions (only in +C optimized mode) +C CALL ML5_0_SET_COUPLINGORDERS_TARGET(SOTARGET) + +C +C Now we can call the matrix element +C + CALL ML5_0_SLOOPMATRIX_THRES(P,MATELEM,-1.0D0,PREC_FOUND + $ ,RETURNCODE) + + CALL ML5_0_COMPUTE_RES_FROM_JAMP(RES,HEL_MULT) +C WRITE(*,*) "ML5_0_COMPUTE_RES_FROM_JAMP", RES(1:3,0) + + + IF (K.EQ.NPSPOINTS) THEN + WRITE (*,*) + WRITE (*,*) ' Phase space point:' + WRITE (*,*) + WRITE (*,*) '---------------------------------' + WRITE (*,*) 'n E px py pz m' + DO I=1,NEXTERNAL + WRITE (*,'(i2,1x,5e15.7)') I, P(0,I),P(1,I),P(2,I),P(3,I) + $ ,DSQRT(DABS(DOT(P(0,I),P(0,I)))) + ENDDO + WRITE (*,*) '---------------------------------' + WRITE (*,*) 'Detailed result for each coupling orders' + $ //' combination.' + + WRITE (*,*) '---------------------------------' + WRITE(*,*) 'All loop contributions are of split orders' + $ //' (QCD=4)' + WRITE (*,*) '---------------------------------' + UNITS=MOD(RETURNCODE,10) + TENS=(MOD(RETURNCODE,100)-UNITS)/10 + HUNDREDS=(RETURNCODE-TENS*10-UNITS)/100 + IF (HUNDREDS.EQ.1) THEN + IF (TENS.EQ.3.OR.TENS.EQ.4) THEN + WRITE(*,*) 'Unknown numerical stability because MadLoop' + $ //' is in the initialization stage.' + ELSE + WRITE(*,*) 'Unknown numerical stability, check CTModeRun' + $ //' value in MadLoopParams.dat.' + ENDIF + ELSEIF (HUNDREDS.EQ.2) THEN + WRITE(*,*) 'Stable kinematic configuration (SPS).' + ELSEIF (HUNDREDS.EQ.3) THEN + WRITE(*,*) 'Unstable kinematic configuration (UPS).' + WRITE(*,*) 'Quadruple precision rescue successful.' + ELSEIF (HUNDREDS.EQ.4) THEN + WRITE(*,*) 'Exceptional kinematic configuration (EPS).' + WRITE(*,*) 'Both double an quadruple precision' + $ //' computations, are unstable.' + ENDIF + IF (TENS.EQ.2.OR.TENS.EQ.4) THEN + WRITE(*,*) 'Quadruple precision computation used.' + ENDIF + IF (HUNDREDS.NE.1) THEN + IF (PREC_FOUND(0).GT.0.0D0) THEN + WRITE(*,'(1x,a23,1x,1e10.2)') 'Relative accuracy =' + $ ,PREC_FOUND(0) + ELSEIF (PREC_FOUND(0).EQ.0.0D0) THEN + WRITE(*,'(1x,a23,1x,1e10.2,1x,a30)') 'Relative accuracy ' + $ //' =',PREC_FOUND(0),'(i.e. beyond double precision)' + ELSE + WRITE(*,*) 'Estimated accuracy could not be computed for' + $ //' an unknown reason.' + ENDIF + ENDIF + WRITE (*,'(1x,a23,3x,i3)') 'MadLoop return code =' + $ ,RETURNCODE + WRITE (*,*) '---------------------------------' + IF (NLOOPCHOSEN.NE.NSQUAREDSO) THEN + WRITE (*,*) 'Selected squared coupling orders combination' + $ //' for the loop summed result below:' + WRITE (*,*) (CHOSEN_LOOP_SO_INDICES(I),I=1,NLOOPCHOSEN) + ENDIF + WRITE (*,*) '---------------------------------' + WRITE (*,*) ' This is a loop induced process, so only the ' + WRITE (*,*) ' unnormalized finite part is output here. Be' + $ //' aware ' + WRITE (*,*) ' that all loops are expected to beUV-finite as' + $ //' no ' + WRITE (*,*) ' renormalization prescription is considered. ' + WRITE (*,*) '---------------------------------' + WRITE (*,*) 'Matrix element finite = ', MATELEM(1,0), + $ ' GeV^',-(2*NEXTERNAL-8) + WRITE (*,*) '---------------------------------' + OPEN(69, FILE='result.dat', ERR=976, ACTION='WRITE') + DO I=1,NEXTERNAL + WRITE (69,'(a2,1x,5ES30.15E3)') 'PS',P(0,I),P(1,I),P(2,I) + $ ,P(3,I) + ENDDO + WRITE (69,'(a3,1x,i3)') 'EXP',-(2*NEXTERNAL-8) + WRITE (69,'(a4,1x,1ES30.15E3)') 'BORN',0.0D0 + WRITE (69,'(a3,1x,1ES30.15E3)') 'FIN',MATELEM(1,0) + WRITE (69,'(a4,1x,1ES30.15E3)') '1EPS',MATELEM(2,0) + WRITE (69,'(a4,1x,1ES30.15E3)') '2EPS',MATELEM(3,0) + WRITE (69,'(a6,1x,1ES30.15E3)') 'ASO2PI',AO2PI + WRITE (69,*) 'Export_Format LoopInduced' + WRITE (69,'(a7,1x,i3)') 'RETCODE',RETURNCODE + WRITE (69,'(a3,1x,1e10.4)') 'ACC',PREC_FOUND(0) + WRITE (69,*) 'Born_kept F' + WRITE (69,*) 'Loop_kept',(CHOSEN_LOOP_SO_CONFIGS(I),I=1 + $ ,NSQUAREDSO) + + WRITE (69,*) 'Split_Orders_Names QCD' + WRITE (69,*) 'Loop_SO_Results 4' + WRITE (69,*) 'SO_Loop ACC ',PREC_FOUND(1) + WRITE (69,*) 'SO_Loop FIN ',MATELEM(1,1) + WRITE (69,*) 'SO_Loop 1EPS ',MATELEM(2,1) + WRITE (69,*) 'SO_Loop 2EPS ',MATELEM(3,1) + + CLOSE(69) + ELSE + WRITE (*,*) 'PS Point #',K,' done.' + ENDIF + ENDDO + +C C +C C Copy down here (or read in) the four momenta as a string. +C C +C C +C buff(1)=" 1 0.5630480E+04 0.0000000E+00 0.0000000E+00 +C 0.5630480E+04" +C buff(2)=" 2 0.5630480E+04 0.0000000E+00 0.0000000E+00 +C -0.5630480E+04" +C buff(3)=" 3 0.5466073E+04 0.4443190E+03 0.2446331E+04 +C -0.4864732E+04" +C buff(4)=" 4 0.8785819E+03 -0.2533886E+03 0.2741971E+03 +C 0.7759741E+03" +C buff(5)=" 5 0.4916306E+04 -0.1909305E+03 -0.2720528E+04 +C 0.4088757E+04" +C C +C C Here the k,E,px,py,pz are read from the string into the +C momenta array. +C C k=1,2 : incoming +C C k=3,nexternal : outgoing +C C +C do i=1,nexternal +C read (buff(i),*) k, P(0,i),P(1,i),P(2,i),P(3,i) +C enddo +C +C C print the momenta out +C +C do i=1,nexternal +C write (*,'(i2,1x,5e15.7)') i, P(0,i),P(1,i),P(2,i),P(3,i), +C .dsqrt(dabs(DOT(p(0,i),p(0,i)))) +C enddo +C +C CALL SLOOPMATRIX(P,MATELEM) +C +C write (*,*) "-------------------------------------------------" +C write (*,*) "Matrix element = ", MATELEM(1), " GeV^", +C &-(2*nexternal-8) +C write (*,*) "-------------------------------------------------" + + END + + + DOUBLE PRECISION FUNCTION DOT(P1,P2) +C ************************************************************* +C 4-Vector Dot product +C ************************************************************* + IMPLICIT NONE + DOUBLE PRECISION P1(0:3),P2(0:3) + DOT=P1(0)*P2(0)-P1(1)*P2(1)-P1(2)*P2(2)-P1(3)*P2(3) + END + + + SUBROUTINE GET_MOMENTA(ENERGY,PMASS,P) +C auxiliary function to change convention between madgraph and +C rambo +C four momenta. + IMPLICIT NONE + INTEGER NEXTERNAL, NINCOMING + PARAMETER (NEXTERNAL=3,NINCOMING=2) +C ARGUMENTS + REAL*8 ENERGY,PMASS(NEXTERNAL),P(0:3,NEXTERNAL),PRAMBO(4,10),WGT +C LOCAL + INTEGER I + REAL*8 ETOT2,MOM,M1,M2,E1,E2 + + ETOT2=ENERGY**2 + M1=PMASS(1) + M2=PMASS(2) + MOM=(ETOT2**2 - 2*ETOT2*M1**2 + M1**4 - 2*ETOT2*M2**2 - 2*M1**2 + $ *M2**2 + M2**4)/(4.*ETOT2) + MOM=DSQRT(MOM) + E1=DSQRT(MOM**2+M1**2) + E2=DSQRT(MOM**2+M2**2) +C write (*,*) e1+e2,mom + + IF(NINCOMING.EQ.2) THEN + + P(0,1)=E1 + P(1,1)=0D0 + P(2,1)=0D0 + P(3,1)=MOM + + P(0,2)=E2 + P(1,2)=0D0 + P(2,2)=0D0 + P(3,2)=-MOM + + CALL RAMBO(NEXTERNAL-2,ENERGY,PMASS(3),PRAMBO,WGT) + DO I=3, NEXTERNAL + P(0,I)=PRAMBO(4,I-2) + P(1,I)=PRAMBO(1,I-2) + P(2,I)=PRAMBO(2,I-2) + P(3,I)=PRAMBO(3,I-2) + ENDDO + + ELSEIF(NINCOMING.EQ.1) THEN + + P(0,1)=ENERGY + P(1,1)=0D0 + P(2,1)=0D0 + P(3,1)=0D0 + + CALL RAMBO(NEXTERNAL-1,ENERGY,PMASS(2),PRAMBO,WGT) + DO I=2, NEXTERNAL + P(0,I)=PRAMBO(4,I-1) + P(1,I)=PRAMBO(1,I-1) + P(2,I)=PRAMBO(2,I-1) + P(3,I)=PRAMBO(3,I-1) + ENDDO + ENDIF + + RETURN + END + + + SUBROUTINE RAMBO(N,ET,XM,P,WT) +C ***************************************************************** +C ***** +C RAMBO * +C RA(NDOM) M(OMENTA) B(EAUTIFULLY) O(RGANIZED) +C * +C * +C A DEMOCRATIC MULTI-PARTICLE PHASE SPACE GENERATOR +C * +C AUTHORS: S.D. ELLIS, R. KLEISS, W.J. STIRLING +C * +C THIS IS VERSION 1.0 - WRITTEN BY R. KLEISS +C * +C -- ADJUSTED BY HANS KUIJF, WEIGHTS ARE LOGARITHMIC (20-08-90) +C * +C * +C N = NUMBER OF PARTICLES +C * +C ET = TOTAL CENTRE-OF-MASS ENERGY +C * +C XM = PARTICLE MASSES ( DIM=NEXTERNAL-nincoming ) +C * +C P = PARTICLE MOMENTA ( DIM=(4,NEXTERNAL-nincoming) ) +C * +C WT = WEIGHT OF THE EVENT +C * +C ***************************************************************** +C ***** + IMPLICIT REAL*8(A-H,O-Z) + INTEGER NEXTERNAL, NINCOMING + PARAMETER (NEXTERNAL=3,NINCOMING=2) + DIMENSION XM(*),P(4,*) + DIMENSION Q(4,NEXTERNAL-NINCOMING),Z(NEXTERNAL-NINCOMING),R(4) + $ ,B(3),P2(NEXTERNAL-NINCOMING),XM2(NEXTERNAL-NINCOMING) + $ ,E(NEXTERNAL-NINCOMING),V(NEXTERNAL-NINCOMING),IWARN(5) + SAVE ACC,ITMAX,IBEGIN,IWARN,TWOPI, Z, PO2LOG + DATA ACC/1.D-14/,ITMAX/6/,IBEGIN/0/,IWARN/5*0/ +C +C INITIALIZATION STEP: FACTORIALS FOR THE PHASE SPACE WEIGHT + IF(IBEGIN.NE.0) GOTO 103 + IBEGIN=1 + TWOPI=8.*DATAN(1.D0) + PO2LOG=LOG(TWOPI/4.) + Z(2)=PO2LOG + DO 101 K=3,(NEXTERNAL-NINCOMING) + 101 Z(K)=Z(K-1)+PO2LOG-2.*LOG(DFLOAT(K-2)) + DO 102 K=3,(NEXTERNAL-NINCOMING) + 102 Z(K)=(Z(K)-LOG(DFLOAT(K-1))) +C +C CHECK ON THE NUMBER OF PARTICLES + 103 IF(N.GT.1.AND.N.LT.101) GOTO 104 + PRINT 1001,N + STOP +C +C CHECK WHETHER TOTAL ENERGY IS SUFFICIENT; COUNT NONZERO MASSES + 104 XMT=0. + NM=0 + DO 105 I=1,N + IF(XM(I).NE.0.D0) NM=NM+1 + 105 XMT=XMT+ABS(XM(I)) + IF(XMT.LE.ET) GOTO 201 + PRINT 1002,XMT,ET + STOP +C +C THE PARAMETER VALUES ARE NOW ACCEPTED +C +C GENERATE N MASSLESS MOMENTA IN INFINITE PHASE SPACE + 201 DO 202 I=1,N + R1=RN(1) + C=2.*R1-1. + S=SQRT(1.-C*C) + F=TWOPI*RN(2) + R1=RN(3) + R2=RN(4) + Q(4,I)=-LOG(R1*R2) + Q(3,I)=Q(4,I)*C + Q(2,I)=Q(4,I)*S*COS(F) + 202 Q(1,I)=Q(4,I)*S*SIN(F) +C +C CALCULATE THE PARAMETERS OF THE CONFORMAL TRANSFORMATION + DO 203 I=1,4 + 203 R(I)=0. + DO 204 I=1,N + DO 204 K=1,4 + 204 R(K)=R(K)+Q(K,I) + RMAS=SQRT(R(4)**2-R(3)**2-R(2)**2-R(1)**2) + DO 205 K=1,3 + 205 B(K)=-R(K)/RMAS + G=R(4)/RMAS + A=1./(1.+G) + X=ET/RMAS +C +C TRANSFORM THE Q'S CONFORMALLY INTO THE P'S + DO 207 I=1,N + BQ=B(1)*Q(1,I)+B(2)*Q(2,I)+B(3)*Q(3,I) + DO 206 K=1,3 + 206 P(K,I)=X*(Q(K,I)+B(K)*(Q(4,I)+A*BQ)) + 207 P(4,I)=X*(G*Q(4,I)+BQ) +C +C CALCULATE WEIGHT AND POSSIBLE WARNINGS + WT=PO2LOG + IF(N.NE.2) WT=(2.*N-4.)*LOG(ET)+Z(N) + IF(WT.GE.-180.D0) GOTO 208 + IF(IWARN(1).LE.5) PRINT 1004,WT + IWARN(1)=IWARN(1)+1 + 208 IF(WT.LE. 174.D0) GOTO 209 + IF(IWARN(2).LE.5) PRINT 1005,WT + IWARN(2)=IWARN(2)+1 +C +C RETURN FOR WEIGHTED MASSLESS MOMENTA + 209 IF(NM.NE.0) GOTO 210 +C RETURN LOG OF WEIGHT + WT=WT + RETURN +C +C MASSIVE PARTICLES: RESCALE THE MOMENTA BY A FACTOR X + 210 XMAX=SQRT(1.-(XMT/ET)**2) + DO 301 I=1,N + XM2(I)=XM(I)**2 + 301 P2(I)=P(4,I)**2 + ITER=0 + X=XMAX + ACCU=ET*ACC + 302 F0=-ET + G0=0. + X2=X*X + DO 303 I=1,N + E(I)=SQRT(XM2(I)+X2*P2(I)) + F0=F0+E(I) + 303 G0=G0+P2(I)/E(I) + IF(ABS(F0).LE.ACCU) GOTO 305 + ITER=ITER+1 + IF(ITER.LE.ITMAX) GOTO 304 + PRINT 1006,ITMAX + GOTO 305 + 304 X=X-F0/(X*G0) + GOTO 302 + 305 DO 307 I=1,N + V(I)=X*P(4,I) + DO 306 K=1,3 + 306 P(K,I)=X*P(K,I) + 307 P(4,I)=E(I) +C +C CALCULATE THE MASS-EFFECT WEIGHT FACTOR + WT2=1. + WT3=0. + DO 308 I=1,N + WT2=WT2*V(I)/E(I) + 308 WT3=WT3+V(I)**2/E(I) + WTM=(2.*N-3.)*LOG(X)+LOG(WT2/WT3*ET) +C +C RETURN FOR WEIGHTED MASSIVE MOMENTA + WT=WT+WTM + IF(WT.GE.-180.D0) GOTO 309 + IF(IWARN(3).LE.5) PRINT 1004,WT + IWARN(3)=IWARN(3)+1 + 309 IF(WT.LE. 174.D0) GOTO 310 + IF(IWARN(4).LE.5) PRINT 1005,WT + IWARN(4)=IWARN(4)+1 +C RETURN LOG OF WEIGHT + 310 WT=WT + RETURN +C + 1001 FORMAT(' RAMBO FAILS: # OF PARTICLES =',I5,' IS NOT ALLOWED') + 1002 FORMAT(' RAMBO FAILS: TOTAL MASS =',D15.6,' IS NOT',' SMALLER' + $ //' THAN TOTAL ENERGY =',D15.6) + 1004 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY UNDERFLOW') + 1005 FORMAT(' RAMBO WARNS: WEIGHT = EXP(',F20.9,') MAY OVERFLOW') + 1006 FORMAT(' RAMBO WARNS:',I3,' ITERATIONS DID NOT GIVE THE', + $ ' DESIRED ACCURACY =',D15.6) + END + + FUNCTION RN(IDUMMY) + REAL*8 RN,RAN + SAVE INIT + DATA INIT /1/ + IF (INIT.EQ.1) THEN + INIT=0 + CALL RMARIN(1802,9373) + END IF +C + 10 CALL RANMAR(RAN) + IF (RAN.LT.1D-16) GOTO 10 + RN=RAN +C + END + + + + SUBROUTINE RANMAR(RVEC) +C ----------------- +C Universal random number generator proposed by Marsaglia and Zaman +C in report FSU-SCRI-87-50 +C In this version RVEC is a double precision variable. + IMPLICIT REAL*8(A-H,O-Z) + COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM + COMMON/ RASET2 / IRANMR,JRANMR + SAVE /RASET1/,/RASET2/ + UNI = RANU(IRANMR) - RANU(JRANMR) + IF(UNI .LT. 0D0) UNI = UNI + 1D0 + RANU(IRANMR) = UNI + IRANMR = IRANMR - 1 + JRANMR = JRANMR - 1 + IF(IRANMR .EQ. 0) IRANMR = 97 + IF(JRANMR .EQ. 0) JRANMR = 97 + RANC = RANC - RANCD + IF(RANC .LT. 0D0) RANC = RANC + RANCM + UNI = UNI - RANC + IF(UNI .LT. 0D0) UNI = UNI + 1D0 + RVEC = UNI + END + + SUBROUTINE RMARIN(IJ,KL) +C ----------------- +C Initializing routine for RANMAR, must be called before generating +C any pseudorandom numbers with RANMAR. The input values should be +C in +C the ranges 0<=ij<=31328 ; 0<=kl<=30081 + IMPLICIT REAL*8(A-H,O-Z) + COMMON/ RASET1 / RANU(97),RANC,RANCD,RANCM + COMMON/ RASET2 / IRANMR,JRANMR + SAVE /RASET1/,/RASET2/ +C This shows correspondence between the simplified input seeds IJ, +C KL +C and the original Marsaglia-Zaman seeds I,J,K,L. +C To get the standard values in the Marsaglia-Zaman paper +C (i=12,j=34 +C k=56,l=78) put ij=1802, kl=9373 + I = MOD( IJ/177 , 177 ) + 2 + J = MOD( IJ , 177 ) + 2 + K = MOD( KL/169 , 178 ) + 1 + L = MOD( KL , 169 ) + DO 300 II = 1 , 97 + S = 0D0 + T = .5D0 + DO 200 JJ = 1 , 24 + M = MOD( MOD(I*J,179)*K , 179 ) + I = J + J = K + K = M + L = MOD( 53*L+1 , 169 ) + IF(MOD(L*M,64) .GE. 32) S = S + T + T = .5D0*T + 200 CONTINUE + RANU(II) = S + 300 CONTINUE + RANC = 362436D0 / 16777216D0 + RANCD = 7654321D0 / 16777216D0 + RANCM = 16777213D0 / 16777216D0 + IRANMR = 97 + JRANMR = 33 + END + + + + + + + diff --git a/tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%loop_matrix.f b/tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%loop_matrix.f new file mode 100644 index 000000000..1e18ff084 --- /dev/null +++ b/tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/loop_induced_standalone_output/%gghLI_IOTest%SubProcesses%P0_gg_h%loop_matrix.f @@ -0,0 +1,3180 @@ +C --=========================================-- +C Main subroutine +C --=========================================-- + + SUBROUTINE ML5_0_SLOOPMATRIX(P_USER,ANS) +C +C Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s +C By the MadGraph5_aMC@NLO Development Team +C Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch +C +C Returns amplitude squared summed/avg over colors +C and helicities for the point in phase space P(0:3,NEXTERNAL) +C and external lines W(0:6,NEXTERNAL) +C +C Process: g g > h WEIGHTED<=4 [ noborn = QCD ] +C +C Modules +C + USE ML5_0_POLYNOMIAL_CONSTANTS + USE ALOHA_OBJECT +C + IMPLICIT NONE +C +C USER CUSTOMIZABLE OPTIONS +C +C The variables below are just used in the context of a JAMP +C consistency check turned off by default. + REAL*8 JAMP_DOUBLECHECK_THRES + PARAMETER (JAMP_DOUBLECHECK_THRES=1.0D-9) + LOGICAL DIRECT_ME_COMPUTATION, ME_COMPUTATION_FROM_JAMP +C Modify the logicals below to chose how the ME must be computed +C DIRECT_ME_COMPUTATION = Each loop amplitude is squared +C individually against all amplitudes with its own color factor. +C ME_COMPUTATION_FROM_JAMP = Amplitudes are first projected onto +C color flows (many less of them) which are then squared to form +C the ME. +C When setting both computation method to .TRUE., their systematic +C comparisons will be printed out. + DATA DIRECT_ME_COMPUTATION/.FALSE./ +C When using this MadLoop output for integration with MadEvent, +C ME_COMPUTATION_FROM_JAMP *must* be set to .True. because it is +C necessary to compute the AMP2 setting up the multichanneling. + DATA ME_COMPUTATION_FROM_JAMP/.TRUE./ +C This parameter is designed for the check timing command of MG5. +C It skips the loop reduction. + LOGICAL SKIPLOOPEVAL + PARAMETER (SKIPLOOPEVAL=.FALSE.) +C For timing checks. Stops the code after having only initialized +C its arrays from the external data files + LOGICAL BOOTANDSTOP + PARAMETER (BOOTANDSTOP=.FALSE.) + INTEGER TIR_CACHE_SIZE +C To change memory foot-print of MadLoop, you can change this +C parameter to be 0,1 or 2 *and recompile*. +C Notice that this will impact MadLoop speed performances in the +C context of stability checks. + INCLUDE 'tir_cache_size.inc' +C +C CONSTANTS +C + CHARACTER*512 PARAMFNAME,HELCONFIGFNAME,LOOPFILTERFNAME + CHARACTER*512 COLORNUMFNAME,COLORDENOMFNAME, HELFILTERFNAME + CHARACTER*512 PROC_PREFIX + PARAMETER ( PARAMFNAME='MadLoopParams.dat') + PARAMETER ( HELCONFIGFNAME='HelConfigs.dat') + PARAMETER ( LOOPFILTERFNAME='LoopFilter.dat') + PARAMETER ( HELFILTERFNAME='HelFilter.dat') + PARAMETER ( COLORNUMFNAME='ColorNumFactors.dat') + PARAMETER ( COLORDENOMFNAME='ColorDenomFactors.dat') + PARAMETER ( PROC_PREFIX='ML5_0_') + + INTEGER NLOOPS, NLOOPGROUPS, NCTAMPS + PARAMETER (NLOOPS=4, NLOOPGROUPS=4, NCTAMPS=2) + INTEGER NLOOPAMPS + PARAMETER (NLOOPAMPS=6) + INTEGER NCOLORROWS + PARAMETER (NCOLORROWS=NLOOPAMPS) + INTEGER NEXTERNAL + PARAMETER (NEXTERNAL=3) + INTEGER NINITIAL + PARAMETER (NINITIAL=2) + INTEGER NWAVEFUNCS,NLOOPWAVEFUNCS + PARAMETER (NWAVEFUNCS=3,NLOOPWAVEFUNCS=12) + INTEGER NCOMB + PARAMETER (NCOMB=4) + REAL*8 ZERO + PARAMETER (ZERO=0D0) + REAL*16 MP__ZERO + PARAMETER (MP__ZERO=0E0_16) + COMPLEX*16 IMAG1 + PARAMETER (IMAG1=(0D0,1D0)) +C These are constants related to the split orders + INTEGER NSQSO_BORN + PARAMETER (NSQSO_BORN=0) + + INTEGER NSO, NSQUAREDSO, NAMPSO + PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=1) + INTEGER ANS_DIMENSION + PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO)) + INTEGER NSQSOXNLG + PARAMETER (NSQSOXNLG=NSQUAREDSO*NLOOPGROUPS) + INTEGER NSQUAREDSOP1 + PARAMETER (NSQUAREDSOP1=NSQUAREDSO+1) +C The total number of loop reduction libraries +C At present, there are only +C CutTools,PJFry++,IREGI,Golem95,Samurai, Ninja and COLLIER + INTEGER NLOOPLIB + PARAMETER (NLOOPLIB=7) +C Only CutTools or possibly Ninja (if installed with qp support) +C provide QP + INTEGER QP_NLOOPLIB + PARAMETER (QP_NLOOPLIB=1) + INTEGER MAXSTABILITYLENGTH + DATA MAXSTABILITYLENGTH/20/ + COMMON/ML5_0_STABILITY_TESTS/MAXSTABILITYLENGTH + +C +C ARGUMENTS +C + REAL*8 P_USER(0:3,NEXTERNAL) +C +C The zeroth component of the second dimension is the result +C summed over all +C contributing split orders. The zeroth component of the first one +C is the Born. +C Notice that the upper bound of the second integer is not number +C of squared orders +C combination for the loops but the maximum between this number +C for the Born +C contributions and the loop ones. There are some cases for which +C the Born contrib. +C has squared split order contributions than the loop does. For +C example +C +C generate u u~ > d d~ QCD^2<=2 QED^2<=99 [virt=QCD] +C +C It is however somehow academical. This is why ANS_DIMENSION is +C not just NSQSO but rather MAX(NSQSO,NSQSO_BORN) +C + REAL*8 ANS(0:3,0:ANS_DIMENSION) +C +C LOCAL VARIABLES +C + INTEGER I,J,K,L,H,HEL_MULT,I_QP_LIB,DUMMY, INDEX_H + + CHARACTER*512 PARAMFN,HELCONFIGFN,LOOPFILTERFN,COLORNUMFN + $ ,COLORDENOMFN,HELFILTERFN + CHARACTER*512 TMP + SAVE PARAMFN + SAVE HELCONFIGFN + SAVE LOOPFILTERFN + SAVE COLORNUMFN + SAVE COLORDENOMFN + SAVE HELFILTERFN + + INTEGER CTMODEINIT_BU + REAL*8 MLSTABTHRES_BU + INTEGER NEWHELREF + LOGICAL HEL_INCONSISTENT + REAL*8 P(0:3,NEXTERNAL) +C DP_RES STORES THE DOUBLE PRECISION RESULT OBTAINED FROM +C DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY. +C THE STAB_STAGE COUNTER I CORRESPONDANCE GOES AS FOLLOWS +C I=1 -> ORIGINAL PS, CTMODE=1 +C I=2 -> ORIGINAL PS, CTMODE=2, (ONLY WITH CTMODERUN=-1) +C I=3 -> PS WITH ROTATION 1, CTMODE=1, (ONLY WITH CTMODERUN=-2) +C I=4 -> PS WITH ROTATION 2, CTMODE=1, (ONLY WITH CTMODERUN=-3) +C I=5 -> POSSIBLY MORE EVALUATION METHODS IN THE FUTURE, MAX IS +C MAXSTABILITYLENGTH +C IF UNSTABLE IT GOES TO THE SAME PATTERN BUT STAB_INDEX IS THEN +C I+20. + LOGICAL EVAL_DONE(MAXSTABILITYLENGTH) + LOGICAL DOING_QP_EVALS + INTEGER STAB_INDEX,BASIC_CT_MODE + +C This is used for loop-induced where the reference scale for +C comparisons is inferred from the first 100 points at most +C (notice that the weight of a given kinematic configuration can +C appear more than once because of the stability tests). +C When changing this parameter, make sure to correspondingly +C update the parameter with the same name in MadLoopCommons.f. + INTEGER MAXNREF_EVALS + PARAMETER (MAXNREF_EVALS=100) + REAL*8 REF_EVALS(MAXNREF_EVALS) + DATA REF_EVALS/MAXNREF_EVALS*ZERO/ + INTEGER NPSPOINTS + DATA NPSPOINTS/0/ + + REAL*8 ACC(0:NSQUAREDSO) + REAL*8 DP_RES(3,0:NSQUAREDSO,MAXSTABILITYLENGTH) +C QP_RES STORES THE QUADRUPLE PRECISION RESULT OBTAINED FROM +C DIFFERENT EVALUATION METHODS IN ORDER TO ASSESS STABILITY. + REAL*8 QP_RES(3,0:NSQUAREDSO,MAXSTABILITYLENGTH) + INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL) + INTEGER NATTEMPTS + DATA NATTEMPTS/0/ + DATA IC/NEXTERNAL*1/ + REAL*8 HELSAVED(3,NCOMB) + INTEGER ITEMP + LOGICAL LTEMP + REAL*8 BORNBUFF(0:NSQSO_BORN),TMPR + REAL*8 BUFFR(3,0:NSQUAREDSO),BUFFR_BIS(3,0:NSQUAREDSO),TEMP(0:3 + $ ,0:NSQUAREDSO),TEMP1(0:NSQUAREDSO) + COMPLEX*16 CTEMP + REAL*8 TEMP2(3) + REAL*8 BUFFRES(0:3,0:NSQUAREDSO) + COMPLEX*16 COEFS(MAXLWFSIZE,0:VERTEXMAXCOEFS-1,MAXLWFSIZE) + COMPLEX*16 CFTOT + LOGICAL FOUNDHELFILTER,FOUNDLOOPFILTER + DATA FOUNDHELFILTER/.TRUE./ + DATA FOUNDLOOPFILTER/.TRUE./ + LOGICAL LOOPFILTERBUFF(NSQUAREDSO,NLOOPGROUPS) + DATA ((LOOPFILTERBUFF(J,I),J=1,NSQUAREDSO),I=1,NLOOPGROUPS) + $ /NSQSOXNLG*.FALSE./ + + LOGICAL AUTOMATIC_CACHE_CLEARING + DATA AUTOMATIC_CACHE_CLEARING/.TRUE./ + COMMON/ML5_0_RUNTIME_OPTIONS/AUTOMATIC_CACHE_CLEARING + + INTEGER IDEN + DATA IDEN/256/ + INTEGER HELAVGFACTOR + DATA HELAVGFACTOR/4/ +C For a 1>N process, them BEAMTWO_HELAVGFACTOR would be set to 1. + INTEGER BEAMS_HELAVGFACTOR(2) + DATA (BEAMS_HELAVGFACTOR(I),I=1,2)/2,2/ + LOGICAL DONEHELDOUBLECHECK + DATA DONEHELDOUBLECHECK/.FALSE./ + INTEGER NEPS + DATA NEPS/0/ +C Below are variables to bypass the checkphase and insure +C stability check to take place + LOGICAL OLD_CHECKPHASE, OLD_HELDOUBLECHECKED + INTEGER OLD_GOODHEL(NCOMB) + LOGICAL OLD_GOODAMP(NSQUAREDSO,NLOOPGROUPS) + LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY + COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY +C +C FUNCTIONS +C + INTEGER ML5_0_TIRCACHE_INDEX + INTEGER ML5_0_ML5SOINDEX_FOR_BORN_AMP + INTEGER ML5_0_ML5SOINDEX_FOR_LOOP_AMP + INTEGER ML5_0_ML5SQSOINDEX + INTEGER ML5_0_ISSAME + LOGICAL ML5_0_ISZERO + LOGICAL ML5_0_IS_HEL_SELECTED + INTEGER SET_RET_CODE_U + REAL*8 MEDIAN +C +C GLOBAL VARIABLES +C + INCLUDE 'process_info.inc' + INCLUDE 'unique_id.inc' + + INCLUDE 'coupl.inc' + INCLUDE 'mp_coupl.inc' + INCLUDE 'MadLoopParams.inc' + + REAL*8 RES_FROM_JAMP(0:3,0:NSQUAREDSO) + COMMON/ML5_0_DOUBLECHECK_JAMP/RES_FROM_JAMP + $ ,DIRECT_ME_COMPUTATION,ME_COMPUTATION_FROM_JAMP + + LOGICAL CHOSEN_SO_CONFIGS(NSQUAREDSO) + DATA CHOSEN_SO_CONFIGS/.TRUE./ + COMMON/ML5_0_CHOSEN_LOOP_SQSO/CHOSEN_SO_CONFIGS + + INTEGER N_DP_EVAL, N_QP_EVAL + DATA N_DP_EVAL/1/ + DATA N_QP_EVAL/1/ + COMMON/ML5_0_N_EVALS/N_DP_EVAL,N_QP_EVAL + + LOGICAL CHECKPHASE + DATA CHECKPHASE/.TRUE./ + LOGICAL HELDOUBLECHECKED + DATA HELDOUBLECHECKED/.FALSE./ + COMMON/ML5_0_INIT/CHECKPHASE, HELDOUBLECHECKED + INTEGER NTRY + DATA NTRY/0/ + REAL*8 REF + DATA REF/0.0D0/ + + LOGICAL MP_DONE + DATA MP_DONE/.FALSE./ + COMMON/ML5_0_MP_DONE/MP_DONE +C A FLAG TO DENOTE WHETHER THE CORRESPONDING LOOPLIBS ARE +C AVAILABLE OR NOT + LOGICAL LOOPLIBS_AVAILABLE(NLOOPLIB) + DATA LOOPLIBS_AVAILABLE/.TRUE.,.FALSE.,.TRUE.,.FALSE.,.FALSE. + $ ,.TRUE.,.TRUE./ + COMMON/ML5_0_LOOPLIBS_AV/ LOOPLIBS_AVAILABLE +C A FLAG TO DENOTE WHETHER THE CORRESPONDING DIRECTION TESTS +C AVAILABLE OR NOT IN THE LOOPLIBS + LOGICAL LOOPLIBS_DIRECTEST(NLOOPLIB) + DATA LOOPLIBS_DIRECTEST /.TRUE.,.TRUE.,.TRUE.,.TRUE.,.TRUE. + $ ,.TRUE.,.TRUE./ +C Specifying for which reduction tool quadruple precision is +C available. +C The index 0 is dummy and simply means that the corresponding +C loop_library is not available +C in which case neither is its quadruple precision version. + LOGICAL LOOPLIBS_QPAVAILABLE(0:7) + DATA LOOPLIBS_QPAVAILABLE /.FALSE.,.TRUE.,.FALSE.,.FALSE. + $ ,.FALSE.,.FALSE.,.FALSE.,.FALSE./ +C PS CAN POSSIBILY BE PASSED THROUGH IMPROVE_PS BUT IS NOT +C MODIFIED FOR THE PURPOSE OF THE STABILITY TEST +C EVEN THOUGH THEY ARE PUT IN COMMON BLOCK, FOR NOW THEY ARE NOT +C USED ANYWHERE ELSE + REAL*8 PS(0:3,NEXTERNAL) + COMMON/ML5_0_PSPOINT/PS +C AGAIN BELOW, MP_PS IS THE FIXED (POSSIBLY IMPROVED) MP PS POINT +C AND MP_P IS THE ONE WHICH CAN BE MODIFIED (I.E. ROTATED ETC.) +C FOR STABILITY PURPOSE + REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL) + COMMON/ML5_0_MP_PSPOINT/MP_PS,MP_P + + REAL*8 LSCALE + INTEGER CTMODE + COMMON/ML5_0_CT/LSCALE,CTMODE + LOGICAL MP_PS_SET + DATA MP_PS_SET/.FALSE./ + +C The parameter below sets the convention for the helicity filter +C For a given helicity, the attached integer 'i' means +C 'i' in ]-inf;-HELOFFSET[ -> Helicity is equal, up to a sign, to +C helicity number abs(i+HELOFFSET) +C 'i' == -HELOFFSET -> Helicity is analytically zero +C 'i' in ]-HELOFFSET,inf[ -> Helicity is contributing with weight +C 'i'. If it is zero, it is skipped. +C Typically, the hel_offset is 10000 + INTEGER HELOFFSET + DATA HELOFFSET/10000/ + INTEGER GOODHEL(NCOMB) + LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS) + COMMON/ML5_0_FILTERS/GOODAMP,GOODHEL,HELOFFSET + + INTEGER HELPICKED + DATA HELPICKED/-1/ + COMMON/ML5_0_HELCHOICE/HELPICKED + INTEGER USERHEL + DATA USERHEL/-1/ + COMMON/ML5_0_USERCHOICE/USERHEL + +C This integer can be accessed by an external user to set its +C target squared split order. +C If set to a value different than -1, the code will try to avoid +C computing anything which +C does not contribute to contributions of squared split orders +C SQSO_TARGET and below. + INTEGER SQSO_TARGET + DATA SQSO_TARGET/-1/ + COMMON/ML5_0_SOCHOICE/SQSO_TARGET +C The following logical are used to broadcast the fact that the +C target 'required' CT and +C loop split orders contributions have been reached already and +C the rest can be skipped. + LOGICAL UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE,CT_REQ_SO_DONE + $ ,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_SO_DONE + $ ,CTCALL_REQ_SO_DONE,FILTER_SO + DATA UVCT_REQ_SO_DONE/.FALSE./ + DATA MP_UVCT_REQ_SO_DONE/.FALSE./ + DATA CT_REQ_SO_DONE/.FALSE./ + DATA MP_CT_REQ_SO_DONE/.FALSE./ + DATA LOOP_REQ_SO_DONE/.FALSE./ + DATA MP_LOOP_REQ_SO_DONE/.FALSE./ + DATA CTCALL_REQ_SO_DONE/.FALSE./ + DATA FILTER_SO/.FALSE./ + COMMON/ML5_0_SO_REQS/UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE + $ ,CT_REQ_SO_DONE,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE + $ ,MP_LOOP_REQ_SO_DONE,CTCALL_REQ_SO_DONE,FILTER_SO + +C Allows to forbid the zero helicity double check, no matter the +C value in MadLoopParams.dat +C This can be accessed with the SET_FORBID_HEL_DOUBLECHECK +C subroutine of MadLoopCommons.dat + LOGICAL FORBID_HEL_DOUBLECHECK + COMMON/FORBID_HEL_DOUBLECHECK/FORBID_HEL_DOUBLECHECK + + INTEGER I_SO + DATA I_SO/1/ + COMMON/ML5_0_I_SO/I_SO + INTEGER I_LIB + DATA I_LIB/1/ + COMMON/ML5_0_I_LIB/I_LIB + LOGICAL QP_TOOLS_AVAILABLE + DATA QP_TOOLS_AVAILABLE/.FALSE./ + INTEGER INDEX_QP_TOOLS(QP_NLOOPLIB+1) + COMMON/ML5_0_LOOP_TOOLS/QP_TOOLS_AVAILABLE,INDEX_QP_TOOLS + + TYPE(ALOHA) W(NWAVEFUNCS) + COMMON/ML5_0_W/W + + TYPE(MP_ALOHA) MPW(NWAVEFUNCS) + COMMON/ML5_0_MP_W/MPW + + COMPLEX*16 WL(MAXLWFSIZE,0:LOOPMAXCOEFS-1,MAXLWFSIZE, + $ -1:NLOOPWAVEFUNCS) + COMPLEX*16 PL(0:3,-1:NLOOPWAVEFUNCS) + COMMON/ML5_0_WL/WL,PL + + COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NLOOPGROUPS) + COMMON/ML5_0_LCOEFS/LOOPCOEFS + +C This flag is used to prevent the re-computation of the OpenLoop +C coefficients when changing the CTMode for the stability test. + LOGICAL SKIP_LOOPNUM_COEFS_CONSTRUCTION + DATA SKIP_LOOPNUM_COEFS_CONSTRUCTION/.FALSE./ + COMMON/ML5_0_SKIP_COEFS/SKIP_LOOPNUM_COEFS_CONSTRUCTION + + LOGICAL TIR_DONE(NLOOPGROUPS) + COMMON/ML5_0_TIRCACHING/TIR_DONE + + COMPLEX*16 AMPL(3,NLOOPAMPS) + COMMON/ML5_0_AMPL/AMPL + + COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS) + LOGICAL S(NSQUAREDSO,NLOOPGROUPS) + COMMON/ML5_0_LOOPRES/LOOPRES,S + + INTEGER CF_D(NCOLORROWS,NLOOPAMPS) + INTEGER CF_N(NCOLORROWS,NLOOPAMPS) + COMMON/ML5_0_CF/CF_D,CF_N + + INTEGER HELC(NEXTERNAL,NCOMB) + COMMON/ML5_0_HELCONFIGS/HELC + + REAL*8 PREC,USER_STAB_PREC + DATA USER_STAB_PREC/-1.0D0/ + COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC + +C Return codes H,T,U correspond to the hundreds, tens and units +C building returncode, i.e. +C RETURNCODE=100*RET_CODE_H+10*RET_CODE_T+RET_CODE_U + + INTEGER RET_CODE_H,RET_CODE_T,RET_CODE_U + REAL*8 ACCURACY(0:NSQUAREDSO) + DATA (ACCURACY(I),I=0,NSQUAREDSO)/NSQUAREDSOP1*1.0D0/ + DATA RET_CODE_H,RET_CODE_T,RET_CODE_U/1,1,0/ + COMMON/ML5_0_ACC/ACCURACY,RET_CODE_H,RET_CODE_T,RET_CODE_U + + LOGICAL MP_DONE_ONCE + DATA MP_DONE_ONCE/.FALSE./ + COMMON/ML5_0_MP_DONE_ONCE/MP_DONE_ONCE + + CHARACTER(512) MLPATH + COMMON/MLPATH/MLPATH + +C This is just so that if the user disabled the computation of +C poles by COLLIER +C using the MadLoop subroutine, we don't overwrite his choice when +C reading the parameters + LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, + $ FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION + LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, + $ COLLIER_IR_POLE_COMPUTATION_CHOICE + DATA FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION + $ ,FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION/.FALSE.,.FALSE./ + COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE + $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, + $ FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION + $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE + $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE + +C This variable controls the general initialization which is +C *common* between all MadLoop SubProcesses. +C For example setting the MadLoopPath or reading the ML runtime +C parameters. + LOGICAL ML_INIT + COMMON/ML_INIT/ML_INIT + +C This variable controls the *local* initialization of this +C particular SubProcess. +C For example, the reading of the filters must be done +C independently by each SubProcess. + LOGICAL LOCAL_ML_INIT + DATA LOCAL_ML_INIT/.TRUE./ + + LOGICAL WARNED_LORENTZ_STAB_TEST_OFF + DATA WARNED_LORENTZ_STAB_TEST_OFF/.FALSE./ + INTEGER NROTATIONS_DP_BU,NROTATIONS_QP_BU + + LOGICAL FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION + DATA FPE_IN_DP_REDUCTION, FPE_IN_QP_REDUCTION/.FALSE.,.FALSE./ + COMMON/ML5_0_FPE_IN_REDUCTION/FPE_IN_DP_REDUCTION, + $ FPE_IN_QP_REDUCTION + +C This array specify potential special requirements on the +C helicities to +C consider. POLARIZATIONS(0,0) is -1 if there is not such +C requirement. + INTEGER POLARIZATIONS(0:NEXTERNAL,0:5) + COMMON/ML5_0_BEAM_POL/POLARIZATIONS + +C ---------- +C BEGIN CODE +C ---------- + + IF(ML_INIT) THEN + ML_INIT = .FALSE. + CALL PRINT_MADLOOP_BANNER() + TMP = 'auto' + CALL SETMADLOOPPATH(TMP) + CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN) + CALL MADLOOPPARAMREADER(PARAMFN,.TRUE.) + IF (FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION) THEN + COLLIERCOMPUTEUVPOLES = COLLIER_UV_POLE_COMPUTATION_CHOICE + ENDIF + IF (FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION) THEN + COLLIERCOMPUTEIRPOLES = COLLIER_IR_POLE_COMPUTATION_CHOICE + ENDIF + IF (FORBID_HEL_DOUBLECHECK) THEN + DOUBLECHECKHELICITYFILTER = .FALSE. + ENDIF + +C Make sure that HELFILTERLEVEL is at most 1 if the beam is +C polarized + IF (POLARIZATIONS(0,0).EQ.0) THEN + IF (HELICITYFILTERLEVEL.GT.1) THEN + WRITE(*,*) '##INFO: When using polarized beam, the' + $ //' helicity filter of MadLoop can be at most 1. Now' + $ //' setting HELICITYFILTERLEVEL to 1.' + HELICITYFILTERLEVEL = 1 + ENDIF + ENDIF + + IF(.NOT.LOOPINITSTARTOVER) THEN + WRITE(*,*) '##INFO: For loop-induced processes it is' + $ //' preferable to always set the parameter' + $ //' LoopInitStartOver to True, so it is hard-set here to' + $ //' True.' + LOOPINITSTARTOVER=.TRUE. + ENDIF + IF(.NOT.HELINITSTARTOVER) THEN + WRITE(*,*) '##INFO: For loop-induced processes it is' + $ //' preferable to always set the parameter HelInitStartOver' + $ //' to True, so it is hard-set here to True.' + HELINITSTARTOVER=.TRUE. + ENDIF + IF (CHECKCYCLE.LT.5) THEN + WRITE(*,*) '##INFO: Due to the dynamic setting of the' + $ //' reference scale for contributions comparisons, it is' + $ //' preferable to set the parameter CheckCycle to a value' + $ //' larger than 4, so it is hard-set here to 5.' + CHECKCYCLE=5 + ENDIF + +C Make sure that NROTATIONS_QP and NROTATIONS_DP are set to zero +C if AUTOMATIC_CACHE_CLEARING is disabled. + IF(.NOT.AUTOMATIC_CACHE_CLEARING) THEN + IF(NROTATIONS_DP.NE.0.OR.NROTATIONS_QP.NE.0) THEN + WRITE(*,*) '##INFO: AUTOMATIC_CACHE_CLEARING is disabled,' + $ //' so MadLoop automatically resets NROTATIONS_DP and' + $ //' NROTATIONS_QP to 0.' + NROTATIONS_QP=0 + NROTATIONS_DP=0 + ENDIF + ENDIF + + ENDIF + + IF (LOCAL_ML_INIT) THEN + LOCAL_ML_INIT = .FALSE. + QP_TOOLS_AVAILABLE=.FALSE. + INDEX_QP_TOOLS(1:QP_NLOOPLIB+1)=0 +C SKIP THE ONES THAT NOT AVAILABLE + J=1 + DO I=1,NLOOPLIB + IF(MLREDUCTIONLIB(J).EQ.0)EXIT + IF(.NOT.LOOPLIBS_AVAILABLE(MLREDUCTIONLIB(J)))THEN + MLREDUCTIONLIB(J:NLOOPLIB-1)=MLREDUCTIONLIB(J+1:NLOOPLIB) + MLREDUCTIONLIB(NLOOPLIB)=0 + ELSE + J=J+1 + ENDIF + ENDDO + IF(MLREDUCTIONLIB(1).EQ.0)THEN + STOP 'No available loop reduction lib is provided. Make sure' + $ //' MLReductionLib is correct.' + ENDIF + J=0 + DO I=1,NLOOPLIB + IF(LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I)))THEN + J=J+1 + IF(.NOT.QP_TOOLS_AVAILABLE) THEN + QP_TOOLS_AVAILABLE=.TRUE. + ENDIF + INDEX_QP_TOOLS(J)=I + ENDIF + ENDDO + +C Setup the file paths + CALL JOINPATH(MLPATH,PARAMFNAME,PARAMFN) + CALL JOINPATH(MLPATH,PROC_PREFIX,TMP) + CALL JOINPATH(TMP,HELCONFIGFNAME,HELCONFIGFN) + CALL JOINPATH(TMP,LOOPFILTERFNAME,LOOPFILTERFN) + CALL JOINPATH(TMP,COLORNUMFNAME,COLORNUMFN) + CALL JOINPATH(TMP,COLORDENOMFNAME,COLORDENOMFN) + CALL JOINPATH(TMP,HELFILTERFNAME,HELFILTERFN) + + CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL) + +C Make sure that the loop filter is disabled when there is +C spin-2 particles for 2>1 or 1>2 processes + IF(MAX_SPIN_EXTERNAL_PARTICLE.GT.3.AND.(NEXTERNAL.LE.3.AND.HELI + $CITYFILTERLEVEL.NE.0)) THEN + WRITE(*,*) '##INFO: Helicity filter deactivated for 2>1' + $ //' processes involving spin 2 particles.' + HELICITYFILTERLEVEL = 0 +C We write a dummy filter for structural reasons here + OPEN(1, FILE=HELFILTERFN, ERR=6116, STATUS='NEW' + $ ,ACTION='WRITE') + DO I=1,NCOMB + WRITE(1,*) 1 + ENDDO + 6116 CONTINUE + CLOSE(1) + ENDIF + + OPEN(1, FILE=COLORNUMFN, ERR=104, STATUS='OLD', + $ ACTION='READ') + DO I=1,NCOLORROWS + READ(1,*,END=105) (CF_N(I,J),J=1,NLOOPAMPS) + ENDDO + GOTO 105 + 104 CONTINUE + STOP 'Color factors could not be initialized from file' + $ //' ML5_0_ColorNumFactors.dat. File not found' + 105 CONTINUE + CLOSE(1) + OPEN(1, FILE=COLORDENOMFN, ERR=106, STATUS='OLD', + $ ACTION='READ') + DO I=1,NCOLORROWS + READ(1,*,END=107) (CF_D(I,J),J=1,NLOOPAMPS) + ENDDO + GOTO 107 + 106 CONTINUE + STOP 'Color factors could not be initialized from file' + $ //' ML5_0_ColorDenomFactors.dat. File not found' + 107 CONTINUE + CLOSE(1) + OPEN(1, FILE=HELCONFIGFN, ERR=108, STATUS='OLD', + $ ACTION='READ') + DO H=1,NCOMB + READ(1,*,END=109) (HELC(I,H),I=1,NEXTERNAL) + ENDDO + GOTO 109 + 108 CONTINUE + STOP 'Color helictiy configurations could not be initialized' + $ //' from file ML5_0_HelConfigs.dat. File not found' + 109 CONTINUE + CLOSE(1) + +C SETUP OF THE COMMON STARTING EXTERNAL LOOP WAVEFUNCTION +C IT IS ALSO PS POINT INDEPENDENT, SO IT CAN BE DONE HERE. +C The index -1 is for the charge-conjugated fermions with +C flipped fermion flow. + DO I=0,3 + PL(I,-1)=DCMPLX(0.0D0,0.0D0) + PL(I,0)=DCMPLX(0.0D0,0.0D0) + ENDDO + DO I=1,MAXLWFSIZE + DO J=0,LOOPMAXCOEFS-1 + DO K=1,MAXLWFSIZE + WL(I,J,K,-1)=(0.0D0,0.0D0) + IF(I.EQ.K.AND.J.EQ.0) THEN + WL(I,J,K,0)=(1.0D0,0.0D0) + ELSE + WL(I,J,K,0)=(0.0D0,0.0D0) + ENDIF + ENDDO + ENDDO + ENDDO + IF(BOOTANDSTOP) THEN + WRITE(*,*) '##Stopped by user request.' + STOP + ENDIF + ENDIF + +C This is the chare conjugate version of the unit 4-currents in +C the canonical cartesian basis. +C This, for now, is only defined for 4-fermionic currents. + WL(1,0,2,-1) = DCMPLX(-1.0D0,0.0D0) + WL(2,0,1,-1) = DCMPLX(1.0D0,0.0D0) + WL(3,0,4,-1) = DCMPLX(1.0D0,0.0D0) + WL(4,0,3,-1) = DCMPLX(-1.0D0,0.0D0) + +C Make sure that lorentz rotation tests are not used if there is +C external loop wavefunction of spin 2 and that one specific +C helicity is asked + NROTATIONS_DP_BU = NROTATIONS_DP + NROTATIONS_QP_BU = NROTATIONS_QP + IF(MAX_SPIN_EXTERNAL_PARTICLE.GT.3.AND.USERHEL.NE.-1) THEN + IF(.NOT.WARNED_LORENTZ_STAB_TEST_OFF) THEN + WRITE(*,*) '##WARNING: Evaluation of a specific helicity was' + $ //' asked for this PS point, and there is a spin-2 (or' + $ //' higher) particle in the external states.' + WRITE(*,*) '##WARNING: As a result, MadLoop disabled the' + $ //' Lorentz rotation test for this phase-space point only.' + WRITE(*,*) '##WARNING: Further warning of that type' + $ //' suppressed.' + WARNED_LORENTZ_STAB_TEST_OFF = .TRUE. + ENDIF + NROTATIONS_QP=0 + NROTATIONS_DP=0 + CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL) + ENDIF + + IF(NTRY.EQ.0) THEN + HELDOUBLECHECKED=(.NOT.DOUBLECHECKHELICITYFILTER) + $ .OR.(HELICITYFILTERLEVEL.EQ.0) + OPEN(1, FILE=LOOPFILTERFN, ERR=100, STATUS='OLD', + $ ACTION='READ') + DO J=1,NLOOPGROUPS + READ(1,*,END=101) (GOODAMP(I,J),I=1,NSQUAREDSO) + ENDDO + GOTO 101 + 100 CONTINUE + FOUNDLOOPFILTER=.FALSE. + DO J=1,NLOOPGROUPS + DO I=1,NSQUAREDSO + GOODAMP(I,J)=(.NOT.USELOOPFILTER) + ENDDO + ENDDO + 101 CONTINUE + CLOSE(1) + + IF (.NOT.USELOOPFILTER) THEN + DO J=1,NLOOPGROUPS + DO I=1,NSQUAREDSO + GOODAMP(I,J)=.TRUE. + ENDDO + ENDDO + ENDIF + + IF (HELICITYFILTERLEVEL.EQ.0) THEN + FOUNDHELFILTER=.TRUE. + DO J=1,NCOMB + GOODHEL(J)=1 + ENDDO + GOTO 122 + ENDIF + OPEN(1, FILE=HELFILTERFN, ERR=102, STATUS='OLD', + $ ACTION='READ') + DO I=1,NCOMB + READ(1,*,END=103) GOODHEL(I) + ENDDO + GOTO 103 + 102 CONTINUE + FOUNDHELFILTER=.FALSE. + DO J=1,NCOMB + GOODHEL(J)=1 + ENDDO + 103 CONTINUE + CLOSE(1) + IF (HELICITYFILTERLEVEL.EQ.1) THEN +C We must make sure to remove the matching-helicity +C optimisation, as requested by the user. + DO J=1,NCOMB + IF ((GOODHEL(J).GT.1).OR.(GOODHEL(J).LT.-HELOFFSET)) THEN + GOODHEL(J)=1 + ENDIF + ENDDO + ENDIF + 122 CONTINUE + ENDIF + +C The born is of course 0 for loop-induced processes. + DO I=0,NSQUAREDSO + ANS(0,I)=0.0D0 + ENDDO + +C For loop-induced, the reference for comparison is set later from +C the total contribution of the previous PS point considered. +C But you can edit here the value to be used for the first PS +C points. + IF (NPSPOINTS.EQ.0) THEN + REF=1.0D-50 + ELSE + IF(NPSPOINTS.GE.MAXNREF_EVALS) THEN + REF=MEDIAN(REF_EVALS,MAXNREF_EVALS) + ELSE + REF=MEDIAN(REF_EVALS,NPSPOINTS) + ENDIF + ENDIF + + MP_DONE=.FALSE. + MP_DONE_ONCE=.FALSE. + MP_PS_SET=.FALSE. + STAB_INDEX=0 + DOING_QP_EVALS=.FALSE. + EVAL_DONE(1)=.TRUE. + DO I=2,MAXSTABILITYLENGTH + EVAL_DONE(I)=.FALSE. + ENDDO + +C For loop-induced processes, we should make sure not to use the +C first points +C to set the filters because it doesn't have a reasonable REF +C scale yet. + IF(.NOT.BYPASS_CHECK.AND.NPSPOINTS.GE.1) THEN + NTRY=NTRY+1 + ENDIF + + IF (USER_STAB_PREC.GT.0.0D0) THEN + MLSTABTHRES_BU=MLSTABTHRES + MLSTABTHRES=USER_STAB_PREC +C In the initialization, I cannot perform stability test and +C therefore guarantee any precision + CTMODEINIT_BU=CTMODEINIT +C So either one choses quad precision directly +C CTMODEINIT=4 +C Or, because this is very slow, we keep the orignal value. The +C accuracy returned is -1 and tells the MC that he should not +C trust the evaluation for checks. + CTMODEINIT=CTMODEINIT_BU + ENDIF + + IF(DONEHELDOUBLECHECK.AND.(.NOT.HELDOUBLECHECKED)) THEN + HELDOUBLECHECKED=.TRUE. + DONEHELDOUBLECHECK=.FALSE. + ENDIF + + CHECKPHASE=(NTRY.LE.CHECKCYCLE).AND.(((.NOT.FOUNDLOOPFILTER) + $ .AND.USELOOPFILTER).OR.(.NOT.FOUNDHELFILTER)) + + IF (WRITEOUTFILTERS) THEN + IF ((HELICITYFILTERLEVEL.NE.0).AND.(.NOT. CHECKPHASE) + $ .AND.(.NOT.FOUNDHELFILTER)) THEN + OPEN(1, FILE=HELFILTERFN, ERR=110, STATUS='NEW' + $ ,ACTION='WRITE') + DO I=1,NCOMB + WRITE(1,*) GOODHEL(I) + ENDDO + 110 CONTINUE + CLOSE(1) + FOUNDHELFILTER=.TRUE. + ENDIF + + IF ((.NOT. CHECKPHASE).AND.(.NOT.FOUNDLOOPFILTER) + $ .AND.USELOOPFILTER) THEN + OPEN(1, FILE=LOOPFILTERFN, ERR=111, STATUS='NEW' + $ ,ACTION='WRITE') + DO J=1,NLOOPGROUPS + WRITE(1,*) (GOODAMP(I,J),I=1,NSQUAREDSO) + ENDDO + 111 CONTINUE + CLOSE(1) + FOUNDLOOPFILTER=.TRUE. + ENDIF + ENDIF + + IF (BYPASS_CHECK) THEN + OLD_CHECKPHASE = CHECKPHASE + OLD_HELDOUBLECHECKED = HELDOUBLECHECKED + CHECKPHASE = .FALSE. + HELDOUBLECHECKED = .TRUE. + DO I=1,NCOMB + OLD_GOODHEL(I)=GOODHEL(I) + GOODHEL(I)=1 + ENDDO + DO I=1,NSQUAREDSO + DO J=1,NLOOPGROUPS + OLD_GOODAMP(I,J)=GOODAMP(I,J) + GOODAMP(I,J)=.TRUE. + ENDDO + ENDDO + ENDIF + + IF(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN + HELPICKED=1 + CTMODE=CTMODEINIT + ELSE + IF (USERHEL.NE.-1) THEN + IF(GOODHEL(USERHEL).EQ.-HELOFFSET) THEN + DO I=0,NSQUAREDSO + ANS(1,I)=0.0D0 + ANS(2,I)=0.0D0 + ANS(3,I)=0.0D0 + ENDDO + GOTO 9999 + ENDIF + ENDIF + HELPICKED=USERHEL + IF (CTMODERUN.NE.-1) THEN + CTMODE=CTMODERUN + ELSE + CTMODE=1 + ENDIF + ENDIF + + DO I=1,NEXTERNAL + DO J=0,3 + PS(J,I)=P_USER(J,I) + ENDDO + ENDDO + +C Make sure we start with empty caches + IF (AUTOMATIC_CACHE_CLEARING) THEN + CALL ML5_0_CLEAR_CACHES() + ENDIF + +C Now make sure to turn on the global COLLIER cache if applicable + CALL ML5_0_SET_COLLIER_GLOBAL_CACHE(.TRUE.) + + IF (IMPROVEPSPOINT.GE.0) THEN +C Make the input PS more precise (exact onshell and +C energy-momentum conservation) + CALL ML5_0_IMPROVE_PS_POINT_PRECISION(PS) + ENDIF + + DO I=1,NEXTERNAL + DO J=0,3 + P(J,I)=PS(J,I) + ENDDO + ENDDO + + DO K=1, 3 + DO I=0,NSQUAREDSO + BUFFR(K,I)=0.0D0 + ENDDO + DO I=1,NLOOPAMPS + AMPL(K,I)=(0.0D0,0.0D0) + ENDDO + ENDDO + +C Start by using the first available loop reduction library and qp +C library. + I_LIB=1 + I_QP_LIB=1 + + GOTO 208 +C MadLoop jumps to this label during stability checks when it +C recomputes a rotated PS point + 200 CONTINUE +C For the computation of a rotated version of this PS point we +C must reset the all MadLoop cache since this changes the +C definition of the loop denominators. +C We don't check for AUTOMATIC_CACHE_CLEARING here because the +C Lorentz test should anyway be disabled if the flag is turned +C off. + CALL ML5_0_CLEAR_CACHES() + 208 CONTINUE + SKIP_LOOPNUM_COEFS_CONSTRUCTION=.FALSE. + GOTO 308 +C MadLoop jumps to this label during stability checks when it +C recomputes the same PS point with a different CTMode + 300 CONTINUE +C Of course the trick of reusing coefficients when reducing at the +C amplitude level only works when computing one helicity at a time + IF (USERHEL.NE.-1) THEN + SKIP_LOOPNUM_COEFS_CONSTRUCTION = .TRUE. + ENDIF + 308 CONTINUE +C We don't want to re-initialized the following quantities when +C checking the helicity filter. (which jumps to label 205 to +C probe each helicity). +C We however want to re-initialize them for each new computation +C part of the stability check (which jumps to label 200) +C This code is therefore placed before 205 and after 200. + CALL ML5_0_REINITIALIZE_CUMULATIVE_ARRAYS() + IF (ME_COMPUTATION_FROM_JAMP) THEN +C If both ME computational methods have been used, then the ME +C computation from color flows was stored in RES_FROM_JAMP and +C we must reset it here. + DO I=0,NSQUAREDSO + DO K=0,3 + RES_FROM_JAMP(K,I)=0.0D0 + ENDDO + ENDDO + ENDIF + IF ((.NOT.DIRECT_ME_COMPUTATION).AND.ME_COMPUTATION_FROM_JAMP) + $ THEN +C When computing the ME with color flows, the Born ME will be +C computed as well, so we reset here the result obtained from +C the smatrix call above. + DO I=0,NSQUAREDSO + ANS(0,I)=0.0D0 + ENDDO + ENDIF + +C Free cache when using IREGI + IF(IREGIRECY.AND.MLREDUCTIONLIB(I_LIB).EQ.3) THEN + CALL IREGI_FREE_PS() + ENDIF + +C Even if the user did ask to turn off the automatic TIR cache +C clearing, we must do it now if the CTModeIndex rolls over the +C size of the TIR cache employed. +C Notice that we must do that only when processing a new CT mode +C as part of the stability test and not when computing a new +C helicity as part of the filtering process. +C This we check that we are not in the initialization phase. +C If we are not in CTModeRun=-1, then we never need to clear the +C cache since the TIR will always be used for a unique +C computation (not stab test). +C Also, it is clear that if we are running OPP when reaching this' +C //' line, then we shouldn't clear the TIR cache as it might +C still be useful later. +C Finally, notice that the conditional statement below should +C never be true except you have TIR library supporting quadruple +C precision or when TIR_CACHE_SIZE<2. + IF((.NOT.CHECKPHASE.AND.(HELDOUBLECHECKED)).AND.CTMODERUN.EQ. + $ -1.AND.(MLREDUCTIONLIB(I_LIB).NE.1.AND.MLREDUCTIONLIB(I_LIB) + $ .NE.5).AND.(ML5_0_TIRCACHE_INDEX(CTMODE).EQ.(TIR_CACHE_SIZE+1))) + $ THEN + CALL ML5_0_CLEAR_TIR_CACHE() + ENDIF + + +C MadLoop jumps to this label during initialization when it goes +C to the computation of the next helicity. + 205 CONTINUE + + IF (.NOT.MP_PS_SET.AND.(CTMODE.EQ.0.OR.CTMODE.GE.4)) THEN + CALL ML5_0_SET_MP_PS(P_USER) + MP_PS_SET = .TRUE. + ENDIF + + LSCALE=DSQRT(ABS((P(0,1)+P(0,2))**2-(P(1,1)+P(1,2))**2-(P(2,1) + $ +P(2,2))**2-(P(3,1)+P(3,2))**2)) + + CTCALL_REQ_SO_DONE=.FALSE. + FILTER_SO = (.NOT.CHECKPHASE) + $ .AND.HELDOUBLECHECKED.AND.(SQSO_TARGET.NE.-1) + + + DO I=1,NLOOPGROUPS + DO J=1,3 + DO K=1,NSQUAREDSO + LOOPRES(J,K,I)=(0.0D0,0.0D0) + ENDDO + ENDDO + ENDDO + + DO K=1,3 + DO I=0,NSQUAREDSO + ANS(K,I)=0.0D0 + ENDDO + ENDDO + +C Check if we directly go to multiple precision + IF (CTMODE.GE.4) THEN + CALL ML5_0_MP_COMPUTE_LOOP_COEFS(MP_P,BUFFR_BIS) + IF ((.NOT.DIRECT_ME_COMPUTATION).AND.ME_COMPUTATION_FROM_JAMP) + $ THEN +C If the ME's are computed from the color flows only, we must +C update the NLO part of ANS from BUFFR_BIS and the Born part +C of ANS using RES_FROM_JAMP(0,*) + DO I=0,NSQUAREDSO + ANS(0,I)=RES_FROM_JAMP(0,I) + DO K=1,3 + ANS(K,I)=BUFFR_BIS(K,I) + ENDDO + ENDDO + ENDIF +C We must skip the double precision computation of both loop +C amplitudes and CT amplitudes because they will all be +C computed in MP_COMPUTE_LOOP_COEFS. + GOTO 301 + ENDIF + + DO H=1,NCOMB + IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1) + $ .AND.(CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED).OR.(GOODHEL(H) + $ .GT.-HELOFFSET.AND.GOODHEL(H).NE.0)))) THEN + +C Handle the possible requirement of specific polarizations + IF ((.NOT.CHECKPHASE) + $ .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0) + $ .EQ.0.AND.(.NOT.ML5_0_IS_HEL_SELECTED(H))) THEN + CYCLE + ENDIF + + DO I=1,NEXTERNAL + NHEL(I)=HELC(I,H) + ENDDO + + UVCT_REQ_SO_DONE=.FALSE. + CT_REQ_SO_DONE=.FALSE. + LOOP_REQ_SO_DONE=.FALSE. + + IF (.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.HELPICKED.EQ.-1) + $ THEN + HEL_MULT=GOODHEL(H) + ELSE + HEL_MULT=1 + ENDIF + + CTCALL_REQ_SO_DONE=.FALSE. + +C The coefficient were already computed previously with +C another CTMode, so we can skip them + IF (SKIP_LOOPNUM_COEFS_CONSTRUCTION) THEN + GOTO 4000 + ENDIF + + DO I=1,NLOOPGROUPS + DO J=0,LOOPMAXCOEFS-1 + LOOPCOEFS(J,I)=(0.0D0,0.0D0) + ENDDO + ENDDO + + DO K=1,3 + DO I=1,NLOOPAMPS + AMPL(K,I)=(0.0D0,0.0D0) + ENDDO + ENDDO + +C Helas calls for the born amplitudes and counterterms +C associated to given loops + CALL ML5_0_HELAS_CALLS_AMPB_1(P,NHEL,H,IC) + 2000 CONTINUE + CT_REQ_SO_DONE=.TRUE. + +C Helas calls for the counterterm of type 'UVtree' in the UFO. +C These are generated irrespectively of the produced loops. +C In general, only wavefunction renormalization counterterms +C (if needed by the loop UFO model) are of this type. +C Quite often and in principle for all loop UFO models from +C FeynRules, there are none of these type of counterterms. + + 3000 CONTINUE + UVCT_REQ_SO_DONE=.TRUE. + + + CALL ML5_0_COEF_CONSTRUCTION_1(P,NHEL,H,IC) + 4000 CONTINUE + LOOP_REQ_SO_DONE=.TRUE. + + IF(SKIPLOOPEVAL.OR.(.NOT.LOOP_REQ_SO_DONE.AND..NOT.MP_LOOP_RE + $Q_SO_DONE)) THEN + GOTO 5000 + ENDIF + DO I=1,NSQUAREDSO + DO J=1,NLOOPGROUPS + S(I,J)=.TRUE. + ENDDO + ENDDO +C We need the dummy argument I_SO for the squared order index +C to conform to the structure that the call to the LOOP* +C subroutine takes for processes with Born diagrams. + I_SO=1 + CALL ML5_0_LOOP_CT_CALLS_1(P,NHEL,H,IC) + 5000 CONTINUE + CTCALL_REQ_SO_DONE=.TRUE. + + IF (DIRECT_ME_COMPUTATION) THEN + DO I=1,NLOOPAMPS + DO J=1,NLOOPAMPS + CFTOT=DCMPLX(CF_N(I,J)/DBLE(ABS(CF_D(I,J))),0.0D0) + IF(CF_D(I,J).LT.0) CFTOT=CFTOT*IMAG1 + ITEMP = + $ ML5_0_ML5SQSOINDEX(ML5_0_ML5SOINDEX_FOR_LOOP_AMP(I) + $ ,ML5_0_ML5SOINDEX_FOR_LOOP_AMP(J)) + TEMP2(1) = HEL_MULT*DBLE(CFTOT*(AMPL(1,I) + $ *DCONJG(AMPL(1,J)))) +C Computing the quantities below is not strictly +C necessary since the result should be finite +C It is however a good cross-check. + TEMP2(2) = HEL_MULT*DBLE(CFTOT*(AMPL(2,I) + $ *DCONJG(AMPL(1,J)) + AMPL(1,I)*DCONJG(AMPL(2,J)))) + TEMP2(3) = HEL_MULT*DBLE(CFTOT*(AMPL(3,I) + $ *DCONJG(AMPL(1,J)) + AMPL(1,I)*DCONJG(AMPL(3,J)) + $ +AMPL(2,I)*DCONJG(AMPL(2,J)))) +C To mimick the structure of the squared amplitude +C reduction, we add here the squared counterterm +C contribution directly to the result ANS() and put the +C loop contributions in the LOOPRES array which will be +C added to ANS later + IF (I.LE.NCTAMPS) THEN + IF (.NOT.FILTER_SO.OR.SQSO_TARGET.EQ.ITEMP) THEN + DO K=1,3 + ANS(K,ITEMP)=ANS(K,ITEMP)+TEMP2(K) + ANS(K,0)=ANS(K,0)+TEMP2(K) + ENDDO + ENDIF + ELSE + DO K=1,3 + LOOPRES(K,ITEMP,I-NCTAMPS)=LOOPRES(K,ITEMP,I + $ -NCTAMPS)+TEMP2(K) +C During the evaluation of the AMPL, we had stored +C the stability in S(1,*) so we now copy over this +C flag to the relevant contributing Squared orders. + S(ITEMP,I-NCTAMPS)=S(1,I-NCTAMPS) + ENDDO + ENDIF + ENDDO + ENDDO + ENDIF + +C We should compute the color flow either if it contributes to +C the final result (i.e. not used just for the filtering), or +C if the computation is only done from the color flows + IF (((.NOT.DIRECT_ME_COMPUTATION) + $ .AND.ME_COMPUTATION_FROM_JAMP) + $ .OR.((H.EQ.USERHEL.OR.USERHEL.EQ.-1).AND.(POLARIZATIONS(0,0) + $ .EQ.-1.OR.ML5_0_IS_HEL_SELECTED(H)))) THEN +C The cumulative quantities must only be computed if that +C helicity contributes according to user request (second +C argument of the subroutine below). + CALL ML5_0_COMPUTE_COLOR_FLOWS(HEL_MULT) + + + IF(ME_COMPUTATION_FROM_JAMP) THEN + CALL ML5_0_COMPUTE_RES_FROM_JAMP(BUFFRES,HEL_MULT) + IF(((.NOT.DIRECT_ME_COMPUTATION) + $ .AND.ME_COMPUTATION_FROM_JAMP)) THEN +C If the computation from the color flow is the only +C form of computation, we directly update the answer. + DO K=0,3 + DO I=0,NSQUAREDSO + ANS(K,I)=ANS(K,I)+BUFFRES(K,I) + ENDDO + ENDDO +C When setting up the loop filter, it is important to +C set the quantitied LOOPRES. +C Notice that you may have a more powerful filter with +C the direct computation mode because it can filter +C vanishing loop contributions for a particular squared +C split order only +C The quantity LOOPRES defined below is not physical,' +C //' but it's ok since it is only intended for the loop +C filtering. + IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN + DO J=1,NLOOPGROUPS + DO I=1,NSQUAREDSO + DO K=1,3 + LOOPRES(K,I,J)=LOOPRES(K,I,J)+AMPL(K,NCTAMPS+J) + ENDDO + ENDDO + ENDDO + ENDIF +C The if statement below is not strictly necessary but +C makes it clear when it is executed. + ELSEIF(H.EQ.USERHEL.OR.USERHEL.EQ.-1) THEN +C Make sure that that no polarization constraint filters +C out this helicity + IF (POLARIZATIONS(0,0).EQ. + $ -1.OR.ML5_0_IS_HEL_SELECTED(H)) THEN +C If both computational method is used, then we must +C just update RES_FROM_JAMP + DO K=0,3 + DO I=0,NSQUAREDSO + RES_FROM_JAMP(K,I)=RES_FROM_JAMP(K,I)+BUFFRES(K + $ ,I) + ENDDO + ENDDO + ENDIF + ENDIF + IF (H.EQ.USERHEL.OR.USERHEL.EQ.-1) THEN +C Make sure that that no polarization constraint filters +C out this helicity + IF (POLARIZATIONS(0,0).EQ. + $ -1.OR.ML5_0_IS_HEL_SELECTED(H)) THEN + CALL + $ ML5_0_COMPUTE_COLOR_FLOWS_DERIVED_QUANTITIES(HEL_MU + $LT) + ENDIF + ENDIF + ENDIF + ENDIF + + ENDIF + ENDDO + + + IF(DIRECT_ME_COMPUTATION) THEN +C Lines below are not necessary when computing the ME from color +C flows + DO I=0,NSQUAREDSO + DO J=1,3 + BUFFR_BIS(J,I)=ANS(J,I) + ENDDO + ENDDO + ENDIF + + +C MadLoop jumps to this label just after having called the +C subroutine ML5_0_MP_COMPUTE_LOOP_COEFS to compute OpenLoop +C coefficients in quadruple precision (and not double precision +C as done above) + 301 CONTINUE + + IF(DIRECT_ME_COMPUTATION) THEN +C Lines below are not necessary when computing the ME from color +C flows + DO I=0,NSQUAREDSO + DO J=1,3 + ANS(J,I)=BUFFR_BIS(J,I) + ENDDO + ENDDO + ENDIF + + IF ((.NOT.DIRECT_ME_COMPUTATION).AND.ME_COMPUTATION_FROM_JAMP) + $ THEN +C We can skip the update of ANS if it was computed from color +C flows + GOTO 1226 + ENDIF + + IF(SKIPLOOPEVAL.OR.(.NOT.LOOP_REQ_SO_DONE.AND..NOT.MP_LOOP_REQ_SO + $_DONE)) THEN + GOTO 1226 + ENDIF + + + DO I=1,NLOOPGROUPS + LTEMP=.TRUE. + DO K=1,NSQUAREDSO + IF (.NOT.FILTER_SO.OR.SQSO_TARGET.EQ.K) THEN + IF (.NOT.S(K,I)) LTEMP=.FALSE. + DO J=1,3 + ANS(J,K)=ANS(J,K)+LOOPRES(J,K,I) + ANS(J,0)=ANS(J,0)+LOOPRES(J,K,I) + ENDDO + ENDIF + ENDDO + IF((CTMODERUN.NE.-1).AND..NOT.CHECKPHASE.AND.(.NOT.LTEMP)) THEN + WRITE(*,*) '##W03 WARNING Contribution ',I,' is unstable.' + ENDIF + ENDDO + +C Make sure that no NaN is present in the result + DO K=1,NSQUAREDSO + DO J=1,3 + IF (.NOT.(ANS(J,K).EQ.ANS(J,K))) THEN + IF (DOING_QP_EVALS) THEN + FPE_IN_QP_REDUCTION = .TRUE. + ELSE + FPE_IN_DP_REDUCTION = .TRUE. + ENDIF + ENDIF + ENDDO + ENDDO + + 1226 CONTINUE + + IF (CHECKPHASE.OR.(.NOT.HELDOUBLECHECKED)) THEN + IF((USERHEL.EQ.-1).OR.(USERHEL.EQ.HELPICKED)) THEN +C Make sure that that no polarization constraint filters out +C this helicity + IF (POLARIZATIONS(0,0).EQ. + $ -1.OR.ML5_0_IS_HEL_SELECTED(HELPICKED)) THEN +C TO KEEP TRACK OF THE FINAL ANSWER TO BE RETURNED DURING +C CHECK PHASE + DO I=0,NSQUAREDSO + DO K=1,3 + BUFFR(K,I)=BUFFR(K,I)+ANS(K,I) + ENDDO + ENDDO + ENDIF + ENDIF +C SAVE RESULT OF EACH INDEPENDENT HELICITY FOR COMPARISON DURING +C THE HELICITY FILTER SETUP + HELSAVED(1,HELPICKED)=0.0D0 + HELSAVED(2,HELPICKED)=0.0D0 + HELSAVED(3,HELPICKED)=0.0D0 + DO I=1,NSQUAREDSO + IF (CHOSEN_SO_CONFIGS(I)) THEN + HELSAVED(1,HELPICKED)=HELSAVED(1,HELPICKED)+ANS(1,I) + HELSAVED(2,HELPICKED)=HELSAVED(2,HELPICKED)+ANS(2,I) + HELSAVED(3,HELPICKED)=HELSAVED(3,HELPICKED)+ANS(3,I) + ENDIF + ENDDO + +C We make sure not to perform any check when NTRY is 0 because +C this means that +C the REF scale has not been set to the result of the evaluation +C for the first PS point. + IF (CHECKPHASE.AND.NTRY.NE.0) THEN +C SET THE HELICITY FILTER + IF(.NOT.FOUNDHELFILTER) THEN + HEL_INCONSISTENT=.FALSE. + IF(ML5_0_ISZERO(DABS(HELSAVED(1,HELPICKED)) + $ +DABS(HELSAVED(2,HELPICKED))+DABS(HELSAVED(3,HELPICKED)) + $ ,REF/DBLE(NCOMB),-1,-1)) THEN + IF(NTRY.EQ.1) THEN + GOODHEL(HELPICKED)=-HELOFFSET + ELSEIF(GOODHEL(HELPICKED).NE.-HELOFFSET) THEN + WRITE(*,*) '##W02A WARNING Inconsistent zero helicity' + $ //' ',HELPICKED + IF(HELINITSTARTOVER) THEN + WRITE(*,*) '##I01 INFO Initialization starting over' + $ //' because of inconsistency in the helicity filter' + $ //' setup.' + NTRY=0 + ELSE + HEL_INCONSISTENT=.TRUE. + ENDIF + ENDIF + ELSEIF(HELICITYFILTERLEVEL.GT.1) THEN + DO H=1,HELPICKED-1 + IF(GOODHEL(H).GT.-HELOFFSET) THEN +C Be looser for helicity check, bring a factor 100 + DUMMY=ML5_0_ISSAME(HELSAVED(1,HELPICKED),HELSAVED(1 + $ ,H),REF,.FALSE.) + IF(DUMMY.NE.0) THEN + IF(NTRY.EQ.1) THEN +C Set the matching helicity to be contributing +C once more + GOODHEL(H)=GOODHEL(H)+DUMMY +C Use an offset to clearly show it is linked to an +C other one and to avoid overlap + GOODHEL(HELPICKED)=-H-HELOFFSET +C Make sure we have paired this hel config to the +C same one last PS point + ELSEIF(GOODHEL(HELPICKED).NE.(-H-HELOFFSET)) THEN + WRITE(*,*) '##W02B WARNING Inconsistent matching' + $ //' helicity ',HELPICKED + IF(HELINITSTARTOVER) THEN + WRITE(*,*) '##I01 INFO Initialization starting' + $ //' over because of inconsistency in the' + $ //' helicity filter setup.' + NTRY=0 + ELSE + HEL_INCONSISTENT=.TRUE. + ENDIF + ENDIF + ENDIF + ENDIF + ENDDO + ENDIF + IF(HEL_INCONSISTENT) THEN +C This helicity has unstable filter so we will always +C compute it by itself. +C We therefore also need to remove it from the +C multiplicative factor of the corresponding helicity. + IF(GOODHEL(HELPICKED).LT.-HELOFFSET) THEN + GOODHEL(-GOODHEL(HELPICKED)-HELOFFSET)=GOODHEL( + $ -GOODHEL(HELPICKED)-HELOFFSET)-1 + ENDIF +C If several helicities were matched to that one, we need +C to chose another one as reference and redirect the +C others to this new one +C Of course if it is one, then we do not need to do +C anything (because with HELINITSTARTOVER=.FALSE. we only +C support exactly identical Hels.) + IF(GOODHEL(HELPICKED).GT. + $ -HELOFFSET.AND.GOODHEL(HELPICKED).NE.1) THEN + NEWHELREF=-1 + DO H=1,NCOMB + IF (GOODHEL(H).EQ.(-HELOFFSET-HELPICKED)) THEN + IF (NEWHELREF.EQ.-1) THEN + NEWHELREF=H + GOODHEL(H)=GOODHEL(HELPICKED)-1 + ELSE + GOODHEL(H)=-NEWHELREF-HELOFFSET + ENDIF + ENDIF + ENDDO + ENDIF +C In all cases, from now on this helicity will be computed +C independantly of the others. +C In particular, it is the only thing to do if the +C helicity was flagged not contributing. + GOODHEL(HELPICKED)=1 + ENDIF + ENDIF + +C SET THE LOOP FILTER + IF(.NOT.FOUNDLOOPFILTER.AND.USELOOPFILTER) THEN + DO I=1,NLOOPGROUPS + DO J=1,NSQUAREDSO + IF(.NOT.ML5_0_ISZERO(ABS(LOOPRES(1,J,I))+ABS(LOOPRES(2 + $ ,J,I))+ABS(LOOPRES(3,J,I)),(REF*1.0D-4),I,J)) THEN + IF(NTRY.EQ.1) THEN + GOODAMP(J,I)=.TRUE. + LOOPFILTERBUFF(J,I)=.TRUE. + ELSEIF(.NOT.LOOPFILTERBUFF(J,I)) THEN + WRITE(*,*) '##W02 WARNING Inconsistent loop amp ' + $ ,I,'.' + IF(LOOPINITSTARTOVER) THEN + WRITE(*,*) '##I01 INFO Initialization starting' + $ //' over because of inconsistency in the loop' + $ //' filter setup.' + NTRY=0 + ELSE + GOODAMP(J,I)=.TRUE. + ENDIF + ENDIF + ENDIF + ENDDO + ENDDO + ENDIF + ELSEIF (.NOT.HELDOUBLECHECKED.AND.NTRY.NE.0)THEN +C DOUBLE CHECK THE HELICITY FILTER + IF (GOODHEL(HELPICKED).EQ.-HELOFFSET) THEN + IF (.NOT.ML5_0_ISZERO(DABS(HELSAVED(1,HELPICKED)) + $ +DABS(HELSAVED(2,HELPICKED))+DABS(HELSAVED(2,HELPICKED)) + $ ,REF/DBLE(NCOMB),-1,-1)) THEN + WRITE(*,*) '##W15 Helicity filter could not be' + $ //' successfully double checked.' + WRITE(*,*) '##One reason for this is that you might have' + $ //' changed sensible parameters which affected what are' + $ //' the zero helicity configurations.' + WRITE(*,*) '##MadLoop will try to reset the Helicity' + $ //' filter with the next PS points it receives.' + NTRY=0 + OPEN(29,FILE=HELFILTERFN,ERR=348) + 348 CONTINUE + CLOSE(29,STATUS='delete') + ENDIF + ENDIF + IF (GOODHEL(HELPICKED).LT.-HELOFFSET.AND.NTRY.NE.0) THEN + IF(ML5_0_ISSAME(HELSAVED(1,HELPICKED),HELSAVED(1 + $ ,ABS(GOODHEL(HELPICKED)+HELOFFSET)),REF,.TRUE.).EQ.0) THEN + WRITE(*,*) '##W15 Helicity filter could not be' + $ //' successfully double checked.' + WRITE(*,*) '##One reason for this is that you might have' + $ //' changed sensible parameters which affected the' + $ //' helicity dependance relations.' + WRITE(*,*) '##MadLoop will try to reset the Helicity' + $ //' filter with the next PS points it receives.' + NTRY=0 + OPEN(30,FILE=HELFILTERFN,ERR=349) + 349 CONTINUE + CLOSE(30,STATUS='delete') + ENDIF + ENDIF +C SET HELDOUBLECHECKED TO .TRUE. WHEN DONE +C even if it failed we do not want to redo the check +C afterwards if HELINITSTARTOVER=.FALSE. + IF (HELPICKED.EQ.NCOMB.AND.(NTRY.NE.0.OR..NOT.HELINITSTARTOVE + $R)) THEN + DONEHELDOUBLECHECK=.TRUE. + ENDIF + ENDIF + +C GOTO NEXT HELICITY OR FINISH + IF(HELPICKED.NE.NCOMB) THEN + HELPICKED=HELPICKED+1 + MP_DONE=.FALSE. + GOTO 205 + ELSE +C Useful printout +C do I=1,NCOMB +C write(*,*) 'HELSAVED(1,',I,')=',HELSAVED(1,I) +C write(*,*) 'HELSAVED(2,',I,')=',HELSAVED(2,I) +C write(*,*) 'HELSAVED(3,',I,')=',HELSAVED(3,I) +C write(*,*) ' GOODHEL(',I,')=',GOODHEL(I) +C ENDDO + DO I=0,NSQUAREDSO + DO K=1,3 + ANS(K,I)=BUFFR(K,I) + ENDDO + ENDDO +C Update of REF_EVALS (only for loop-induced processes). + TMPR = ABS(ANS(1,0)) + ABS(ANS(2,0)) + ABS(ANS(3,0)) +C We add one here to the number of PS points used for building +C the reference scale for comparisons. +C It might be that when asking for specific helicities, the +C user started with a vanishing helicity +C not filtered yet. In this case, the new ref would remain +C zero. So we want to check for this +C and wait for a point for which the evaluation isn't be zero. + IF(TMPR.NE.0.0D0) THEN + REF_EVALS(MOD(NPSPOINTS,MAXNREF_EVALS)+1) = TMPR + NPSPOINTS = NPSPOINTS+1 + ENDIF + IF(NTRY.EQ.0) THEN + NATTEMPTS=NATTEMPTS+1 + IF(NATTEMPTS.EQ.MAXATTEMPTS) THEN + WRITE(*,*) '##E01 ERROR Could not initialize the filters' + $ //' in ',MAXATTEMPTS,' trials' + STOP 1 + ENDIF + ENDIF + ENDIF + ELSE +C When not in checking mode, update the ref for the first +C MAXNREF_EVALS points (The ref. scale could still be used +C after this stage if BYPASS_CHECK was set to true at some +C point using SLOOPMATRIX_THRES). +C Is is possible to simply remove the if statement below to have +C a running reference scale which always depends on the +C MAXNREF_EVALS *last* evaluations. + IF (NPSPOINTS.LE.MAXNREF_EVALS) THEN + TMPR = ABS(ANS(1,0)) + ABS(ANS(2,0)) + ABS(ANS(3,0)) + IF (TMPR.NE.0.0D0) THEN + REF_EVALS(MOD(NPSPOINTS,MAXNREF_EVALS)+1) = TMPR + NPSPOINTS = NPSPOINTS+1 + ENDIF + ENDIF + ENDIF + +C When computing the ME from the color flows, we also compute the +C born ME from them, so we must apply the normalization factors +C to the born ME as well. + IF(((.NOT.DIRECT_ME_COMPUTATION).AND.ME_COMPUTATION_FROM_JAMP)) + $ THEN + ITEMP=0 + ELSE + ITEMP=1 + ENDIF + DO K=ITEMP,3 + DO I=0,NSQUAREDSO + ANS(K,I)=ANS(K,I)/DBLE(IDEN) + IF (USERHEL.NE.-1) THEN + ANS(K,I)=ANS(K,I)*HELAVGFACTOR + ELSE + DO J=1,NINITIAL + IF (POLARIZATIONS(J,0).NE.-1) THEN + ANS(K,I)=ANS(K,I)*BEAMS_HELAVGFACTOR(J) + ANS(K,I)=ANS(K,I)/POLARIZATIONS(J,0) + ENDIF + ENDDO + ENDIF + ENDDO + ENDDO + + IF (DIRECT_ME_COMPUTATION.AND.ME_COMPUTATION_FROM_JAMP) THEN + WRITE(*,*) ' ================================= ' + WRITE(*,*) ' === JAMP double-checking test === ' + WRITE(*,*) ' ================================= ' + CALL ML5_0_WRITE_MOM(P) + DO J=1,NSQUAREDSO+1 +C We should finish by the summed orders + I = MOD(J,NSQUAREDSO+1) + IF (I.EQ.0) THEN + WRITE(*,*) ' > Checking the sum of all chosen squared' + $ //' split orders' + ELSE + WRITE(*,*) ' > Checking squared split order #',I + ENDIF + DO K=1,1 + RES_FROM_JAMP(K,I)=RES_FROM_JAMP(K,I)/DBLE(IDEN) + IF (USERHEL.NE.-1) THEN + RES_FROM_JAMP(K,I)=RES_FROM_JAMP(K,I)*HELAVGFACTOR + ELSE + DO L=1,NINITIAL + IF (POLARIZATIONS(L,0).NE.-1) THEN + RES_FROM_JAMP(K,I)=RES_FROM_JAMP(K,I) + $ *BEAMS_HELAVGFACTOR(L) + RES_FROM_JAMP(K,I)=RES_FROM_JAMP(K,I) + $ /POLARIZATIONS(L,0) + ENDIF + ENDDO + ENDIF + IF (K.EQ.0) WRITE(*,*) ' || Born :' + IF (K.EQ.1) WRITE(*,*) ' || Finite part :' + IF (K.EQ.2) WRITE(*,*) ' || Single pole residue :' + IF (K.EQ.3) WRITE(*,*) ' || Double pole residue :' + WRITE(*,*) ' --> Direct result =',ANS(K,I) + WRITE(*,*) ' --> Computed from JAMPS =',RES_FROM_JAMP(K,I) + IF((RES_FROM_JAMP(K,I)+ANS(K,I)).EQ.0.0D0) THEN + TMPR = ABS(RES_FROM_JAMP(K,I)-ANS(K,I)) + ELSE + TMPR = ABS((ANS(K,I)-RES_FROM_JAMP(K,I))/((ANS(K,I) + $ +RES_FROM_JAMP(K,I))/2.0D0)) + ENDIF + WRITE(*,*) ' --> Relative diff. =',TMPR + IF(TMPR.GT.JAMP_DOUBLECHECK_THRES) THEN + STOP 'Consistency cross-check of JAMPS failed.' + ENDIF + ENDDO + ENDDO + WRITE(*,*) ' ================================= ' + ENDIF + + IF(.NOT.CHECKPHASE.AND.HELDOUBLECHECKED.AND.(CTMODERUN.EQ.-1)) + $ THEN + STAB_INDEX=STAB_INDEX+1 + IF(DOING_QP_EVALS.AND.LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I_LIB) + $ )) THEN +C Only run over the reduction algorithms which support +C quadruple precision + DO I=0,NSQUAREDSO + DO K=1,3 + QP_RES(K,I,STAB_INDEX)=ANS(K,I) + ENDDO + ENDDO + ELSE + DO I=0,NSQUAREDSO + DO K=1,3 + DP_RES(K,I,STAB_INDEX)=ANS(K,I) + ENDDO + ENDDO + ENDIF + + IF(DOING_QP_EVALS.AND.LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I_LIB) + $ )) THEN + BASIC_CT_MODE=4 + ELSE + BASIC_CT_MODE=1 + ENDIF + +C BEGINNING OF THE DEFINITIONS OF THE DIFFERENT EVALUATION +C METHODS + + IF(.NOT.EVAL_DONE(2)) THEN + EVAL_DONE(2)=.TRUE. + IF(LOOPLIBS_DIRECTEST(MLREDUCTIONLIB(I_LIB)))THEN + CTMODE=BASIC_CT_MODE+1 + GOTO 300 + ELSE +C If some TIR library would not support the loop direction +C test (they all do for now), then we would just copy the +C answer from mode 1 and carry on. + STAB_INDEX=STAB_INDEX+1 + IF(DOING_QP_EVALS)THEN + DO I=0,NSQUAREDSO + DO K=1,3 + QP_RES(K,I,STAB_INDEX)=ANS(K,I) + ENDDO + ENDDO + ELSE + DO I=0,NSQUAREDSO + DO K=1,3 + DP_RES(K,I,STAB_INDEX)=ANS(K,I) + ENDDO + ENDDO + ENDIF + ENDIF + ENDIF + + CTMODE=BASIC_CT_MODE + + IF(.NOT.EVAL_DONE(3).AND. + $ ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.1) + $ .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.1)) ) THEN + EVAL_DONE(3)=.TRUE. + CALL ML5_0_ROTATE_PS(PS,P,1) + IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,1) + GOTO 200 + ENDIF + + IF(.NOT.EVAL_DONE(4).AND. + $ ((DOING_QP_EVALS.AND.NROTATIONS_QP.GE.2) + $ .OR.((.NOT.DOING_QP_EVALS).AND.NROTATIONS_DP.GE.2)) ) THEN + EVAL_DONE(4)=.TRUE. + CALL ML5_0_ROTATE_PS(PS,P,2) + IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,2) + GOTO 200 + ENDIF + + CALL ML5_0_ROTATE_PS(PS,P,0) + IF (DOING_QP_EVALS) CALL ML5_0_MP_ROTATE_PS(MP_PS,MP_P,0) + +C END OF THE DEFINITIONS OF THE DIFFERENT EVALUATION METHODS + + IF(DOING_QP_EVALS.AND.LOOPLIBS_QPAVAILABLE(MLREDUCTIONLIB(I_LIB) + $ )) THEN + CALL ML5_0_COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS) +C If a floating point exception was encountered during the +C reduction, +C the result cannot be trusted at all and we hardset all +C accuracies to 1.0 + IF(FPE_IN_QP_REDUCTION) THEN + DO I=0,NSQUAREDSO + ACC(I)=1.0D0 + ENDDO + ENDIF + DO I=0,NSQUAREDSO + ACCURACY(I)=ACC(I) + ENDDO + RET_CODE_H=3 + RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.TRUE. + $ ,.TRUE.) + IF(MAXVAL(ACC).GE.MLSTABTHRES) THEN + I_QP_LIB=I_QP_LIB+1 + IF(I_QP_LIB.GT.QP_NLOOPLIB.OR.INDEX_QP_TOOLS(I_QP_LIB) + $ .EQ.0)THEN + RET_CODE_H=4 + RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.TRUE. + $ ,.FALSE.) + NEPS=NEPS+1 + CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,TEMP1,TEMP) + CALL ML5_0_COMPUTE_ACCURACY(QP_RES,N_QP_EVAL,ACC,ANS) + IF(NEPS.LE.10) THEN + WRITE(*,*) '##W03 WARNING An unstable PS point was', + $ ' detected.' + IF(FPE_IN_QP_REDUCTION) THEN + WRITE(*,*) '## The last QP reduction was deemed' + $ //' unstable because a floating point exception was' + $ //' encountered.' + ENDIF + IF (NSQUAREDSO.NE.1) THEN + WRITE(*,*) '##Accuracies for each split order,' + $ //' starting with the summed case' + WRITE(*,*) '##DP accuracies (for each split order):' + $ //' ',(TEMP1(I),I=0,NSQUAREDSO) + WRITE(*,*) '##QP accuracies (for each split order):' + $ //' ',(ACC(I),I=0,NSQUAREDSO) + ELSE + WRITE(*,*) '##DP accuracy: ',TEMP1(1) + WRITE(*,*) '##QP accuracy: ',ACC(1) + ENDIF + DO J=0,NSQUAREDSO + IF (NSQUAREDSO.NE.1.OR.J.NE.0) THEN + IF (J.EQ.0) THEN + WRITE(*,*) 'Details for all split orders summed' + $ //' :' + ELSE + WRITE(*,*) 'Details for split order index : ',J + ENDIF + WRITE(*,*) 'Best estimate (fin,1eps,2eps):',(ANS(I + $ ,J),I=1,3) + WRITE(*,*) 'Finite double precision evaluations :' + $ ,(DP_RES(1,J,I),I=1,N_DP_EVAL) + WRITE(*,*) 'Finite quad precision evaluations :' + $ ,(QP_RES(1,J,I),I=1,N_QP_EVAL) + ENDIF + ENDDO + WRITE(*,*) 'PS point specification :' + WRITE(*,*) 'Renormalization scale MU_R=',MU_R + DO I=1,NEXTERNAL + WRITE (*,'(i2,1x,4e27.17)') I, P(0,I),P(1,I),P(2,I) + $ ,P(3,I) + ENDDO + ENDIF + IF(NEPS.EQ.10) THEN + WRITE(*,*) 'Further output of the details of these' + $ //' unstable PS points will now be suppressed.' + ENDIF + ELSE +C A new reduction tool will be used. Reinitialize the FPE +C flags. + FPE_IN_DP_REDUCTION=.FALSE. + FPE_IN_QP_REDUCTION=.FALSE. + I_LIB=INDEX_QP_TOOLS(I_QP_LIB) + EVAL_DONE(1)=.TRUE. + DO I=2,MAXSTABILITYLENGTH + EVAL_DONE(I)=.FALSE. + ENDDO + STAB_INDEX=0 + IF(NROTATIONS_QP.GE.1)THEN + GOTO 200 + ELSE + GOTO 300 + ENDIF + ENDIF + ENDIF + ELSEIF(.NOT.DOING_QP_EVALS)THEN + CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS) +C If a floating point exception was encountered during the +C reduction, +C the result cannot be trusted at all and we hardset all +C accuracies to 1.0 + IF(FPE_IN_DP_REDUCTION) THEN + DO I=0,NSQUAREDSO + ACC(I)=1.0D0 + ENDDO + ENDIF + IF(MAXVAL(ACC).GE.MLSTABTHRES) THEN + I_LIB=I_LIB+1 + IF((I_LIB.GT.NLOOPLIB.OR.MLREDUCTIONLIB(I_LIB).EQ.0) + $ .AND.QP_TOOLS_AVAILABLE)THEN + I_LIB=INDEX_QP_TOOLS(1) +C A new reduction tool will be used. Reinitialize the FPE +C flags. + FPE_IN_DP_REDUCTION=.FALSE. + FPE_IN_QP_REDUCTION=.FALSE. + I_QP_LIB=1 + DOING_QP_EVALS=.TRUE. + EVAL_DONE(1)=.TRUE. + DO I=2,MAXSTABILITYLENGTH + EVAL_DONE(I)=.FALSE. + ENDDO + STAB_INDEX=0 + CTMODE=4 + GOTO 200 + ELSEIF(I_LIB.LE.NLOOPLIB.AND.MLREDUCTIONLIB(I_LIB).GT.0) + $ THEN +C A new reduction tool will be used. Reinitialize the FPE +C flags. + FPE_IN_DP_REDUCTION=.FALSE. + FPE_IN_QP_REDUCTION=.FALSE. + EVAL_DONE(1)=.TRUE. + DO I=2,MAXSTABILITYLENGTH + EVAL_DONE(I)=.FALSE. + ENDDO + STAB_INDEX=0 + IF(NROTATIONS_DP.GE.1)THEN + GOTO 200 + ELSE + GOTO 300 + ENDIF + ELSE + DO I=0,NSQUAREDSO + ACCURACY(I)=ACC(I) + ENDDO + RET_CODE_H=4 + RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE. + $ ,.FALSE.) + NEPS=NEPS+1 + IF(NEPS.LE.10) THEN + WRITE(*,*) '##W03 WARNING An unstable PS point was', + $ ' detected.' + WRITE(*,*) '##W03 WARNING No quadruple precision will' + $ //' be used.' + IF(FPE_IN_DP_REDUCTION) THEN + WRITE(*,*) '## The last DP reduction was deemed' + $ //' unstable because a floating point exception was' + $ //' encountered.' + ENDIF + CALL ML5_0_COMPUTE_ACCURACY(DP_RES,N_DP_EVAL,ACC,ANS) + IF (NSQUAREDSO.NE.1) THEN + WRITE(*,*) 'Accuracies for each split order,' + $ //' starting with the summed case' + WRITE(*,*) 'DP accuracies (for each split order): ' + $ ,(ACC(I),I=0,NSQUAREDSO) + ELSE + WRITE(*,*) 'DP accuracy: ',ACC(1) + ENDIF + DO J=0,NSQUAREDSO + IF (NSQUAREDSO.NE.1.OR.J.NE.0) THEN + IF (J.EQ.0) THEN + WRITE(*,*) 'Details for all split orders summed' + $ //' :' + ELSE + WRITE(*,*) 'Details for split order index : ',J + ENDIF + WRITE(*,*) 'Best estimate (fin,1eps,2eps):',(ANS(I + $ ,J),I=1,3) + WRITE(*,*) 'Finite double precision evaluations :' + $ ,(DP_RES(1,J,I),I=1,N_DP_EVAL) + ENDIF + ENDDO + WRITE(*,*) 'PS point specification :' + WRITE(*,*) 'Renormalization scale MU_R=',MU_R + DO I=1,NEXTERNAL + WRITE (*,'(i2,1x,4e27.17)') I, P(0,I),P(1,I),P(2,I) + $ ,P(3,I) + ENDDO + ENDIF + IF(NEPS.EQ.10) THEN + WRITE(*,*) 'Further output of the details of these' + $ //' unstable PS points will now be suppressed.' + ENDIF + ENDIF + ELSE + DO I=0,NSQUAREDSO + ACCURACY(I)=ACC(I) + ENDDO + RET_CODE_H=2 + RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE. + $ ,.TRUE.) + ENDIF + ENDIF + ELSE + RET_CODE_H=1 + DO I=0,NSQUAREDSO + ACCURACY(I)=-1.0D0 + ENDDO + RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE. + $ ,.FALSE.) + ENDIF + + 9999 CONTINUE + +C Finalize the return code + IF (MP_DONE_ONCE) THEN + RET_CODE_T=2 + ELSE + RET_CODE_T=1 + ENDIF + IF(CHECKPHASE.OR..NOT.HELDOUBLECHECKED) THEN + RET_CODE_H=1 + RET_CODE_U=SET_RET_CODE_U(MLREDUCTIONLIB(I_LIB),.FALSE. + $ ,.FALSE.) + RET_CODE_T=RET_CODE_T+2 + DO I=0,NSQUAREDSO + ACCURACY(I)=-1.0D0 + ENDDO + ENDIF + +C Finally for the summed result in ANS(1:3,0), make sure to only +C consider the squared order asked for by the user. +C Notice that this filtering using CHOSEN_SO_CONFIGS happens +C here only while everywhere else one always considers the sum. + DO J=1,3 + ANS(J,0)=0.0D0 + ENDDO + DO I=1,NSQUAREDSO + IF (CHOSEN_SO_CONFIGS(I)) THEN + DO J=1,3 + ANS(J,0)=ANS(J,0)+ANS(J,I) + ENDDO + ENDIF + ENDDO + +C Reinitialize the default threshold if it was specified by the +C user + IF (USER_STAB_PREC.GT.0.0D0) THEN + MLSTABTHRES=MLSTABTHRES_BU + CTMODEINIT=CTMODEINIT_BU + ENDIF + +C Reinitialize the Lorentz test if it had been disabled because +C spin-2 particles are in the external states. + NROTATIONS_DP = NROTATIONS_DP_BU + NROTATIONS_QP = NROTATIONS_QP_BU + +C Reinitialize the check phase logicals and the filters if check +C bypassed + IF (BYPASS_CHECK) THEN + CHECKPHASE = OLD_CHECKPHASE + HELDOUBLECHECKED = OLD_HELDOUBLECHECKED + DO I=1,NCOMB + GOODHEL(I)=OLD_GOODHEL(I) + ENDDO + DO I=1,NSQUAREDSO + DO J=1,NLOOPGROUPS + GOODAMP(I,J)=OLD_GOODAMP(I,J) + ENDDO + ENDDO + ENDIF + +C Make sure that we finish by emptying caches + IF (AUTOMATIC_CACHE_CLEARING) THEN + CALL ML5_0_CLEAR_CACHES() + ENDIF + +C Now make sure to turn off the global COLLIER cache if applicable + CALL ML5_0_SET_COLLIER_GLOBAL_CACHE(.FALSE.) + + END + + SUBROUTINE ML5_0_CLEAR_CACHES() +C +C This routine can be called directly from the user if +C AUTOMATIC_CACHE_CLEARING is set to False. It must then be called +C after +C ech event +C + CALL ML5_0_CLEAR_TIR_CACHE() + CALL NINJA_CLEAR_INTEGRAL_CACHE() + CALL ML5_0_CLEAR_COLLIER_CACHE() + END + +C --=========================================-- +C General Helper functions and subroutine +C for the main sloopmatrix subroutine +C --=========================================-- + + LOGICAL FUNCTION ML5_0_IS_HEL_SELECTED(HELID) + IMPLICIT NONE +C +C CONSTANTS +C + INTEGER NEXTERNAL + PARAMETER (NEXTERNAL=3) + INTEGER NCOMB + PARAMETER (NCOMB=4) +C +C ARGUMENTS +C + INTEGER HELID +C +C LOCALS +C + INTEGER I,J + LOGICAL FOUNDIT +C +C GLOBALS +C + INTEGER HELC(NEXTERNAL,NCOMB) + COMMON/ML5_0_HELCONFIGS/HELC + + INTEGER POLARIZATIONS(0:NEXTERNAL,0:5) + COMMON/ML5_0_BEAM_POL/POLARIZATIONS +C ---------- +C BEGIN CODE +C ---------- + + ML5_0_IS_HEL_SELECTED = .TRUE. + IF (POLARIZATIONS(0,0).EQ.-1) THEN + RETURN + ENDIF + + DO I=1,NEXTERNAL + IF (POLARIZATIONS(I,0).EQ.-1) THEN + CYCLE + ENDIF + FOUNDIT = .FALSE. + DO J=1,POLARIZATIONS(I,0) + IF (HELC(I,HELID).EQ.POLARIZATIONS(I,J)) THEN + FOUNDIT = .TRUE. + EXIT + ENDIF + ENDDO + IF(.NOT.FOUNDIT) THEN + ML5_0_IS_HEL_SELECTED = .FALSE. + RETURN + ENDIF + ENDDO + RETURN + + END + + LOGICAL FUNCTION ML5_0_ISZERO(TOTEST, REFERENCE_VALUE, LOOP, + $ SOINDEX) + IMPLICIT NONE +C +C CONSTANTS +C + INTEGER NLOOPGROUPS + PARAMETER (NLOOPGROUPS=4) + INTEGER NSQUAREDSO + PARAMETER (NSQUAREDSO=1) +C +C ARGUMENTS +C + REAL*8 TOTEST, REFERENCE_VALUE + INTEGER LOOP, SOINDEX +C +C GLOBAL +C + INCLUDE 'MadLoopParams.inc' + COMPLEX*16 LOOPRES(3,NSQUAREDSO,NLOOPGROUPS) + LOGICAL S(NSQUAREDSO,NLOOPGROUPS) + COMMON/ML5_0_LOOPRES/LOOPRES,S +C ---------- +C BEGIN CODE +C ---------- + IF(ABS(REFERENCE_VALUE).EQ.0.0D0) THEN + ML5_0_ISZERO=.FALSE. + WRITE(*,*) '##E02 ERRROR Reference value for comparison is' + $ //' zero.' + STOP 1 + ELSE + ML5_0_ISZERO=((ABS(TOTEST)/ABS(REFERENCE_VALUE)).LT.ZEROTHRES) + ENDIF + + IF(LOOP.NE.-1) THEN + IF((.NOT.ML5_0_ISZERO).AND.(.NOT.S(SOINDEX,LOOP))) THEN + WRITE(*,*) '##W01 WARNING Contribution ',LOOP,' of split' + $ //' order ',SOINDEX,' is detected as contributing with CR=' + $ ,(ABS(TOTEST)/ABS(REFERENCE_VALUE)),' but is unstable.' + ENDIF + ENDIF + + END + + INTEGER FUNCTION ML5_0_ISSAME(RESA,RESB,REF,USEMAX) + IMPLICIT NONE +C This function compares the result from two different helicity +C configuration A and B +C It returns 0 if they are not related and (+/-wgt) if +C A=(+/-wgt)*B. +C For now, the only wgt implemented is the integer 1 or -1. +C If useMax is .TRUE., it uses all implemented weights no matter +C what is HELINITSTARTOVER +C +C CONSTANTS +C + INTEGER MAX_WGT_TO_TRY + PARAMETER (MAX_WGT_TO_TRY=2) +C +C ARGUMENTS +C + REAL*8 RESA(3), RESB(3) + REAL*8 REF + LOGICAL USEMAX +C +C LOCAL VARIABLES +C + LOGICAL ML5_0_ISZERO + INTEGER I,J + INTEGER N_WGT_TO_TRY + INTEGER WGT_TO_TRY(MAX_WGT_TO_TRY) + DATA WGT_TO_TRY/1,-1/ +C +C INCLUDES +C + INCLUDE 'MadLoopParams.inc' +C ---------- +C BEGIN CODE +C ---------- + ML5_0_ISSAME=0 + +C If the helicity can be constructed progressively while allowing +C inconsistency, then we only allow for weight one comparisons. + IF (.NOT.HELINITSTARTOVER.AND..NOT.USEMAX) THEN + N_WGT_TO_TRY=1 + ELSE + N_WGT_TO_TRY=MAX_WGT_TO_TRY + ENDIF + + DO I=1,N_WGT_TO_TRY + DO J=1,3 + IF (ML5_0_ISZERO(ABS(RESB(J)),REF,-1,-1)) THEN + IF(.NOT.ML5_0_ISZERO(ABS(RESB(J))+ABS(RESA(J)),REF,-1,-1)) + $ THEN + GOTO 1231 + ENDIF +C Be looser for helicity comparison, so bring a factor 100 + ELSEIF(.NOT.ML5_0_ISZERO(ABS((RESA(J)/RESB(J)) + $ -DBLE(WGT_TO_TRY(I))),1.0D0,-1,-1)) THEN + GOTO 1231 + ENDIF + ENDDO + ML5_0_ISSAME = WGT_TO_TRY(I) + RETURN + 1231 CONTINUE + ENDDO + END + + SUBROUTINE ML5_0_COMPUTE_ACCURACY(FULLLIST, LENGTH, ACC, + $ ESTIMATE) + IMPLICIT NONE +C +C PARAMETERS +C + INTEGER MAXSTABILITYLENGTH + COMMON/ML5_0_STABILITY_TESTS/MAXSTABILITYLENGTH + INTEGER NSQUAREDSO + PARAMETER (NSQUAREDSO=1) +C +C ARGUMENTS +C + REAL*8 FULLLIST(3,0:NSQUAREDSO,MAXSTABILITYLENGTH) + INTEGER LENGTH + REAL*8 ACC(0:NSQUAREDSO), ESTIMATE(0:3,0:NSQUAREDSO) +C +C LOCAL VARIABLES +C + LOGICAL MASK(MAXSTABILITYLENGTH) + LOGICAL MASK3(3) + DATA MASK3/.TRUE.,.TRUE.,.TRUE./ + INTEGER I,J,K + REAL*8 AVG + REAL*8 DIFF + REAL*8 ACCURACIES(3) + REAL*8 LIST(MAXSTABILITYLENGTH) + +C +C GLOBAL VARIABLES +C + LOGICAL CHOSEN_SO_CONFIGS(NSQUAREDSO) + COMMON/ML5_0_CHOSEN_LOOP_SQSO/CHOSEN_SO_CONFIGS + INTEGER I_LIB + COMMON/ML5_0_I_LIB/I_LIB + INCLUDE 'MadLoopParams.inc' + +C ---------- +C BEGIN CODE +C ---------- + DO I=1,LENGTH + MASK(I)=.TRUE. + ENDDO + DO I=LENGTH+1,MAXSTABILITYLENGTH + MASK(I)=.FALSE. +C For some architectures, it is necessary to initialize all the +C elements of fulllist(i,j) +C Beware that if the length provided is incorrect, then this can +C corrup the fulllist given in argument. + DO J=0,NSQUAREDSO + DO K=1,3 + FULLLIST(K,J,I)=0.0D0 + ENDDO + ENDDO + ENDDO + + DO K=0,NSQUAREDSO + + DO I=1,3 + DO J=1,MAXSTABILITYLENGTH + LIST(J)=FULLLIST(I,K,J) + ENDDO + DIFF=MAXVAL(LIST,1,MASK)-MINVAL(LIST,1,MASK) + AVG=(MAXVAL(LIST,1,MASK)+MINVAL(LIST,1,MASK))/2.0D0 + ESTIMATE(I,K)=AVG + IF (AVG.EQ.0.0D0) THEN + ACCURACIES(I)=DIFF + ELSE + ACCURACIES(I)=DIFF/ABS(AVG) + ENDIF + ENDDO + +C The technique below is too sensitive, typically to +C unstablities in very small poles +C acc(k)=MAXVAL(ACCURACIES,1,MASK3) +C The following is used instead + ACC(K) = 0.0D0 + AVG = 0.0D0 + DO I=1,3 + ACC(K) = ACC(K) + ACCURACIES(I)*ABS(ESTIMATE(I,K)) + AVG = AVG + ESTIMATE(I,K) + ENDDO + IF (AVG.NE.0.0D0) THEN + ACC(K) = ACC(K) / ( ABS(AVG) / 3.0D0) + ENDIF + +C When using COLLIER with the internal stability test, the first +C evaluation is typically more reliable so we do not want to +C use the average but rather the first evaluation. + IF (MLREDUCTIONLIB(I_LIB) + $ .EQ.7.AND.COLLIERUSEINTERNALSTABILITYTEST) THEN + DO I=1,3 + ESTIMATE(I,K) = FULLLIST(I,K,1) + ENDDO + ENDIF + +C Make sure to hard-set to zero accuracies of coupling orders +C not included + IF (K.NE.0) THEN + IF (.NOT.CHOSEN_SO_CONFIGS(K)) THEN + ACC(K) = 0.0D0 + ENDIF + ENDIF + +C If NaN are present in the evaluation, automatically set the +C accuracy to 1.0d99. + DO I=1,3 + DO J=1,MAXSTABILITYLENGTH + IF (ISNAN(FULLLIST(I,K,J))) THEN + ACC(K) = 1.0D99 + ENDIF + ENDDO + ENDDO + + ENDDO + + END + + SUBROUTINE ML5_0_SET_N_EVALS(N_DP_EVALS,N_QP_EVALS) + + IMPLICIT NONE + INTEGER N_DP_EVALS, N_QP_EVALS + + INCLUDE 'MadLoopParams.inc' + + IF(CTMODERUN.LE.-1) THEN + N_DP_EVALS=2+NROTATIONS_DP + N_QP_EVALS=2+NROTATIONS_QP + ELSE + N_DP_EVALS=1 + N_QP_EVALS=1 + ENDIF + + IF(N_DP_EVALS.GT.20.OR.N_QP_EVALS.GT.20) THEN + WRITE(*,*) 'ERROR:: Increase hardcoded maxstabilitylength.' + STOP 1 + ENDIF + + END + +C THIS SUBROUTINE SIMPLY SET THE GLOBAL PS CONFIGURATION GLOBAL +C VARIABLES FROM A GIVEN VARIABLE IN DOUBLE PRECISION + SUBROUTINE ML5_0_SET_MP_PS(P) + + INTEGER NEXTERNAL + PARAMETER (NEXTERNAL=3) + REAL*16 MP_PS(0:3,NEXTERNAL),MP_P(0:3,NEXTERNAL) + COMMON/ML5_0_MP_PSPOINT/MP_PS,MP_P + REAL*8 P(0:3,NEXTERNAL) + + DO I=1,NEXTERNAL + DO J=0,3 + MP_PS(J,I)=P(J,I) + ENDDO + ENDDO + CALL ML5_0_MP_IMPROVE_PS_POINT_PRECISION(MP_PS) + DO I=1,NEXTERNAL + DO J=0,3 + MP_P(J,I)=MP_PS(J,I) + ENDDO + ENDDO + + END + +C --=========================================-- +C Functions for dealing with the ordering +C and indexing of split order contributions +C --=========================================-- + + SUBROUTINE ML5_0_GET_NSQSO_LOOP(NSQSO) +C +C Simple subroutine returning the number of squared split order +C contributions returned in ANS when calling sloopmatrix +C + INTEGER NSQUAREDSO + PARAMETER (NSQUAREDSO=1) + + INTEGER NSQSO + + NSQSO=NSQUAREDSO + + END + + SUBROUTINE ML5_0_GET_ANSWER_DIMENSION(ANS_DIM) +C +C MadLoop subroutines return an array of dimension +C ANS(0:3,0:ANS_DIM) +C In order for the user program to be able to correctly declare +C this +C array when calling MadLoop, this subroutine returns its dimension +C + INTEGER NSQUAREDSO + PARAMETER (NSQUAREDSO=1) + INTEGER ANS_DIM + + INTEGER NSQSO_BORN + PARAMETER (NSQSO_BORN=0) + + + ANS_DIM=MAX(NSQSO_BORN,NSQUAREDSO) + + END + + INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS(ORDERS) +C +C This functions returns the integer index identifying the split +C orders list passed in argument which correspond to the values +C of the following list of couplings (and in this order): +C ['QCD'] +C +C CONSTANTS +C + INTEGER NSO, NSQSO + PARAMETER (NSO=1, NSQSO=1) +C +C ARGUMENTS +C + INTEGER ORDERS(NSO) +C +C LOCAL VARIABLES +C + INTEGER I,J + INTEGER SQPLITORDERS(NSQSO,NSO) + DATA (SQPLITORDERS( 1,I),I= 1, 1) / 4/ + COMMON/ML5_0_ML5SQPLITORDERS/SQPLITORDERS +C +C BEGIN CODE +C + DO I=1,NSQSO + DO J=1,NSO + IF (ORDERS(J).NE.SQPLITORDERS(I,J)) GOTO 1009 + ENDDO + ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS = I + RETURN + 1009 CONTINUE + ENDDO + + WRITE(*,*) 'ERROR:: Stopping function' + $ //' ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS' + WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO) + STOP + + END + + INTEGER FUNCTION ML5_0_GETORDPOWFROMINDEX_ML5(IORDER, INDX) +C +C Return the power of the IORDER-th order appearing at position +C INDX +C in the split-orders output +C +C ['QCD'] +C +C CONSTANTS +C + INTEGER NSO, NSQSO + PARAMETER (NSO=1, NSQSO=1) +C +C ARGUMENTS +C + INTEGER ORDERS(NSO) +C +C LOCAL VARIABLES +C + INTEGER I,J + INTEGER SQPLITORDERS(NSQSO,NSO) + DATA (SQPLITORDERS( 1,I),I= 1, 1) / 4/ +C +C BEGIN CODE +C + IF (IORDER.GT.NSO.OR.IORDER.LT.1) THEN + WRITE(*,*) 'INVALID IORDER ML5', IORDER + WRITE(*,*) 'SHOULD BE BETWEEN 1 AND ', NSO + STOP + ENDIF + + IF (INDX.GT.NSQSO.OR.INDX.LT.1) THEN + WRITE(*,*) 'INVALID INDX ML5', INDX + WRITE(*,*) 'SHOULD BE BETWEEN 1 AND ', NSQSO + STOP + ENDIF + + ML5_0_GETORDPOWFROMINDEX_ML5=SQPLITORDERS(INDX, IORDER) + + END + + INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_BORN_AMP(AMPID) +C +C For a given born amplitude number, it returns the ID of the +C split orders it has +C +C CONSTANTS +C + INTEGER NBORNAMPS + PARAMETER (NBORNAMPS=0) +C +C ARGUMENTS +C + INTEGER AMPID +C +C LOCAL VARIABLES +C + INTEGER BORNAMPORDERS(NBORNAMPS) + +C ----------- +C BEGIN CODE +C ----------- + IF (AMPID.GT.NBORNAMPS) THEN + WRITE(*,*) 'ERROR:: Born amplitude ID ',AMPID,' above the' + $ //' maximum ',NBORNAMPS + ENDIF + ML5_0_ML5SOINDEX_FOR_BORN_AMP = BORNAMPORDERS(AMPID) + + END + + INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_LOOP_AMP(AMPID) +C +C For a given loop amplitude number, it returns the ID of the +C split orders it has +C +C CONSTANTS +C + INTEGER NLOOPAMPS + PARAMETER (NLOOPAMPS=6) +C +C ARGUMENTS +C + INTEGER AMPID +C +C LOCAL VARIABLES +C + INTEGER LOOPAMPORDERS(NLOOPAMPS) + DATA (LOOPAMPORDERS(I),I= 1, 5) / 1, 1, 1, 1, 1/ + DATA (LOOPAMPORDERS(I),I= 6, 6) / 1/ +C ----------- +C BEGIN CODE +C ----------- + IF (AMPID.GT.NLOOPAMPS) THEN + WRITE(*,*) 'ERROR:: Loop amplitude ID ',AMPID,' above the' + $ //' maximum ',NLOOPAMPS + ENDIF + ML5_0_ML5SOINDEX_FOR_LOOP_AMP = LOOPAMPORDERS(AMPID) + + END + + + INTEGER FUNCTION ML5_0_ML5SQSOINDEX(ORDERINDEXA, ORDERINDEXB) +C +C This functions plays the role of the interference matrix. It can +C be hardcoded or +C made more elegant using hashtables if its execution speed ever +C becomes a relevant +C factor. From two split order indices, it return the +C corresponding index in the squared +C order canonical ordering. +C +C CONSTANTS +C + INTEGER NSO, NSQUAREDSO, NAMPSO + PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=1) +C +C ARGUMENTS +C + INTEGER ORDERINDEXA, ORDERINDEXB +C +C LOCAL VARIABLES +C + INTEGER I, SQORDERS(NSO) + INTEGER AMPSPLITORDERS(NAMPSO,NSO) + DATA (AMPSPLITORDERS( 1,I),I= 1, 1) / 2/ + COMMON/ML5_0_ML5AMPSPLITORDERS/AMPSPLITORDERS +C +C FUNCTION +C + INTEGER ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS +C +C BEGIN CODE +C + DO I=1,NSO + SQORDERS(I)=AMPSPLITORDERS(ORDERINDEXA,I) + $ +AMPSPLITORDERS(ORDERINDEXB,I) + ENDDO + ML5_0_ML5SQSOINDEX=ML5_0_ML5SOINDEX_FOR_SQUARED_ORDERS(SQORDERS) + END + +C This is the inverse subroutine of ML5SOINDEX_FOR_SQUARED_ORDERS. +C Not directly useful, but provided nonetheless. + SUBROUTINE ML5_0_ML5GET_SQUARED_ORDERS_FOR_SOINDEX(SOINDEX + $ ,ORDERS) +C +C This functions returns the orders identified by the squared +C split order index in argument. Order values correspond to +C following list of couplings (and in this order): +C ['QCD'] +C +C CONSTANTS +C + INTEGER NSO, NSQSO + PARAMETER (NSO=1, NSQSO=1) +C +C ARGUMENTS +C + INTEGER SOINDEX, ORDERS(NSO) +C +C LOCAL VARIABLES +C + INTEGER I + INTEGER SQPLITORDERS(NSQSO,NSO) + COMMON/ML5_0_ML5SQPLITORDERS/SQPLITORDERS +C +C BEGIN CODE +C + IF (SOINDEX.GT.0.AND.SOINDEX.LE.NSQSO) THEN + DO I=1,NSO + ORDERS(I) = SQPLITORDERS(SOINDEX,I) + ENDDO + RETURN + ENDIF + + WRITE(*,*) 'ERROR:: Stopping function' + $ //' ML5_0_ML5GET_SQUARED_ORDERS_FOR_SOINDEX' + WRITE(*,*) 'Could not find squared orders index ',SOINDEX + STOP + + END SUBROUTINE + +C This is the inverse subroutine of getting amplitude SO orders. +C Not directly useful, but provided nonetheless. + SUBROUTINE ML5_0_ML5GET_ORDERS_FOR_AMPSOINDEX(SOINDEX,ORDERS) +C +C This functions returns the orders identified by the split order +C index in argument. Order values correspond to following list of +C couplings (and in this order): +C ['QCD'] +C +C CONSTANTS +C + INTEGER NSO, NAMPSO + PARAMETER (NSO=1, NAMPSO=1) +C +C ARGUMENTS +C + INTEGER SOINDEX, ORDERS(NSO) +C +C LOCAL VARIABLES +C + INTEGER I + INTEGER AMPSPLITORDERS(NAMPSO,NSO) + COMMON/ML5_0_ML5AMPSPLITORDERS/AMPSPLITORDERS +C +C BEGIN CODE +C + IF (SOINDEX.GT.0.AND.SOINDEX.LE.NAMPSO) THEN + DO I=1,NSO + ORDERS(I) = AMPSPLITORDERS(SOINDEX,I) + ENDDO + RETURN + ENDIF + + WRITE(*,*) 'ERROR:: Stopping function' + $ //' ML5_0_ML5GET_ORDERS_FOR_AMPSOINDEX' + WRITE(*,*) 'Could not find amplitude split orders index ',SOINDEX + STOP + + END SUBROUTINE + + +C This function is not directly useful, but included for +C completeness + INTEGER FUNCTION ML5_0_ML5SOINDEX_FOR_AMPORDERS(ORDERS) +C +C This functions returns the integer index identifying the +C amplitude split orders passed in argument which correspond to +C the values of the following list of couplings (and in this +C order): +C ['QCD'] +C +C CONSTANTS +C + INTEGER NSO, NAMPSO + PARAMETER (NSO=1, NAMPSO=1) +C +C ARGUMENTS +C + INTEGER ORDERS(NSO) +C +C LOCAL VARIABLES +C + INTEGER I,J + INTEGER AMPSPLITORDERS(NAMPSO,NSO) + COMMON/ML5_0_ML5AMPSPLITORDERS/AMPSPLITORDERS +C +C BEGIN CODE +C + DO I=1,NAMPSO + DO J=1,NSO + IF (ORDERS(J).NE.AMPSPLITORDERS(I,J)) GOTO 1009 + ENDDO + ML5_0_ML5SOINDEX_FOR_AMPORDERS = I + RETURN + 1009 CONTINUE + ENDDO + + WRITE(*,*) 'ERROR:: Stopping function' + $ //' ML5_0_ML5SOINDEX_FOR_AMPORDERS' + WRITE(*,*) 'Could not find squared orders ',(ORDERS(I),I=1,NSO) + STOP + + END + +C --=========================================-- +C Definition of additional access routines +C --=========================================-- + + SUBROUTINE ML5_0_COLLIER_COMPUTE_UV_POLES(ONOFF) +C +C This function can be called by the MadLoop user so as to chose +C to have COLLIER +C compute the UV pole or not (it costs more time). +C + LOGICAL ONOFF + + INCLUDE 'MadLoopParams.inc' + + LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, + $ FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION + LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, + $ COLLIER_IR_POLE_COMPUTATION_CHOICE + COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE + $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, + $ FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION + $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE + $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE + + COLLIERCOMPUTEUVPOLES = ONOFF +C This is just so that if we read the param again, we don't +C overwrite the choice made here + FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION = .TRUE. + COLLIER_UV_POLE_COMPUTATION_CHOICE = ONOFF + + END SUBROUTINE + + SUBROUTINE ML5_0_COLLIER_COMPUTE_IR_POLES(ONOFF) +C +C This function can be called by the MadLoop user so as to chose +C to have COLLIER +C compute the IR pole or not (it costs more time). +C + LOGICAL ONOFF + + INCLUDE 'MadLoopParams.inc' + + LOGICAL FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, + $ FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION + LOGICAL COLLIER_UV_POLE_COMPUTATION_CHOICE, + $ COLLIER_IR_POLE_COMPUTATION_CHOICE + COMMON/ML5_0_COLLIERPOLESFORCEDCHOICE + $ /FORCED_CHOICE_OF_COLLIER_UV_POLE_COMPUTATION, + $ FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION + $ ,COLLIER_UV_POLE_COMPUTATION_CHOICE + $ ,COLLIER_IR_POLE_COMPUTATION_CHOICE + + COLLIERCOMPUTEIRPOLES = ONOFF +C This is just so that if we read the param again, we don't +C overwrite the choice made here + FORCED_CHOICE_OF_COLLIER_IR_POLE_COMPUTATION = .TRUE. + COLLIER_IR_POLE_COMPUTATION_CHOICE = ONOFF + + END SUBROUTINE + + SUBROUTINE ML5_0_FORCE_STABILITY_CHECK(ONOFF) +C +C This function can be called by the MadLoop user so as to always +C have stability +C checked, even during initialisation, when calling the *_thres +C routines. +C + LOGICAL ONOFF + + LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY + DATA BYPASS_CHECK, ALWAYS_TEST_STABILITY /.FALSE.,.FALSE./ + COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY + + ALWAYS_TEST_STABILITY = ONOFF + + END SUBROUTINE + + SUBROUTINE ML5_0_SET_AUTOMATIC_CACHE_CLEARING(ONOFF) +C +C This function can be called by the MadLoop user so as to +C manually chose when +C to reset the TIR cache. +C + IMPLICIT NONE + + INCLUDE 'MadLoopParams.inc' + + LOGICAL ONOFF + + LOGICAL AUTOMATIC_CACHE_CLEARING + DATA AUTOMATIC_CACHE_CLEARING/.TRUE./ + COMMON/ML5_0_RUNTIME_OPTIONS/AUTOMATIC_CACHE_CLEARING + + INTEGER N_DP_EVAL, N_QP_EVAL + COMMON/ML5_0_N_EVALS/N_DP_EVAL,N_QP_EVAL + + + AUTOMATIC_CACHE_CLEARING = ONOFF + + IF (NROTATIONS_DP.NE.0.OR.NROTATIONS_QP.NE.0) THEN + WRITE(*,*) 'Warning: One cannot remove the TIR cache automatic' + $ //' clearing while at the same time keeping Lorentz rotations' + $ //' for stability tests.' + WRITE(*,*) 'MadLoop will therefore automatically set' + $ //' NRotations_DP and NRotations_QP to 0.' + NROTATIONS_DP = 0 + NROTATIONS_QP = 0 + CALL ML5_0_SET_N_EVALS(N_DP_EVAL,N_QP_EVAL) + ENDIF + END SUBROUTINE + + SUBROUTINE ML5_0_SET_COUPLINGORDERS_TARGET(SOTARGET) + IMPLICIT NONE +C +C This routine can be accessed by an external user to set the +C squared split order target. +C If set to a value different than -1, the code will try to avoid +C computing anything which +C does not contribute to contributions of squared split orders +C SQSO_TARGET and below. +C This can considerably speed up the code. However, keep in mind +C that any contribution of +C 'squared order index' larger than SQSO_TARGET cannot be trust. +C +C ARGUMENTS +C + INTEGER SOTARGET +C +C GLOBAL +C + INTEGER SQSO_TARGET + COMMON/ML5_0_SOCHOICE/SQSO_TARGET +C ---------- +C BEGIN CODE +C ---------- + SQSO_TARGET = SOTARGET + END + + SUBROUTINE ML5_0_SET_LEG_POLARIZATION(LEG_ID, LEG_POLARIZATION) + IMPLICIT NONE +C +C ARGUMENTS +C + INTEGER LEG_ID + INTEGER LEG_POLARIZATION +C +C LOCALS +C + INTEGER I + INTEGER LEG_POLARIZATIONS(0:5) +C ---------- +C BEGIN CODE +C ---------- + + IF (LEG_POLARIZATION.EQ.-10000) THEN + LEG_POLARIZATIONS(0)=-1 + DO I=1,5 + LEG_POLARIZATIONS(I)=-10000 + ENDDO + ELSE + LEG_POLARIZATIONS(0)=1 + LEG_POLARIZATIONS(1)=LEG_POLARIZATION + DO I=2,5 + LEG_POLARIZATIONS(I)=-10000 + ENDDO + ENDIF + CALL ML5_0_SET_LEG_POLARIZATIONS(LEG_ID,LEG_POLARIZATIONS) + + END + + SUBROUTINE ML5_0_SET_LEG_POLARIZATIONS(LEG_ID, LEG_POLARIZATIONS) + IMPLICIT NONE +C +C CONSTANTS +C + INTEGER NEXTERNAL + PARAMETER (NEXTERNAL=3) + INTEGER NPOLENTRIES + PARAMETER (NPOLENTRIES=(NEXTERNAL+1)*6) + INTEGER NCOMB + PARAMETER (NCOMB=4) +C +C ARGUMENTS +C + INTEGER LEG_ID + INTEGER LEG_POLARIZATIONS(0:5) +C +C LOCALS +C + INTEGER I,J + LOGICAL ALL_SUMMED_OVER +C +C GLOBALS +C +C Entry 0 of the first dimension is all -1 if there is no +C polarization requirement. +C Then for each leg with ID legID, it is either summed over if +C POLARIZATIONS(legID,0) is -1, or the list of helicity considered +C for that +C leg is POLARIZATIONS(legID,1: POLARIZATIONS(legID,0) ). + INTEGER POLARIZATIONS(0:NEXTERNAL,0:5) + DATA ((POLARIZATIONS(I,J),I=0,NEXTERNAL),J=0,5)/NPOLENTRIES*-1/ + COMMON/ML5_0_BEAM_POL/POLARIZATIONS + + INTEGER BORN_POLARIZATIONS(0:NEXTERNAL,0:5) + COMMON/ML5_0_BORN_BEAM_POL/BORN_POLARIZATIONS + +C ---------- +C BEGIN CODE +C ---------- + + IF (LEG_POLARIZATIONS(0).EQ.-1) THEN + DO I=0,5 + POLARIZATIONS(LEG_ID,I)=-1 + ENDDO + ELSE + DO I=0,LEG_POLARIZATIONS(0) + POLARIZATIONS(LEG_ID,I)=LEG_POLARIZATIONS(I) + ENDDO + DO I=LEG_POLARIZATIONS(0)+1,5 + POLARIZATIONS(LEG_ID,I)=-10000 + ENDDO + ENDIF + + ALL_SUMMED_OVER = .TRUE. + DO I=1,NEXTERNAL + IF (POLARIZATIONS(I,0).NE.-1) THEN + ALL_SUMMED_OVER = .FALSE. + EXIT + ENDIF + ENDDO + IF (ALL_SUMMED_OVER) THEN + DO I=0,5 + POLARIZATIONS(0,I)=-1 + ENDDO + ELSE + DO I=0,5 + POLARIZATIONS(0,I)=0 + ENDDO + ENDIF + + DO I=0,NEXTERNAL + DO J=0,5 + BORN_POLARIZATIONS(I,J) = POLARIZATIONS(I,J) + ENDDO + ENDDO + + + RETURN + + END + + SUBROUTINE ML5_0_SLOOPMATRIXHEL(P,HEL,ANS) + IMPLICIT NONE +C +C CONSTANTS +C + INTEGER NEXTERNAL + PARAMETER (NEXTERNAL=3) + INTEGER NSQUAREDSO + PARAMETER (NSQUAREDSO=1) + INTEGER NSQSO_BORN + PARAMETER (NSQSO_BORN=0) + +C +C ARGUMENTS +C + REAL*8 P(0:3,NEXTERNAL) + INTEGER ANS_DIMENSION + PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO)) + REAL*8 ANS(0:3,0:ANS_DIMENSION) + INTEGER HEL, USERHEL + COMMON/ML5_0_USERCHOICE/USERHEL +C ---------- +C BEGIN CODE +C ---------- + USERHEL=HEL + CALL ML5_0_SLOOPMATRIX(P,ANS) + END + + SUBROUTINE ML5_0_SLOOPMATRIXHEL_THRES(P,HEL,ANS,PREC_ASKED + $ ,PREC_FOUND,RET_CODE) + IMPLICIT NONE +C +C CONSTANTS +C + INTEGER NEXTERNAL + PARAMETER (NEXTERNAL=3) + INTEGER NSQUAREDSO + PARAMETER (NSQUAREDSO=1) +C +C ARGUMENTS +C + REAL*8 P(0:3,NEXTERNAL) + INTEGER NSQSO_BORN + PARAMETER (NSQSO_BORN=0) + + INTEGER ANS_DIMENSION + PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO)) + REAL*8 ANS(0:3,0:ANS_DIMENSION) + INTEGER HEL, RET_CODE + REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO) +C +C LOCAL VARIABLES +C + INTEGER I +C +C GLOBAL VARIABLES +C + REAL*8 USER_STAB_PREC + COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC + + INTEGER H,T,U + REAL*8 ACCURACY(0:NSQUAREDSO) + COMMON/ML5_0_ACC/ACCURACY,H,T,U + + LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY + COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY + +C ---------- +C BEGIN CODE +C ---------- + USER_STAB_PREC = PREC_ASKED + + CALL ML5_0_SLOOPMATRIXHEL(P,HEL,ANS) + IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0)) + $ THEN + BYPASS_CHECK = .TRUE. + CALL ML5_0_SLOOPMATRIXHEL(P,HEL,ANS) + BYPASS_CHECK = .FALSE. +C Make sure we correctly return an initialization-type T code + IF (T.EQ.2) T=4 + IF (T.EQ.1) T=3 + ENDIF + +C Reset it to default value not to affect next runs + USER_STAB_PREC = -1.0D0 + + DO I=0,NSQUAREDSO + PREC_FOUND(I)=ACCURACY(I) + ENDDO + RET_CODE=100*H+10*T+U + + END + + SUBROUTINE ML5_0_SLOOPMATRIX_THRES(P,ANS,PREC_ASKED,PREC_FOUND + $ ,RET_CODE) +C +C Inputs are: +C P(0:3, Nexternal) double :: Kinematic configuration +C (E,px,py,pz) +C PEC_ASKED double :: Target relative accuracy, -1 for +C default +C +C Outputs are: +C ANS(3) double :: Result (finite, single pole, +C double pole) +C PREC_FOUND double :: Relative accuracy estimated for +C the result +C Returns -1 if no stab test could be performed. +C RET_CODE integer :: Return code. See below for details +C +C Return code conventions: RET_CODE = H*100 + T*10 + U +C +C H == 1 +C Stability unknown. +C H == 2 +C Stable PS (SPS) point. +C No stability rescue was necessary. +C H == 3 +C Unstable PS (UPS) point. +C Stability rescue necessary, and successful. +C H == 4 +C Exceptional PS (EPS) point. +C Stability rescue attempted, but unsuccessful. +C +C T == 1 +C Default computation (double prec.) was performed. +C T == 2 +C Quadruple precision was used for this PS point. +C T == 3 +C MadLoop in initialization phase. Only double precision used. +C T == 4 +C MadLoop in initialization phase. Quadruple precision used. +C +C U == 0 +C Not stable. +C U == 1 +C Stable with CutTools in double precision. +C U == 2 +C Stable with PJFry++. +C U == 3 +C Stable with IREGI. +C U == 4 +C Stable with Golem95 +C U == 5 +C Stable with Samurai +C U == 6 +C Stable with Ninja in double precision +C U == 8 +C Stable with Ninja in quadruple precision +C U == 9 +C Stable with CutTools in quadruple precision. +C + IMPLICIT NONE +C +C CONSTANTS +C + INTEGER NEXTERNAL + PARAMETER (NEXTERNAL=3) + INTEGER NSQUAREDSO + PARAMETER (NSQUAREDSO=1) +C +C ARGUMENTS +C + REAL*8 P(0:3,NEXTERNAL) + INTEGER NSQSO_BORN + PARAMETER (NSQSO_BORN=0) + + INTEGER ANS_DIMENSION + PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO)) + REAL*8 ANS(0:3,0:ANS_DIMENSION) + REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO) + INTEGER RET_CODE +C +C LOCAL VARIABLES +C + INTEGER I +C +C GLOBAL VARIABLES +C + REAL*8 USER_STAB_PREC + COMMON/ML5_0_USER_STAB_PREC/USER_STAB_PREC + + INTEGER H,T,U + REAL*8 ACCURACY(0:NSQUAREDSO) + COMMON/ML5_0_ACC/ACCURACY,H,T,U + + LOGICAL BYPASS_CHECK, ALWAYS_TEST_STABILITY + COMMON/ML5_0_BYPASS_CHECK/BYPASS_CHECK, ALWAYS_TEST_STABILITY + +C ---------- +C BEGIN CODE +C ---------- + USER_STAB_PREC = PREC_ASKED + CALL ML5_0_SLOOPMATRIX(P,ANS) + IF(ALWAYS_TEST_STABILITY.AND.(H.EQ.1.OR.ACCURACY(0).LT.0.0D0)) + $ THEN + BYPASS_CHECK = .TRUE. + CALL ML5_0_SLOOPMATRIX(P,ANS) + BYPASS_CHECK = .FALSE. +C Make sure we correctly return an initialization-type T code + IF (T.EQ.2) T=4 + IF (T.EQ.1) T=3 + ENDIF + +C Reset it to default value not to affect next runs + USER_STAB_PREC = -1.0D0 + DO I=0,NSQUAREDSO + PREC_FOUND(I)=ACCURACY(I) + ENDDO + RET_CODE=100*H+10*T+U + + END + +C The subroutine below perform clean-up duties for MadLoop like +C de-allocating +C arrays + SUBROUTINE ML5_0_EXIT_MADLOOP() + CALL ML5_0_DEALLOCATE_COLOR_FLOWS() + CONTINUE + END + diff --git a/tests/unit_tests/loop/test_loop_induced_output.py b/tests/unit_tests/loop/test_loop_induced_output.py new file mode 100644 index 000000000..a9d33822e --- /dev/null +++ b/tests/unit_tests/loop/test_loop_induced_output.py @@ -0,0 +1,335 @@ +################################################################################ +# +# Copyright (c) 2009 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 +# +################################################################################ + +"""Which output formats a loop-induced ([noborn=...]) process may be sent to. + +A loop-induced process does not reach the exporters the way a [virt=...] one +does. master_interface borrows the MadLoop interface only long enough to +validate the model, then switches *back* to 'MadGraph' and calls +create_loop_induced -- so the process is exported by the ordinary tree-level +output machinery (madgraph_interface.do_output, then ExportV4Factory with +output_type='default', or ExportCPPFactory) even though the amplitude is a +LoopAmplitude and the matrix element a LoopHelasMatrixElement. + +Every format used to hand that loop matrix element to a tree-level exporter and +die deep inside it -- 'output standalone' with an IndexError in write_check_sa, +'output matrix' with "wavefunction_rank has not been computed", 'output mg7' +with a KeyError on the first loop leg, which the mg7 exporter's edge-name map +does not contain. Two things were wrong: + + * `HelasMatrixElement._flavor_enumeration_context` counted *every* motherless + wavefunction as an external leg. A LoopHelasMatrixElement's L-cut + wavefunctions are motherless too, so `get_external_flavors()` returned + tuples of length nexternal+2 -- which is what write_check_sa tripped over; + * the exporter itself was the tree-level one. 'standalone' is now routed to + the MadLoop standalone exporter (LOOP_INDUCED_FORMATS), joining 'madevent' + which has always had the LoopInducedExporterME* exporters. + +The formats with no MadLoop backend at all still cannot serve the process, and +refuse it up front pointing at [sqrvirt=...], which stays in the MadLoop +interface and yields the very same matrix element. + +The [sqrvirt=]/[virt=] tests at the end are important to keep green: the +refusal keys off the amplitude being a LoopAmplitude, and those amplitudes are +LoopAmplitudes too -- they are spared only because they never reach these two +factories. If that ever stops being true, the escape hatch the error message +recommends would be refused along with everything else. +""" + +from __future__ import absolute_import + +import os +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.interface.master_interface as MGCmd +import madgraph.iolibs.export_v4 as export_v4 +import madgraph.loop.loop_diagram_generation as loop_diagram_generation +import madgraph.loop.loop_exporters as loop_exporters +import madgraph.loop.loop_helas_objects as loop_helas_objects + +from madgraph import InvalidCmd, MG5DIR + +pjoin = os.path.join + +# The cheapest loop-induced process there is: no tree-level diagram exists for +# g g > h in loop_sm, so [noborn=QCD] gives the quark-loop amplitude alone. +LOOP_INDUCED_PROCESS = 'g g > h [noborn=QCD]' +# A 2 -> 2 one, to check the leg count is not accidentally right. +LOOP_INDUCED_PROCESS_2TO2 = 'g g > z z [noborn=QCD]' +# The same matrix element, generated the way the refusal message recommends. +SQRVIRT_PROCESS = 'g g > h [sqrvirt=QCD]' +# An ordinary virtual correction, for good measure. +VIRTUAL_PROCESS = 'u u~ > d d~ [virt=QCD]' + + + +def get_interface(process=LOOP_INDUCED_PROCESS, model='loop_sm'): + """A fresh MasterCmd with 'process' generated, ready to be output. + + Not cached: 'output' populates _curr_matrix_elements, and a second output + of a different format on the same interface would reuse them. + """ + + interface = MGCmd.MasterCmd() + interface.no_notification() + if model: + interface.exec_cmd('import model %s' % model, printcmd=False, + precmd=True) + interface.exec_cmd('generate %s' % process, printcmd=False, precmd=True) + return interface + + +# building these is the slow part and they are only read, so cache them +_matrix_elements = {} + + +def get_matrix_element(process): + """The LoopHelasMatrixElement of a loop-induced process.""" + + if process not in _matrix_elements: + # compute_loop_nc as group_subprocs builds them: building a loop ME + # without it leaves colour data that a later output picks up (a + # pre-existing cross-talk, reproducible on a clean tree) + _matrix_elements[process] = loop_helas_objects.LoopHelasProcess( + get_interface(process)._curr_amps, optimized_output=True, + compute_loop_nc=True).get_matrix_elements()[0] + return _matrix_elements[process] + + +#=============================================================================== +# TestLoopInducedExternalFlavors +#=============================================================================== +class TestLoopInducedExternalFlavors(unittest.TestCase): + """get_external_flavors() must describe the external legs only. + + The L-cut wavefunctions of a loop matrix element are motherless like the + external ones and carry number_external = nexternal+1 / nexternal+2, so they + used to be counted as two extra external legs. + """ + + def test_loop_induced_external_flavors(self): + """g g > h [noborn=QCD] has 3 external legs, not 5.""" + + me = get_matrix_element(LOOP_INDUCED_PROCESS) + self.assertEqual(me.get_nexternal_ninitial(), (3, 2)) + flavors, pdgs = me.get_external_flavors(return_pdgs=True) + self.assertEqual([tuple(f) for f in flavors], [(1, 1, 1)]) + self.assertEqual([list(p) for p in pdgs], [[21, 21, 25]]) + + def test_loop_induced_external_flavors_4legs(self): + """Same for a 2 -> 2 loop-induced process.""" + + me = get_matrix_element(LOOP_INDUCED_PROCESS_2TO2) + self.assertEqual(me.get_nexternal_ninitial(), (4, 2)) + flavors, pdgs = me.get_external_flavors(return_pdgs=True) + self.assertEqual([tuple(f) for f in flavors], [(1, 1, 1, 1)]) + self.assertEqual([list(p) for p in pdgs], [[21, 21, 23, 23]]) + + +#=============================================================================== +# TestLoopInducedOutput +#=============================================================================== +class TestLoopInducedOutput(unittest.TestCase): + """The output formats a loop-induced process may and may not be sent to.""" + + def setUp(self): + # 'output' chdir's into the directory it writes; come back before it is + # removed, or the next MasterCmd() cannot even call os.getcwd() + os.chdir(MG5DIR) + self.tmpdir = tempfile.mkdtemp(prefix='loop_induced_output') + + def tearDown(self): + os.chdir(MG5DIR) + shutil.rmtree(self.tmpdir, ignore_errors=True) + + #=========================================================================== + # helpers + #=========================================================================== + def get_exporter(self, interface, format='standalone'): + """The exporter ExportV4Factory picks for 'format'.""" + + interface._export_format = format + interface._export_dir = pjoin(self.tmpdir, 'factory') + return export_v4.ExportV4Factory(interface, True, + group_subprocesses=False) + + def assert_output_refused(self, format): + """'output ' must raise InvalidCmd, and say what to do instead. + + Anything that is not an InvalidCmd -- an IndexError, a KeyError, a + MadGraph5Error from inside a tree-level exporter -- means the loop + matrix element reached an exporter that cannot write it. + """ + + interface = get_interface() + out_dir = pjoin(self.tmpdir, format.replace(' ', '_')) + try: + interface.exec_cmd('output %s %s -f' % (format, out_dir), + printcmd=False, precmd=True) + except InvalidCmd as error: + self.assertTrue('sqrvirt' in str(error), + "'output %s' refused %s without pointing at " + "[sqrvirt=]: %s" % (format, LOOP_INDUCED_PROCESS, + error)) + return + except Exception as error: + raise AssertionError( + "'output %s' crashed on %s with %s: %s" + % (format, LOOP_INDUCED_PROCESS, type(error).__name__, error)) + raise AssertionError( + "'output %s' silently accepted %s; no exporter for that format can " + "write a LoopHelasMatrixElement" % (format, LOOP_INDUCED_PROCESS)) + + def assert_output_succeeds(self, process, format='standalone'): + """'output ' must go through, and write a MadLoop directory.""" + + interface = get_interface(process) + out_dir = pjoin(self.tmpdir, 'ok') + interface.exec_cmd('output %s %s -f' % (format, out_dir), + printcmd=False, precmd=True) + # A MadLoop output, not a tree-level one: this file only exists when a + # loop exporter ran. + self.assertTrue( + os.path.exists(pjoin(out_dir, 'Cards', 'MadLoopParams.dat')), + '%s did not produce a MadLoop output' % process) + + #=========================================================================== + # 'standalone' is served by the MadLoop standalone exporter + #=========================================================================== + def test_standalone_factory_uses_the_madloop_exporter(self): + """ExportV4Factory must not hand a loop-induced process to the + tree-level standalone exporter.""" + + interface = get_interface() + self.assertTrue(isinstance(interface._curr_amps[0], + loop_diagram_generation.LoopAmplitude), + "%s did not produce a LoopAmplitude" + % LOOP_INDUCED_PROCESS) + self.assertIsInstance(self.get_exporter(interface), + loop_exporters.LoopProcessExporterFortranSA) + + def test_standalone_factory_unchanged_for_tree(self): + """A tree process must still get the tree-level exporter.""" + + interface = get_interface('g g > t t~', model='sm') + exporter = self.get_exporter(interface) + self.assertIsInstance(exporter, export_v4.ProcessExporterFortranSA) + self.assertFalse(isinstance( + exporter, loop_exporters.LoopProcessExporterFortranSA)) + + def test_output_standalone_accepts_loop_induced(self): + """'output standalone' used to die with an IndexError in + write_check_sa; it now writes a MadLoop standalone directory.""" + + self.assert_output_succeeds(LOOP_INDUCED_PROCESS) + + def test_no_model_imported(self): + """validate_model is called before create_loop_induced's check_add, so + it has to bootstrap the model itself.""" + + interface = MGCmd.MasterCmd() + interface.no_notification() + self.assertFalse(interface._curr_model) + interface.exec_cmd('generate %s' % LOOP_INDUCED_PROCESS, + printcmd=False, precmd=True) + self.assertEqual(interface._curr_model.get('name'), 'loop_sm') + + #=========================================================================== + # ... the formats with no MadLoop backend refuse it + #=========================================================================== + def test_loop_induced_formats_are_matched_exactly(self): + """'standalone' is a prefix of formats that have no loop backend; the + allow-list must never be tested with startswith.""" + + for format in ['standalone_cpp', 'standalone_mg7', 'standalone_msP', + 'standalone_msF', 'standalone_rw']: + self.assertNotIn(format, export_v4.LOOP_INDUCED_FORMATS) + + def test_output_matrix_refuses_loop_induced(self): + """'output matrix' shares the standalone branch of the factory; it used + to die with "wavefunction_rank has not been computed".""" + + self.assert_output_refused('matrix') + + def test_output_mg7_refuses_loop_induced(self): + """'output mg7' -- the default format -- used to die with a KeyError on + the first loop leg. mg7/madmatrix has no MadLoop backend at all.""" + + self.assert_output_refused('mg7') + + def test_output_standalone_cpp_refuses_loop_induced(self): + """The prefix trap: standalone_cpp must not ride on 'standalone'.""" + + self.assert_output_refused('standalone_cpp') + + def test_output_standalone_msP_refuses_loop_induced(self): + """Same for the MadSpin standalone variants.""" + + self.assert_output_refused('standalone_msP') + + def test_refusal_leaves_an_existing_directory_alone(self): + """The refusal comes before the rmtree that cleans an existing output + directory; an already existing one must survive it.""" + + out_dir = pjoin(self.tmpdir, 'existing') + os.mkdir(out_dir) + sentinel = pjoin(out_dir, 'sentinel.txt') + open(sentinel, 'w').write('do not delete me\n') + + interface = get_interface() + self.assertRaises(InvalidCmd, interface.exec_cmd, + 'output mg7 %s -f' % out_dir) + self.assertTrue(os.path.exists(sentinel), + 'the refused output wiped the existing directory') + + #=========================================================================== + # ... and the routes that already worked must keep working + #=========================================================================== + def test_output_madevent_still_accepts_loop_induced(self): + """madevent is the format loop-induced processes are *for*. + + It has the LoopInducedExporterMEGroup / ...MENoGroup exporters, so it + must go straight through the refusal above. + """ + + interface = get_interface() + out_dir = pjoin(self.tmpdir, 'me') + interface.exec_cmd('output madevent %s -f' % out_dir, + printcmd=False, precmd=True) + self.assertTrue( + os.path.isdir(pjoin(out_dir, 'SubProcesses', 'MadLoop5_resources')), + 'madevent did not produce a loop-induced output') + + def test_sqrvirt_standalone_output_still_works(self): + """The escape hatch the refusal message recommends must actually work. + + [sqrvirt=] gives a LoopAmplitude just like [noborn=] does; it is spared + the refusal only because master_interface keeps it in the MadLoop + interface, which passes output_type='madloop'. + """ + + self.assert_output_succeeds(SQRVIRT_PROCESS) + + def test_virt_standalone_output_still_works(self): + """Same for an ordinary [virt=] process.""" + + self.assert_output_succeeds(VIRTUAL_PROCESS) From 7fc627c09dd4f0c12b9368b33be529734a4f2a7a Mon Sep 17 00:00:00 2001 From: Olivier Mattelaer Date: Sat, 22 Aug 2026 00:10:55 +0200 Subject: [PATCH 6/6] refuse a loop-induced --me_exporter= before the output dir is cleaned The guard in do_output tested only self._export_format, so a refusable format asked for as the *second* matrix-element exporter fell through to the factory backstop -- which runs after the shutil.rmtree that cleans an existing output directory: output madevent --me_exporter=mg7 -f # g g > h [noborn=QCD] -> InvalidCmd from ExportCPPFactory, and already deleted Not a regression (a clean tree wipes it too, then crashes deeper), but holding that directory is the entire reason the check sits above the rmtree. Check both names in the same loop. A plugin me_exporter is named by its plugin, so it is refused like any other unlisted name -- which is what the backstops already do with it. Also trim the comments back to the house rule (they outnumbered the product code): drop the block comment on LOOP_INDUCED_FORMATS restated by the docstring under it, keeping the exact-membership warning, which is the load-bearing part. The bare `except Exception` around perturbation_couplings is gone with it: both callers pass a Process or ProcessDefinition, which always has that key. And note why the 'matchbox' factory branch carries no backstop, unlike its neighbours. Co-Authored-By: Claude Opus 5 --- madgraph/interface/madgraph_interface.py | 17 +++++--- madgraph/iolibs/export_cpp.py | 7 +--- madgraph/iolibs/export_v4.py | 40 +++++-------------- .../loop/test_loop_induced_output.py | 27 ++++++++++--- 4 files changed, 44 insertions(+), 47 deletions(-) diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py index 5e5ff0c5a..ea8873ff2 100755 --- a/madgraph/interface/madgraph_interface.py +++ b/madgraph/interface/madgraph_interface.py @@ -9810,13 +9810,18 @@ def do_output(self, line): # create_loop_induced), but only the formats in LOOP_INDUCED_FORMATS # have a loop backend to route it to. Refuse the others here, ahead of # the directory cleaning just below, so that a guaranteed refusal never - # deletes an existing output directory first. The exporter factories - # carry the same check as a backstop. - if self._export_format not in export_v4.LOOP_INDUCED_FORMATS and \ - self._curr_amps and isinstance(self._curr_amps[0], + # deletes an existing output directory first. The factories carry the + # same check, but they run after that cleaning. + if self._curr_amps and isinstance(self._curr_amps[0], loop_diagram_generation.LoopAmplitude): - raise self.InvalidCmd(export_v4.loop_induced_not_supported_msg( - self._export_format, self._curr_amps[0].get('process'))) + # --me_exporter= writes into the same directory, so it has to be + # checked here too + for format in [self._export_format, + options['me_exporter'].get('name')]: + if format and format not in export_v4.LOOP_INDUCED_FORMATS: + raise self.InvalidCmd( + export_v4.loop_induced_not_supported_msg( + format, self._curr_amps[0].get('process'))) # check if os.path.realpath(self._export_dir) == os.getcwd(): diff --git a/madgraph/iolibs/export_cpp.py b/madgraph/iolibs/export_cpp.py index 7561e961c..076944005 100755 --- a/madgraph/iolibs/export_cpp.py +++ b/madgraph/iolibs/export_cpp.py @@ -3562,11 +3562,8 @@ def ExportCPPFactory(cmd, group_subprocesses=False, cmd_options={}): opt['output_options'] = cmd_options cformat = cmd._export_format - # None of the C++ exporters below has a MadLoop backend, so a loop-induced - # process would reach them as a LoopHelasMatrixElement whose loop legs they - # cannot even index (the mg7 exporter builds its edge names from the - # external legs alone). Refuse it here instead. Plugins are left alone: - # they are free to implement their own loop support. + # No C++ exporter has a MadLoop backend (the mg7 one cannot even index the + # loop legs: it builds its edge names from the external legs alone). if cformat not in export_v4.LOOP_INDUCED_FORMATS and cmd._curr_amps and \ isinstance(cmd._curr_amps[0], loop_diagram_generation.LoopAmplitude): raise InvalidCmd(export_v4.loop_induced_not_supported_msg( diff --git a/madgraph/iolibs/export_v4.py b/madgraph/iolibs/export_v4.py index b86ff5e1d..a1561f7e0 100755 --- a/madgraph/iolibs/export_v4.py +++ b/madgraph/iolibs/export_v4.py @@ -11503,34 +11503,16 @@ def create_param_card(self, write_special=True): mssm_convert=True, write_special=write_special) -# The output formats that can serve a loop-induced ([noborn=]) process coming -# through the tree-level do_output: 'madevent' has the LoopInducedExporterME* -# exporters, 'standalone' is routed to the MadLoop standalone exporter, and a -# plugin is free to bring its own. Every other format sends the -# LoopHelasMatrixElement to a tree-level exporter that cannot write it. -# Membership must be tested exactly: 'standalone' is a prefix of -# standalone_cpp / _mg7 / _msP / _msF / _rw, none of which has a loop backend. +# Output formats with a loop backend for a loop-induced ([noborn=]) process +# coming through the tree-level do_output. Test membership EXACTLY: 'standalone' +# is a prefix of standalone_cpp / _mg7 / _msP / _msF / _rw, which have none. LOOP_INDUCED_FORMATS = ['madevent', 'plugin', 'standalone'] def loop_induced_not_supported_msg(format, process=None): - """Error text for an output format that has no MadLoop backend. + """Refusal text for a format that cannot write a LoopHelasMatrixElement.""" - A loop-induced ([noborn=]) process is exported by the *tree-level* output - machinery: master_interface only borrows the MadLoop interface to validate - the model, then switches back to 'MadGraph' and calls create_loop_induced. - So a format whose exporter cannot write a LoopHelasMatrixElement has to say - so here rather than let the tree-level exporter fail deep inside. - - The same matrix element is available through [sqrvirt=], which does stay in - the MadLoop interface and therefore reaches the MadLoop exporters. - """ - - orders = 'QCD' - if process: - try: - orders = ' '.join(process.get('perturbation_couplings')) or orders - except Exception: - pass + orders = (' '.join(process.get('perturbation_couplings')) if process + else '') or 'QCD' return """The '%(format)s' output format does not support loop-induced processes. Generate the process with [sqrvirt=%(orders)s] rather than [noborn=%(orders)s] to obtain the @@ -11686,14 +11668,10 @@ def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True opt['madanalysis5'] = cmd.options['madanalysis5_path'] if format == 'matrix' or format.startswith('standalone'): - # A loop-induced ([noborn=]) process reaches this factory through the - # MadGraph interface (master_interface switches back to it after - # generation), but ProcessExporterFortranSA cannot write a - # LoopHelasMatrixElement. if cmd._curr_amps and isinstance( cmd._curr_amps[0], loop_diagram_generation.LoopAmplitude): - # Only plain 'standalone' has a MadLoop backend to route to, as - # the madevent branches below have for their own format. + # of the formats sharing this branch only 'standalone' has a + # MadLoop backend; ProcessExporterFortranSA has none if format not in LOOP_INDUCED_FORMATS: raise InvalidCmd( loop_induced_not_supported_msg(format, curr_proc)) @@ -11731,6 +11709,8 @@ def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True else: return ProcessExporterFortranME(cmd._export_dir,opt) elif format in ['matchbox']: + # no loop-induced backstop needed: do_output refuses 'matchbox' + # before any factory runs, and loop_interface never comes here return ProcessExporterFortranMatchBox(cmd._export_dir,opt) elif cmd._export_format in ['madweight'] and group_subprocesses: diff --git a/tests/unit_tests/loop/test_loop_induced_output.py b/tests/unit_tests/loop/test_loop_induced_output.py index a9d33822e..fafe88ad0 100644 --- a/tests/unit_tests/loop/test_loop_induced_output.py +++ b/tests/unit_tests/loop/test_loop_induced_output.py @@ -286,9 +286,12 @@ def test_output_standalone_msP_refuses_loop_induced(self): self.assert_output_refused('standalone_msP') - def test_refusal_leaves_an_existing_directory_alone(self): - """The refusal comes before the rmtree that cleans an existing output - directory; an already existing one must survive it.""" + def assert_refusal_keeps_directory(self, line): + """A refused 'output' must not delete the directory it would write to. + + This is why do_output checks ahead of its own rmtree: the factories + refuse too, but they only run once the directory is already gone. + """ out_dir = pjoin(self.tmpdir, 'existing') os.mkdir(out_dir) @@ -296,10 +299,22 @@ def test_refusal_leaves_an_existing_directory_alone(self): open(sentinel, 'w').write('do not delete me\n') interface = get_interface() - self.assertRaises(InvalidCmd, interface.exec_cmd, - 'output mg7 %s -f' % out_dir) + self.assertRaises(InvalidCmd, interface.exec_cmd, line % out_dir) self.assertTrue(os.path.exists(sentinel), - 'the refused output wiped the existing directory') + '%s wiped the existing directory' % (line % out_dir)) + + def test_refusal_leaves_an_existing_directory_alone(self): + """The refused format is the one asked for.""" + + self.assert_refusal_keeps_directory('output mg7 %s -f') + + def test_me_exporter_refusal_leaves_an_existing_directory_alone(self): + """--me_exporter= writes into the same directory, so it has to be + refused by the same early check; testing _export_format alone let it + fall through to the factory, which runs after the rmtree.""" + + self.assert_refusal_keeps_directory( + 'output madevent %s --me_exporter=mg7 -f') #=========================================================================== # ... and the routes that already worked must keep working