C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                  ---------------------------
C                  beowulfsubs.f  Version 1.0
C                  ---------------------------
C
C23456789012345678901234567890123456789012345678901234567890123456789012
      SUBROUTINE VERSION
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      WRITE(NOUT,*)'beowulf 1.0 subroutines 30 September 1999'
      RETURN
      END
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C Subroutines for numerical integration of jet cross sections in
C electron-positron annihilation. -- D. E. Soper
C
C First code: 29 November 1992.
C Latest revision: see Version subroutine above.
C Special note: modified 8 December 1999 to change tabs to spaces
C and to correct the header in line 4 above.
C
C The main program and subroutines that a user might want to modify
C are contained in the companion package, beowulf.f. In particular, a
C user may very well want to modify parameter settings in the main
C program and to change the observables calculated in the subroutine
C HROTHGAR and in the functions CALStype(NCUT,KCUT,index). Subroutines
C that can be modified only at extreme peril to the reliability of
C the results are in this package, beowulfsubs.f.
C
C There are two parallel calculations. Program beowulf calculates a
C sample integral, which by default is the average value of 
C (1 - thrust)^2. These are summed in the variable INTEGRAL and
C reported upon completion of the program. The program also computes
C a simple check integral in order to check on the jacobians etc.
C In the meantime, for each point in loop space and each final
C state cut, the program reports the corresponding point in the space
C of final state momenta along with the corresponding weight (Feynman
C diagram times jacobian factors) to the subroutine HROTHGAR, which
C multiplies by the measurement functions CALS corresponding to the
C measurements desired and accumulates the results.
C 
C In order to control roundoff errors, a point in loop space is rejected
C if the point is too near a singularity or if there is too much
C cancellation in the contribution from that point to INTEGRAL.
C 
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C                        PROGRAM STRUCTURE
C
C * denotes routines found in beowulf.f. 
C Other routines are in beowulfsubs.f
C
C PROGRAM BEOWULF (*) <setup and final reporting>
C    SUBROUTINE MAKECUTINFO
C       SUBROUTINE NEWGRAPH
C          SUBROUTINE NEWCHOICE
C          SUBROUTINE NEXTCHOICE
C          FUNCTION ONEPI
C          SUBROUTINE CHECK
C             SUBROUTINE CHECKOUT
C             SUBROUTINE EXCHANGE
C       SUBROUTINE NEWCUT
C       SUBROUTINE FINDQS
C          FUNCTION DEPENDENT
C             SUBROUTINE FINDA <see below>
C          SUBROUTINE SOFTHELP
C             SUBROUTINE NEWCUT
C          SUBROUTINE FINDA <see below>
C             FUNCTION PROPSIGN
C          SUBROUTINE NEWCUT
C    SUBROUTINE GETCOUNTS (*)
C    SUBROUTINE HROTHGAR (*) <see below>
C    SUBROUTINE RANDOMINIT
C       SUBROUTINE NEWRAN
C    SUBROUTINE DAYTIME (*)
C    SUBROUTINE VERSION
C    SUBROUTINE RENO <see below>
C    SUBROUTINE DIAGNOSTIC
C     SUBROUTINE NEWGRAPH ...
C        SUBROUTINE FINDQS ...
C       SUBROUTINE CHECKPOINT
C       SUBROUTINE CALCULATE <see below>
C end
C
C SUBROUTINE RENO <choses graphs and points>
C    SUBROUTINE TIMING (*)
C    SUBROUTINE HROTHGAR (*) <see below>
C    SUBROUTINE NEWGRAPH ...
C    SUBROUTINE FINDQS ...
C    SUBROUTINE FINDA <see below>
C       FUNCTION PROPSIGN
C    SUBROUTINE NEWPOINT
C       SUBROUTINE AXES
C       FUNCTION RANDOM
C          SUBROUTINE NEWRAN
C    SUBROUTINE CHECKPOINT
C    SUBROUTINE CALCULATE <see below>
C return
C
C SUBROUTINE CALCULATE <the main calculation routines.
C    SUBROUTINE REFLECT <not used by default>
C       SUBROUTINE REFLECTSPEC
C    FUNCTION DENSITY
C    FUNCTION CALS0 (*)
C       FUNCTION THRUST (*)
C    SUBROUTINE GETCUTINFO
C    SUBROUTINE DEFORM
C    SUBROUTINE FINDA <see below>
C    SUBROUTINE CHECKDEFORM
C    FUNCTION RNUMERATOR
C       SUBROUTINE QPROPR
C       SUBROUTINE GPROPR
C       SUBROUTINE QPROP
C          SUBROUTINE CHECKDEFORM2
C       FUNCTION RQQP3AQ
C       FUNCTION RQQP3BQ
C       FUNCTION RQQP3CQ
C       FUNCTION RQQG3AG
C       FUNCTION RQQG3BG
C       FUNCTION RQQG3CG
C       FUNCTION RQQG3AQ
C       FUNCTION RQQG3BQ
C       FUNCTION RQQG3CQ
C    FUNCTION NUMERATOR
C       SUBROUTINE QPROP
C          SUBROUTINE CHECKDEFORM2
C       SUBROUTINE GPROP
C    FUNCTION SMEAR
C    SUBROUTINE CHECKCALC
C    SUBROUTINE HROTHGAR (*) <see below>
C return
C 
C SUBROUTINE FINDA <finds matrix to get all k's from the loop l's>
C    FUNCTION PROPSIGN
C return
C 
C SUBROUTINE HROTHGAR (*) 
C  Beowulf serves Hothgar, who accepts the points in the space of final
C  state momenta with the corresponding weights, multiplies by
C  desired measurement functions CALS, and accumulates results.
C    FUNCTION CALSTHRUST (*)
C       FUNCTION THRUST (*)
C       FUNCTION KN (*)
C    FUNCTION CALS3JET (*)
C       SUBROUTINE COMBINEJETS (*)
C       FUNCTION BETHKE (*)
C    FUNCTION KN (*)
C    FUNCTION BETHKE (*)
C return
C
C Simple functions called from routines above, with calls
C not listed above:
C 
C FUNCTION XXREAL
C FUNCTION XXIMAG
C FUNCTION COMPLEXSQRT
C FUNCTION FACTORIAL
C FUNCTION SINHINV
C FUNCTION DELTA
C 
C23456789012345678901234567890123456789012345678901234567890123456789012
C
C A brief introduction to the variables used:
C
C Size of the calculation:
C NLOOPS = number of loops (in cut photon self energy graph).
C NPROPS = number of propagators in graph, = 3 * NLOOPS - 1.
C NVERTS = number of vertices in graph, = 2 * NLOOPS.
C CUTMAX = NLOOPS + 1
C        = maximum number of cut propagators;
C        = number of independent loop momenta needed to determine the
C          propagator momenta, counting the virtual photon momentum.
C The current program is restricted to 0 and 1 virtual loops.
C
C Labels:
C L    = index of loop momenta, L = 0,1,...,NLOOPS.
C        L = 0 normally denontes the virtual photon momentum.
C P    = index of propagator, P = 0,1,...,NPROPS.
C        P = 0 denotes the virtual photon momentum.
C Q(L) = index P of propagator carrying the Lth loop momentum.
C V    = index of vertices, V = 1,...,NVERTS
C
C Momentum variables (MU = 0,1,2,3):
C K(P,MU)       = Momentum of Pth propagator.
C                 For P = 0, this is the virtual photon momentum:
C                 K(0,MU) = 0 for MU = 1, 2, 3 while K(0,0) = RTS.
C ABSK(P)       = Square of the three momentum of Pth propagator.
C KINLOOP(J,MU) = K(LOOPINDEX(J),MU) = momenta of loop propagators.
C KCUT(I,MU)    = K(CUTINDEX(I),MU) = momenta of cut propagators.
C K(Q(L),MU)    = Lth loop momentum, L = 0,...,NLOOPS;
C KC(P,MU)      = complex propagator momenta.
C A(P,L)        = Matrix relating propagator momenta to loop momenta.
C                 K(P,MU) = SUM_{L=0}^{NLOOPS} A(P,L) K(Q(L),MU)
C
C Variables from NEWGRAPH:
C VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of
C              of propagator P. Specifies the supergraph.
C PROP(V,I) = Index of Ith propagator attached to vertex V, I = 1,2,3.
C              Also specifies the supergraph.
C SELFPROP(P) = True if propagator P is part of a one loop self-energy
C               subgraph or attaches to a such a subgraph.
C
C Variables associated with NEWPOINT and FINDQS:
C NMAPS = Number of different maps from random x's to momenta.
C MAPNMUMBER = Number labelling a certain map.
C QS(MAPNUMBER,II) = Label of the IIth propagator that is special
C                   in map number MAPNUMBER.
C MAPTYPES(MAPNUMBER) = IRFULL, IRPART, MEDIUM, or UVHARD.
C A1S(MAPNUMBER,L),A2S(MAPNUMBER,L) specify the momenta of the special
C    propagators needed for an IRFULL map:
C           P1(mu) = sum_{L=2}^N A1(L) K(Q(L),mu)
C           P2(mu) = sum_{L=2}^N A2(L) K(Q(L),mu)
C
C JACNEWPOINT =1/DENSITY(GRAPHNUMBER,K,ABSK,QS,A1S,A2S,MAPTYPES,NMAPS)
C             = Jacobian for loop momenta L.
C
C Variables from NEWCUT:
C NEWCUTINIT: .TRUE. tells NEWCUT to initialize itself.
C NCUT = Number of cut propagators.
C ISIGN(P) = +1 if propagator P is left of cut, -1 if right, 0 if cut.
C CUTINDEX(I) = Index P of cut propagator I, I = 1,...,CUTMAX.
C CUTSIGN(I) = Sign of cut propagator I I = 1,...,CUTMAX.
C              (+1 if K(P,0) >0 for cut propagator.)
C LEFTLOOP = True iff there is a virtual loop to the left of the cut.
C RIGHTLOOP = True iff there is a virtual loop to the right of the cut.
C NINLOOP = Number of propagators in loop.
C LOOPINDEX(NP) = Index P of NPth propagator around the loop.
C LOOPSIGN(NP) = 1 if propagator direction is same as loop direction.
C               -1 if  direction is opposite to loop direction.
C NP = JCUT: Propagator cut by loopcut.
C CUTFOUND: .TRUE. if NEWCUT found a new cut.
C
C In RENO we use CUTINDEX to define CUT(P) = True if propagator
C    P is cut.
C
C Solving for the propagator energies:
C For NCUT = CUTMAX, cut propagators are P = CUTINDEX(I).
C   with direction of positive energy given by CUTSIGN(I).
C For NCUT = CUTMAX - 1, we define a "loopcut" on the propagator
C   numbered JCUT in order around the loop, 1.LE.JCUT.LE.NINLOOP:
C         CUTINDEX(CUTMAX) = LOOPINDEX(JCUT) and
C         CUTSIGN(CUTMAX) = LOOPSIGN(JCUT).
C Energies of cut propagators are
C   E(I-1) = K(CUTINDEX(I),0) for I = 1,...,CUTMAX.
C and are determined from
C   E(I-1) = CUTSIGN(I) * SQRT( Sum_J [  K(CUTINDEX(I),J)**2 ] ).
C This gives energies E(L) for L = 0,...,NLOOPS. We consider the
C propagators designated by QE(L) = CUTINDEX(L+1) as independent
C and generate the matrix AE(P,L) that gives the propagator energies
C in terms of these independent momenta. This gives the propagator
C energies.
C
C Contour deformation:
C   NEWKINLOOP(MU) =  addition to the momentum going around the loop
C                     caused by deforming the contour. We have
C   Im[ KC(LOOPINDEX(J,MU)) ] = LOOPSIGN(LOOPINDEX(J)) 
C                              * Im[ NEWKINLOOP(J,MU) ] for MU = 1,2,3.
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE RENO(
     >          SUMR,ERRORR,SUMI,ERRORI,
     >          SUMBIS,ERRORBIS,
     >          SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI,FLUCT,
     >          INCLUDED,EXTRACOUNT,OMITTED,
     >          NVALPT,VALPTMAX,KBAD,BADGRAPHNUMBER,BADMAPNUMBER,
     >          NRENO,CPUTIME)
C Array sizes:
      INTEGER SIZE,MAXGRAPHS,MAXMAPS
      PARAMETER (SIZE = 3)
      PARAMETER (MAXGRAPHS = 10)
      PARAMETER (MAXMAPS = 64)
C Out:
      REAL*8 SUMR,ERRORR,SUMI,ERRORI
      REAL*8 SUMBIS,ERRORBIS
      REAL*8 SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI
      REAL*8 FLUCT(MAXGRAPHS,MAXMAPS)
      INTEGER*8 INCLUDED,EXTRACOUNT,OMITTED
      INTEGER NVALPT(-9:6)
      REAL*8 VALPTMAX
      REAL*8 KBAD(0:3*SIZE-1,0:3)
      INTEGER BADGRAPHNUMBER,BADMAPNUMBER
      INTEGER NRENO
      REAL*8 CPUTIME
C
C Computes the cross section integral by Monte Carlo integration.
C
C Latest revision 11 August 1999
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      REAL*8 ENERGYSCALE
      COMMON /MSCALE/ ENERGYSCALE
      REAL*8 BADNESSLIMIT,CANCELLIMIT,THRUSTCUT
      COMMON /LIMITS/ BADNESSLIMIT,CANCELLIMIT,THRUSTCUT
      REAL*8 TIMELIMIT
      COMMON /MAXTIME/ TIMELIMIT
C Graphs to include
      LOGICAL USEGRAPH(MAXGRAPHS)
      COMMON /WHICHGRAPHS/ USEGRAPH
C How many graphs and how many cuts and maps for each:
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C Momenta:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
C Matrices:
      INTEGER A(0:3*SIZE-1,0:SIZE)
C NEWGRAPH variables:
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
      LOGICAL SELFPROP(3*SIZE-1)
      LOGICAL GRAPHFOUND
      INTEGER GRAPHNUMBER
C FINDA variable:
      LOGICAL QOK
C MAP variables:
      CHARACTER*6 MAPTYPE
      INTEGER NMAPS,MAPNUMBER
      INTEGER QS(256,0:SIZE)
      INTEGER Q(0:SIZE)
      INTEGER A1S(256,2:SIZE),A2S(256,2:SIZE)
      INTEGER A1(2:SIZE),A2(2:SIZE)
      CHARACTER*6 MAPTYPES(256)
C Variable from CHECKPOINT:
      REAL*8 BADNESS
C Logical variables to tell how to treat point:
      LOGICAL XTRAPOINTQ, BADPOINTQ
C Functions:
      REAL*8 XXREAL,XXIMAG
C Index variables:
      INTEGER L,P,MU
C Hrothgar dummy variables:
      REAL*8 KCUT0(SIZE+1,0:3)
C Reno size and counting variables:
      INTEGER GROUPSIZE(MAXGRAPHS,MAXMAPS)
      INTEGER GROUPSIZEGRAPH(MAXGRAPHS)
      INTEGER GROUPSIZETOTAL
      COMMON /MONTECARLO/GROUPSIZE,GROUPSIZEGRAPH,GROUPSIZETOTAL
      INTEGER POINT
C Reno results variables:
      REAL*8 SQRSUMR,SQRSUMCHKR
      REAL*8 SQRSUMI,SQRSUMCHKI
      REAL*8 SQRSUMBIS
      COMPLEX*16 INTEGRAL,INTEGRALCHK
      REAL*8 INTEGRALBIS
C Calculate variables:
      COMPLEX*16 VALUE,VALUECHK
      REAL*8 MAXPART
      REAL*8 VALPT,LOGVALPT
      LOGICAL REPORT,DETAILS
      COMMON /CALCULOOK/ REPORT,DETAILS
C Timing variables
      REAL*8 DELTATIME
C
C------------------------------ Begin ----------------------------------
C
C Dummy variables for Hrothgar.
C
      DO L = 1,SIZE+1
        DO MU = 0,3
         KCUT0(L,MU) = 1.0D0
        ENDDO
      ENDDO
C
C Initialize CPUTIME and NRENO. Call to TIMING starts the clock.
C
      CPUTIME = 0.0
      NRENO = 0
      CALL TIMING(DELTATIME)
C
C Initialize sums for loop over groups of Reno points. The sums
C will be updated for each group. Within a group, the quantities
C corresponding to SUMxxR + i SUMxxI are complex variables called
C INTEGRALxx.
C
      SUMR = 0.0D0
      SUMI = 0.0D0
      SUMBIS = 0.0D0
      SUMCHKR = 0.0D0
      SUMCHKI = 0.0D0
C
      SQRSUMR = 0.0D0
      SQRSUMI = 0.0D0
      SQRSUMBIS = 0.0D0
      SQRSUMCHKR = 0.0D0
      SQRSUMCHKI = 0.0D0
C
      DO GRAPHNUMBER = 1,NUMBEROFGRAPHS
      DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER)
        FLUCT(GRAPHNUMBER,MAPNUMBER) = 0.0D0
      ENDDO
      ENDDO
C
      DO L = -9,6
        NVALPT(L) = 0
      ENDDO
      VALPTMAX = 0.0D0
      INCLUDED = 0
      EXTRACOUNT = 0
      OMITTED = 0
C
C Tell CALCULATE not to report its findings for each calculation
C
      REPORT = .FALSE.
C
C Initialize integrals for first group.
C
      INTEGRAL = (0.0D0,0.0D0)
      INTEGRALBIS = 0.0D0
      INTEGRALCHK = (0.0D0,0.0D0)
C
C Loop over groups of points.
C
      DO WHILE (CPUTIME.LT.TIMELIMIT)
      NRENO = NRENO + 1
C
C Call Hrothgar to tell him to that we are starting a new group.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'STARTGROUP')
C
C Get a new graph.
C
      GRAPHFOUND = .TRUE.
      GRAPHNUMBER = 0
C
      DO WHILE (GRAPHFOUND)
      CALL NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND)
      IF (GRAPHFOUND) THEN
      GRAPHNUMBER = GRAPHNUMBER + 1
C
C Calculate number of maps NMAPS, index arrays QS,
C types MAPTYPES, and matrices A1S and A2S
C associated with the maps.
C
      CALL FINDQS(VRTX,PROP,NMAPS,QS,A1S,A2S,MAPTYPES)
      NUMBEROFMAPS(GRAPHNUMBER) = NMAPS
C
C Check if we were supposed to use this graph (USEGRAPH is
C set in the main program.)
C
      IF (USEGRAPH(GRAPHNUMBER)) THEN
C
C Loop over choices of maps from x's to loop momenta.
C
      DO MAPNUMBER = 1,NMAPS
C
      MAPTYPE = MAPTYPES(MAPNUMBER)
      DO L = 0,NLOOPS
        Q(L) = QS(MAPNUMBER,L)
      ENDDO
      DO L = 2,NLOOPS
        A1(L) = A1S(MAPNUMBER,L)
        A2(L) = A2S(MAPNUMBER,L)
      ENDDO
C
      CALL FINDA(VRTX,Q,NLOOPS,A,QOK)
C
C Loop over Reno points within a group.
C
      DO POINT = 1,GROUPSIZE(GRAPHNUMBER,MAPNUMBER)
C
C Call Hrothgar to tell him that we are starting a new point.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'STARTPOINT')
C
C Get a new point. Check on its badness. If it is too bad, we omit
C the point after notifying Hrothgar. Also, if NEWPOINT runs into
C trouble, it can set BADPOINTQ to .TRUE..
C
      BADPOINTQ = .FALSE.
      XTRAPOINTQ = .FALSE.
      CALL NEWPOINT(Q,A,VRTX,PROP,MAPTYPE,A1,A2,K,ABSK,BADPOINTQ)
      IF (BADPOINTQ) THEN
        CALL HROTHGAR(1,KCUT0,1.0D0,1,'BADPOINT  ')
      ELSE
        CALL CHECKPOINT(K,ABSK,PROP,BADNESS)
        IF (BADNESS.GT.100*BADNESSLIMIT) THEN
          CALL HROTHGAR(1,KCUT0,1.0D0,1,'BADPOINT  ')
          BADPOINTQ = .TRUE.
        ELSEIF (BADNESS.GT.BADNESSLIMIT) THEN
          CALL HROTHGAR(1,KCUT0,1.0D0,1,'XTRAPOINT ')
          XTRAPOINTQ = .TRUE.
        ENDIF
      ENDIF
C
C If the point is not too bad, we can call CALCULATE.
C The final state momenta found, KCUT, along with the corresponding
C weights, are reported to Hrothgar by CACULATE.
C Then call Hrothgar to tell him that we are done with this point.
C
      IF (.NOT.BADPOINTQ) THEN
        CALL CALCULATE(VRTX,SELFPROP,GRAPHNUMBER,K,ABSK,
     >               QS,A1S,A2S,MAPTYPES,NMAPS,VALUE,MAXPART,VALUECHK)
      ENDIF
C
C Add contribution from this point to integral.
C We count the point if Maxvalue/|Value| < Cancellimit.
C
      IF (.NOT.BADPOINTQ) THEN
        IF ( MAXPART.GT. 100*CANCELLIMIT*ABS(XXREAL(VALUE)) ) THEN
          CALL HROTHGAR(1,KCUT0,1.0D0,1,'BADPOINT  ')
          BADPOINTQ = .TRUE.
        ELSEIF ( MAXPART.GT. CANCELLIMIT*ABS(XXREAL(VALUE)) ) THEN
          CALL HROTHGAR(1,KCUT0,1.0D0,1,'XTRAPOINT ')
          XTRAPOINTQ = .TRUE.
        ENDIF
      ENDIF
C      
      IF ( (.NOT.BADPOINTQ).AND.(.NOT.XTRAPOINTQ) ) THEN
        INTEGRAL = INTEGRAL + VALUE
        FLUCT(GRAPHNUMBER,MAPNUMBER) = FLUCT(GRAPHNUMBER,MAPNUMBER)
     >    + XXREAL(VALUE)**2/GROUPSIZE(GRAPHNUMBER,MAPNUMBER)
        INTEGRALCHK = INTEGRALCHK + VALUECHK
        INCLUDED = INCLUDED + 1
C
C For diagnostic purposes, we need VALPT, the contribution to
C the integral being calculated from this point, normalized such
C that the integral is the sum over all points chosen of VALPT
C divided by the total number of points, NRENO * GROUPSIZETOTAL.
C
        VALPT = ABS(XXREAL(VALUE))*GROUPSIZETOTAL
        LOGVALPT = LOG10(VALPT)
        DO L = -9,6
          IF((LOGVALPT.GE.L).AND.(LOGVALPT.LT.(L+1))) THEN
            NVALPT(L) = NVALPT(L) + 1
          ENDIF
        ENDDO
        IF (VALPT.GT.VALPTMAX) THEN
          VALPTMAX = VALPT
          DO P = 1,NPROPS
          DO MU = 1,3
            KBAD(P,MU) = K(P,MU)
          ENDDO
          ENDDO
          BADGRAPHNUMBER = GRAPHNUMBER
          BADMAPNUMBER = MAPNUMBER
        ENDIF
      ELSE IF ((.NOT.BADPOINTQ).AND.(XTRAPOINTQ) ) THEN
C
C For points that are 'extra', we include the value of
C the integrand in the INTEGRALBIS, which will provide an estimate
C or the effect of the cutoffs.
C
        INTEGRALBIS = INTEGRALBIS + XXREAL(VALUE)
        EXTRACOUNT = EXTRACOUNT + 1
C
      ELSE
        OMITTED = OMITTED + 1
      ENDIF
C
C End of loop over POINT.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'POINTDONE ')
      ENDDO
C
C End of loop over MAPNUMBER.
C
      ENDDO
C
C End for IF (USEGRAPH(GRAPHNUMBER)) THEN
C
      ENDIF
C
C End of loop DO WHILE (GRAPHFOUND)/ IF (GRAPHFOUND).
C
      ENDIF
      ENDDO
C
C Call Hrothgar to tell him that we are done with this group.
C
      CALL HROTHGAR(1,KCUT0,1.0D0,1,'GROUPDONE ')
C
C Add results from this group to the SUM variables. 
C
      SUMR = SUMR + XXREAL(INTEGRAL)
      SUMI = SUMI + XXIMAG(INTEGRAL)
      SUMBIS = SUMBIS + INTEGRALBIS
      SUMCHKR = SUMCHKR + XXREAL(INTEGRALCHK)
      SUMCHKI = SUMCHKI + XXIMAG(INTEGRALCHK)
C
      SQRSUMR = SQRSUMR + XXREAL(INTEGRAL)**2
      SQRSUMI = SQRSUMI + XXIMAG(INTEGRAL)**2
      SQRSUMBIS = SQRSUMBIS + INTEGRALBIS**2
      SQRSUMCHKR = SQRSUMCHKR + XXREAL(INTEGRALCHK)**2
      SQRSUMCHKI = SQRSUMCHKI + XXIMAG(INTEGRALCHK)**2
C
C Reset the INTEGRAL variables for the next group.
C
      INTEGRAL = (0.0D0,0.0D0)
      INTEGRALBIS = 0.0D0
      INTEGRALCHK = (0.0D0,0.0D0)
C
C End of loop DO WHILE (CPUTIME.LT.TIMELIMIT)
C
      CALL TIMING(DELTATIME)
      CPUTIME = CPUTIME + DELTATIME
      ENDDO
C
C Calculate the SUM results.
C
      SUMR = SUMR/NRENO
      SUMI = SUMI/NRENO
      SUMBIS = SUMBIS/NRENO
      SUMCHKR = SUMCHKR/NRENO
      SUMCHKI = SUMCHKI/NRENO
C
      SQRSUMR = SQRSUMR/NRENO
      SQRSUMI = SQRSUMI/NRENO
      SQRSUMBIS = SQRSUMBIS/NRENO
      SQRSUMCHKR = SQRSUMCHKR/NRENO
      SQRSUMCHKI = SQRSUMCHKI/NRENO
C
      ERRORR = SQRT((SQRSUMR - SUMR**2)/(NRENO - 1))
      ERRORI = SQRT((SQRSUMI - SUMI**2)/(NRENO - 1))
      ERRORBIS = SQRT((SQRSUMBIS - SUMBIS**2)/(NRENO - 1))
      ERRORCHKR = SQRT((SQRSUMCHKR - SUMCHKR**2)/(NRENO-1))
      ERRORCHKI = SQRT((SQRSUMCHKI - SUMCHKI**2)/(NRENO-1))
C
      DO GRAPHNUMBER = 1,NUMBEROFGRAPHS
      DO MAPNUMBER = 1,NUMBEROFMAPS(GRAPHNUMBER)
        FLUCT(GRAPHNUMBER,MAPNUMBER) =
     >         FLUCT(GRAPHNUMBER,MAPNUMBER)/NRENO
      ENDDO
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C                 Subroutines associated with NEWGRAPH                 C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C Out:
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
      LOGICAL SELFPROP(3*SIZE-1)
      LOGICAL GRAPHFOUND
C
C 8 November 1992 Home fixup of bugs.
C 28 November 1992 Add check that we get each graph only once.
C 13 July 1994
C 13 April 1996
C  1 January 1998 Add output variable SELFPROP. Omit NPERMS as output.
C----------
C Varibles:
C VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of
C             of propagator P. Specifies the supergraph for output.
C PROP(V,I) = Index of Ith propagator attached to vertex V, I = 1,2,3.
C              Also specifies the supergraph for output.
C SELFPROP(P) = True if propagator P is part of a one loop self-energy
C               subgraph or attaches to a such a subgraph.
C C(V,I) = Index of Ith vertex connected to vertex V.
C        V = 1,...,NVERTS; I =1,2,3; C(V,I) = 1,...,NVERTS and -1,-2.
C        Here C(V,1).LE.C(V,2).LE.C(V,3).
C        This is the fundamental specification of the supergraph.
C N = Number of permutations of the vertices that give same graph.
C GRAPHFOUND = True when the subroutine finds a new graph.
C COUNT(V) = Number of vertices connected to vertex V.
C Vertex 1 is automatically connected to the photon "-1":C(1,1) = -1.
C Vertex 2 is automatically connected to the photon "-2":C(2,1) = -2.
C The freedom to renumber the vertices 3,...,NVERTS is used to choose
C a standard numbering:
C  We choose the numbering with the smallest value of C(1,1);
C  For numberings with equal values of C(1,1) we choose the numbering
C    with the smallest value of C(1,2);
C  For numberings with equal values of C(1,2) we choose the numbering
C    with the smallest value of C(1,3);
C  For numberings with equal values of C(1,3) we choose the numbering
C    with the smallest value of C(2,1); et cetera.
C
C The connections are generated starting with vertex 1.  We make
C a choice of connections for vertex V, then move on to make a choice
C for connections to vertex V + 1.  When we are out of choices for
C connections to vertex V, we step back and try the next choice for
C vertex V - 1.
C
C Connections to the external boson:
C In C(V,I) we assign the first connection of vertex 1 to be vertex "-1"
C while the first connection of vertex 2 is vertex "-2." This numbering
C is convenient for working out C(V,I). In reporting the results, 
C however, we label the external boson with propagator 0, so that
C PROP(1,1) = PROP(2,1) = 0. Then propagator 0 attaches to vertices
C 1 and 2: VERT(0,1) = 2, VERT(0,2) = 1.
C----------
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER C(2*SIZE,3),COUNT(2*SIZE)
      INTEGER NUSED(2*SIZE),VA,VB
      INTEGER V,VV,I,P,NPERMS
      LOGICAL ONEPI,OK
      LOGICAL FAIL,NEWSTART,UP
      DATA NEWSTART/.TRUE./
      SAVE
C
C Initializations.
C
      IF (NEWSTART) THEN
        DO  VV = 1,NVERTS
          COUNT(VV) = 0
          DO I = 1,3
            C(VV,I) = 0
          ENDDO
        ENDDO
        C(1,1) = -1
        COUNT(1) = 1
        C(2,1) = -2
        COUNT(2) = 1
        V = 1
        UP = .TRUE.
      ENDIF
C
C Move from level to level in tree structure of choices. When UP
C is true, we have moved to a higher V; when UP is false, we have
C moved to a smaller V.
C
      DO WHILE (.TRUE.)
C
      IF (UP) THEN
        CALL NEWCHOICE(C,COUNT,V,FAIL)
      ELSE
        CALL NEXTCHOICE(C,COUNT,V,FAIL)
      ENDIF
      IF (FAIL) THEN
C
C If we couldn't find connectections for vertex V, then we should
C step back and look for the next connections for vertex V-1.  But if
C V is currently 1, then we can't step back, so we have found all
C the graphs.
C
        IF (V.GT.1) THEN
          V = V - 1
          UP = .FALSE.
        ELSE
          NEWSTART = .TRUE.
          GRAPHFOUND = .FALSE.
          DO P = 0,NPROPS
           DO I = 1,2
            VRTX(P,I) = 0
           ENDDO
          ENDDO
          RETURN
        ENDIF
C
C If we did find connections for vertex V, then we should step onward
C and look for new connections for vertex V+1.  But if V is currently
C equal to NVERTS, then we must have found a graph.  We check for
C validity. If it is valid, we exit with the results, setting V and UP
C so that the next time the subroutine is called we will start looking
C for the next connections for vertex V-1.  If our graph is not valid
C (eg. one particle reducible) then we step back to look for new
C connections for vertex V-1 right away.
C
      ELSE
        IF (V.LT.NVERTS) THEN
          V = V + 1
          UP = .TRUE.
        ELSE
          V = V - 1
          UP = .FALSE.
          IF (ONEPI(C)) THEN
          CALL CHECK(C,NPERMS,OK)
            IF (OK) THEN
              NEWSTART = .FALSE.
              GRAPHFOUND = .TRUE.
C
C Exit. We translate the results for C(V,I) into VRTX(P,I), I = 1,2,
C and PROP(P,I), I = 1,2,3. Here NUSED(V) denotes how many propagators
C we have so far assigned connecting to vertex V.
C
              DO VV = 1,NVERTS
                NUSED(VV) = 0
              ENDDO
              VRTX(0,1) = 2
              VRTX(0,2) = 1
              PROP(1,1) = 0
              NUSED(1)  = 1
              PROP(2,1) = 0
              NUSED(2)  = 1
              P = 1
              DO VV = 1,NVERTS
              DO I = 1,3
                IF (C(VV,I).GT.VV) THEN
                  VA = VV
                  VB = C(VV,I)
                  VRTX(P,1) = VA
                  NUSED(VA) = NUSED(VA) + 1
                  PROP(VA,NUSED(VA)) = P
                  VRTX(P,2) = VB
                  NUSED(VB) = NUSED(VB) + 1
                  PROP(VB,NUSED(VB)) = P
                  P = P+1
                ENDIF
              ENDDO
              ENDDO
              IF (P.NE.NPROPS+1) THEN
                WRITE(NOUT,*)'SNAFU in NEWGRAPH',P-1,NPROPS
                STOP
              ENDIF
              DO VV = 1,NVERTS
                IF (NUSED(VV).NE.3) THEN
                  WRITE(NOUT,*)'Problem in NEWWGRAPH',VV,NUSED(VV)
                  STOP
                ENDIF
              ENDDO
C
C We also need to report which propagators P connect to a vertex
C that is part of a one loop self-energy subgraph.
C
              DO P = 1,NPROPS
                SELFPROP(P) = .FALSE.
              ENDDO
              DO VV = 1,NVERTS
                IF ((C(VV,1).EQ.C(VV,2)).OR.(C(VV,2).EQ.C(VV,3))) THEN
                  DO I = 1,3
                    SELFPROP(PROP(VV,I)) = .TRUE.
                  ENDDO
                ENDIF
              ENDDO
C
C OK. We are ready to return.
C
              RETURN
C
            ENDIF
          ENDIF
        ENDIF
      ENDIF
C
C End main loop "DO WHILE (.TRUE.)"
C
      ENDDO
C
      END
C
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEWCHOICE(C,COUNT,V,FAIL)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C
      INTEGER C(2*SIZE,3),COUNT(2*SIZE)
      INTEGER V
      LOGICAL FAIL
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER VV,K
      LOGICAL FOUND
      SAVE
C
C If COUNT(V) = 3, then we don't need any more connections to this
C vertex. If COUNT(V) = 0, then we appear to be starting to make a
C vacuum graph after having completed a graph with too few loops, so
C we should just quit.
C
      IF (COUNT(V).EQ.3) THEN
          FAIL = .FALSE.
          RETURN
      ELSE IF (COUNT(V).EQ.0) THEN
          FAIL = .TRUE.
          RETURN
      ENDIF
C
C Generate starting choice for new vertices to connect to V. We connect
C to the vertices with the smallest possible indices.
C
      FAIL = .FALSE.
      VV = V + 1
      DO K = (COUNT(V) + 1),3
        FOUND = .FALSE.
        DO WHILE (.NOT.FOUND)
          IF (VV.GT.NVERTS) THEN
              FAIL = .TRUE.
              RETURN
          ENDIF
          IF ( COUNT(VV).LT.3) THEN
            COUNT(V) = COUNT(V) + 1
            C(V,K) = VV
            COUNT(VV) = COUNT(VV) + 1
            C(VV,COUNT(VV)) = V
            FOUND = .TRUE.
          ELSE
            VV = VV + 1
          ENDIF
        ENDDO
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEXTCHOICE(C,COUNT,V,FAIL)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C
      INTEGER C(2*SIZE,3),COUNT(2*SIZE)
      INTEGER V
      LOGICAL FAIL
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER VV,VVV,V2,V3,I
      LOGICAL FOUND
      SAVE
C
C First, erase any connections among higher index vertices.
C
      DO VV = V+1,NVERTS
      DO VVV = V+1,NVERTS
      DO I = 1,3
        IF (C(VV,I).EQ.VVV) THEN
          C(VV,I) = 0
          COUNT(VV) = COUNT(VV) - 1
        ENDIF
      ENDDO
      ENDDO
      ENDDO
C
C Next, get the next connection set for vertex V.
C First, we try to find a new third connection for V.
C
      V3 = C(V,3)
C If third connection was to a lower index vertex, we can't change it.
      IF (V3.LE.V) THEN
        FAIL = .TRUE.
        RETURN
      ENDIF
C Erase third connection:
      C(V,3) = 0
      C(V3,COUNT(V3)) = 0
      COUNT(V) = COUNT(V) - 1
      COUNT(V3) = COUNT(V3) - 1
C Look for a new one:
      DO WHILE (V3.LT.NVERTS)
        V3 = V3 + 1
        IF ((COUNT(V3-1).GT.0).AND.(COUNT(V3).LT.3)) THEN
          COUNT(V) = COUNT(V) + 1
          COUNT(V3) = COUNT(V3) + 1
          C(V,3) = V3
          C(V3,COUNT(V3)) = V
          FAIL = .FALSE.
          RETURN
        ENDIF
      ENDDO
C
C We have failed to find a new third connection for V, so
C try for a second connection.  
C
      V2 = C(V,2)
C If second connection was to a lower index vertex, we can't change it.
      IF (V2.LE.V) THEN
        FAIL = .TRUE.
        RETURN
      ENDIF
C Erase second connection:
      C(V,2) = 0
      C(V2,COUNT(V2)) = 0
      COUNT(V) = COUNT(V) - 1
      COUNT(V2) = COUNT(V2) - 1
C Look for a new one:
      DO WHILE (V2.LT.NVERTS)
        V2 = V2 + 1
        IF ((COUNT(V2-1).GT.0).AND.(COUNT(V2).LT.3)) THEN
          COUNT(V) = COUNT(V) + 1
          COUNT(V2) = COUNT(V2) + 1
          C(V,2) = V2
          C(V2,COUNT(V2)) = V
C We found a new second connection.  Now get a third connection.
C---
          V3 = V2
          FOUND = .FALSE.
          DO WHILE (.NOT.FOUND)
            IF ( COUNT(V3).LT.3) THEN
              COUNT(V) = COUNT(V) + 1
              COUNT(V3) = COUNT(V3) + 1
              C(V,3) = V3
              C(V3,COUNT(V3)) = V
              FOUND = .TRUE.
            ELSE
              V3 = V3 + 1
              IF (V3.GT.NVERTS) THEN
                FAIL = .TRUE.
                RETURN
              ENDIF
            ENDIF
          ENDDO
C---
C We have found a good third connection also, so we are done!
          FAIL = .FALSE.
          RETURN
        ENDIF
      ENDDO
C We couldn't find a second connection
C
      FAIL = .TRUE.
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      LOGICAL FUNCTION ONEPI(CIN)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER CIN(2*SIZE,3)
C
C Checks that the graph is connected and 1 particle irreducible.
C Modified 26 July 1994.
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      LOGICAL LEFT(2*SIZE),CHANGE
      INTEGER C(2*SIZE,3)
      INTEGER V,I,V1,V2,I1,I2
      SAVE
C
C Initialize
C
      ONEPI = .TRUE.
C
      DO V = 1,NVERTS
      DO I = 1,3
        C(V,I) = CIN(V,I)
      ENDDO
      ENDDO
C
C Set up loops to successively erase each propagator.
C
      DO V1 = 1,NVERTS
      DO I1 = 1,3
      V2 = C(V1,I1)
      IF (V2.GT.V1) THEN
      DO I2 = 1,3
      IF (C(V2,I2).EQ.V1) THEN
        C(V1,I1) = 0
        C(V2,I2) = 0
C--We have now erased the propagator from V1 to V2.  Let's see if
C  the remaining graph is connected.
        DO V = 1,NVERTS
         LEFT(V) = .FALSE.
        ENDDO
C Construct Left set.
        LEFT(1) = .TRUE.
        CHANGE = .TRUE.
        DO WHILE (CHANGE)
          CHANGE = .FALSE.
          DO V = 1,NVERTS
          DO I = 1,3
            IF ( (1.LE.C(V,I)).AND.(C(V,I).LE.NVERTS) ) THEN
             IF ( LEFT(V) .AND. (.NOT.LEFT(C(V,I))) ) THEN
              CHANGE = .TRUE.
              LEFT(C(V,I)) = .TRUE.
             ENDIF
            ENDIF
          ENDDO
          ENDDO
        ENDDO
C Check for connectedness
        DO V = 1,NVERTS
         IF ( .NOT.LEFT(V) ) THEN
            ONEPI = .FALSE.
            RETURN
         ENDIF
        ENDDO
C--OK, that remaining graph was OK.  Restore the graph.
        C(V1,I1) = V2
        C(V2,I2) = V1
      ENDIF
      ENDDO
      ENDIF
      ENDDO
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHECK(CIN,NPERMS,OK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER CIN(2*SIZE,3),NPERMS
      LOGICAL OK
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      INTEGER C(2*SIZE,3),V(2*SIZE)
      INTEGER L,I,VV
C
      DO VV = 1,NVERTS
      DO I = 1,3
       C(VV,I) = CIN(VV,I)
      ENDDO
      ENDDO
      L = NVERTS
      NPERMS = 0
      OK = .TRUE.
C
C "CALL PERMUTATIONS(L,C)"
C
C-----
C "SUBROUTINE PERMUTATIONS(L,C)"
C Mock subroutine that generates each element of the permutation
C group S_(L-2), applies it to C, and calls CHECKOUT(C,CIN,N,OK).
C If OK = False is returned, the graph C was no good and we exit
C from CHECK immediately.
C
  1   CONTINUE
      IF (L.EQ.4) THEN
         CALL CHECKOUT(C,CIN,NPERMS,OK)
         IF (.NOT.OK) RETURN
         CALL EXCHANGE(3,4,C)
         CALL CHECKOUT(C,CIN,NPERMS,OK)
         IF (.NOT.OK) RETURN
         CALL EXCHANGE(3,4,C)
C "RETURN"
         GO TO 2
      ENDIF
C "DO V(L) = L,3,-1"
      V(L) = L
  3   CONTINUE
      CALL EXCHANGE(V(L),L,C)
      L = L - 1
C "CALL PERMUTATIONS(L,C)"
      GO TO 1
C Return from mock subroutine comes here:
  2   CONTINUE
      L = L + 1
      CALL EXCHANGE(V(L),L,C)
      V(L) = V(L) - 1
C "ENDDO"
      IF (V(L).GE.3) THEN
        GO TO 3
      ENDIF
C "RETURN" Executed from level L as long as L < NVERTS, else we are done.
      IF (L.LT.NVERTS) THEN
        GO TO 2
      ENDIF
C-----
C Come to here if graph CIN was OK, with NPERMS = number of permutations
C that gave a distinct numbering of the vertices.
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHECKOUT(C,CIN,NPERMS,OK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER C(2*SIZE,3),CIN(2*SIZE,3),NPERMS
      LOGICAL OK
C
C Test if graph C (with vertices permuted) is "less than," or
C "greater than," or equal to the original graph CIN using
C the standard ordering of graphs.  If C > CIN we leave unchanged the
C count NPERMS of how many vertex interchanges give the same graph and
C return OK = True.  If C = CIN, we add one to NPERMS, and we
C still return OK = True.  If C < CIN, then we should not have
C generated CIN, so we return OK = False.
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      INTEGER V,I
C
      DO V = 1,NVERTS
      DO I = 1,3
       IF (C(V,I).LT.CIN(V,I)) THEN
         OK = .FALSE.
         RETURN
       ELSE IF (C(V,I).GT.CIN(V,I)) THEN
         OK = .TRUE.
         RETURN
       ENDIF
      ENDDO
      ENDDO
C
C Come to here if the new graph C is the same as CIN.
C
      OK = .TRUE.
      NPERMS = NPERMS + 1
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE EXCHANGE(V1,V2,C)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
      INTEGER C(2*SIZE,3)
      INTEGER V1,V2
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      INTEGER TEMP1,TEMP2,I,V
      LOGICAL CHANGE
C
      DO I = 1,3
        TEMP1 = C(V1,I)
        TEMP2 = C(V2,I)
        C(V1,I) = TEMP2
        C(V2,I) = TEMP1
      ENDDO
C
      DO V = 1,NVERTS
        CHANGE = .FALSE.
        DO I = 1,3
          IF (C(V,I).EQ.V1) THEN
            C(V,I) = V2
            CHANGE = .TRUE.
          ELSEIF (C(V,I).EQ.V2) THEN
            C(V,I) = V1
            CHANGE = .TRUE.
          ENDIF
        ENDDO
        IF (CHANGE) THEN
C
C Put vertices connected to vertex V in order
C--
          IF (C(V,2).LT.C(V,1)) THEN
            TEMP1 = C(V,1)
            TEMP2 = C(V,2)
            C(V,1) = TEMP2
            C(V,2) = TEMP1
          ENDIF
          IF (C(V,3).LT.C(V,1)) THEN
            TEMP1 = C(V,1)
            TEMP2 = C(V,3)
            C(V,1) = TEMP2
            C(V,3) = TEMP1
          ELSE IF (C(V,3).LT.C(V,2)) THEN
            TEMP1 = C(V,2)
            TEMP2 = C(V,3)
            C(V,2) = TEMP2
            C(V,3) = TEMP1
          ENDIF
C--
        ENDIF
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C            End of subroutines associated with NEWGRAPH               C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE FINDA(VRTX,Q,NQ,A,QOK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2),Q(0:SIZE),NQ
C Out:
      INTEGER A(0:3*SIZE-1,0:SIZE)
      LOGICAL QOK
C
C Finds matrix A relating propagator momenta to loop momenta.
C
C VRTX(P,N)     specifies the graph considered
C Q(L)           specifies the propagators to be considered independent
C NQ             specifies how many entries of Q should be considered
C   NQ = NLOOPS  all the entries in Q should be considered.
C                If Q(0),Q(1),...,Q(NLOOPS) are independent then
C                FINDA generates the matrix A and sets QOK = .TRUE.
C                Otherwise the generation of A fails and QOK = .FALSE.
C   NQ < NLOOPS  only first NQ entries in Q should be considered.
C                If Q(0),Q(1),...,Q(NQ) are independent then
C                FINDA sets QOK = .TRUE.
C                Otherwise  QOK = .FALSE.
C                In either case, a complete A is not generated.
C
C L              index of loop momenta, L = 0,1,...,NLOOPS.
C                L = 0 normally denontes the virtual photon momentum.
C P              index of propagator, P = 0,1,...,NPROPS.
C                P = 0 denotes the virtual photon momentum.
C V              index of vertices, V = 1,...,NVERTS
C A(P,L)         matrix relating propagator momenta to loop momenta.
C                K(P) = Sum_L A(P,L) K(Q(L)).
C VRTX(P,1) = V means that the vertex connected to the tail of
C                propagator P is V.
C VRTX(P,2) = V means that the vertex connected to the head of
C                propagator P is V.
C Q(L) = P       means that we consider the Lth loop momentum to
C                be that carried by propagator P.
C CONNECTED(V,J) = P means that the Jth propagator connected to
C                vertex V is P.
C FIXED(P) = True means that we have determined the momentum carried
C                by propagator P.
C FINISHED(V) = True means that we have determined the momenta carried
C                by all the propagators connected to vertex V.
C PROPSIGN(VRTX,P,V) is a function that returns +1 if the head of
C                propagator P is at V, -1 if the tail is at V.
C COUNT          is the number of propagators connected to the vertex
C                under consideration such that FIXED(P) = True. If
C                COUNT = 2, then we can fix another propagator momentum.
C 3 July 1994
C 19 December 1995
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER L,P,V,J,L1,L2
      INTEGER CONNECTED(2*SIZE,3)
      LOGICAL FIXED(0:3*SIZE-1),FINISHED(2*SIZE)
      LOGICAL CHANGE
      INTEGER PROPSIGN,SIGN
      INTEGER SUM(0:SIZE)
      INTEGER COUNT
      INTEGER PTOFIX
C
      IF((NQ.LT.1).OR.(NQ.GT.NLOOPS)) THEN
        WRITE(NOUT,*)'NQ out of range in FINDA'
      ENDIF
C
C First check to see that the same propagator hasn't been
C assigned to two loop variables.
C
      DO L1 = 0,NQ-1
       DO L2 = L1+1,NQ
        IF (Q(L1).EQ.Q(L2)) THEN
          QOK = .FALSE.
          RETURN
        ENDIF
       ENDDO
      ENDDO
C
C Initialization.
C
      QOK = .FALSE.
C
      DO V = 1,NVERTS
        J = 1
        DO P = 0,NPROPS
         IF( (VRTX(P,1).EQ.V).OR.(VRTX(P,2).EQ.V) ) THEN
          CONNECTED(V,J) = P
          J = J+1
         ENDIF
        ENDDO
      ENDDO
C
      DO P = 0,NPROPS
        DO L = 0,NLOOPS
          A(P,L) = 0
        ENDDO
      ENDDO
      DO L = 0,NQ
        A(Q(L),L) = 1
      ENDDO
C
      DO P = 0,NPROPS
        FIXED(P) = .FALSE.
      ENDDO
      DO L = 0,NQ
        FIXED(Q(L)) = .TRUE.
      ENDDO
C
      DO V = 1,NVERTS
        FINISHED(V) = .FALSE.
      ENDDO
C
      CHANGE = .TRUE.
C
C Start.
C
      DO WHILE (CHANGE)
      CHANGE = .FALSE.
C
      DO V = 1,NVERTS
      IF (.NOT.FINISHED(V)) THEN
        COUNT = 0
        DO J = 1,3
         P = CONNECTED(V,J)
         IF ( FIXED(P) ) THEN
           COUNT = COUNT + 1
         ENDIF
        ENDDO
C
C There are 3 already fixed propagators conencted to this vertex, so
C we must check to see if the momenta coming into the vertex sum to
C zero.
C
        IF (COUNT.EQ.3) THEN
          DO L = 0,NQ
            SUM(L) = 0
          ENDDO
          DO J = 1,3
            P = CONNECTED(V,J)
            SIGN = PROPSIGN(VRTX,P,V)
            DO L = 0,NQ
              SUM(L) = SUM(L) + SIGN * A(P,L)
            ENDDO
          ENDDO
          DO L = 0,NQ
C
C Dependent propagators given to FINDA.
C
            IF (SUM(L).NE.0) THEN
              QOK = .FALSE.
              RETURN
C
            ENDIF
          ENDDO
          FINISHED(V) = .TRUE.
          CHANGE = .TRUE.
C
C There are two already fixed propagators connected to this vertex,
C so we should determine the momentum carried by the remaining,
C unfixed, propagator.
C
        ELSEIF (COUNT.EQ.2) THEN
          DO L = 0,NQ
            SUM(L) = 0
          ENDDO
          DO J = 1,3
            P = CONNECTED(V,J)
            IF ( FIXED(P) ) THEN
              SIGN = PROPSIGN(VRTX,P,V)
              DO L = 0,NQ
                SUM(L) = SUM(L) + SIGN * A(P,L)
              ENDDO
            ELSE
              PTOFIX = P
            ENDIF
          ENDDO
          SIGN = PROPSIGN(VRTX,PTOFIX,V)
          DO L = 0,NQ
            A(PTOFIX,L) =  - SIGN * SUM(L)
          ENDDO
          FIXED(PTOFIX) = .TRUE.
          FINISHED(V) = .TRUE.
          CHANGE = .TRUE.
        ENDIF
C
C Close loop DO V = 1,NVERTS ; IF (.NOT.FINISHED(V)) THEN.
C
      ENDIF
      ENDDO
C
C Close loop DO WHILE (CHANGE)
C
      ENDDO
C
C At this point, we have not found a contradiction with momentum
C conservation, so the Q's must have been OK:
C
      QOK = .TRUE.
C
C If we had been given a complete set of Q's, then we should have
C fixed each propagator at each vertex. Check just to make sure.
C
      IF (NQ.EQ.NLOOPS) THEN
      DO V = 1,NVERTS
        IF (.NOT.FINISHED(V) ) THEN
          WRITE(NOUT,*)'SNAFU in FINDA'
          write(nout,*)'v = ',v,'  nq =',nq
          write(nout,*)'q =',q(0),q(1),q(2),q(3)
          STOP
        ENDIF
      ENDDO
      ENDIF
C
      RETURN
C
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      INTEGER FUNCTION PROPSIGN(VRTX,P,V)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2)
      INTEGER P,V
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      IF ( VRTX(P,1).EQ.V ) THEN
          PROPSIGN = -1
          RETURN
      ELSEIF ( VRTX(P,2).EQ.V ) THEN
          PROPSIGN = 1
          RETURN
      ELSE
          WRITE(NOUT,*)'PROPSIGN called for P not connected to V.'
          STOP
      ENDIF
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE FINDQS(VRTX,PROP,NMAPS,QS,A1S,A2S,MAPTYPES)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
C Out:
      INTEGER NMAPS,QS(256,0:SIZE)
      INTEGER A1S(256,2:SIZE),A2S(256,2:SIZE)
      CHARACTER*6 MAPTYPES(256)
C
C Finds labels Q of 'special' propagators for map MAPNUMBER.
C The idea is that there should be nested softness of the loop
C momenta, with K(Q(1)) < K(Q(1)) < ... < K(Q(N)). There is also
C a nested collinearity to the direction of K(Q(N)).
C Generates the maptypes,
C IRFULL for infrared map with Q(1) points concentrated in a plane.
C IRPART for infrared map with Q(1) points concentrated in a line.
C MEDIUM for an IR map that with a mild Q(1) soft singularity.
C UVHARD for concentration of points for Q(N) in the UV.
C For maptype IRFULL, gives propagator identifiers P1, P2 and a SIGN
C for use in constructing the map.
C Also reports the total number of maps, NMAPS.
C  1 July 1993
C 25 June 1994
C 10 July 1994
C 19 December 1995
C  4 May 1996
C 29 January 1998
C  2 February 1999
C
C This subroutine explores a tree structure of choices for
C (Q(1),...,Q(NLOOPS)). At the base level, there are choices for Q(1).
C Then for each choice of Q(1), there are choices for Q(2), etc.
C At each level, the choice must pass certain tests.
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER L,LL,J,P,P0
      INTEGER Q(0:SIZE)
      INTEGER NUSED(SIZE),QUSED(SIZE,3*SIZE-1)
      LOGICAL DEPENDENT
C
      INTEGER Q1,V1,V2
C
C FINDA variables
C
      INTEGER A(0:3*SIZE-1,0:SIZE)
      LOGICAL QOK
C
C SOFTHELP variables
C
      INTEGER P1,P2,SIGN1,SIGN2
      LOGICAL YES
C
      LOGICAL MORENEEDED
C
C NEWCUT variables:
C
      INTEGER CUTINDEX(SIZE+1),CUTSIGN(SIZE+1),NCUT
      INTEGER ISIGN(3*SIZE-1)
      LOGICAL LEFTLOOP,RIGHTLOOP
      INTEGER LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1),NINLOOP
      LOGICAL NEWCUTINIT,CUTFOUND
C
      DO L = 1,NLOOPS
        NUSED(L) = 0
      ENDDO
C
      NMAPS = 0
      L = 1
      P0 = 1
      Q(0) = 0
C
C Enter this loop searching for Q(L) for loopindex L, starting at
C propagator P0. We try successive choices Q(L) = P.
C
1     CONTINUE
      DO P = P0,NPROPS
        Q(L) = P
C
C Test 1. (Q(1),...,Q(L)) must be independent. If this fails, GOTO 2
C to try another P.
C
        IF (DEPENDENT(VRTX,Q,L)) THEN
          GOTO 2
        ENDIF
C
C Test 2. (Q(1),...,Q(L)) must be inequivalent to previous choices
C (Q(1),...,Q(L)'). That is, (Q(1),...,Q(L),Q(L)') must be independent.
C There have been NUSED(L) previous choices for Q(L), recorded in
C QUSED(L,J) for J = 1,...,NUSED(L). If this test fails, GOTO 2
C to try another P.
C
        DO J = 1,NUSED(L)
          Q(L+1) = QUSED(L,J)
          IF (DEPENDENT(VRTX,Q,L+1)) THEN
            GOTO 2
          ENDIF
        ENDDO
C
C Success. (Q(1),...,Q(L)) is good.
C If L = NLOOPS, we have a complete choice.
C If L < NLOOPS, move to higher loopindex and look for Q(L+1).
C
        IF (L.EQ.NLOOPS) THEN
C
C We have a complete choice. Record it.
C
          NMAPS = NMAPS + 1
          IF (NMAPS.GT.256) THEN
            WRITE(NOUT,*)'NMAPS.GT.256'
            STOP
          ENDIF
          DO LL = 0,NLOOPS
            QS(NMAPS,LL) = Q(LL)
          ENDDO
C
C We also need to record the correct MAPTYPE. It is a soft map.
C But what kind of a soft map do we want? The choices are 
C IRFULL for infrared map with Q(1) points concentrated in a plane.
C IRPART for infrared map with Q(1) points concentrated in a line.
C MEDIUM for an IR map that with a mild Q(1) soft singularity.
C UVHARD for concentration of points for Q(N) in the UV.
C 1) The MEDIUM case occurs when the propagator Q(1) connects
C    to vertex 1 or vertex 2. Then we need only a slight concentration
C    of points in the soft region.
C 2) In some cases we want an IRFULL map. This occurs if there is a 
C    cut for which propagator Q(1) is part of a virtual loop and 
C    connects to cut lines. In this case subrountine SOFTHELP returns 
C    YES = .TRUE. The IRFULL case requires a concentration of points 
C    near theta = 0 in a certain coordinate frame in order to partially 
C    cancel the scattering singularities. In this case, we need to
C    record the propagator indices P1 and P2 and the sign SIGN for
C    use in constructing the map.
C 3) For these cases, we also want an IRPART map, which does not
C    have a special concentration of points in a plane, but does have
C    a special concentration of points in a line. Thus we also create
C    an IRPART entry with the same Q(J).
C 4) In this case subrountine SOFTHELP returns YES = .FALSE., we want
C    only an IRPART entry.
C 5) The UVHARD case is for maps with a in which loop momentum 1
C    gives a hard virtual loop. (The Qs for these maps are generated
C    further on in this subroutine.)
C
          Q1 = Q(1)
          V1 = VRTX(Q1,1)
          V2 = VRTX(Q1,2)
          MAPTYPES(NMAPS) = '      '
          DO L = 2,NLOOPS
            A1S(NMAPS,L) = 0
            A2S(NMAPS,L) = 0
          ENDDO
          IF ((V1.LE.2).OR.(V2.LE.2)) THEN
            MAPTYPES(NMAPS) = 'MEDIUM'
          ELSE
C
C Check SOFTHELP for whether to record this choice just once
C as an IRPART map or twice as an IRPART map and an IRFULL map.
C
            CALL SOFTHELP(Q1,VRTX,PROP,P1,P2,SIGN1,SIGN2,YES)
            IF (YES) THEN
              MAPTYPES(NMAPS) = 'IRPART'
C
C We want to record this choice again as an IRFULL map.
C
              NMAPS = NMAPS + 1
              IF (NMAPS.GT.256) THEN
                WRITE(NOUT,*)'NMAPS.GT.256'
                STOP
              ENDIF
              DO LL = 0,NLOOPS
                QS(NMAPS,LL) = Q(LL)
              ENDDO            
              MAPTYPES(NMAPS) = 'IRFULL'
              CALL FINDA(VRTX,Q,NLOOPS,A,QOK)
              IF(.NOT.QOK) THEN
                WRITE(NOUT,*) 'Problem from FINDA in FINDQS'
                STOP
              ENDIF
              DO L = 2,NLOOPS
                A1S(NMAPS,L) = SIGN1 * A(P1,L)
                A2S(NMAPS,L) = SIGN2 * A(P2,L)
              ENDDO
            ELSE
C
C We want to record this just once as an IRPART map.
C           
              MAPTYPES(NMAPS) = 'IRPART'
            ENDIF
          ENDIF
C
C Then drop back to L = NLOOPS -1 and pick up the search. 
C [Note that we want only one successful choice of Q(NLOOPS) 
C for each possible choice of (Q(1),...,Q(NLOOPS - 1)).]
C
          L = NLOOPS -1
          P0 = Q(L) + 1
          GOTO 1
C         
        ELSE IF (L.LT.NLOOPS) THEN
C
C Move to higher loopindex and look for Q(L+1). First, we record our
C position at this value of L, then we go on to the higher value of L.
C
          NUSED(L) = NUSED(L) + 1
          QUSED(L,NUSED(L)) = Q(L)
          L = L + 1
          P0 = 1
          GOTO 1
        ELSE
          WRITE(NOUT,*)'SNAFU in FINDQS'
          STOP
        ENDIF
C
C Our choice P for Q(L) failed. We must try another.
C
2     CONTINUE
      ENDDO
C
C We have tried all of the choices for this L.
C If L = 1, we are done.
C If L > 1, we drop back to a lower value of L and continue our search.
C
      IF (L.GT.1) THEN
        NUSED(L) = 0
        L = L - 1
        P0 = Q(L) + 1
        GOTO 1      
      ENDIF
C
C Come to here when we have generated all of the soft maps.
C Now we want the UV maps.  We use NEWCUT for this. We want those cuts
C for which NCUT.EQ.(CUTMAX-1) , which indicates a virtual loop.
C We run through all the cuts of the present graph, until NEWCUT
C runs out of cuts with NCUT.EQ.(CUTMAX-1) or returns CUTFOUND = .FALSE.
C (and thereby sets its internal variable NEWCUTINIT to .TRUE.).
C
C This code has the "feature" that if a virtual loop can appear both
C on the left and the right of the cut, then two corresponding Q(i)
C sets will be generated.
C
C This code is a simplified version for NLOOPS = 2 or 3 only. But in
C the case NLOOPS = 4 there are in any case a lot of things to change.
C For NLOOPS = 3, we want to allow any of the final state particles
C to be soft (or for all three to be collinear). It is the particle
C with label Q(1) that will be allowed to be soft by NEWPOINT.
C
      MORENEEDED = .TRUE.
      NEWCUTINIT = .TRUE.
      DO WHILE (MORENEEDED)
      CALL NEWCUT(VRTX,NEWCUTINIT,NCUT,ISIGN,
     >            CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP,
     >            NINLOOP,LOOPINDEX,LOOPSIGN,CUTFOUND)
        IF ((CUTFOUND).AND.(NCUT.EQ.(CUTMAX-1))) THEN
        DO L = 1,NLOOPS
          NMAPS = NMAPS + 1
          MAPTYPES(NMAPS) = 'UVHARD'
          DO LL = 2,NLOOPS
            A1S(NMAPS,LL) = 0
            A2S(NMAPS,LL) = 0
          ENDDO
          QS(NMAPS,0) = 0
          QS(NMAPS,NLOOPS) = LOOPINDEX(1)
          QS(NMAPS,1) = CUTINDEX(L)
          DO J = 2,NLOOPS-1
            LL = L + J - 1
            IF (LL.GT.NLOOPS) THEN
              LL = LL - NLOOPS
            ENDIF
            QS(NMAPS,J) = CUTINDEX(LL)
          ENDDO  
        ENDDO
        ELSE
          MORENEEDED = .FALSE.
        ENDIF
      ENDDO
C
C We now have the Q vectors for both the IR maps and the UV maps,
C so we are done.
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      LOGICAL FUNCTION DEPENDENT(VRTX,Q,NQ)
C
C  Checks whether the list of propagators Q(0),Q(1),...,Q(NQ) could
C  be independent loop momenta. If so, DEPENDENT = F. 
C  20 December 1995
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2),Q(0:SIZE),NQ
C
      INTEGER DUMMY(0:3*SIZE-1,0:SIZE)
      LOGICAL QOK
C
      CALL FINDA(VRTX,Q,NQ,DUMMY,QOK)
      IF (QOK) THEN
        DEPENDENT = .FALSE.
      ELSE
        DEPENDENT = .TRUE.
      ENDIF
      RETURN
C
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE SOFTHELP(Q1,VRTX,PROP,P1,P2,SIGN1,SIGN2,YES)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER Q1
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
C Out:
      INTEGER P1,P2,SIGN1,SIGN2
      LOGICAL YES
C
C Given a propagator index Q1, this subroutine examines whether 
C there is a cut graph such that
C  1) Q1 is part of a virtual loop 
C  2) one of the vertices to which Q1 attaches is attached to a cut
C     propagator.
C  3) the other vertex to which Q1 attaches is also attached to a cut
C     propagator.
C If not, YES = .FALSE. is returned. If there is at least one such
C case then the subroutine returns YES = .TRUE. along with P1,P2,SIGN.
C 1) P1 is the index of one of the propagators that attaches to the
C    vertex V1 at the tail of propagator Q1.
C    [Of the two possibilities, other than Q1, P1 is the smaller.]
C 2) P2 is the index of one of the propagators that attaches to the
C    vertex V2 at the head of propagator Q1.
C    [Of the two possibilities, other than Q1, P2 is the smaller.]
C The SIGNs are determined as follows. Suppose that we look at a point
C in loop momentum space such that k_{Q1} = 0. Consider one of the
C cases mentioned above. There are two possibilities.
C 1) Q1 is to the left of the cut:
C     Let pa be the momentum leaving V1 and crossing the cut.
C     Let pb be the momentum leaving V2 and crossing the cut. 
C 2) Q1 is to the right of the cut:
C     Let pa be the momentum crossing the cut and entering V1.
C     Let pb be the momentum crossing the cut and entering V2.
C Now define SIGN1 by k_{P1} = SIGN1 * pa 
C and define SIGN2 by k_{P2} = SIGN2 * pb.
C I think that one gets the same value of SIGNs for all applicable cuts,
C but I am not sure, so the program simply checks.
C
C 21 April 1996
C  2 May 1996
C
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      REAL*8 V1,V2,P1B,P2B
      LOGICAL MORENEEDED
C NEWCUT input:
      LOGICAL NEWCUTINIT
C NEWCUT output:
      INTEGER CUTINDEX(SIZE+1),CUTSIGN(SIZE+1),NCUT
      INTEGER ISIGN(3*SIZE-1)
      LOGICAL LEFTLOOP,RIGHTLOOP
      INTEGER LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1),NINLOOP
      LOGICAL CUTFOUND
C
      INTEGER P,I
      LOGICAL CUTDIRECTION(0:3*SIZE-1),LOOP(0:3*SIZE-1)
      INTEGER TRIAL1,TRIAL2
C
C First, we need V1, V2, P1, P2, and the other propagators connecting
C to V1 and V2, which we call P1B and P2B respectively.
C
      V1 = VRTX(Q1,1)
      V2 = VRTX(Q1,2)
      P1  = PROP(V1,1)
      P1B = PROP(V1,2)
      IF (P1.EQ.Q1) THEN
        P1  = PROP(V1,2)
        P1B = PROP(V1,3)
      ELSEIF (P1B.EQ.Q1) THEN
        P1B = PROP(V1,3)
      ENDIF
      P2  = PROP(V2,1)
      P2B = PROP(V2,2)
      IF (P2.EQ.Q1) THEN
        P2  = PROP(V2,2)
        P2B = PROP(V2,3)
      ELSEIF (P2B.EQ.Q1) THEN
        P2B = PROP(V2,3)
      ENDIF
      SIGN1 = 0
      SIGN2 = 0
      TRIAL1 = 0
      TRIAL2 = 0
C
C Now we cycle through the cuts that leave a virtual loop.
C
      MORENEEDED = .TRUE.
      NEWCUTINIT = .TRUE.
      DO WHILE (MORENEEDED)
      CALL NEWCUT(VRTX,NEWCUTINIT,NCUT,ISIGN,
     >            CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP,
     >            NINLOOP,LOOPINDEX,LOOPSIGN,CUTFOUND)
      IF ((CUTFOUND).AND.(NCUT.EQ.(CUTMAX-1))) THEN
C
C We have a cut. Define a logical variables LOOP(P) that is true
C if propagator P is in the loop. Define an integer variable
C CUTDIRECTION(P) that is 0 if propagator P is not cut and is
C +1 if the propagator is cut and its momentum flows through the
C cut in the positive direction, -1 if the propagator is cut and 
C its momentum flows through the cut in the negative direction.
C
        DO P = 0,NPROPS
          CUTDIRECTION(P) = 0
          LOOP(P) = .FALSE.
        ENDDO
        DO I = 1,NCUT
          CUTDIRECTION(CUTINDEX(I)) = CUTSIGN(I)
        ENDDO
        DO I = 1,NINLOOP
          LOOP(LOOPINDEX(I)) = .TRUE.
        ENDDO
C
C We look further only if Q1 is in the loop.
C
        IF (LOOP(Q1)) THEN
C
C First sign.
C
          IF (CUTDIRECTION(P1).NE.0) THEN
            TRIAL1 = CUTDIRECTION(P1)
          ELSEIF (CUTDIRECTION(P1B).NE.0) THEN
            IF     ( (VRTX(P1,1).EQ.V1).AND.(VRTX(P1B,1).EQ.V1) ) THEN
              TRIAL1 = - CUTDIRECTION(P1B)
            ELSEIF ( (VRTX(P1,2).EQ.V1).AND.(VRTX(P1B,2).EQ.V1) ) THEN
              TRIAL1 = - CUTDIRECTION(P1B)
            ELSEIF ( (VRTX(P1,1).EQ.V1).AND.(VRTX(P1B,2).EQ.V1) ) THEN
              TRIAL1 = CUTDIRECTION(P1B)
            ELSEIF ( (VRTX(P1,2).EQ.V1).AND.(VRTX(P1B,1).EQ.V1) ) THEN
              TRIAL1 = CUTDIRECTION(P1B)
            ELSE
              WRITE(NOUT,*) 'Failure in SOFTHELP'
              STOP
            ENDIF
          ELSE
            TRIAL1 = 0
          ENDIF
C
C Second sign.
C
          IF (CUTDIRECTION(P2).NE.0) THEN
            TRIAL2 = CUTDIRECTION(P2)
          ELSEIF (CUTDIRECTION(P2B).NE.0) THEN
            IF     ( (VRTX(P2,1).EQ.V2).AND.(VRTX(P2B,1).EQ.V2) ) THEN
              TRIAL2 = - CUTDIRECTION(P2B)
            ELSEIF ( (VRTX(P2,2).EQ.V2).AND.(VRTX(P2B,2).EQ.V2) ) THEN
              TRIAL2 = - CUTDIRECTION(P2B)
            ELSEIF ( (VRTX(P2,1).EQ.V2).AND.(VRTX(P2B,2).EQ.V2) ) THEN
              TRIAL2 = CUTDIRECTION(P2B)
            ELSEIF ( (VRTX(P2,2).EQ.V2).AND.(VRTX(P2B,1).EQ.V2) ) THEN
              TRIAL2 = CUTDIRECTION(P2B)
            ELSE
              WRITE(NOUT,*) 'Failure in SOFTHELP'
              STOP
            ENDIF
          ELSE
            TRIAL2 = 0
          ENDIF
C
C Net sign. Check whether we have found an applicable case and,
C if so, whether it is consistent with previous cases found.
C
          IF ((TRIAL1*TRIAL2).NE.0) THEN
            IF ((SIGN1*SIGN2).EQ.0) THEN
              SIGN1 = TRIAL1
              SIGN2 = TRIAL2
            ELSE IF (TRIAL1.NE.SIGN1) THEN
                WRITE(NOUT,*)'Ambiguity discovered in SOFTHELP'
                STOP
            ELSE IF (TRIAL2.NE.SIGN2) THEN
                WRITE(NOUT,*)'Ambiguity discovered in SOFTHELP'
                STOP
            ENDIF
          ENDIF
C
C End of IF (LOOP(Q1)) THEN
C
        ENDIF
C
C End of IF ((CUTFOUND).AND.(NCUT.EQ.(CUTMAX-1))) THEN
C
      ELSE
        MORENEEDED = .FALSE.
      ENDIF
C
C End of DO WHILE (MORENEEDED)
C
      ENDDO
C
      IF (SIGN1*SIGN2.NE.0) THEN
        YES = .TRUE.
      ELSE
        YES = .FALSE.
      ENDIF
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE NEWPOINT(Q,A,VRTX,PROP,MAPTYPE,A1,A2,K,ABSK,BADPOINTQ)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER Q(0:SIZE)
      INTEGER A(0:3*SIZE-1,0:SIZE)
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
      CHARACTER*6 MAPTYPE
      INTEGER A1(2:SIZE),A2(2:SIZE)
C Out:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
      LOGICAL BADPOINTQ
C
C Chooses a new Monte Carlo point in the space of loop 3-momenta.
C 4 March 1993
C 12 July 1993
C 17 July 1994
C  2 May 1996
C  5 February 1997
C  4 February 1999
C 10 March 1999
C  9 April 1999
C 20 August 1999
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      REAL*8 ENERGYSCALE
      COMMON /MSCALE/ ENERGYSCALE
      REAL*8 AS,BS,AM,BM,AH,BH,AUV,BUV,CTHETA,AC,BC,FC
      COMMON /POINTS/ AS,BS,AM,BM,AH,BH,AUV,BUV,CTHETA,AC,BC,FC
C
      REAL*8 PI
      PARAMETER (PI = 3.141592653589793239D0)
C
      REAL*8 MAGK,EPSILON,GAMMA
      REAL*8 TEMP,TEMP1,TEMP2
      REAL*8 K1(3),K1HAT(3),K1SQ,K1ABS
      REAL*8 K2(3),K2HAT(3),K2SQ,K2ABS
      REAL*8 NORMA
      REAL*8 EA(3),EB(3),EC(3)
      REAL*8 ELL(SIZE,3)
      REAL*8 COSTHETA,SINTHETA,PHI
      REAL*8 N(3),V(3),U(3)
      REAL*8 Z,ZM,ZP
      REAL*8 X,RANDOM
      REAL*8 KSQ,PREVIOUSK,PREVIOUSKT
      INTEGER L,P,MU
      REAL*8 C,EPS
      REAL*8 SINHINV
C
C------------
C
      IF (NLOOPS.LT.1) THEN
        WRITE(NOUT,*) 'NLOOPS less than 1 in NEWPOINT'
        STOP
      ENDIF
C
C Logical structure below is
C IF (MAPTYPE.NE.'UVHARD') THEN
C     IF (MAPTYPE.EQ.'MEDIUM') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRPART') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRFULL') THEN ...
C     ELSE  <error>
C     ENDIF
C  ELSE   <maptype must have been UVHARD> ...
C  ENDIF
C
      IF (MAPTYPE.NE.'UVHARD') THEN
C
C We want a soft map. To start, we use a "not-soft" but "not-UV"
C map for k_{Q(N)}. This map normally has AH = 3, BH = 3.
C We also make a unit vector N(mu) in the direction of k_{Q(N)}.
C
      X = RANDOM(1)
      MAGK = ENERGYSCALE * ( X**(-AH/BH) - 1.0D0 )**(1.0D0/AH)
      X = RANDOM(1)
      COSTHETA = 2.0D0*X - 1.0D0
      SINTHETA = SQRT(1.0D0 - COSTHETA**2)
      X = RANDOM(1)
      PHI = 2.0D0 * PI * X
      N(1) = SINTHETA * COS(PHI)
      ELL(NLOOPS,1) = MAGK * N(1)
      N(2) = SINTHETA * SIN(PHI)
      ELL(NLOOPS,2) = MAGK * N(2)
      N(3) = COSTHETA
      ELL(NLOOPS,3) = MAGK * N(3)
C
C Now we have to complete the map for ELL(L,mu), L < Nloops.
C What kind of map do we want? The choices correspond to the MAPTYPE:
C MEDIUM for an IR map that with a mild soft singularity.
C IRPART for infrared map with Q(1) points concentrated in a line.
C IRFULL for infrared map with Q(1) points concentrated in a plane.
C For maptype IRFULL, we use identifiers SOFTP1,SOFTP2,SOFTSIGN
C in constructing the map.
C
      IF (MAPTYPE.EQ.'MEDIUM') THEN
C
C For the other k_{Q(i)}, we use a nested "mediumsoft" map
C with AM = 2, BM = 2 normally.  
C For NLOOPS = 3, this evaluates ELL(2,mu), ELL(1,mu).
C
        PREVIOUSK = ENERGYSCALE
        DO L = NLOOPS-1,1,-1
          X = RANDOM(1)
          MAGK = PREVIOUSK * ( X**(-AM/BM) - 1.0D0 )**(1.0D0/AM)
          PREVIOUSK = MAGK
          X = RANDOM(1)
          COSTHETA = 2.0D0*X - 1.0D0
          SINTHETA = SQRT(1.0D0 - COSTHETA**2)
          X = RANDOM(1)
          PHI = 2.0D0 * PI * X
          ELL(L,1) = MAGK * SINTHETA * COS(PHI)
          ELL(L,2) = MAGK * SINTHETA * SIN(PHI)
          ELL(L,3) = MAGK * COSTHETA
        ENDDO
C
      ELSE IF (MAPTYPE.EQ.'IRPART') THEN
C
C We generate the other k_{Q(i)}.
C We use a nested "soft" map with AS = 1, BS = 2 normally.
C For NLOOPS = 3, this evaluates ELL(2,mu), ELL(1,MU).
C For these maps, we concentrate the points in the direction of N(mu).
C
        PREVIOUSK = ENERGYSCALE
        PREVIOUSKT = ENERGYSCALE
        DO L = NLOOPS-1,1,-1
C
C First, we generate the magnitude of our vector.
C
          X = RANDOM(1)
          MAGK = PREVIOUSK * ( X**(-AS/BS) - 1.0D0 )**(1.0D0/AS)
          PREVIOUSK = MAGK
C
C Now we want an X that varies from -1 to 1.
C
          X = RANDOM(1)
          X = 2.0D0*X - 1.0D0
C
C Given X, we compute Z = Cos(theta), where theta is the angle
C between ell(mu) and N(mu). Normally, AC = 1, BC = 1/2,
C FC = 1/2.
C
C The parameter EPS, with 0 < EPS < 1 (or EPS = 1), provides
C a favored scale for small angles:
C
          IF (L.EQ.(NLOOPS - 1)) THEN
            EPS = 1.0D0
          ELSE
            EPS = PREVIOUSKT/MAX(MAGK,PREVIOUSKT)
          ENDIF
C
C There are two functions, depending on how big X is.
C Let ZM = 1 - |Z| and ZP = 1 + |Z|. Then Cos(theta) is Z
C and Sin(theta) is Sqrt(ZP*ZM), which is an accurate way to
C do it if ZM << 1. Note that Fortran cleverly provides  
C SIGN(xx,X) to give xx with the sign of X.
C
          IF( ABS(X) .GT. (1.0D0 - EPS)*FC ) THEN
            C = (1.0D0 - FC + EPS*FC)*(1.0D0 + EPS**AC)**(BC/AC)
            TEMP = ((1.0D0 - ABS(X))/C)**(-AC/BC)
            ZM   = EPS**2 * (TEMP - 1.0D0)**(-1.0D0/AC)
            ZP = 2.0D0 - ZM
            COSTHETA = SIGN(1.0D0 - ZM,X)
            SINTHETA = SQRT(ZP*ZM)
          ELSE
            Z = X/FC
            COSTHETA = Z
            SINTHETA = SQRT(1.0D0 - Z**2)
          ENDIF
C
C Generate axes with first axis in the direction of N(mu).
C
          CALL AXES(N,EB,EC)
C
C Generate ell_L(mu) using our unit vectors.
C
          X = RANDOM(1)
          PHI = 2.0D0 * PI * X
          DO MU = 1,3
           ELL(L,MU) = MAGK * ( COSTHETA * N(MU)  + SINTHETA *
     >          ( COS(PHI) * EB(MU) + SIN(PHI) * EC(MU) )  )                       
          ENDDO  
C
C The magnitude of the part transverse to N(mu) is
C
          PREVIOUSKT = MAGK * SINTHETA
C
C Close DO L = NLOOPS-1,1,-1
C
        ENDDO
C
      ELSE IF (MAPTYPE.EQ.'IRFULL') THEN
C
C Here we choose the same map as above for IRPART for K_{Q(2)}, or,
C more generally, for K_{Q(Nloops - 1)},...,K_{Q(2)}. Then, we
C will do something special for K_{Q(1)}.
C -----
C Begin copy of IRPART code with "DO L = NLOOPS-1,1,-1" changed
C to "DO L = NLOOPS-1,2,-1":
C
        PREVIOUSK = ENERGYSCALE
        PREVIOUSKT = ENERGYSCALE
        DO L = NLOOPS-1,2,-1
          X = RANDOM(1)
          MAGK = PREVIOUSK * ( X**(-AS/BS) - 1.0D0 )**(1.0D0/AS)
          PREVIOUSK = MAGK
          X = RANDOM(1)
          X = 2.0D0*X - 1.0D0
          IF (L.EQ.(NLOOPS - 1)) THEN
            EPS = 1.0D0
          ELSE
            EPS = PREVIOUSKT/MAX(MAGK,PREVIOUSKT)
          ENDIF
          IF( ABS(X) .GT. (1.0D0 - EPS)*FC ) THEN
            C = (1.0D0 - FC + EPS*FC)*(1.0D0 + EPS**AC)**(BC/AC)
            TEMP = ((1.0D0 - ABS(X))/C)**(-AC/BC)
            ZM   = EPS**2 * (TEMP - 1.0D0)**(-1.0D0/AC)
            ZP = 2.0D0 - ZM
            COSTHETA = SIGN(1.0D0 - ZM,X)
            SINTHETA = SQRT(ZP*ZM)
          ELSE
            Z = X/FC
            COSTHETA = Z
            SINTHETA = SQRT(1.0D0 - Z**2)
          ENDIF
          CALL AXES(N,EB,EC)
          X = RANDOM(1)
          PHI = 2.0D0 * PI * X
          DO MU = 1,3
           ELL(L,MU) = MAGK * ( COSTHETA * N(MU)  + SINTHETA *
     >          ( COS(PHI) * EB(MU) + SIN(PHI) * EC(MU) )  )                       
          ENDDO  
          PREVIOUSKT = MAGK * SINTHETA
        ENDDO
C
C End copy of IRPART code with "DO L = NLOOPS-1,1,-1" changed
C to "DO L = NLOOPS-1,2,-1".
C -----
C Now we generate K_{Q(1)}. We use the AS, BS map, 
C nested with the previous maps. BUT this
C time we need a special concentration of points near the
C plane theta = 0 in a certain coordinate system.
C
C Step 1:
C We need the appropriate unit vectors for the map.
C
      K1SQ = 0.0D0
      K2SQ = 0.0D0
      DO MU = 1,3
        TEMP1 = 0.0D0
        TEMP2 = 0.0D0
        DO L = 2,NLOOPS
          TEMP1 = TEMP1 + A1(L) * ELL(L,MU)
          TEMP2 = TEMP2 + A2(L) * ELL(L,MU)
        ENDDO
        K1(MU) = TEMP1
        K1SQ = K1SQ + TEMP1**2
        K2(MU) = TEMP2
        K2SQ = K2SQ + TEMP2**2
      ENDDO
      IF ((K1SQ.LT.1.0D-30).OR.(K2SQ.LT.1.0D-30)) THEN
        WRITE(NOUT,*)'K1SQ or K2SQ too small in NEWPOINT'
        STOP
      ENDIF  
      K1ABS = SQRT(K1SQ)
      K2ABS = SQRT(K2SQ)
      DO MU = 1,3
        K1HAT(MU) = K1(MU)/K1ABS
        K2HAT(MU) = K2(MU)/K2ABS
      ENDDO
C
      NORMA = 0.0D0
      DO MU = 1,3
        NORMA = NORMA + (K1HAT(MU) - K2HAT(MU))**2
      ENDDO
      NORMA = SQRT(NORMA)
      IF ((NORMA.LT.1.0D-10)) THEN
        WRITE(NOUT,*)'NORMA too small in NEWPOINT'
        BADPOINTQ = .TRUE.
        RETURN
      ENDIF
      DO MU = 1,3
        EA(MU) = (K1HAT(MU) - K2HAT(MU))/NORMA
      ENDDO
      CALL AXES(EA,EB,EC)
C
C Step 2:
C We want a soft map concentrating points near cos(theta) = 0.
C
      X = RANDOM(1)
      MAGK = PREVIOUSK * ( X**(-AS/BS) - 1.0D0 )**(1.0D0/AS)
      EPSILON = CTHETA * MAGK/PREVIOUSK
      GAMMA = SINHINV(1.0D0/EPSILON)    
      X = RANDOM(1)
      COSTHETA = EPSILON * SINH( (2.0D0*X - 1.0D0)*GAMMA )
      SINTHETA = SQRT(1.0D0 - COSTHETA**2)
      X = RANDOM(1)
      PHI = 2.0D0 * PI * X
C
C Step 3: Put this together using our unit vectors.
C
      DO MU = 1,3
        ELL(1,MU) = MAGK * ( COSTHETA * EA(MU)  + SINTHETA *
     >          ( COS(PHI) * EB(MU) + SIN(PHI) * EC(MU) )  )                       
      ENDDO        
C
C End IF (MAPTYPE.EQ.'MEDIUM') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRPART') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRFULL') THEN ...
C
      ELSE
        WRITE(NOUT,*)'This cannot happen in NEWPOINT'
        STOP
      ENDIF
C
C This completes the IF (MAPTYPE.NE.'UVHARD') THEN ...
C
      ELSE
C
C Now we generate the maps for nearly UV divergent loops.
C
C The AUV = 3, BUV  = 1 (normally) map for loop momentum NLOOPS.
C
        X = RANDOM(1)
        MAGK = ENERGYSCALE*( X**(-AUV/BUV) - 1.0D0 )**(1.0D0/AUV)
        X = RANDOM(1)
        COSTHETA = 2.0D0*X - 1.0D0
        SINTHETA = SQRT(1.0D0 - COSTHETA**2)
        X = RANDOM(1)
        PHI = 2.0D0 * PI * X
        ELL(NLOOPS,1) = MAGK * SINTHETA * COS(PHI)
        ELL(NLOOPS,2) = MAGK * SINTHETA * SIN(PHI)
        ELL(NLOOPS,3) = MAGK * COSTHETA
C
C The AH = 3, BH = 3 (normally) map for loop momentum NLOOPS-1
C
        L = NLOOPS-1
        X = RANDOM(1)
        MAGK = ENERGYSCALE * ( X**(-AH/BH) - 1.0D0 )**(1.0D0/AH)
        X = RANDOM(1)
        COSTHETA = 2.0D0*X - 1.0D0
        SINTHETA = SQRT(1.0D0 - COSTHETA**2)
        X = RANDOM(1)
        PHI = 2.0D0 * PI * X
        N(1) = SINTHETA * COS(PHI)
        ELL(L,1) = MAGK * N(1)
        N(2) = SINTHETA * SIN(PHI)
        ELL(L,2) = MAGK * N(2)
        N(3) = COSTHETA
        ELL(L,3) = MAGK * N(3)
C
C For NLOOPS = 3, we now generate an AS = 1, BS = 2 (normally) 
C map for loop momentum 1. This map includes a concentration
C of soft points and of points in the direction of ELL(NLOOPS-1). 
C (In case NLOOPS = 2, we do nothing. In case NLOOPS > 3, we simply
C repeat the same map.)
C
        DO L = NLOOPS-2,1,-1
C
C First, we generate the magnitude of our vector.
C
          X = RANDOM(1)
          MAGK = ENERGYSCALE * ( X**(-AS/BS) - 1.0D0 )**(1.0D0/AS)
          DO MU = 1,3
            ELL(L,MU) = MAGK * U(MU)
          ENDDO
C
C Now, we need a unit vector V(mu) randomly distributed on the 
C unit sphere.
C
          X = RANDOM(1)
          COSTHETA = 2.0D0*X - 1.0D0
          SINTHETA = SQRT(1.0D0 - COSTHETA**2)
          X = RANDOM(1)
          PHI = 2.0D0 * PI * X
          V(1) = SINTHETA * COS(PHI)
          V(2) = SINTHETA * SIN(PHI)
          V(3) = COSTHETA
C
C Now we want X = V * N. (Note that we have been using X for a 
C random number. This X is essentially a random number with a uniform
C distribution between -1 and 1, so we can use the same name).
C
          X = 0.0D0
          DO MU = 1,3
            X = X + V(MU)*N(MU)
          ENDDO
C
C Here we compute the dot product Z that we want to have between
C our new unit vector U(mu) and N(mu). Normally, AC = 1, BC = 1/2,
C FC = 1/2.
C
          EPS = 1.0D0
          C = (1.0D0 - FC + EPS*FC)*(1.0D0 + EPS**AC)**(BC/AC)
          Z = ((1.0D0 - ABS(X))/C)**(-AC/BC)
          Z = EPS**2 * (Z - 1.0D0)**(-1.0D0/AC)
          Z = 1.0D0 - Z
          Z = SIGN(Z,X)
C
C Now we need a new unit vector U(mu) with U * N = Z.
C
          TEMP = SQRT( (1.0D0 - Z**2)/(1.0D0 - X**2) )
          DO MU = 1,3
            U(MU) = TEMP * ( V(MU) - X * N(MU) ) + Z * N(MU)
          ENDDO
C
C Our vector is MAGK times the unit vector U(mu).
C
          DO MU = 1,3
            ELL(L,MU) = MAGK * U(MU)
          ENDDO
C
C Close DO L = NLOOPS-2,1,-1.
C
        ENDDO
C
C End of IF (MAPTYPE.NE.'UVHARD') THEN ... ELSE ...
C
      ENDIF
C
C Now we have ELL(L,MU) and we need to translate to the propagator
C momenta K(P,MU).
C
      DO P = 1,NPROPS
        KSQ = 0.0D0
        DO MU = 1,3
          TEMP = 0.0D0
          DO L = 1,NLOOPS
            TEMP = TEMP + A(P,L) * ELL(L,MU)
          ENDDO
          K(P,MU) = TEMP
          KSQ = KSQ + TEMP**2
      ENDDO
        ABSK(P) = SQRT(KSQ)
        K(P,0) = 0.0D0
      ENDDO
      DO MU = 0,3
        K(0,MU) = 0.0D0
      ENDDO
      ABSK(0) = 0.0D0
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CHECKPOINT(K,ABSK,PROP,BADNESS)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
      INTEGER PROP(2*SIZE,3)
C Out:
      REAL*8 BADNESS
C
C Calculates the BADNESS of a point chosen by NEWPOINT. If there
C are very collinear particles meeting at a vertex or of there is a
C very soft particle, then the badness is big. Specifically, for
C each vertex V we order the momenta entering the vertex Kmin, Kmid
C Kmax in order of |K|. Then
C
C   Kmin (Kmin + Kmid - Kmax )/Kmax^2
C
C is the 1/badness^2 for that vertex. The badness for the point is the 
C largest of the badness values of all the vertices.
C 
C Revised 13 may 1998
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      REAL*8 SMALLNESSV,SMALLNESS
      INTEGER V
      REAL*8 KMIN,KMID,KMAX,K1,K2,K3
C
      SMALLNESS = 1.0D0
      DO V = 3,NVERTS
       K1 = ABSK(PROP(V,1))
       K2 = ABSK(PROP(V,2))
       K3 = ABSK(PROP(V,3))
       IF (K1.LT.K2) THEN
         KMIN = K1
         KMAX = K2
       ELSE  
         KMIN = K2
         KMAX = K1
       ENDIF
       IF (K3.LT.KMIN) THEN
         KMID = KMIN
         KMIN = K3
       ELSE IF (K3.GT.KMAX) THEN
         KMID = KMAX
         KMAX = K3
       ELSE
         KMID = K3
       ENDIF
       SMALLNESSV = KMIN * (KMIN + KMID - KMAX) /KMAX**2
       IF( SMALLNESSV .LT. SMALLNESS ) THEN
         SMALLNESS = SMALLNESSV
       ENDIF
       ENDDO
       IF (SMALLNESS.LT.1.0D-30) THEN
         BADNESS = 1.0D15
       ELSE
         BADNESS = SQRT(1.0D0/SMALLNESS)
       ENDIF
       RETURN
       END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE AXES(EA,EB,EC)
C
C In:
      REAL*8 EA(3)
C Out:
      REAL*8 EB(3),EC(3)
C
C Given a unit vector E_a(mu), generates unit vectors E_b(mu) and
C E_c(mu) such that E_a*E_b = E_b*E_c = E_c*E_a = 0.
C
C The vector E_b will lie in the plane formed by the z-axis and
C E_a unless E_a itself is nearly aligned along the z-axis, in which
C case E_b will lie in the plane formed by the x-axis and E_a.
C
C 18 April 1996
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C
      REAL*8 COSTHETASQ,SINTHETAINV
C
C For check
C
      INTEGER MU
      REAL*8 DOTAA,DOTBB,DOTCC,DOTAB,DOTAC,DOTBC
C
      COSTHETASQ = EA(3)**2
      IF(COSTHETASQ.LT.0.9D0) THEN
        SINTHETAINV = 1.0D0/SQRT(1.0D0 - COSTHETASQ)
        EC(1) = - EA(2)*SINTHETAINV
        EC(2) =   EA(1)*SINTHETAINV
        EC(3) =   0.0D0
      ELSE
        COSTHETASQ = EA(1)**2
        SINTHETAINV = 1.0D0/SQRT(1.0D0 - COSTHETASQ)
        EC(1) =   0.0D0
        EC(2) = - EA(3)*SINTHETAINV
        EC(3) =   EA(2)*SINTHETAINV
      ENDIF
      EB(1) = EC(2)*EA(3) - EC(3)*EA(2)
      EB(2) = EC(3)*EA(1) - EC(1)*EA(3)
      EB(3) = EC(1)*EA(2) - EC(2)*EA(1) 
C     
C Check:
C
      DOTAA = 0.0D0
      DOTBB = 0.0D0
      DOTCC = 0.0D0
      DOTAB = 0.0D0
      DOTAC = 0.0D0
      DOTBC = 0.0D0
      DO MU = 1,3
        DOTAA = DOTAA + EA(MU)*EA(MU)
        DOTBB = DOTBB + EB(MU)*EB(MU)
        DOTCC = DOTCC + EC(MU)*EC(MU)
        DOTAB = DOTAB + EA(MU)*EB(MU)
        DOTAC = DOTAC + EA(MU)*EC(MU)
        DOTBC = DOTBC + EB(MU)*EC(MU)
      ENDDO
      IF (ABS(DOTAA - 1.0D0).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTAA messed up in AXES'
        STOP
      ELSE IF (ABS(DOTBB - 1.0D0).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTBB messed up in AXES'
        STOP
      ELSE IF (ABS(DOTCC - 1.0D0).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTCC messed up in AXES'
        STOP
      ELSE IF (ABS(DOTAB).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTAB messed up in AXES'
        STOP
      ELSE IF (ABS(DOTAC).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTAC messed up in AXES'
        STOP
      ELSE IF (ABS(DOTBC).GT.1.0D20) THEN
        WRITE(NOUT,*)'DOTBC messed up in AXES'
        STOP
      ENDIF
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C                 Subroutine to calculate integrand                    C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE CALCULATE(VRTX,SELFPROP,GRAPHNUMBER,KIN,ABSKIN,
     >               QS,A1S,A2S,MAPTYPES,NMAPS,VALUE,MAXPART,VALUECHK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2)
      LOGICAL SELFPROP(3*SIZE-1)
      INTEGER GRAPHNUMBER
      REAL*8 KIN(0:3*SIZE-1,0:3),ABSKIN(0:3*SIZE-1)
      INTEGER QS(256,0:SIZE)
      INTEGER A1S(256,2:SIZE),A2S(256,2:SIZE)
      CHARACTER*6 MAPTYPES(256)
      INTEGER NMAPS
C Out:
      COMPLEX*16 VALUE,VALUECHK
      REAL*8 MAXPART
C
C Calculates the value of the graph specified by VRTX at the point K,
C returning result in VALUE, which includes the division by the density
C of points and the jacobian for deforming the contour. Also reports
C MAXPART, the biggest absolute value of the contributions to Re(VALUE).
C This helps us to keep track of cancellations and thus to abort the 
C calculation if too much cancellation among terms will be required.
C
C*********************
C
C Max number of graphs, for array size:
      INTEGER MAXGRAPHS
      PARAMETER (MAXGRAPHS = 10)
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
      REAL*8 MUOVERRTS
      COMMON /RENORMALIZE/ MUOVERRTS
      LOGICAL REPORT,DETAILS
      COMMON /CALCULOOK/ REPORT,DETAILS
      REAL*8 NC,NF
      COMMON /COLORFACTORS/ NC,NF
C How many graphs and how many cuts and maps for each:
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C Labels:
      INTEGER QE(0:SIZE)
C Momenta:
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
      REAL*8 KINLOOP(SIZE+1,0:3)
      REAL*8 KCUT(SIZE+1,0:3)
      COMPLEX*16 NEWKINLOOP(0:3)
      COMPLEX*16 KC(0:3*SIZE-1,0:3)
      COMPLEX*16 ELLSQ,ELL
      REAL*8 E(0:SIZE),RTS
C Renormalization:
      REAL*8 MUMSBAR
C Matrices:
      INTEGER AE(0:3*SIZE-1,0:SIZE)
C FINDA variable:
      LOGICAL QOK
C Whether to reflect across collinear singularities:
      LOGICAL REFLECTQ
      COMMON /WANTREFLECT/ REFLECTQ
C REFLECT variables:
      REAL*8 KREFLECTED(0:3*SIZE-1,0:3),ABSKREFLECTED(0:3*SIZE-1)
      INTEGER NREFLECT,JREFLECT
C DENSITY variables:
      REAL*8 JACNEWPOINT,DENSITY,RHO
C NEWCUT variables:
      INTEGER CUTINDEX(SIZE+1),CUTSIGN(SIZE+1),NCUT
      INTEGER ISIGN(3*SIZE-1)
      LOGICAL LEFTLOOP,RIGHTLOOP
      INTEGER LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1),NINLOOP
C Loopcut variables:
      LOGICAL CALCMORE
      INTEGER JCUT
C DEFORM variables:
      COMPLEX*16 JACDEFORM
C Functions:
      REAL*8 CALS0,SMEAR
      REAL*8 XXREAL,XXIMAG
      COMPLEX*16 COMPLEXSQRT
C Index variables:
      INTEGER P,MU,I,J,CUTNUMBER
C Propagator properties and momenta
      LOGICAL INLOOP(3*SIZE-1),CUT(3*SIZE-1)
      LOGICAL LOOPCUT(3*SIZE-1)
      INTEGER CUTSIGNP(3*SIZE-1)
      INTEGER LOOPLABEL(3*SIZE-1)
      REAL*8 KSQ
      COMPLEX*16 KCSQ
C Results variables:
      REAL*8 CALSVAL
      REAL*8 WEIGHT,MAXWEIGHT
      COMPLEX*16 NUMERATOR,RNUMERATOR,FEYNMAN
      COMPLEX*16 LOOPDENOM
      REAL*8 PLAINDENOM
      REAL*8 PREFACTOR
      COMPLEX*16 INTEGRAND
      COMPLEX*16 CHECK
C Useful constants:
      REAL*8 METRIC(0:3)
      DATA METRIC /+1.0D0,-1.0D0,-1.0D0,-1.0D0/
      REAL*8 PI
      DATA PI /3.1415926535898D0/
C      
C-----------------------------------------------------------------------
C Latest revision: 11 May 1996
C                  24 October 1996 (call to CHECKDEFORM)
C                  15 November 1996 (remove finite 'i epsilon')
C                  18 November 1996 (add CHECKVALUE)
C                  22 November 1996 Bug fixed.
C                  27 November 1996 (complex checkvalue)
C                  29 November 1996 (branchcut check; better checkvalue)
C                  27 February 1997 renormalization; reporting
C                  25 July     1997 renormalization; self-energy graphs
C                  17 September 1997 more renormalization & self-energy
C                  21 September 1997 finish DEFORM
C                  24 September 1997 fix bugs
C                  19 October 1997 fix cutsign bug
C                  22 October 1997 fix renormalization sign bug
C                   6 November 1997 improvements for deformation
C                  28 November 1997 more work on deformation
C                   2 December 1997 more precision in "report" numbers
C                   4 January  1998 revisions for self-energy graphs
C                  11 January 1998 renormalizaion for self-energy graphs
C                  27 February 1998 use Hrothgar for output
C                   5 March 1998 integrate Hrothgar
C                  14 March 1998 restore checks of deformation direction
C                  24 July 1998 use countfactor(graphnumber)
C                   4 August 1998 better CHECKDEFORM
C                   5 August 1998 change to groupsize(graphnumber)
C                  22 August 1998 add color factors
C                  22 December 1998 precalculate cut structure in RENO
C                  26 April 1999 omit REFLECT except as option
C----------------------------------
C
C We do not want to change the value of KIN and ABSKIN, even though
C K and ABSK get changed by the reflection feature of the subroutine.
C
      DO P = 1,NPROPS
        ABSK(P) = ABSKIN(P)
        DO MU = 0,3
          K(P,MU) = KIN(P,MU)
        ENDDO
      ENDDO
C
C Initialize contribution to integral from this point. Also initialize
C BIGGEST, which will be the biggest absolute value of the contributions
C to VALUE. This helps us to keep track of cancellations and thus to
C abort the calculation if too much cancellation among terms will be
C required.
C
      MAXPART = 0.0D0
      VALUE = (0.0D0,0.0D0)
      VALUECHK = (0.0D0,0.0D0)
C
C Do we want a point reflected across the nearest collinear singularity?
C Probalby not, but for special choices of the measurement function
C the collinear cancellation may reduce 1/k_T^2 to just k_T*v/k_T^2,
C where v is a vector that is fixed in the collinear limit. Then
C a collinear reflection will eliminate the net 1/K_T singularity.
C Since we normally don't need this, the default value of REFLECTQ
C is false.
C
      IF (REFLECTQ) THEN
C
C Normally, REFLECT returns NREFLECT = 2. In exceptional
C cases, the reflection transformation applied twice does not bring
C us back to the original point. In this case, REFLECT returns
C NREFLECT = 1 and we do not use a reflected point.
C
      CALL REFLECT(VRTX,ABSK,K,ABSKREFLECTED,KREFLECTED,NREFLECT)
C
      ELSE
C
C In the default case that REFLECTQ is false, we set NREFLECT to 1,
C so that we do not use a reflected point.
C
        NREFLECT = 1
C
      ENDIF
C
C Calculate jacobian.
C
      RHO = DENSITY(GRAPHNUMBER,K,ABSK,QS,A1S,A2S,MAPTYPES,NMAPS)
      IF (NREFLECT.EQ.2) THEN
        RHO = RHO + DENSITY(GRAPHNUMBER,KREFLECTED,ABSKREFLECTED,
     >                      QS,A1S,A2S,MAPTYPES,NMAPS)
      ENDIF
      JACNEWPOINT = 1.0D0/RHO
C
C Loop over both the main point K and the reflected point (if there
C is one).
C
      DO JREFLECT = 1,NREFLECT
C
      IF (JREFLECT.EQ.2) THEN
        DO P = 0,NPROPS
          ABSK(P) = ABSKREFLECTED(P)
          DO MU = 0,3
            K(P,MU) = KREFLECTED(P,MU)
          ENDDO
        ENDDO
      ENDIF
C
C Loop over cuts.
C
      DO CUTNUMBER = 1,NUMBEROFCUTS(GRAPHNUMBER)
      CALL GETCUTINFO(GRAPHNUMBER,CUTNUMBER,NCUT,ISIGN,
     >            CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP,
     >            NINLOOP,LOOPINDEX,LOOPSIGN)
C....
      IF (REPORT) THEN
        WRITE(NOUT,301)NCUT,CUTINDEX(1),CUTINDEX(2),
     >                    CUTINDEX(3),CUTINDEX(4)
301     FORMAT('Ncut =',I2,' CUTINDEX =',4I2)
      ENDIF
C''''
C
C Calculate Sqrt(s) and the renormalization scale MUMSBAR.
C
      RTS = 0.0D0
      DO J=1,NCUT
        RTS = RTS + ABSK(CUTINDEX(J))
      ENDDO
      MUMSBAR = MUOVERRTS * RTS
C
C Calculate final state momenta.
C Then we can also calculate CALSVAL and the PREFACTOR.
C
      DO I = 1,NCUT
        KCUT(I,0) = ABSK(CUTINDEX(I))
        DO MU = 1,3
          KCUT(I,MU) = CUTSIGN(I) * K(CUTINDEX(I),MU)
        ENDDO
      ENDDO
      CALSVAL = CALS0(NCUT,KCUT)
      PREFACTOR = 1.0D0 / (NC * RTS**2 * (2.0D0 * PI)**NLOOPS )
C
C Calculate momenta around loop (if any). In case NINLOOP = 0, this
C DO loop is skipped.
C
      DO J = 1,NINLOOP
        DO MU = 1,3
         KINLOOP(J,MU) = LOOPSIGN(J) * K(LOOPINDEX(J),MU)
        ENDDO
      ENDDO
C
C Please note that at this point the energy in the loop, KINLOOP(J,0),
C is not calculated. We have to wait until we have a loop cut to
C do this.
C
C Now KINLOOP(J,MU) gets an imaginary part for MU = 1,2,3.
C DEFORM calculates NEWKINLOOP and the associated jacobian, JACDEFORM.
C In case NINLOOP = 0, this subroutine just returns NEWKINLOOP(MU) = 0
C and JACDEFORM = 1.
C
      CALL DEFORM(VRTX,LOOPINDEX,RTS,LEFTLOOP,RIGHTLOOP,
     >      NINLOOP,KINLOOP,NEWKINLOOP,JACDEFORM)
C
C If there is a loop, we need to go around the loop and generate
C a "loopcut."  There are three cases.
C 1) NINLOOP = 0, with NCUT = CUTMAX.  
C Then we are ready to proceed, and we should calculate only once 
C before going back to NEWCUT. Therefore we set CALCMORE to .FALSE. 
C so that we do not enter this code again.
C 2) NINLOOP = 2, with NCUT = CUTMAX - 1.
C Then the loop is a self-energy subgraph and, with our dispersive
C treatment of these graphs, there is one term, with JCUT = 1.
C 3) NINLOOP > 2, with NCUT = CUTMAX - 1  
C Then we should loop over JCUT = 1,2,...,NINLOOP and
C set CUTINDEX(CUTMAX) = LOOPINDEX(JCUT).  When we are done with this
C we set CALCMORE to .FALSE. .
C
C If we need to renormalize this loop, we will do it
C as part of the JCUT = 1 calculation.
C
C We initialize the weight, then add to it for the renormalization
C counterterm and for each loopcut.
C
      WEIGHT = 0.0D0
      MAXWEIGHT = 0.0D0
C
      JCUT = 0
      CALCMORE = .TRUE.
      DO WHILE (CALCMORE)
      IF (NINLOOP.EQ.0) THEN
         CALCMORE = .FALSE.
      ELSE
         JCUT = JCUT + 1
         CUTINDEX(CUTMAX) = LOOPINDEX(JCUT)
         CUTSIGN(CUTMAX) = LOOPSIGN(JCUT)
         IF ((NINLOOP.EQ.2).OR.(JCUT.EQ.NINLOOP)) THEN
             CALCMORE = .FALSE.
         ENDIF
      ENDIF
C
C Calculate matrix AE(P,I) relating propagator energies to energies of 
C cut lines. NOTE that the index I here is displaced by 1.
C
      DO I = 0,NLOOPS
        QE(I) = CUTINDEX(I+1)
      ENDDO
      CALL FINDA(VRTX,QE,NLOOPS,AE,QOK)
      IF (.NOT.QOK) THEN
         WRITE(NOUT,*)'AE not found'
         STOP
      ENDIF
C
C Find which propagators are which. A propagator can be exactly
C on shell even if it isn't cut if it is linked by a self-energy
C correction to a propagator that is cut. The matrix AE(P,I) will
C tell us.
C
C Define logical and sign variables:
C CUT(P) = .TRUE. if propagator P crosses the final state cut.
C CUTSIGNP(P) = CUTSIGN(I) for P = CUTINDEX(I).
C LOOPCUT(P) = .TRUE. propagator P crosses the loopcut.
C INLOOP(P) = .TRUE. if it is in a virtual loop.
C LOOPLABEL(P) = label 1,2,... counting around loop for propagator P in loop.
C SELFPROP(P) = .TRUE. if it is part of a one loop self-energy diagram
C                      or attaches to such a diagram.
C
      DO P = 1,NPROPS
         CUT(P) = .FALSE.
         LOOPCUT(P) = .FALSE.
         INLOOP(P) = .FALSE.
         CUTSIGNP(P) = 0
         LOOPLABEL(P) = 0
      ENDDO
      DO I = 1,CUTMAX
         CUT(CUTINDEX(I)) = .TRUE.
         CUTSIGNP(CUTINDEX(I)) = CUTSIGN(I)
      ENDDO
      IF (NINLOOP.GT.0) THEN
         CUT(CUTINDEX(CUTMAX)) = .FALSE.
         LOOPCUT(CUTINDEX(CUTMAX)) = .TRUE.
      ENDIF
      DO J = 1,NINLOOP
         INLOOP(LOOPINDEX(J)) = .TRUE.
         LOOPLABEL(LOOPINDEX(J)) = J
      ENDDO
C
C Calculate part of the energies of cut lines corresponding to the
C real part of the loop three-momenta. NOTE that I is displaced by 1
C in order to work with the matrix AE(P,I).
C
      DO I = 0,NLOOPS
        E(I) = CUTSIGN(I+1) * ABSK(CUTINDEX(I+1))
      ENDDO
C
C Calculate part of the propagator energies corresponding to the
C real part of the loop three-momenta.
C
      DO P = 0,NPROPS
         K(P,0) = 0.0D0
         DO I = 0,NLOOPS
           K(P,0) = K(P,0) + AE(P,I) * E(I)
         ENDDO
      ENDDO
      IF ( ABS(RTS - K(0,0)).GT.1.0D-8 ) THEN
         WRITE(NOUT,*)'Oops, the calculation of RTS did not work'
         STOP
      ENDIF
C
C Calculate the added complex loop energy. Check that we do not
C cross the cut of Sqrt(ELLSQ) by using COMPLEXSQRT(ELLSQ).
C
      IF (NINLOOP.GT.0) THEN
        KINLOOP(JCUT,0) =  LOOPSIGN(JCUT) * K(LOOPINDEX(JCUT),0)
        ELLSQ = (0.0D0,0.0D0)
        DO MU = 1,3
          ELLSQ = ELLSQ + ( KINLOOP(JCUT,MU) + NEWKINLOOP(MU) )**2
        ENDDO
        ELL = COMPLEXSQRT(ELLSQ)
        NEWKINLOOP(0) = ELL - KINLOOP(JCUT,0)
      ELSE
        NEWKINLOOP(0) = (0.0D0,0.0D0)
      ENDIF
C....
      IF (REPORT) THEN
        IF( DETAILS .AND. (NINLOOP.GT.0) ) THEN
          WRITE(NOUT,340)NEWKINLOOP(0),XXIMAG(NEWKINLOOP(1)),
     >                 XXIMAG(NEWKINLOOP(2)),XXIMAG(NEWKINLOOP(3))
340       FORMAT('NEWKINLOOP =',2(1P G12.3),' AND',3(1P G12.3))
        ENDIF
      ENDIF
C''''
C Now we calculate the complex propagator momenta.
C
      DO P = 0,NPROPS
        DO MU = 0,3
         KC(P,MU) = K(P,MU)
        ENDDO
      ENDDO
C
      DO J = 1,NINLOOP
       DO MU = 0,3
         KC(LOOPINDEX(J),MU) = KC(LOOPINDEX(J),MU)
     >       + LOOPSIGN(J) * NEWKINLOOP(MU)
        ENDDO
      ENDDO
C
C Calculate denominator.
C
C We calculate two denominators: the part from the loop propagators
C (LOOPDENOM) and the part from the other propagators (PLAINDENOM).
C The renormalization counterterm uses only PLAINDENOM.
C
C Propagators that are part of a one loop self-energy subgraph, or
C attached to such a subgraph, do not contribute to the denominator
C factor at all. The functions QPROP and GPROP, called by NUMERATOR,
C take care of the factors associated with these propagators.
C
      PLAINDENOM = 1.0D0
      LOOPDENOM = (1.0D0,0.0D0)
      DO P = 1,NPROPS
C
      IF (.NOT.SELFPROP(P)) THEN
        IF (INLOOP(P)) THEN
C
C P is in the loop:
C
          IF(LOOPCUT(P)) THEN
            LOOPDENOM   = LOOPDENOM * 2.0D0 * CUTSIGNP(P) * KC(P,0)
          ELSE
            KCSQ = 0.0D0
            DO MU = 0,3
              KCSQ = KCSQ + METRIC(MU) * KC(P,MU)**2
            ENDDO
            CALL CHECKDEFORM(KCSQ,LEFTLOOP,RIGHTLOOP,LOOPLABEL(P),JCUT,
     >                       GRAPHNUMBER,CUT,LOOPCUT)
            LOOPDENOM = LOOPDENOM * KCSQ
          ENDIF
C
        ELSE
C
C P is not in the loop:
C
          IF (CUT(P)) THEN
            PLAINDENOM = PLAINDENOM * 2.0D0 * CUTSIGNP(P) * K(P,0)
          ELSE
            KSQ = 0.0D0
            DO MU = 0,3
              KSQ = KSQ + METRIC(MU) * K(P,MU)**2
            ENDDO
            PLAINDENOM = PLAINDENOM * KSQ
          ENDIF
C
C End   IF (INLOOP(P)) ... ELSE ...
C End IF (.NOT.SELFPROP(P)) ...
C
        ENDIF
      ENDIF
C........ 
        IF (REPORT.AND.DETAILS) THEN
        IF (.NOT.SELFPROP(P)) THEN
          IF (INLOOP(P)) THEN
            IF(LOOPCUT(P)) THEN
              WRITE(NOUT,350)P,CUTSIGNP(P) * KC(P,0)
350           FORMAT('Loopcut propagator',I3,' Energy =',2(1P G12.3))
            ELSE
              WRITE(NOUT,351)P,KCSQ
351           FORMAT('Loop propagator',I3,' KCSQ =',2(1P G12.3))
            ENDIF
          ELSE
            IF (CUT(P)) THEN
              WRITE(NOUT,352)P,CUTSIGNP(P) * K(P,0)
352           FORMAT('Cut propagator',I3,' Energy =',(1P G12.3))
            ELSE
              WRITE(NOUT,353)P,KSQ
353           FORMAT('Tree propagator',I3,' KSQ =',(1P G12.3))
            ENDIF
          ENDIF
        ENDIF
        ENDIF
C'''''''
C
C End DO P = 1,NPROPS
C
      ENDDO
C
C Calculate graph.
C
C Add to contribution for this point.
C
C If we have a virtual loop with 2 or 3 lines, then we need the
C renormalization counter term. We are in a loop over JCUT. We
C will include the counter term when JCUT = 1.
C
C There is a minus sign because this is the counter term and we
C want to subtract it.
C      
      IF (((NINLOOP.EQ.2).OR.(NINLOOP.EQ.3)).AND.(JCUT.EQ.1)) THEN
C
        FEYNMAN = RNUMERATOR(GRAPHNUMBER,KC,MUMSBAR,CUT)/PLAINDENOM
C
        INTEGRAND = - PREFACTOR * JACNEWPOINT * JACDEFORM 
     >              * FEYNMAN * SMEAR(RTS)
        MAXWEIGHT = MAX(MAXWEIGHT,ABS(XXREAL(INTEGRAND)))
        WEIGHT = WEIGHT + XXREAL(INTEGRAND)
C
        INTEGRAND = INTEGRAND * CALSVAL
        MAXPART = MAX(MAXPART,ABS(XXREAL(INTEGRAND)))
        VALUE = VALUE + INTEGRAND
C....
        IF (REPORT) THEN
         IF (DETAILS) THEN
          WRITE(NOUT,360)
360       FORMAT('PREFACTOR * JACNEWPOINT * (JACDEFORM-R JACDEFORM-I)',
     >           ' (FEYNMAN-R FEYNMAN-I) * CALSVAL * SMEAR(RTS)')
          WRITE(NOUT,361)PREFACTOR,JACNEWPOINT,JACDEFORM,
     >                 FEYNMAN,CALSVAL,SMEAR(RTS)
361       FORMAT(8(1P G12.3))
         ENDIF
         WRITE(NOUT,362)INTEGRAND
362      FORMAT('Contribution (CT):',2(1P G18.10))
         IF (DETAILS) THEN
           WRITE(NOUT,*)' '
         ENDIF
        ENDIF
C''''
      ENDIF
C
C Done with the counter term (if any), now we do the main term.
C
      FEYNMAN = NUMERATOR(GRAPHNUMBER,KC,CUT)/LOOPDENOM /PLAINDENOM
      INTEGRAND = PREFACTOR * JACNEWPOINT * JACDEFORM 
     >            * FEYNMAN * SMEAR(RTS)
      MAXWEIGHT = MAX(MAXWEIGHT,ABS(XXREAL(INTEGRAND)))
      WEIGHT = WEIGHT + XXREAL(INTEGRAND)
C
      INTEGRAND = INTEGRAND * CALSVAL
      MAXPART = MAX(MAXPART,ABS(XXREAL(INTEGRAND)))
      VALUE = VALUE + INTEGRAND
C....
      IF (REPORT) THEN
        IF (DETAILS) THEN
          WRITE(NOUT,370)
370       FORMAT('PREFACTOR * JACNEWPOINT * (JACDEFORM-R JACDEFORM-I)',
     >           ' (FEYNMAN-R FEYNMAN-I) * CALSVAL * SMEAR(RTS)')
          WRITE(NOUT,371)PREFACTOR,JACNEWPOINT,JACDEFORM,
     >                 FEYNMAN,CALSVAL,SMEAR(RTS)
371       FORMAT(8(1P G12.3))
        ENDIF
        IF (NINLOOP.GT.0) THEN
          WRITE(NOUT,373),LOOPINDEX(JCUT),INTEGRAND
373       FORMAT(I3,'  Contribution:',2(1P G18.10))
        ELSE
          WRITE(NOUT,374)INTEGRAND
374       FORMAT('     Contribution:',2(1P G18.10))        
        ENDIF
        IF (DETAILS) THEN
          WRITE(NOUT,*)' '
        ENDIF
      ENDIF
C''''
C
C Compute a known integral to see if we have it right.
C Subroutine CHECKCALC calculates CHECK.
C
      CALL 
     >  CHECKCALC(GRAPHNUMBER,CUTINDEX,KC,JACNEWPOINT,JACDEFORM,CHECK)
      VALUECHK = VALUECHK + CHECK
C
C End of loop DO WHILE (CALCMORE) that runs over loopcuts.
C
      ENDDO
C
C We are ready to call Hrothgar to process the result for this cut.
C
      CALL HROTHGAR(NCUT,KCUT,WEIGHT,1,'NEWRESULT ')
C
C Close loop DO CUTNUMBER = 1,NUMBEROFCUTS(GRAPHNUMBER)
C
      ENDDO
C
C End of loop  DO JREFLECT = 1,NREFLECT 
C     
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C             End of subroutine to calculate integrand                 C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE 
     >  CHECKCALC(GRAPHNUMBER,CUTINDEX,KC,JACNEWPOINT,JACDEFORM,CHECK)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER GRAPHNUMBER,CUTINDEX(SIZE+1)
      COMPLEX*16 KC(0:3*SIZE-1,0:3)
      REAL*8 JACNEWPOINT
      COMPLEX*16 JACDEFORM
C Out:
      COMPLEX*16 CHECK
C
C Compute a known integral to see if we have it right.
C This subroutine calculates the integrand.
C The check is based on 
C   Int d^3 p [p^2 + M^2]^(-3) = Pi^2/ (4 M^3).
C   Int d^3 p [p^2 (p^2 + M^2)]^(-1) = 2 Pi^2 /M
C Note that we look at just one term in the sum over cuts
C and loopcuts: 
C For graph 10, we take Cutindex = (7,5,4,1);
C For graph 8, we take Cutindex = (8,6,4,1), etc.
C
C Latest modification: 24 December 1998
C
C Max number of graphs, cuts, maps for array sizes:
      INTEGER MAXGRAPHS,MAXMAPS
      PARAMETER (MAXGRAPHS = 10)
      PARAMETER (MAXMAPS = 64)
C Input and output units.
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C How many graphs and how many cuts and maps for each:
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C Reno size and counting variables:
      INTEGER GROUPSIZE(MAXGRAPHS,MAXMAPS)
      INTEGER GROUPSIZEGRAPH(MAXGRAPHS)
      INTEGER GROUPSIZETOTAL
      COMMON /MONTECARLO/GROUPSIZE,GROUPSIZEGRAPH,GROUPSIZETOTAL
C
      REAL*8 MM
      DATA MM /3.0D-1/
      REAL*8 PI
      DATA PI /3.1415926535898D0/
      COMPLEX*16 TEMP1,TEMP2,TEMP3
      INTEGER MU
C
C If it is not the right graph and the right cut, this default
C value will be returned.
C
      CHECK = (0.0D0,0.0D0)
C
      TEMP1 = 0.0D0
      TEMP2 = 0.0D0
      TEMP3 = 0.0D0
C
      IF (GRAPHNUMBER.EQ.10) THEN
C
        IF (   (CUTINDEX(1).EQ.7).AND.(CUTINDEX(2).EQ.5)
     >    .AND.(CUTINDEX(3).EQ.4).AND.(CUTINDEX(4).EQ.1) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(6,MU)*KC(6,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSEIF (GRAPHNUMBER.EQ.9) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.7)
     >    .AND.(CUTINDEX(3).EQ.3).AND.(CUTINDEX(4).EQ.5) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(8,MU)*KC(8,MU)
            TEMP2 = TEMP2 + KC(6,MU)*KC(6,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C      
      ELSEIF (GRAPHNUMBER.EQ.8) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.6)
     >    .AND.(CUTINDEX(3).EQ.4).AND.(CUTINDEX(4).EQ.1) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(8,MU)*KC(8,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSEIF (GRAPHNUMBER.EQ.7) THEN
C
        IF (   (CUTINDEX(1).EQ.5).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.6) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(7,MU)*KC(7,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSEIF (GRAPHNUMBER.EQ.6) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.5)
     >    .AND.(CUTINDEX(3).EQ.3).AND.(CUTINDEX(4).EQ.6) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(7,MU)*KC(7,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSEIF (GRAPHNUMBER.EQ.5) THEN
C
        IF (   (CUTINDEX(1).EQ.8).AND.(CUTINDEX(2).EQ.7)
     >    .AND.(CUTINDEX(3).EQ.3).AND.(CUTINDEX(4).EQ.1) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(5,MU)*KC(5,MU)
            TEMP2 = TEMP2 + KC(8,MU)*KC(8,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSEIF (GRAPHNUMBER.EQ.4) THEN
C
        IF (   (CUTINDEX(1).EQ.6).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.7) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSEIF (GRAPHNUMBER.EQ.3) THEN
C
        IF (   (CUTINDEX(1).EQ.5).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.6) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSEIF (GRAPHNUMBER.EQ.2) THEN
C
        IF (   (CUTINDEX(1).EQ.7).AND.(CUTINDEX(2).EQ.6)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.4) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C
      ELSEIF (GRAPHNUMBER.EQ.1) THEN
C
        IF (   (CUTINDEX(1).EQ.5).AND.(CUTINDEX(2).EQ.4)
     >    .AND.(CUTINDEX(3).EQ.1).AND.(CUTINDEX(4).EQ.7) ) THEN
          DO MU = 1,3
            TEMP1 = TEMP1 + KC(7,MU)*KC(7,MU)
            TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU)
            TEMP3 = TEMP3 + KC(1,MU)*KC(1,MU) 
          ENDDO
        ELSE
          RETURN
        ENDIF
C 
      ELSE
        WRITE(NOUT,*)'Problem with graph number in CHECKCALC'
        STOP
      ENDIF
C
      CHECK =         TEMP1 * (TEMP1 + MM**2)
      CHECK = CHECK * TEMP2 * (TEMP2 + MM**2)
      CHECK = CHECK * (TEMP3 + MM**2)**3
      CHECK = (MM**5/PI**6) /CHECK
      CHECK = JACDEFORM * JACNEWPOINT * CHECK
C
C Weight according to the number of points devoted to the current
C graph.
C
      CHECK = CHECK * GROUPSIZEGRAPH(GRAPHNUMBER)/GROUPSIZETOTAL
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE REFLECT(VRTX,ABSK,K,NEWABSK,NEWK,NREFLECT)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2)
      REAL*8 ABSK(0:3*SIZE-1)
      REAL*8 K(0:3*SIZE-1,0:3)
C Out:
      REAL*8 NEWABSK(0:3*SIZE-1)
      REAL*8 NEWK(0:3*SIZE-1,0:3)
      INTEGER NREFLECT
C
C For the given set of propagator momenta K, find the 'most collinear'
C vertex, VSTAR. Reverse the parts of the momenta at this vertex that
C are transverse to the collinear direction, while keeping a certain
C set of independent propagator momenta fixed. This gives new 
C propagator momenta NEWK.
C
C The work of finding which propagators to use to specify the
C the transformation is done by REFLECTSPEC.
C
C The transformation K --> NEWK is 'valid' if the choice of propagators
C used in the transformation is the same as the choice corresponding
C to the reflected point NEWK, as given by REFLECTSPEC.
C For a valid transformation, NREFLECT = 2.
C Otherwise, NREFLECT = 1.
C
C Note that the jacobian associated with this transformation is 1.
C
C 19 July 1994
C 26 December 1995
C 20 November 1998
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER P,L,MU
      INTEGER Q(0:SIZE),NEWQ(0:SIZE)
      INTEGER A(0:3*SIZE-1,0:SIZE)
      REAL*8 TEMP,DOT12
      REAL*8 FACTOR
      REAL*8 KSQUARE
      LOGICAL QOK
C
C NREFLECT set to 1.
C
      NREFLECT = 1
C
      IF (NLOOPS.LT.2) RETURN
C
C Get specification for the reflection transformation in the form of
C labels Q(0),...,Q(nloops) of independent loop momenta.
C
      CALL REFLECTSPEC(VRTX,ABSK,K,Q)
C
C All energies are zero, incoming momentum is zero.
C
      DO P = 0,NPROPS
        NEWK(P,0) = 0.0D0
      ENDDO
      DO MU = 1,3
        NEWK(0,MU) = 0.0D0
      ENDDO
C
C The Q(L) label our choice of independent loop momenta. Start with
C the new loop momenta equal to the old.
C
      DO L = 1,NLOOPS
      DO MU = 1,3
        NEWK(Q(L),MU) = K(Q(L),MU)
      ENDDO
      ENDDO
C
C Now the reflection transformation: we reflect K(Q(2),mu) across
C a line in the direction of K(Q(1),mu).
C
      DOT12 = 0.0D0
      DO MU = 1,3
        DOT12 = DOT12 + K(Q(1),MU)*K(Q(2),MU)
      ENDDO
      FACTOR = 2.0D0/ABSK(Q(1))**2

      DO MU = 1,3
        NEWK(Q(2),MU) = FACTOR * DOT12 * K(Q(1),MU) - K(Q(2),MU)
      ENDDO
C
C This gives the loop momenta. Now get the corresponding matrix A so
C that we can construct the other propagator momenta.
C
      CALL FINDA(VRTX,Q,NLOOPS,A,QOK)
      IF (.NOT.QOK) THEN
        WRITE(NOUT,*) 'Could not find A in REFLECT'
      ENDIF
C
C Use this matrix A to construct the propagator momenta. Note that
C when P is one of the loop momenta, say P = Q(L*), then
C A(P,L) =  delta(P,L*) so K(P,MU) doesn't change.
C
C We also need the absolute values of the new propagator momenta.
C
      DO P = 1,NPROPS
        KSQUARE = 0.0D0
        DO MU = 1,3
          TEMP = 0.0D0
          DO L = 1,NLOOPS
           TEMP = TEMP + A(P,L) * NEWK(Q(L),MU)
          ENDDO
          NEWK(P,MU) = TEMP
          KSQUARE = KSQUARE + TEMP**2
        ENDDO
        NEWABSK(P) = SQRT(KSQUARE)
      ENDDO
C
C -----        Check on validity of the transformation.           -----
C
C The transformation L --> NEWL is 'valid' if the set of Q(i) for the
C transformed point is the same as the set for the untransformed point.
C For a valid transformation, return NREFLECT = 2. For an invalid
C transformation, return NREFLECT = 1.
C 
      CALL REFLECTSPEC(VRTX,NEWABSK,NEWK,NEWQ)
C
      NREFLECT = 2
      DO L = 1,NLOOPS
        IF (NEWQ(L).NE.Q(L)) THEN
          NREFLECT = 1
        ENDIF
      ENDDO
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE REFLECTSPEC(VRTX,ABSK,K,Q)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER VRTX(0:3*SIZE-1,2)
      REAL*8 ABSK(0:3*SIZE-1)
      REAL*8 K(0:3*SIZE-1,0:3)
C Out:
      INTEGER Q(0:SIZE)
C
C Given graph specified by VRTX with momenta specified by ABSK and K,
C this subroutine returns a specification for the reflection
C transformation in the form of a list of propagator indices
C (Q(0),Q(1),...,Q(NLOOPS)) transformation.
C
C First, Q(0) = 0.
C Then, find the most collinear vertex, labelled by VSTAR:
C   At vertex V, the propagators are labelled P1,P2,P3 with
C   ABSK(P1) > ABSK(P2) > ABSK(P3). The degree of collinearity
C   is specified by Cos(Theta_23). The vertex with the largest
C   value of Cos(Theta_23) is the `most collinear.' (There is a
C   small modification of this rule, as explained later below.)
C Now, choose Q(1) = P1 and Q(2) = P3 corresponding the the most
C   collinear vertex.
C Finally, choose Q(3),...,Q(nloops) so that the Q's specify an
C independent set of loop momenta and the ABSK(Q(3)),...,ABSK(Q(nloops))
C are as small as possible.
C
C 18 July 1994 (Subroutine MOSTCOLLINEAR)
C 19 Dec 1995
C 20 Nov 1998 (New version named REFLECTSPEC)
C 30 Nov 1998
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C
      INTEGER V,L,P,P1,P2,P3
      INTEGER VSTAR,P1STAR,P2STAR,P3STAR
      INTEGER MU
      LOGICAL NONESOFAR
      INTEGER PROPSIGN
      REAL*8 COS23,COSMAX,MAXSOFAR,MINSOFAR
      REAL*8 TINY
      PARAMETER (TINY = 1.0D-10)
      REAL*8 ALMOST
      PARAMETER (ALMOST = 0.9999999999D0)
      LOGICAL LOOKMORE,DEPENDENT
      REAL*8 ABSKBEST
      INTEGER PBEST
C      
      VSTAR = 0
      COSMAX = -1.0D0
C
      DO V = 3,NVERTS
C
C For each vertex, find the indices P1, P2, P3 of the propagators
C that connect to it. Choose the labels such that
C ABSK(P1) > ABSK(P2) > ABSK(P3). The first P found is assigned
C tentatively to P1 and P3. Then when the second and third P's are
C found we proceed according to which of the following three conditions
C are met:
C
C Minsofar .LE. Maxsofar .LT. k
C Minsofar .LE.     k    .LE. Maxsofar
C        k .LT. Minsofar .LE. Maxsofar
C
      NONESOFAR = .TRUE.
      DO P = 1,NPROPS
        IF ((VRTX(P,1).EQ.V).OR.(VRTX(P,2).EQ.V)) THEN
          IF (NONESOFAR) THEN
            P1 = P
            MAXSOFAR = ABSK(P)
            P3 = P
            MINSOFAR = ABSK(P)
            NONESOFAR = .FALSE.
          ELSE
            IF (ABSK(P).GE.MINSOFAR) THEN
              IF (ABSK(P).GT.MAXSOFAR) THEN
                P2 = P1
                P1 = P
                MAXSOFAR = ABSK(P)
              ELSE
                P2 = P
              ENDIF
            ELSE
              P2 = P3
              P3 = P
              MINSOFAR = ABSK(P)
            ENDIF
          ENDIF
        ENDIF
      ENDDO
C
C Now we test to see if vertex V is 'more collinear' than VSTAR.
C We calculate cos( theta_{23} ). This quantity is nearly 1 for
C a highly collinear splitting. 
C
      COS23 = 0.0D0
      DO MU = 0,3
        COS23 = COS23 + K(P2,MU)*K(P3,MU)
      ENDDO
      COS23 = COS23/(ABSK(P2)*ABSK(P3))
      COS23 = PROPSIGN(VRTX,P2,V) * PROPSIGN(VRTX,P3,V) * COS23
C
C Compare to the previous highest value of COS23, corresponding to 
C vertex VSTAR. To take care of the possibility that vertex V is
C equivalent to VSTAR, and thus that COS23 equals COSMAX up to 
C roundoff errors, we demand that COS23 be bigger than COSMAX 
C plus a tiny amount.
C
C In the case that COS23 .GT. COSMAX + TINY, we have a new 'most
C collinear' vertex.
C
      IF (COS23.GT.COSMAX+TINY) THEN
C
        VSTAR = V
        P1STAR = P1
        P2STAR = P2
        P3STAR = P3
        COSMAX = COS23
C
C Close IF (COS23.GT.COSMAX).
C
      ENDIF
C
C End DO V = 3,NVERTS.
C
      ENDDO
C
C Now we have the most collinear vertex VSTAR, attached to
C propagators P1,P2,P3. We set Q(0)= 0, Q(1) = P1STAR, 
C and Q(2) = P3STAR, so that, in subroutine REFLECT, we 
C will reflect K(P3STAR,mu) across a line in the direction
C of K(P1STAR,mu).
C
      Q(0) = 0
      Q(1) = P1STAR
      Q(2) = P3STAR
C
C Now we look for more independent propagators with
C labels Q(3),...,Q(n) in the case of n loops.
C For each L from L = 3 to L = NLOOPS, we find Q(L) that labels
C the propagator with the smallest momentum, ABSK(Q(L)), provided
C that (Q(1), Q(2),...,Q(L)) is an independent set of propagators.
C In subroutine REFLECT, we will hold the corresponding propagator
C momenta fixed.
C
      DO L = 3,NLOOPS
C
C First, we need to find at least one possible choice for the
C best propagator index, PBEST.
C
        LOOKMORE = .TRUE.
        P=0
        DO WHILE (LOOKMORE)
          P = P+1
          IF(P.GT.NPROPS) THEN
            WRITE(NOUT,*)'Failure in REFLECTDATA'
            STOP
          ENDIF
          Q(L) = P
          IF (.NOT.DEPENDENT(VRTX,Q,L)) THEN
            ABSKBEST = ABSK(P)
            PBEST = P
            LOOKMORE = .FALSE.
          ENDIF
        ENDDO
C
C Now, we need to get the best possible choice, PBEST. The comparison
C uses ALMOST = 0.999... to force a choice based on propagator
C labelling rather than roundoff errors in the case of two propagators
C that, in principle, carry the same momentum.  
C
        LOOKMORE = .TRUE.
        DO WHILE (LOOKMORE)
          P = P+1
          IF(P.GT.NPROPS) THEN
            LOOKMORE = .FALSE.
          ELSE
            IF (ABSK(P).LT.ALMOST*ABSKBEST) THEN
              Q(L) = P
              IF (.NOT.DEPENDENT(VRTX,Q,L)) THEN
                PBEST = P
                ABSKBEST = ABSK(P)
                LOOKMORE = .FALSE.
              ENDIF
            ENDIF
          ENDIF
        ENDDO
C
C We have the best choice for this L. Record Q(L).
C
      Q(L) = PBEST
C
C Close loop over L = 3,NLOOPS.
C
      ENDDO
C
C Now we have (Q(1),...,Q(NLOOPS)), so we are done.
C
      RETURN
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      REAL*8 FUNCTION
     >    DENSITY(GRAPHNUMBER,K,ABSK,QS,A1S,A2S,MAPTYPES,NMAPS)
C
      INTEGER SIZE
      PARAMETER (SIZE = 3)
C In:
      INTEGER GRAPHNUMBER
      REAL*8 K(0:3*SIZE-1,0:3),ABSK(0:3*SIZE-1)
      INTEGER QS(256,0:SIZE),A1S(256,2:SIZE),A2S(256,2:SIZE)
      CHARACTER*6 MAPTYPES(256)
      INTEGER NMAPS
C
C Density of Monte Carlo points as a function of |K(p)|'s.
C We calculate the jacobian J of associated with the each map,
C then take DENSITY = 1/J, and finally sum DENSITY over maps.
C 
C 29 June 1993
C 12 July 1993
C 17 July 1994
C  4 May 1996
C 21 November 1996
C  5 December 1996
C  5 February 1997
C 15 December 1998
C 23 December 1998
C  9 February 1999
C 10 March    1999
C 20 August   1999
C
C Max number of graphs and maps for array sizes:
      INTEGER MAXGRAPHS,MAXMAPS
      PARAMETER (MAXGRAPHS = 10)
      PARAMETER (MAXMAPS = 64)
C
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C Parameters governing the choice of points
      REAL*8 ENERGYSCALE
      COMMON /MSCALE/ ENERGYSCALE
      REAL*8 AS,BS,AM,BM,AH,BH,AUV,BUV,CTHETA,AC,BC,FC
      COMMON /POINTS/ AS,BS,AM,BM,AH,BH,AUV,BUV,CTHETA,AC,BC,FC
      INTEGER GROUPSIZE(MAXGRAPHS,MAXMAPS)
      INTEGER GROUPSIZEGRAPH(MAXGRAPHS)
      INTEGER GROUPSIZETOTAL
      COMMON /MONTECARLO/GROUPSIZE,GROUPSIZEGRAPH,GROUPSIZETOTAL
C How many graphs and how many cuts and maps for each:
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C
      INTEGER Q(0:SIZE),A1(2:SIZE),A2(2:SIZE)
      CHARACTER*6 MAPTYPE
      REAL*8 MAGK,K0,PREVIOUSK,PREVIOUSKT
      REAL*8 K1SQ,K2SQ,K1ABS,K2ABS,TEMP1,TEMP2,NORMA
      REAL*8 K1(3),K2(3),K1HAT(3),K2HAT(3),EA(3)
      REAL*8 EPSILON,COSTHETA,GAMMA
      REAL*8 N(3),EPS,TEMP,RHOCOLL,ONEMINUSZ,ONEPLUSZ
      INTEGER MAPNUMBER,L,MU
      REAL*8 J
      REAL*8 KRATIO
      REAL*8 FOURPI
      REAL*8 SINHINV
      PARAMETER (FOURPI = 12.5663706143591730D0)
C
      IF (NLOOPS.LT.1) THEN
        WRITE(NOUT,*) 'NLOOPS less than 1 in DENSITY'
        STOP
      ENDIF
C
C We construct the density as a sum.
C
      DENSITY = 0.0D0
      DO MAPNUMBER = 1,NMAPS
C
      MAPTYPE = MAPTYPES(MAPNUMBER)
      DO L = 0,NLOOPS
       Q(L) = QS(MAPNUMBER,L)
      ENDDO
      DO L = 2,NLOOPS
        A1(L) = A1S(MAPNUMBER,L)
        A2(L) = A2S(MAPNUMBER,L)
      ENDDO
C
C We construct the jacobian J as a product.
C
      J = 1.0D0
C
C Logical structure below is
C IF (MAPTYPE.NE.'UVHARD') THEN
C     IF (MAPTYPE.EQ.'MEDIUM') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRPART') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRFULL') THEN ...
C     ELSE  <error>
C     ENDIF
C ELSE   <maptype must have been UVHARD> ...
C ENDIF
C
      IF (MAPTYPE.NE.'UVHARD') THEN
C
C We have a soft map. To start, we have a "not-soft" but "not-UV"
C map for k_{Q(N)}. This map has AH = 3, BH = 3 normally.
C We also make a unit vector N(mu) in the direction of k_{Q(N)}. 
C
      MAGK = ABSK(Q(NLOOPS))
      K0 = ENERGYSCALE
      KRATIO = MAGK/K0
      J = J * FOURPI /BH * K0**3 * KRATIO**(3.0D0 - AH)
     >      * ( 1.0D0 + KRATIO**AH )**(1.0D0 + BH/AH)
      DO MU = 1,3
        N(MU) = K(Q(NLOOPS),MU)/MAGK
      ENDDO
C
C Now we have to complete the jacobian. What kind
C of map do we have? The choices correspond to the MAPTYPE:
C MEDIUM for an IR map that with a mild soft singularity.
C IRPART for infrared map with Q(1) points concentrated in a line.
C IRFULL for infrared map with Q(1) points concentrated in a plane.
C For maptype IRFULL, we use identifiers SOFTP1,SOFTP2,SOFTSIGN
C in constructing the map.
C
      IF (MAPTYPE.EQ.'MEDIUM') THEN
C
C For the other k_{Q(i)}, we use a nested "mediumsoft" map
C with AM = 2, BM = 2 normally.  
C For Nloops = 3, this evaluates K_{Q(2)}, K_{Q(1)}.
C
        PREVIOUSK = ENERGYSCALE
        DO L = NLOOPS-1,1,-1
          MAGK = ABSK(Q(L))
          K0 = PREVIOUSK
          KRATIO = MAGK/K0
          J = J * FOURPI /BM * K0**3 * KRATIO**(3.0D0 - AM)
     >          * ( 1.0D0 + KRATIO**AM )**(1.0D0 + BM/AM)
          PREVIOUSK = MAGK
        ENDDO
C
      ELSE IF (MAPTYPE.EQ.'IRPART') THEN
C
C For the other k_{Q(i)} except K_{Q(1)}, we use a nested "soft" map
C with AS = 1, BS = 2 normally. For these maps, we have
C concentrated the points in the direction of N(mu).
C For Nloops = 3, this evaluates K_{Q(2)}, K_{Q(1)}.
C
C First, we generate the jacobian for the map with a uniform angular
C distribution. We also calculate epsilon, the parameter
C that is used below in the angular distribution.
C
        PREVIOUSK = ENERGYSCALE
        PREVIOUSKT = ENERGYSCALE
        DO L = NLOOPS-1,1,-1
          MAGK = ABSK(Q(L))
          IF (L.EQ.(NLOOPS - 1)) THEN
            EPS = 1.0D0
          ELSE
            EPS = PREVIOUSKT/MAX(MAGK,PREVIOUSKT)
          ENDIF
          K0 = PREVIOUSK
          KRATIO = MAGK/K0
          J = J * FOURPI /BS * K0**3 * KRATIO**(3.0D0 - AS)
     >      * ( 1.0D0 + KRATIO**AS )**(1.0D0 + BS/AS)
C
C Modify the jacobian to account for our mapping toward the direction
C collinear to N(mu).
C
C First, we need Z = cosine of the angle between K(Q(L),mu) and N(mu),
C ONEMINUSZ = 1 - Z and ONEPLUSZ = 1 + Z.
C
          ONEMINUSZ = 0.0D0
          ONEPLUSZ  = 0.0D0
          DO MU = 1,3
           ONEMINUSZ = ONEMINUSZ + 0.5D0*(K(Q(L),MU)/MAGK - N(MU))**2
           ONEPLUSZ  = ONEPLUSZ  + 0.5D0*(K(Q(L),MU)/MAGK + N(MU))**2
          ENDDO
C
C Now we have two possibilities for the density of points in angle
C depending on the value of 1 - |Z| = MIN(ONEMINUSZ,ONEPLUSZ).
C
          IF ( MIN(ONEMINUSZ,ONEPLUSZ).LT.EPS ) THEN
            TEMP = EPS**2/MIN(ONEMINUSZ,ONEPLUSZ)
            RHOCOLL = (1.0D0 - FC + EPS*FC)*(1 + EPS**AC)**(BC/AC)
            RHOCOLL = RHOCOLL * BC /EPS**2
            RHOCOLL = RHOCOLL * TEMP**(1.0D0 + AC)
            RHOCOLL = RHOCOLL * (1.0D0 + TEMP**AC)**(-(1.0D0 + BC/AC))
          ELSE
            RHOCOLL = FC
          ENDIF
          J = J/RHOCOLL
C
          PREVIOUSKT = MAGK * SQRT(ONEPLUSZ*ONEMINUSZ)
          PREVIOUSK = MAGK
C
        ENDDO
C
      ELSE IF (MAPTYPE.EQ.'IRFULL') THEN
C
C Here we have choosen the same map as above for IRPART for K_{Q(2)},
C or, more generally, for K_{Q(Nloops - 1)},...,K_{Q(2)}. Then, we
C have done something special for K_{Q(1)}.
C -----
C Begin copy of IRPART code with "DO L = NLOOPS-1,1,-1" changed
C to "DO L = NLOOPS-1,2,-1":
C
        PREVIOUSK = ENERGYSCALE
        PREVIOUSKT = ENERGYSCALE
        DO L = NLOOPS-1,2,-1
          MAGK = ABSK(Q(L))
          IF (L.EQ.(NLOOPS - 1)) THEN
            EPS = 1.0D0
          ELSE
            EPS = PREVIOUSKT/MAX(MAGK,PREVIOUSKT)
          ENDIF
          K0 = PREVIOUSK
          KRATIO = MAGK/K0
          J = J * FOURPI /BS * K0**3 * KRATIO**(3.0D0 - AS)
     >      * ( 1.0D0 + KRATIO**AS )**(1.0D0 + BS/AS)
          ONEMINUSZ = 0.0D0
          ONEPLUSZ  = 0.0D0
          DO MU = 1,3
           ONEMINUSZ = ONEMINUSZ + 0.5D0*(K(Q(L),MU)/MAGK - N(MU))**2
           ONEPLUSZ  = ONEPLUSZ  + 0.5D0*(K(Q(L),MU)/MAGK + N(MU))**2
          ENDDO
          IF ( MIN(ONEMINUSZ,ONEPLUSZ).LT.EPS ) THEN
            TEMP = EPS**2/MIN(ONEMINUSZ,ONEPLUSZ)
            RHOCOLL = (1.0D0 - FC + EPS*FC)*(1 + EPS**AC)**(BC/AC)
            RHOCOLL = RHOCOLL * BC /EPS**2
            RHOCOLL = RHOCOLL * TEMP**(1.0D0 + AC)
            RHOCOLL = RHOCOLL * (1.0D0 + TEMP**AC)**(-(1.0D0 + BC/AC))
          ELSE
            RHOCOLL = FC
          ENDIF
          J = J/RHOCOLL
          PREVIOUSKT = MAGK * SQRT(ONEPLUSZ*ONEMINUSZ)
          PREVIOUSK = MAGK
        ENDDO
C
C End copy of IRPART code with "DO L = NLOOPS-1,1,-1" changed
C to "DO L = NLOOPS-1,2,-1".
C -----
C Next, we have generated K_{Q(1)}. We used the AS, BS map, 
C nested with the previous maps. BUT this
C time we created a special concentration of points near the
C plane theta = 0 in a certain coordinate system.
C
C Step 1:
C We need the appropriate unit vectors for the map.
C
      K1SQ = 0.0D0
      K2SQ = 0.0D0
      DO MU = 1,3
        TEMP1 = 0.0D0
        TEMP2 = 0.0D0
        DO L = 2,NLOOPS
          TEMP1 = TEMP1 + A1(L) * K(Q(L),MU)
          TEMP2 = TEMP2 + A2(L) * K(Q(L),MU)
        ENDDO
        K1(MU) = TEMP1
        K1SQ = K1SQ + TEMP1**2
        K2(MU) = TEMP2
        K2SQ = K2SQ + TEMP2**2
      ENDDO
      IF ((K1SQ.LT.1.0D-30).OR.(K2SQ.LT.1.0D-30)) THEN
        WRITE(NOUT,*)'K1SQ or K2SQ too small in DENSITY'
        STOP
      ENDIF  
      K1ABS = SQRT(K1SQ)
      K2ABS = SQRT(K2SQ)
      DO MU = 1,3
        K1HAT(MU) = K1(MU)/K1ABS
        K2HAT(MU) = K2(MU)/K2ABS
      ENDDO
C
      NORMA = 0.0D0
      DO MU = 1,3
        NORMA = NORMA + (K1HAT(MU) - K2HAT(MU))**2
      ENDDO
      NORMA = SQRT(NORMA)
      IF ((NORMA.LT.0.99D-10)) THEN
        WRITE(NOUT,*)'NORMA too small in DENSITY'
        STOP
      ENDIF
      DO MU = 1,3
        EA(MU) = (K1HAT(MU) - K2HAT(MU))/NORMA
      ENDDO
C
C Step 2:
C We want a soft map with a contentration of points near cos(theta) = 0.
C First, we have the jacobian for the AS = 1, BS = 2 (normally) map, 
C nested with the previous maps.
C
        MAGK = ABSK(Q(1))
        K0 = PREVIOUSK
        KRATIO = MAGK/K0
        J = J * FOURPI /BS * K0**3 * KRATIO**(3.0D0 - AS)
     >        * ( 1.0D0 + KRATIO**AS )**(1.0D0 + BS/AS)
C
C Step 3:
C Construct the "extra" jacobian for the angular map.
C
      EPSILON = CTHETA * MAGK/K0
      COSTHETA = 0.0D0
      DO MU = 1,3
        COSTHETA = COSTHETA + K(Q(1),MU) * EA(MU)
      ENDDO
      COSTHETA = COSTHETA/MAGK
      GAMMA = SINHINV(1.0D0/EPSILON)
      J = J * GAMMA * SQRT(COSTHETA**2 + EPSILON**2)
C
C End IF (MAPTYPE.EQ.'MEDIUM') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRPART') THEN ...
C     ELSE IF (MAPTYPE.EQ.'IRFULL') THEN ...
C
      ELSE
        WRITE(NOUT,*)'This cannot happen in DENSITY'
        STOP
      ENDIF
C
C This completes the IF (MAPTYPE.NE.'UVHARD') THEN ...
C
      ELSE
C
C Now we generate the jacobian for nearly UV divergent loops. There
C are 'medium hard' maps for k_{Q(i)}, i = 1,...,NLOOPS -1 and
C a 'hard' map for the possible virutal loop momentum k_{Q(NLOOPS)}.
C
C The AUV = 3, BUV = 1 (normally) map for loop momentum NLOOPS.
C
        K0 = ENERGYSCALE
        MAGK = ABSK(Q(NLOOPS))
        KRATIO = MAGK/K0
        J = J * FOURPI /BUV * K0**3 * KRATIO**(3.0D0 - AUV)
     >        * ( 1.0D0 + KRATIO**AUV )**(1.0D0 + BUV/AUV)
C
C
C The AH = 3, BH = 3 (normally) map for loop momentum NLOOPS-1.
C We generate the appropriate jacobian and also calculate the
C vector N(mu) used to define a collinear direction.
C
        L = NLOOPS - 1
        K0 = ENERGYSCALE
        MAGK = ABSK(Q(L))
        KRATIO = MAGK/K0
        J = J * FOURPI /BH * K0**3 * KRATIO**(3.0D0 - AH)
     >        * ( 1.0D0 + KRATIO**AH )**(1.0D0 + BH/AH)
        DO MU = 1,3
          N(MU) = K(Q(L),MU)/MAGK
        ENDDO
C
C
C For NLOOPS = 3, we have generated an AS = 1, BS = 2 (normally) 
C map for loop momentum 1. This map includes a concentration
C of soft points and of points in the direction of ELL(NLOOPS-1). 
C (In case NLOOPS = 2, we do nothing. In case NLOOPS > 3, we simply
C repeat the calculation for the same map.)
C
        DO L = NLOOPS-2,1,-1
C
        K0 = ENERGYSCALE
        MAGK = ABSK(Q(L))
        KRATIO = MAGK/K0
        J = J * FOURPI /BS * K0**3 * KRATIO**(3.0D0 - AS)
     >        * ( 1.0D0 + KRATIO**AS )**(1.0D0 + BS/AS)
C
C Modify the jacobian to account for our mapping toward the direction
C collinear to N(mu).
C
C First, we need Z = cosine of the angle between K(Q(L),mu) and N(mu),
C ONEMINUSZ = 1 - Z and ONEPLUSZ = 1 + Z.
C
          ONEMINUSZ = 0.0D0
          ONEPLUSZ  = 0.0D0
          DO MU = 1,3
           ONEMINUSZ = ONEMINUSZ + 0.5D0*(K(Q(L),MU)/MAGK - N(MU))**2
           ONEPLUSZ  = ONEPLUSZ  + 0.5D0*(K(Q(L),MU)/MAGK + N(MU))**2
          ENDDO
C
C Now we evaluate the density of points in angle.
C
          TEMP = 1.0D0/MIN(ONEMINUSZ,ONEPLUSZ)
          RHOCOLL = BC * (1.0D0 + 1.0D0**AC)**(BC/AC)
          RHOCOLL = RHOCOLL * TEMP**(1.0D0 + AC)
          RHOCOLL = RHOCOLL * (1.0D0 + TEMP**AC)**(-(1.0D0 + BC/AC))
          J = J/RHOCOLL
C
C End of loop DO L = NLOOPS-2,1,-1.
C
        ENDDO

C End of IF (MAPTYPE.NE.'UVHARD') THEN ... ELSE ...
C
      ENDIF
C
      DENSITY = DENSITY + GROUPSIZE(GRAPHNUMBER,MAPNUMBER) /J
C 
C End of DO MAPNUMBER = 1,NMAPS
C
      ENDDO
      RETURN
C
      END
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C                  Subroutines associated with NEWCUT                  2
C23456789012345678901234567890123456789012345678901234567890123456789012
C
      SUBROUTINE MAKECUTINFO
C
C Input: none
C Output: none
C
C This subroutine creates the information about the various cuts
C of each graph and stores it in the common block CUTINFORMATION.
C Also counts the number of graphs, the number of cuts for each
C graph, and the number of maps for each cut and stores this in the
C common block GRAPHCOUNTS.
C
C In early versions, beowulf called NEWCUT to generate the cut
C information each time a new cut was started, but this happens so
C often that some thirty percent of the computer time was devoted to
C NEWCUT. 
C
C Latest revision: 5 January 1999.
C
C Array sizes. (We check MAXGRAPHS,MAXCUTS,MAXMAP here.):
      INTEGER SIZE,MAXGRAPHS,MAXCUTS,MAXMAPS
      PARAMETER (SIZE = 3)
      PARAMETER (MAXGRAPHS = 10)
      PARAMETER (MAXCUTS = 9)
      PARAMETER (MAXMAPS = 64)
C Input and output units.      
      INTEGER NIN,NOUT
      COMMON /IOUNIT/ NIN,NOUT
C Graph size variables.
      INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX
      COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX
C Information on cut structure:
      INTEGER NCUTINFO(MAXGRAPHS,MAXCUTS)
      INTEGER ISIGNINFO(MAXGRAPHS,MAXCUTS,3*SIZE + 1)
      INTEGER CUTINDEXINFO(MAXGRAPHS,MAXCUTS,SIZE + 1)
      INTEGER CUTSIGNINFO(MAXGRAPHS,MAXCUTS,SIZE + 1)
      LOGICAL LEFTLOOPINFO(MAXGRAPHS,MAXCUTS)
      LOGICAL RIGHTLOOPINFO(MAXGRAPHS,MAXCUTS)
      INTEGER NINLOOPINFO(MAXGRAPHS,MAXCUTS)
      INTEGER LOOPINDEXINFO(MAXGRAPHS,MAXCUTS,SIZE+1)
      INTEGER LOOPSIGNINFO(MAXGRAPHS,MAXCUTS,SIZE+1)
      COMMON /CUTINFORMATION/ NCUTINFO,ISIGNINFO,
     >   CUTINDEXINFO,CUTSIGNINFO,LEFTLOOPINFO,RIGHTLOOPINFO,
     >   NINLOOPINFO,LOOPINDEXINFO,LOOPSIGNINFO
C
      INTEGER NUMBEROFGRAPHS
      INTEGER NUMBEROFCUTS(MAXGRAPHS)
      INTEGER NUMBEROFMAPS(MAXGRAPHS)
      COMMON /GRAPHCOUNTS/ NUMBEROFGRAPHS,NUMBEROFCUTS,NUMBEROFMAPS
C
C NEWGRAPH variables:
      INTEGER VRTX(0:3*SIZE-1,2),PROP(2*SIZE,3)
      LOGICAL SELFPROP(3*SIZE-1)
      LOGICAL GRAPHFOUND
      INTEGER GRAPHNUMBER
C NEWCUT variables:
      INTEGER CUTINDEX(SIZE+1),CUTSIGN(SIZE+1),NCUT
      INTEGER ISIGN(3*SIZE-1)
      LOGICAL LEFTLOOP,RIGHTLOOP
      INTEGER LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1),NINLOOP
C FINDQS variables:
      INTEGER NMAPS,QS(256,0:SIZE)
      INTEGER A1S(256,2:SIZE),A2S(256,2:SIZE)
      CHARACTER*6 MAPTYPES(256)
C
      LOGICAL NEWCUTINIT,CUTFOUND
C
      INTEGER P,I,NP
      INTEGER CUTNUMBER
C
C---------
C Get a new graph.
C
      GRAPHFOUND = .TRUE.
      GRAPHNUMBER = 0
C
      DO WHILE (GRAPHFOUND)
      CALL NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND)
      IF (GRAPHFOUND) THEN
      GRAPHNUMBER = GRAPHNUMBER + 1
C
C Get a new cut.
C
      CUTFOUND = .TRUE.
      NEWCUTINIT = .TRUE.
      CUTNUMBER = 0
      DO WHILE (CUTFOUND)
      CALL NEWCUT(VRTX,NEWCUTINIT,NCUT,ISIGN,
     >            CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP,
     >            NINLOOP,LOOPINDEX,LOOPSI