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) 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']: diff --git a/madgraph/interface/madgraph_interface.py b/madgraph/interface/madgraph_interface.py index e7cd08fcf..ea8873ff2 100755 --- a/madgraph/interface/madgraph_interface.py +++ b/madgraph/interface/madgraph_interface.py @@ -9806,6 +9806,23 @@ 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 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): + # --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(): if len(args) == 0: diff --git a/madgraph/iolibs/export_cpp.py b/madgraph/iolibs/export_cpp.py index 1b83b6f67..076944005 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,14 @@ def ExportCPPFactory(cmd, group_subprocesses=False, cmd_options={}): opt = dict(cmd.options) opt['output_options'] = cmd_options cformat = cmd._export_format - + + # 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( + 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 459e43a1d..a1561f7e0 100755 --- a/madgraph/iolibs/export_v4.py +++ b/madgraph/iolibs/export_v4.py @@ -11502,7 +11502,23 @@ def create_param_card(self, write_special=True): rule_card_path=rule_card, mssm_convert=True, write_special=write_special) - + +# 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): + """Refusal text for a format that cannot write a LoopHelasMatrixElement.""" + + 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 +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. @@ -11652,8 +11668,28 @@ def ExportV4Factory(cmd, noclean, output_type='default', group_subprocesses=True opt['madanalysis5'] = cmd.options['madanalysis5_path'] if format == 'matrix' or format.startswith('standalone'): + if cmd._curr_amps and isinstance( + cmd._curr_amps[0], loop_diagram_generation.LoopAmplitude): + # 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)) + 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): @@ -11673,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/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, 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..fafe88ad0 --- /dev/null +++ b/tests/unit_tests/loop/test_loop_induced_output.py @@ -0,0 +1,350 @@ +################################################################################ +# +# 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 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) + sentinel = pjoin(out_dir, 'sentinel.txt') + open(sentinel, 'w').write('do not delete me\n') + + interface = get_interface() + self.assertRaises(InvalidCmd, interface.exec_cmd, line % out_dir) + self.assertTrue(os.path.exists(sentinel), + '%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 + #=========================================================================== + 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)