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 (*) 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 C SUBROUTINE SOFTHELP C SUBROUTINE NEWCUT C SUBROUTINE FINDA C FUNCTION PROPSIGN C SUBROUTINE NEWCUT C SUBROUTINE GETCOUNTS (*) C SUBROUTINE HROTHGAR (*) C SUBROUTINE RANDOMINIT C SUBROUTINE NEWRAN C SUBROUTINE DAYTIME (*) C SUBROUTINE VERSION C SUBROUTINE RENO C SUBROUTINE DIAGNOSTIC C SUBROUTINE NEWGRAPH ... C SUBROUTINE FINDQS ... C SUBROUTINE CHECKPOINT C SUBROUTINE CALCULATE C end C C SUBROUTINE RENO C SUBROUTINE TIMING (*) C SUBROUTINE HROTHGAR (*) C SUBROUTINE NEWGRAPH ... C SUBROUTINE FINDQS ... C SUBROUTINE FINDA C FUNCTION PROPSIGN C SUBROUTINE NEWPOINT C SUBROUTINE AXES C FUNCTION RANDOM C SUBROUTINE NEWRAN C SUBROUTINE CHECKPOINT C SUBROUTINE CALCULATE C return C C SUBROUTINE CALCULATE C SUBROUTINE REFLECTSPEC C FUNCTION DENSITY C FUNCTION CALS0 (*) C FUNCTION THRUST (*) C SUBROUTINE GETCUTINFO C SUBROUTINE DEFORM C SUBROUTINE FINDA 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 (*) C return C C SUBROUTINE FINDA 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 C ENDIF C ELSE ... 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 C ENDIF C ELSE ... 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,LOOPSIGN,CUTFOUND) IF (CUTFOUND) THEN CUTNUMBER = CUTNUMBER + 1 C NCUTINFO(GRAPHNUMBER,CUTNUMBER) = NCUT DO P = 1,NPROPS ISIGNINFO(GRAPHNUMBER,CUTNUMBER,P) = ISIGN(P) ENDDO DO I = 1,CUTMAX CUTINDEXINFO(GRAPHNUMBER,CUTNUMBER,I) = CUTINDEX(I) CUTSIGNINFO(GRAPHNUMBER,CUTNUMBER,I) = CUTSIGN(I) ENDDO LEFTLOOPINFO(GRAPHNUMBER,CUTNUMBER) = LEFTLOOP RIGHTLOOPINFO(GRAPHNUMBER,CUTNUMBER) = RIGHTLOOP NINLOOPINFO(GRAPHNUMBER,CUTNUMBER) = NINLOOP DO NP = 1,CUTMAX LOOPINDEXINFO(GRAPHNUMBER,CUTNUMBER,NP) = LOOPINDEX(NP) LOOPSIGNINFO(GRAPHNUMBER,CUTNUMBER,NP) = LOOPSIGN(NP) ENDDO C C Close loop DO WHILE (CUTFOUND), IF (CUTFOUND) THEN C ENDIF ENDDO C IF (CUTNUMBER.GT.MAXCUTS) THEN WRITE(NOUT,*)'More cuts than I thought.' STOP ENDIF NUMBEROFCUTS(GRAPHNUMBER) = CUTNUMBER 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 IF (NMAPS.GT.MAXMAPS) THEN WRITE(NOUT,*)'Ooops, more maps than we anticipated.' STOP ENDIF C C Close loop DO WHILE (GRAPHFOUND), IF (GRAPHFOUND) THEN C ENDIF ENDDO C IF (GRAPHNUMBER.GT.MAXGRAPHS) THEN WRITE(NOUT,*)'More graphs than I thought.' STOP ENDIF NUMBEROFGRAPHS = GRAPHNUMBER C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE GETCUTINFO(GRAPHNUMBER,CUTNUMBER,NCUT,ISIGN, > CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP, > NINLOOP,LOOPINDEX,LOOPSIGN) C INTEGER SIZE PARAMETER (SIZE = 3) C Input: INTEGER GRAPHNUMBER,CUTNUMBER C Output: INTEGER NCUT,ISIGN(3*SIZE-1) INTEGER CUTINDEX(SIZE+1),CUTSIGN(SIZE+1) LOGICAL LEFTLOOP,RIGHTLOOP INTEGER NINLOOP,LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1) C C This subroutine reads from the information recorded in the common C CUTINFORMATION and returns the information relevant for C the current graph,specified by GRAPHNUMBER and the current cut, C specified by CUTNUMBER. See the subroutine NEWCUT for definition C of the variables returned. C C Latest revision: 5 January 1999. C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX C C Information on cut structure: INTEGER MAXGRAPHS,MAXCUTS PARAMETER (MAXGRAPHS = 10) PARAMETER (MAXCUTS = 9) 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 P,I,NP C C-------- C NCUT = NCUTINFO(GRAPHNUMBER,CUTNUMBER) DO P = 1,NPROPS ISIGN(P) = ISIGNINFO(GRAPHNUMBER,CUTNUMBER,P) ENDDO DO I = 1,CUTMAX CUTINDEX(I) = CUTINDEXINFO(GRAPHNUMBER,CUTNUMBER,I) CUTSIGN(I) = CUTSIGNINFO(GRAPHNUMBER,CUTNUMBER,I) ENDDO LEFTLOOP = LEFTLOOPINFO(GRAPHNUMBER,CUTNUMBER) RIGHTLOOP = RIGHTLOOPINFO(GRAPHNUMBER,CUTNUMBER) NINLOOP = NINLOOPINFO(GRAPHNUMBER,CUTNUMBER) DO NP = 1,CUTMAX LOOPINDEX(NP) = LOOPINDEXINFO(GRAPHNUMBER,CUTNUMBER,NP) LOOPSIGN(NP) = LOOPSIGNINFO(GRAPHNUMBER,CUTNUMBER,NP) ENDDO C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE NEWCUT(XVRTX,NEWCUTINIT,XNCUT,XISIGN, > XCUTINDEX,XCUTSIGN,XLEFTLOOP,XRIGHTLOOP, > XNINLOOP,XLOOPINDEX,XLOOPSIGN,CUTFOUND) C INTEGER SIZE PARAMETER (SIZE = 3) C Input: INTEGER XVRTX(0:3*SIZE-1,2) C Input and output: LOGICAL NEWCUTINIT C Output: INTEGER XCUTINDEX(SIZE+1),XCUTSIGN(SIZE+1),XNCUT INTEGER XISIGN(3*SIZE-1) LOGICAL XLEFTLOOP,XRIGHTLOOP INTEGER XLOOPINDEX(SIZE+1),XLOOPSIGN(SIZE+1),XNINLOOP LOGICAL CUTFOUND C C This subroutine generates valid cuts for a given graph. C In its present form, it generates cuts with CUTMAX lines cut C and with (CUTMAX - 1) lines cut. In the case that (CUTMAX - 1) C lines are cut, it also finds the (single) virtual loop. C C The action of the subroutine depends on its state when called. It C has two possible states. If NEWCUTINIT is true, NEWCUT is ready C start generating cuts for a new graph. This is the state when the C subroutine is called for the first time. If NEWCUTINIT is false, the C subroutine is ready to generate a new cut for the current graph. C When NEWCUT is called with NEWCUTINIT = False but it cannot find C a new cut, it exits with NEWCUTINIT = True and the output variable C CUTFOUND = False. This tells the mainprogram to produce a new graph. C C Notation: [X variables are interchanged with the external world.] C VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of C of propagator P. Specifies the supergraph. 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,...,NCUT. C CUTSIGN(I) = Sign of cut propagator I (+1 if from Left to Right). C If NCUT = CUTMAX - 1 then there is a virtual loop and C we define CUTINDEX(CUTMAX) = 0. C But in subroutine LOOPCUT we will let I = CUTMAX designate C the loopcut: C CUTINDEX(CUTMAX) = LOOPINDEX(JCUT) C CUTSIGN(CUTMAX) = LOOPSIGN(CUTMAX) 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 The loop begins (with NP = 1) at the starting vertex that is C defined as the current vertex if the the loop includes the C current vertex or the vertex in the loop that is attached to the C uncut propagator of the lowest index P. C LOOPSIGN(NP) = 1 if propagator direction is same as loop direction. C -1 if direction is opposite to loop direction. C CUTFOUND = False if we can't find a next cut. C C CUT(P) = True iff propagator P is cut. C VALIDCUT = True iff the cut is OK. C LEFT(V) = True iff vertex V is to the left of the cut. C Vertex 1, the left hand current, is always in LEFT. C NLEFT = Number of vertices to the left of the cut. C RIGHT(V) = True iff vertex V is to the right of the cut. C Vertex 2, the right hand current, is always in RIGHT. C NRIGHT = Number of vertices to the right of the cut. C C LOOPVERTEX(V) = True if vertex V is in the loop. C LOOPPROP(P) = True if propagator P is in the loop. C NCONNECTED = Number of propagators connected to a vertex. C STARTVERTEX = Starting vertex in a loop. C HOTVERTEX = Vertex to which next loop propagator should be added. C C Logical state variables: C NEWCUTINIT = True if NEWCUT called for first time with a new graph, C else False. C C ----- Outline ----- C C Output variables -> default values. C IF (NEWCUTINIT) THEN C Initialize, including NEXTCUT = True. C ENDIF C START = .FALSE. C VALIDCUT = .FALSE. C DO WHILE (.NOT.VALIDCUT) C Generate next cut, specified by CUTINDEX(I). C Here CUTINDEX(CUTMAX) = 0 indicates a virtual loop. C In this case, set NCUT = CUTMAX - 1. Else NCUT = CUTMAX. C Check cut, setting VALIDCUT to True if cut is OK. C ENDDO C IF (NCUT.EQ.CUTMAX) THEN C Set LEFTLOOP and RIGHTLOOP to False C Normal Return. C ELSE set LEFTLOOP or RIGHTLOOP to True as appropriate. C ENDIF C Find the loop, determining LOOPINDEX(NP) and LOOPSIGN(NP), C where NP = 1 designates the propagator starting at the C starting vertex defined above. C NINLOOP = number of propagators in loop. C Normal Return. C C Normal Return: C Set output variables to values of internal variables. C C 2 November 1992 C 26 January 1992 C 20 June 1993 C 21 August 1993 C 26 June 1994 C 11 July 1994 C 20 March 1996 C C ------------------- C C----------------------------------------------------------------------- C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX C INTEGER V,P INTEGER I,II,IP,V1,V2 LOGICAL LEFT(2*SIZE),RIGHT(2*SIZE) LOGICAL CHANGE LOGICAL VALIDCUT LOGICAL LOOPVERTEX(2*SIZE),LOOPPROP(3*SIZE-1) INTEGER NLEFT,NRIGHT,NCONNECTED LOGICAL CUT(3*SIZE-1),LEFTLOOP,RIGHTLOOP INTEGER CUTSIGN(SIZE+1),NCUT INTEGER LOOPINDEX(SIZE+1),LOOPSIGN(SIZE+1) INTEGER NINLOOP LOGICAL LOOKMORE INTEGER STARTVERTEX,HOTVERTEX,PREVIOUSPROP C C Internal state variables to be saved C INTEGER VRTX(0:3*SIZE-1,2) SAVE VRTX INTEGER CUTINDEX(SIZE+1) SAVE CUTINDEX C C Default values for the output variables C XNCUT = 0 XLEFTLOOP = .FALSE. XRIGHTLOOP = .FALSE. XNINLOOP = 0 DO I = 1,CUTMAX XCUTINDEX(I) = 0 XCUTSIGN(I) = 0 XLOOPINDEX(I) = 0 XLOOPSIGN(I) = 0 ENDDO CUTFOUND = .FALSE. C C If we should start generating cuts anew, we initialize. Else C we do some checking (just in case). C IF (NEWCUTINIT) THEN CUTINDEX(1) = CUTMAX - 2 DO I = 2, CUTMAX CUTINDEX(I) = CUTMAX - I ENDDO DO P = 0,NPROPS DO I = 1,2 VRTX(P,I) = XVRTX(P,I) ENDDO ENDDO ELSE DO P = 0,NPROPS DO I = 1,2 IF (.NOT.(VRTX(P,I).EQ.XVRTX(P,I))) THEN WRITE(NOUT,*)'SNAFU in NEWCUT. VRTX changed.' STOP ENDIF ENDDO ENDDO ENDIF C NEWCUTINIT = .FALSE. C C Initialization complete. C C Loop to find new valid cut. C VALIDCUT = .FALSE. DO WHILE (.NOT.VALIDCUT) C C Generate a new choice of the CUTINDEX(I). C We choose CUTINDEX(1),CUTINDEX(2),...,CUTINDEX(CUTMAX) C with 0 .LE. CUTINDEX(I) .LE. NPROPS and CUTINDEX(I) > CUTINDEX(I+1). C CUTINDEX(CUTMAX) = 0 indicates that there is a virtual loop. C Example, for NLOOPS = 3: CUTINDEX(I) is initialized to (2,2,1,0). C From this, we generate first CUTINDEX(I) = (3,2,1,0). On successive C runs through this code, we generate (4,2,1,0), (5,2,1,0), (6,2,1,0), C (7,2,1,0), (8,2,1,0), (4,3,1,0), (5,3,1,0),..., (8,3,1,0), (5,4,1,0), C (6,4,1,0),..., (8,7,6,0), (4,3,2,1), (5,3,2,1),..., (8,7,6,5). C In case the choice previously analyzed was the last one, then C there are no more cuts, we set CUTFOUND to False to abort further C analysis in the main program, and set NEWCUTINIT to True so that we C start over again next time this subroutine is called, and then return. C I = 1 77 IF( CUTINDEX(I).LT.NPROPS + 1 - I) THEN CUTINDEX(I) = CUTINDEX(I) + 1 DO II = 1,I - 1 CUTINDEX(II) = CUTINDEX(I) + I - II ENDDO ELSE I = I + 1 IF (I.LE.CUTMAX) THEN GO TO 77 ELSE CUTFOUND = .FALSE. NEWCUTINIT = .TRUE. RETURN ENDIF ENDIF C C Now that we have the CUTINDEX(I), set CUT(CUTINDEX(I)) = True. C CUTINDEX(CUTMAX) = 0 indicates that CUTMAX - 1 propagators were cut. C DO P= 1,NPROPS CUT(P) = .FALSE. ENDDO DO I = 1, CUTMAX - 1 CUT(CUTINDEX(I)) = .TRUE. ENDDO IF (CUTINDEX(CUTMAX).EQ.0) THEN NCUT = CUTMAX - 1 ELSE NCUT = CUTMAX CUT(CUTINDEX(CUTMAX)) = .TRUE. ENDIF C C Construct Left and Right sets. Any vertex that is connected to a Left C vertex by an uncut propagator is in Left. Similarly for Right. C Start with vertex 1 in Left and vertex 2 in Right. C DO V = 1,NVERTS LEFT(V) = .FALSE. RIGHT(V) = .FALSE. ENDDO LEFT(1) = .TRUE. RIGHT(2) = .TRUE. NLEFT = 1 NRIGHT = 1 C C Now add vertices that are connected to the Left vertices. C CHANGE = .TRUE. DO WHILE (CHANGE) CHANGE = .FALSE. DO P = 1,NPROPS IF(.NOT.CUT(P)) THEN C IF(LEFT(VRTX(P,1)).AND.(.NOT.LEFT(VRTX(P,2)))) THEN LEFT(VRTX(P,2)) = .TRUE. CHANGE = .TRUE. NLEFT = NLEFT + 1 ELSEIF(LEFT(VRTX(P,2)).AND.(.NOT.LEFT(VRTX(P,1)))) THEN LEFT(VRTX(P,1)) = .TRUE. CHANGE = .TRUE. NLEFT = NLEFT + 1 ENDIF C ENDIF ENDDO ENDDO C C Now add vertices that are connected to the Right vertices. C CHANGE = .TRUE. DO WHILE (CHANGE) CHANGE = .FALSE. DO P = 1,NPROPS IF(.NOT.CUT(P)) THEN C IF(RIGHT(VRTX(P,1)).AND.(.NOT.RIGHT(VRTX(P,2)))) THEN RIGHT(VRTX(P,2)) = .TRUE. CHANGE = .TRUE. NRIGHT = NRIGHT + 1 ELSEIF(RIGHT(VRTX(P,2)).AND.(.NOT.RIGHT(VRTX(P,1)))) THEN RIGHT(VRTX(P,1)) = .TRUE. CHANGE = .TRUE. NRIGHT = NRIGHT + 1 ENDIF C ENDIF ENDDO ENDDO C C Check for validity of the cut. Cut is not valid unless each vertex is C in Left or Right but not both. C VALIDCUT = .TRUE. DO V = 1,NVERTS IF (.NOT.(LEFT(V).XOR.RIGHT(V))) THEN VALIDCUT = .FALSE. ENDIF ENDDO C C Check that each cut propagator divides the Left set from the Right set. C DO I = 1,NCUT IF(LEFT(VRTX(CUTINDEX(I),1)).AND. > RIGHT(VRTX(CUTINDEX(I),2)) ) THEN CUTSIGN(I) = 1 ELSE IF(LEFT(VRTX(CUTINDEX(I),2)).AND. > RIGHT(VRTX(CUTINDEX(I),1)) ) THEN CUTSIGN(I) = -1 ELSE VALIDCUT = .FALSE. ENDIF ENDDO C C End of loop to generate a new cut C ENDDO C C Are there virtual loops? If not, just return (GOTO 1), C IF (NCUT.EQ.CUTMAX) THEN LEFTLOOP = .FALSE. RIGHTLOOP = .FALSE. NINLOOP = 0 DO I = 1,CUTMAX LOOPINDEX(I) = 0 LOOPSIGN(I) = 0 ENDDO GOTO 1 ELSE IF (NCUT.EQ.(CUTMAX - 1)) THEN IF ((NLEFT - NRIGHT).EQ.2) THEN LEFTLOOP = .TRUE. RIGHTLOOP = .FALSE. ELSEIF ((NRIGHT - NLEFT).EQ.2) THEN RIGHTLOOP = .TRUE. LEFTLOOP = .FALSE. ELSE WRITE(NOUT,*)'NRIGHT,NLEFT out of bounds' STOP ENDIF ELSE WRITE(NOUT,*)'NCUT out of bounds' STOP ENDIF C C Find the virtual loops. C C First, initialize all vertices to the left of the cut to be candidate C vertices for the left-loop. Alternatively, if we have a right-loop, C initialize all vertices to the right of the cut to be candidate C vertices for the right-loop C IF (LEFTLOOP) THEN DO V = 1,NVERTS IF (LEFT(V)) THEN LOOPVERTEX(V) = .TRUE. ELSE LOOPVERTEX(V) = .FALSE. ENDIF ENDDO ELSEIF (RIGHTLOOP) THEN DO V = 1,NVERTS IF (RIGHT(V)) THEN LOOPVERTEX(V) = .TRUE. ELSE LOOPVERTEX(V) = .FALSE. ENDIF ENDDO ENDIF C C Now we will iteratively remove candidate vertices for the loop if they C are not properly connected. C CHANGE = .TRUE. DO WHILE (CHANGE) CHANGE = .FALSE. C C A loop propagator is one that joins two loop vertices. C DO P = 1,NPROPS IF(LOOPVERTEX(VRTX(P,1)).AND.LOOPVERTEX(VRTX(P,2)))THEN LOOPPROP(P) = .TRUE. ELSE LOOPPROP(P) = .FALSE. ENDIF ENDDO C C Now a loop vertex is one that connects at least two loop propagators. C If a vertex fails this test, remove it from the loop list. C DO V = 1,NVERTS IF (LOOPVERTEX(V)) THEN C NCONNECTED = 0 DO P = 1,NPROPS IF( LOOPPROP(P).AND. > ((VRTX(P,1).EQ.V).OR.(VRTX(P,2).EQ.V)) ) THEN NCONNECTED = NCONNECTED + 1 ENDIF ENDDO C IF(NCONNECTED.LT.2) THEN LOOPVERTEX(V) = .FALSE. CHANGE = .TRUE. ENDIF ENDIF ENDDO C C Close loop over removal of loop vertices C ENDDO C C We now know which propagators are in the loop. Next we need to C find the starting vertex in the loop. The STARTVERTEX is either C vertex 1 or vertex 2 (connected to the external lines) or it is C a LOOPVERTEX connected connected to a propagator that is not in C the loop and not cut. We take the first one that we find. C IF (LOOPVERTEX(1)) THEN STARTVERTEX = 1 ELSE IF (LOOPVERTEX(2)) THEN STARTVERTEX = 2 ELSE P = 1 LOOKMORE = .TRUE. DO WHILE (LOOKMORE) IF (P.GT.NPROPS) THEN WRITE(NOUT,*) 'SNAFU 1 in NEWCUT, P too big' STOP ENDIF IF ( (.NOT.LOOPPROP(P)).AND.(.NOT.CUT(P)) ) THEN IF ( LOOPVERTEX(VRTX(P,1)) ) THEN STARTVERTEX = VRTX(P,1) LOOKMORE = .FALSE. ELSE IF ( LOOPVERTEX(VRTX(P,2)) ) THEN STARTVERTEX = VRTX(P,2) LOOKMORE = .FALSE. ENDIF ENDIF P = P + 1 ENDDO ENDIF C C Now add first propagator in the loop. C P = 1 LOOKMORE = .TRUE. DO WHILE (LOOKMORE) IF (P.GT.NPROPS) THEN WRITE(NOUT,*) 'SNAFU 2 in NEWCUT, P too big' STOP ENDIF IF(LOOPPROP(P)) THEN IF((VRTX(P,1).EQ.STARTVERTEX)) THEN IP = 1 LOOPINDEX(IP) = P PREVIOUSPROP = P LOOPSIGN(IP) = 1 HOTVERTEX = VRTX(P,2) LOOKMORE = .FALSE. ELSEIF((VRTX(P,2).EQ.STARTVERTEX)) THEN IP = 1 LOOPINDEX(IP) = P PREVIOUSPROP = P LOOPSIGN(IP) = -1 HOTVERTEX = VRTX(P,1) LOOKMORE = .FALSE. ENDIF ENDIF P = P+1 ENDDO C C Now add propagators around the loop. C DO WHILE (HOTVERTEX.NE.STARTVERTEX) P = 1 LOOKMORE = .TRUE. DO WHILE (LOOKMORE) IF (P.GT.NPROPS) THEN WRITE(NOUT,*) 'SNAFU 3 in NEWCUT, P too big' STOP ENDIF IF(LOOPPROP(P).AND.(.NOT.(PREVIOUSPROP.EQ.P))) THEN IF((VRTX(P,1).EQ.HOTVERTEX)) THEN IP = IP + 1 LOOPINDEX(IP) = P PREVIOUSPROP = P LOOPSIGN(IP) = 1 HOTVERTEX = VRTX(P,2) LOOKMORE = .FALSE. ELSEIF((VRTX(P,2).EQ.HOTVERTEX)) THEN IP = IP + 1 LOOPINDEX(IP) = P PREVIOUSPROP = P LOOPSIGN(IP) = -1 HOTVERTEX = VRTX(P,1) LOOKMORE = .FALSE. ENDIF ENDIF P = P + 1 ENDDO ENDDO NINLOOP = IP C C Come to here for a normal return. C 1 CUTFOUND = .TRUE. XLEFTLOOP = LEFTLOOP XRIGHTLOOP = RIGHTLOOP XNINLOOP = NINLOOP XNCUT = NCUT C DO I = 1,CUTMAX XCUTINDEX(I) = CUTINDEX(I) XCUTSIGN(I) = CUTSIGN(I) ENDDO C DO I = 1,NINLOOP XLOOPINDEX(I) = LOOPINDEX(I) XLOOPSIGN(I) = LOOPSIGN(I) ENDDO C DO P = 1,NPROPS V1 = VRTX(P,1) V2 = VRTX(P,2) IF (LEFT(V1).AND.LEFT(V2)) THEN XISIGN(P) = 1 ELSE IF (RIGHT(V1).AND.RIGHT(V2)) THEN XISIGN(P) = -1 ELSE XISIGN(P) = 0 ENDIF ENDDO C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C End of subroutines associated with NEWCUT 2 C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE DEFORM(VRTX,LOOPINDEX,RTS,LEFTLOOP,RIGHTLOOP, > NINLOOP,KINLOOP,NEWKINLOOP,JACDEFORM) C INTEGER SIZE PARAMETER (SIZE = 3) C In: INTEGER VRTX(0:3*SIZE-1,2) INTEGER LOOPINDEX(SIZE+1) REAL*8 RTS LOGICAL LEFTLOOP,RIGHTLOOP INTEGER NINLOOP REAL*8 KINLOOP(SIZE+1,0:3) C Out: COMPLEX*16 NEWKINLOOP(0:3) COMPLEX*16 JACDEFORM C C Contour deformation. Note that this simple algorithm should C work for NINLOOP = 3 and for NINLOOP = 4 if the sum of the C three 3-momenta exiting the loop vanishes. C C In variables: C VRTX(P,I) = Index of vertex at beginning (i= 1) and end (I = 2) of C of propagator P. Specifies the supergraph. C LOOPINDEX(NP) = Index P of NPth propagator around the loop. C RTS = energy of final state. C LEFTLOOP = T if there is a loop to the left of the cut. C RIGHTLOOP = T if there is a loop to the right of the cut. C NINLOOP = number of propagators in the loop. C KINLOOP(J,MU) = momentum of Jth propagator in loop (real part) C Out variables: C NEWKINLOOP(MU) = added part of loop momentum. C (purely imaginary for MU = 1,2,3) C JACDEFORM = jacobian associated with contour deformation. C C Our notation is C vec Q(j) = vec L(j) - vec L(j+1) j = 1,...,Ninloop - 1 C L(j) = |L(j)| j = 1,...,Ninloop - 1 C Q(j) = |Q(j)| j = 1,...,Ninloop - 1 C C 27 October 1992 first DEFORM C 1 February 1998 latest version C 23 February 1998 minor revision to rename deform variables C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX REAL*8 DEFORMALPHA,DEFORMBETA,DEFORMGAMMA COMMON /DEFORMSCALES/DEFORMALPHA,DEFORMBETA,DEFORMGAMMA C REAL*8 S INTEGER SIGN REAL*8 L(SIZE+1,3) C REAL*8 Q(SIZE,3),QSQ(SIZE),QABS(SIZE) REAL*8 LHAT(SIZE,3),LSQ(SIZE),LABS(SIZE) REAL*8 W(SIZE,3),WSQ(SIZE),WABS(SIZE) REAL*8 ACRIT2,ACRIT3,A2,A3 REAL*8 DELTA REAL*8 M1(3,3),M2(3,3),M3(3,3) REAL*8 D1,D2,D3,DSQ,GRADDSQ(3) REAL*8 FRACTION,GRADF(3) REAL*8 G2,G3,DG2DA2,DG3DA3 REAL*8 C,DLNCDDSQ REAL*8 TERMC,TERMF,TERMG2,TERMG3,TERMW2,TERMW3,TERMS COMPLEX*16 A(3,3) C LOGICAL CONNECTSTOCURRENT REAL*8 TEMP,TEMP1,TEMP2,TEMP3 INTEGER J,MU,NU C C Calculate s. C S = RTS**2 C C Initialize with default value. C DO MU = 0,3 NEWKINLOOP(MU) = (0.0D0,0.0D0) ENDDO C JACDEFORM = (1.0D0,0.0D0) C C Check to see if we should actually do anything C IF (NINLOOP.LT.2) THEN RETURN ENDIF C C Set C SIGN = +1 and L(J,MU) = KINLOOP(J,MU) for a left loop, C SIGN = -1 and L(J,MU) = KINLOOP(NINLOOP-J+1,MU) for a right loop. C IF (LEFTLOOP) THEN SIGN = + 1 DO J = 1,NINLOOP DO MU = 1,3 L(J,MU) = KINLOOP(J,MU) ENDDO ENDDO ELSE IF (RIGHTLOOP) THEN SIGN = - 1 DO J = 1,NINLOOP DO MU = 1,3 L(J,MU) = KINLOOP(NINLOOP-J+1,MU) ENDDO ENDDO ELSE WRITE(NOUT,*) 'Snafu in DEFORM' STOP ENDIF C C Two particles in the loop. C IF (NINLOOP.EQ.2) THEN C C Calculate Q(3,mu), |L(j)|^2, |Q(3)|^2, |L(j)|, |Q(3)|, hat L(j,mu). C QSQ(3) = 0.0D0 DO MU = 1,3 Q(3,MU) = L(1,MU) - L(2,MU) QSQ(3) = QSQ(3) + Q(3,MU)**2 ENDDO QABS(3) = SQRT(QSQ(3)) DO J = 1,2 LSQ(J) = 0.0D0 DO MU = 1,3 LSQ(J) = LSQ(J) + L(J,MU)**2 ENDDO LABS(J) = SQRT(LSQ(J)) DO MU = 1,3 LHAT(J,MU) = L(J,MU)/LABS(J) ENDDO ENDDO C C Calculate the vector W(3,mu), along with the corresponding C normalization factor. C WSQ(3) = 0.0D0 DO MU = 1,3 W(3,MU) = LHAT(1,MU) + LHAT(2,MU) WSQ(3) = WSQ(3) + W(3,MU)**2 ENDDO WABS(3) = SQRT(WSQ(3)) C C The size of the critical ellipse. C ACRIT3 = RTS - 2.0D0*QABS(3) C C The size of the ellipse at point L. C A3 = LABS(1) + LABS(2) - QABS(3) C C Calculate the matrix d W(3,mu)/ d L(mu) = M3(mu,nu). C DO MU = 1,3 DO NU = 1,3 TEMP = DELTA(MU,NU) TEMP1 = (TEMP - LHAT(1,MU)*LHAT(1,NU))/LABS(1) TEMP2 = (TEMP - LHAT(2,MU)*LHAT(2,NU))/LABS(2) M3(MU,NU) = TEMP1 + TEMP2 ENDDO ENDDO C C The "distance" to the collinear line. C D3 = LABS(1)*LABS(2)*WABS(3)/QABS(3) C C The square of the distance its gradient. C DSQ = D3**2 DO NU = 1,3 TEMP = 0.0D0 DO MU = 1,3 TEMP = TEMP + W(3,MU)*M3(MU,NU) ENDDO TEMP = TEMP/WSQ(3) TEMP = TEMP + L(1,NU)/LSQ(1) + L(2,NU)/LSQ(2) GRADDSQ(NU) = 2.0D0*DSQ*TEMP ENDDO C C The function G3 and its derivative. C G3 = 1.0D0/(ACRIT3 + DEFORMGAMMA*A3) DG3DA3 = - DEFORMGAMMA/(ACRIT3 + DEFORMGAMMA*A3)**2 C C Calculate the function C(DSQ) and its derivative. Note that we change C the sign of C in the case of a loop to the right of the cut. C C We effectively make DEFORMALPHA smaller by a factor 10 for the two C point function so as to avoid crossing branch cut of SQRT C C = SIGN*DEFORMALPHA*DSQ /(1.0D0 + DEFORMBETA*DSQ/QSQ(3)) C = C * ACRIT3/RTS DLNCDDSQ = 1.0D0/DSQ/(1.0D0 + DEFORMBETA*DSQ/QSQ(3)) C C Calculate the imaginary part of the loop momentum L(mu). C DO MU = 1,3 NEWKINLOOP(MU) = (0.0D0,-1.0D0) * C * G3 * W(3,MU) ENDDO C C Calculate the jacobian. C First, we need the comlex matrix A(mu,nu), the derivative C of ComplexL(mu) with respecdt to L(nu). C DO MU = 1,3 DO NU = 1,3 TERMC = G3*W(3,MU)*DLNCDDSQ*GRADDSQ(NU) TERMG3 = DG3DA3*W(3,MU)*W(3,NU) TERMW3 = G3*M3(MU,NU) TERMS = TERMC + TERMG3 + TERMW3 A(MU,NU) = DELTA(MU,NU) + (0.0D0,-1.0D0) * C * TERMS ENDDO ENDDO C C Finally, the jacobian is the determinant of A C JACDEFORM = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) > + A(1,2) * ( A(2,3) * A(3,1) - A(2,1) * A(3,3) ) > + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) C C End of Ninloop = 2 calculation C C Three particles in the loop. C ELSE IF (NINLOOP.EQ.3) THEN C C First we need to determine if our loop connects to the current vertex. C IF (LEFTLOOP) THEN IF ((VRTX(LOOPINDEX(1),1).EQ.1)) THEN CONNECTSTOCURRENT = .TRUE. ELSE CONNECTSTOCURRENT = .FALSE. ENDIF ELSE IF ((VRTX(LOOPINDEX(1),1).EQ.2)) THEN CONNECTSTOCURRENT = .TRUE. ELSE CONNECTSTOCURRENT = .FALSE. ENDIF ENDIF IF (CONNECTSTOCURRENT) THEN C C Calculation for a three particle loop that connects to the current. C C Calculate Q(3,mu), |L(j)|^2, |Q(3)|^2, |L(j)|, |Q(3)|, hat L(j,mu). C QSQ(3) = 0.0D0 DO MU = 1,3 Q(3,MU) = L(1,MU) - L(2,MU) QSQ(3) = QSQ(3) + Q(3,MU)**2 ENDDO QABS(3) = SQRT(QSQ(3)) DO J = 1,2 LSQ(J) = 0.0D0 DO MU = 1,3 LSQ(J) = LSQ(J) + L(J,MU)**2 ENDDO LABS(J) = SQRT(LSQ(J)) DO MU = 1,3 LHAT(J,MU) = L(J,MU)/LABS(J) ENDDO ENDDO C C Calculate the vector W(3,mu), along with the corresponding C normalization factor. C WSQ(3) = 0.0D0 DO MU = 1,3 W(3,MU) = LHAT(1,MU) + LHAT(2,MU) WSQ(3) = WSQ(3) + W(3,MU)**2 ENDDO WABS(3) = SQRT(WSQ(3)) C C The size of the critical ellipse. C ACRIT3 = RTS - 2.0D0*QABS(3) C C The size of the ellipse at point L. C A3 = LABS(1) + LABS(2) - QABS(3) C C Calculate the matrix d W(3,mu)/ d L(mu) = M3(mu,nu). C DO MU = 1,3 DO NU = 1,3 TEMP = DELTA(MU,NU) TEMP1 = (TEMP - LHAT(1,MU)*LHAT(1,NU))/LABS(1) TEMP2 = (TEMP - LHAT(2,MU)*LHAT(2,NU))/LABS(2) M3(MU,NU) = TEMP1 + TEMP2 ENDDO ENDDO C C The "distance" to the collinear line. C D3 = LABS(1)*LABS(2)*WABS(3)/QABS(3) C C The square of the distance its gradient. C DSQ = D3**2 DO NU = 1,3 TEMP = 0.0D0 DO MU = 1,3 TEMP = TEMP + W(3,MU)*M3(MU,NU) ENDDO TEMP = TEMP/WSQ(3) TEMP = TEMP + L(1,NU)/LSQ(1) + L(2,NU)/LSQ(2) GRADDSQ(NU) = 2.0D0*DSQ*TEMP ENDDO C C The function G3 and its derivative. C G3 = 1.0D0/(ACRIT3 + DEFORMGAMMA*A3) DG3DA3 = - DEFORMGAMMA/(ACRIT3 + DEFORMGAMMA*A3)**2 C C Calculate the function C(DSQ) and its derivative. Note that we change C the sign of C in the case of a loop to the right of the cut. C C = SIGN*DEFORMALPHA*DSQ /(1.0D0 + DEFORMBETA*DSQ/QSQ(3)) DLNCDDSQ = 1.0D0/DSQ/(1.0D0 + DEFORMBETA*DSQ/QSQ(3)) C C Calculate the imaginary part of the loop momentum L(mu). C DO MU = 1,3 NEWKINLOOP(MU) = (0.0D0,-1.0D0) * C * G3 * W(3,MU) ENDDO C C Calculate the jacobian. C First, we need the comlex matrix A(mu,nu), the derivative C of ComplexL(mu) with respecdt to L(nu). C DO MU = 1,3 DO NU = 1,3 TERMC = G3*W(3,MU)*DLNCDDSQ*GRADDSQ(NU) TERMG3 = DG3DA3*W(3,MU)*W(3,NU) TERMW3 = G3*M3(MU,NU) TERMS = TERMC + TERMG3 + TERMW3 A(MU,NU) = DELTA(MU,NU) + (0.0D0,-1.0D0) * C * TERMS ENDDO ENDDO C C Finally, the jacobian is the determinant of A C JACDEFORM = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) > + A(1,2) * ( A(2,3) * A(3,1) - A(2,1) * A(3,3) ) > + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) C C End of calculation for Ninloop = 3 for a loop connecting to the C current ( IF(CONNECTSTOCURRENT) ). C ELSE C C Calculation for a three particle loop that does not connect C to the current. C C Calculate Q(j,mu), |L(j)|^2, |Q(j)|^2, |L(j)|, |Q(j)|, hat L(j,mu). C DO MU = 1,3 Q(1,MU) = L(2,MU) - L(3,MU) Q(2,MU) = L(3,MU) - L(1,MU) Q(3,MU) = L(1,MU) - L(2,MU) ENDDO DO J = 1,3 LSQ(J) = 0.0D0 QSQ(J) = 0.0D0 DO MU = 1,3 LSQ(J) = LSQ(J) + L(J,MU)**2 QSQ(J) = QSQ(J) + Q(J,MU)**2 ENDDO QABS(J) = SQRT(QSQ(J)) LABS(J) = SQRT(LSQ(J)) DO MU = 1,3 LHAT(J,MU) = L(J,MU)/LABS(J) ENDDO ENDDO C C The vectors W(j,mu) and their squares and their absolute values. C DO MU = 1,3 W(1,MU) = LHAT(2,MU) + LHAT(3,MU) W(2,MU) = LHAT(3,MU) + LHAT(1,MU) W(3,MU) = LHAT(1,MU) + LHAT(2,MU) ENDDO DO J = 1,3 WSQ(J) = 0.0D0 DO MU = 1,3 WSQ(J) = WSQ(J) + W(J,MU)**2 ENDDO WABS(J) = SQRT(WSQ(J)) ENDDO C C The size of the critical ellipses. C ACRIT2 = RTS - 2.0D0*QABS(2) ACRIT3 = RTS - 2.0D0*QABS(3) C C The sizes of the ellipses at point L. C A2 = LABS(3) + LABS(1) - QABS(2) A3 = LABS(1) + LABS(2) - QABS(3) C C Calculate the matrix d W(j,mu)/ d L(mu) = Mj(mu,nu). C DO MU = 1,3 DO NU = 1,3 TEMP = DELTA(MU,NU) TEMP1 = (TEMP - LHAT(1,MU)*LHAT(1,NU))/LABS(1) TEMP2 = (TEMP - LHAT(2,MU)*LHAT(2,NU))/LABS(2) TEMP3 = (TEMP - LHAT(3,MU)*LHAT(3,NU))/LABS(3) M1(MU,NU) = TEMP2 + TEMP3 M2(MU,NU) = TEMP3 + TEMP1 M3(MU,NU) = TEMP1 + TEMP2 ENDDO ENDDO C C The "distances" to the collinear lines. In this case we do not need D2. C D1 = LABS(2)*LABS(3)*WABS(1)/QABS(1) D3 = LABS(1)*LABS(2)*WABS(3)/QABS(3) C C The square of the smaller of D1 and D3 and its gradient. C IF (D1.LT.D3) THEN C C D1 is the smaller distance. C DSQ = D1**2 DO NU = 1,3 TEMP = 0.0D0 DO MU = 1,3 TEMP = TEMP + W(1,MU)*M1(MU,NU) ENDDO TEMP = TEMP/WSQ(1) TEMP = TEMP + L(2,NU)/LSQ(2) + L(3,NU)/LSQ(3) GRADDSQ(NU) = 2.0D0*DSQ*TEMP ENDDO C ELSE C C D3 is the smaller distance. C DSQ = D3**2 DO NU = 1,3 TEMP = 0.0D0 DO MU = 1,3 TEMP = TEMP + W(3,MU)*M3(MU,NU) ENDDO TEMP = TEMP/WSQ(3) TEMP = TEMP + L(1,NU)/LSQ(1) + L(2,NU)/LSQ(2) GRADDSQ(NU) = 2.0D0*DSQ*TEMP ENDDO C ENDIF C C The function G2 and its derivative. C G2 = 1.0D0/(ACRIT2 + DEFORMGAMMA*A2) DG2DA2 = - DEFORMGAMMA/(ACRIT2 + DEFORMGAMMA*A2)**2 C C Calculate the function C(DSQ) and its derivative. Note that we change C the sign of C in the case of a loop to the right of the cut. C C = SIGN*DEFORMALPHA*DSQ /(1.0D0 + 4.0D0*DEFORMBETA*DSQ/S) DLNCDDSQ = 1.0D0/DSQ/(1.0D0 + 4.0D0*DEFORMBETA*DSQ/S) C C Calculate the imaginary part of the loop momentum L(mu). C DO MU = 1,3 NEWKINLOOP(MU) = (0.0D0,-1.0D0) * C * G2 * W(2,MU) ENDDO C C Calculate the jacobian. C First, we need the comlex matrix A(mu,nu), the derivative C of ComplexL(mu) with respecdt to L(nu). C DO MU = 1,3 DO NU = 1,3 TERMC = G2*W(2,MU)*DLNCDDSQ*GRADDSQ(NU) TERMG2 = DG2DA2*W(2,MU)*W(2,NU) TERMW2 = G2*M2(MU,NU) TERMS = TERMC + TERMG2 + TERMW2 A(MU,NU) = DELTA(MU,NU) + (0.0D0,-1.0D0) * C * TERMS ENDDO ENDDO C C Finally, the jacobian is the determinant of A C JACDEFORM = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) > + A(1,2) * ( A(2,3) * A(3,1) - A(2,1) * A(3,3) ) > + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) C C End of calculation for Ninloop = 3 for a loop not connecting to the C current ( IF(CONNECTSTOCURRENT) ... ELSE ...). C ENDIF C C End of Ninloop = 3 calculation C C Four particles in the loop.------------ C ELSE IF (NINLOOP.EQ.4) THEN C C Calculate Q(j,mu), |L(j)|^2, |Q(j)|^2, |L(j)|, |Q(j)|, hat L(j,mu). C DO MU = 1,3 Q(1,MU) = L(2,MU) - L(3,MU) Q(2,MU) = L(3,MU) - L(1,MU) Q(3,MU) = L(1,MU) - L(2,MU) ENDDO DO J = 1,3 LSQ(J) = 0.0D0 QSQ(J) = 0.0D0 DO MU = 1,3 LSQ(J) = LSQ(J) + L(J,MU)**2 QSQ(J) = QSQ(J) + Q(J,MU)**2 ENDDO QABS(J) = SQRT(QSQ(J)) LABS(J) = SQRT(LSQ(J)) DO MU = 1,3 LHAT(J,MU) = L(J,MU)/LABS(J) ENDDO ENDDO C C The vectors W(j,mu) and their squares and their absolute values. C DO MU = 1,3 W(1,MU) = LHAT(2,MU) + LHAT(3,MU) W(2,MU) = LHAT(3,MU) + LHAT(1,MU) W(3,MU) = LHAT(1,MU) + LHAT(2,MU) ENDDO DO J = 1,3 WSQ(J) = 0.0D0 DO MU = 1,3 WSQ(J) = WSQ(J) + W(J,MU)**2 ENDDO WABS(J) = SQRT(WSQ(J)) ENDDO C C The size of the critical ellipses. C ACRIT2 = RTS - 2.0D0*QABS(2) ACRIT3 = RTS - 2.0D0*QABS(3) C C The sizes of the ellipses at point L. C A2 = LABS(3) + LABS(1) - QABS(2) A3 = LABS(1) + LABS(2) - QABS(3) C C Calculate the matrix d W(j,mu)/ d L(mu) = Mj(mu,nu). C DO MU = 1,3 DO NU = 1,3 TEMP = DELTA(MU,NU) TEMP1 = (TEMP - LHAT(1,MU)*LHAT(1,NU))/LABS(1) TEMP2 = (TEMP - LHAT(2,MU)*LHAT(2,NU))/LABS(2) TEMP3 = (TEMP - LHAT(3,MU)*LHAT(3,NU))/LABS(3) M1(MU,NU) = TEMP2 + TEMP3 M2(MU,NU) = TEMP3 + TEMP1 M3(MU,NU) = TEMP1 + TEMP2 ENDDO ENDDO C C The "distances" to the collinear lines. C D1 = LABS(2)*LABS(3)*WABS(1)/QABS(1) D2 = LABS(3)*LABS(1)*WABS(2)/QABS(2) D3 = LABS(1)*LABS(2)*WABS(3)/QABS(3) C C The square of the smallest distance its gradient. C IF ((D1.LT.D2).AND.(D1.LT.D3)) THEN C C D1 is the smallest distance C DSQ = D1**2 DO NU = 1,3 TEMP = 0.0D0 DO MU = 1,3 TEMP = TEMP + W(1,MU)*M1(MU,NU) ENDDO TEMP = TEMP/WSQ(1) TEMP = TEMP + L(2,NU)/LSQ(2) + L(3,NU)/LSQ(3) GRADDSQ(NU) = 2.0D0*DSQ*TEMP ENDDO C ELSE IF (D2.LT.D3) THEN C C D2 is the smallest distance C DSQ = D2**2 DO NU = 1,3 TEMP = 0.0D0 DO MU = 1,3 TEMP = TEMP + W(2,MU)*M2(MU,NU) ENDDO TEMP = TEMP/WSQ(2) TEMP = TEMP + L(3,NU)/LSQ(3) + L(1,NU)/LSQ(1) GRADDSQ(NU) = 2.0D0*DSQ*TEMP ENDDO C ELSE C C D3 is the smallest distance. C DSQ = D3**2 DO NU = 1,3 TEMP = 0.0D0 DO MU = 1,3 TEMP = TEMP + W(3,MU)*M3(MU,NU) ENDDO TEMP = TEMP/WSQ(3) TEMP = TEMP + L(1,NU)/LSQ(1) + L(2,NU)/LSQ(2) GRADDSQ(NU) = 2.0D0*DSQ*TEMP ENDDO C ENDIF C C The mixing fraction FRACTION and its gradient. C FRACTION = LSQ(3) /(LSQ(2) + LSQ(3)) DO NU = 1,3 TEMP = LSQ(2)*L(3,NU) - LSQ(3)*L(2,NU) GRADF(NU) = 2.0D0*TEMP/(LSQ(2) + LSQ(3))**2 ENDDO C C The functions G2 and G3 and their derivatives. C G2 = 1.0D0/(ACRIT2 + DEFORMGAMMA*A2) G3 = 1.0D0/(ACRIT3 + DEFORMGAMMA*A3) DG2DA2 = - DEFORMGAMMA/(ACRIT2 + DEFORMGAMMA*A2)**2 DG3DA3 = - DEFORMGAMMA/(ACRIT3 + DEFORMGAMMA*A3)**2 C C Calculate the function C(DSQ) and its derivative. Note that we change C the sign of C in the case of a loop to the right of the cut. C C = SIGN*DEFORMALPHA*DSQ /(1.0D0 + 4.0D0*DEFORMBETA*DSQ/S) DLNCDDSQ = 1.0D0/DSQ/(1.0D0 + 4.0D0*DEFORMBETA*DSQ/S) C C Calculate the imaginary part of the loop momentum L(mu). C DO MU = 1,3 NEWKINLOOP(MU) = (0.0D0,-1.0D0) * C > * (FRACTION*G2*W(2,MU) + (1.0D0 - FRACTION)*G3*W(3,MU)) ENDDO C C Calculate the jacobian. C First, we need the comlex matrix A(mu,nu), the derivative C of ComplexL(mu) with respecdt to L(nu). C DO MU = 1,3 DO NU = 1,3 TERMC = FRACTION*G2*W(2,MU) + (1.0D0 - FRACTION)*G3*W(3,MU) TERMC = TERMC * DLNCDDSQ*GRADDSQ(NU) TERMF = ( G2*W(2,MU) - G3*W(3,MU) )*GRADF(NU) TERMG2 = FRACTION *DG2DA2*W(2,MU)*W(2,NU) TERMG3 = (1.0D0 - FRACTION)*DG3DA3*W(3,MU)*W(3,NU) TERMW2 = FRACTION *G2*M2(MU,NU) TERMW3 = (1.0D0 - FRACTION)*G3*M3(MU,NU) TERMS = TERMC + TERMF + TERMG2 + TERMG3 + TERMW2 + TERMW3 A(MU,NU) = DELTA(MU,NU) + (0.0D0,-1.0D0) * C * TERMS ENDDO ENDDO C C Finally, the jacobian is the determinant of A C JACDEFORM = A(1,1) * ( A(2,2) * A(3,3) - A(2,3) * A(3,2) ) > + A(1,2) * ( A(2,3) * A(3,1) - A(2,1) * A(3,3) ) > + A(1,3) * ( A(2,1) * A(3,2) - A(2,2) * A(3,1) ) C C End of Ninloop = 4 calculation C ELSE WRITE(NOUT,*) 'Not programed for NINLOOP > 4 yet.' STOP ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION DELTA(MU,NU) C C In: INTEGER MU,NU C C Kroneker delta. C IF (MU.EQ.NU) THEN DELTA = 1.0D0 ELSE DELTA = 0.0D0 ENDIF RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE CHECKDEFORM(KCSQ,LEFTLOOP,RIGHTLOOP,JPROBED,JCUT, > GRAPHNUMBER,CUT,LOOPCUT) C INTEGER SIZE PARAMETER (SIZE = 3) C In: COMPLEX*16 KCSQ LOGICAL LEFTLOOP,RIGHTLOOP INTEGER JPROBED,JCUT INTEGER GRAPHNUMBER LOGICAL CUT(3*SIZE-1),LOOPCUT(3*SIZE-1) C Out: None C C This subroutine just checks that the contour deformation was in C the right direction. Here JPROBED is the label 1,2,3... counting C around of the loop of the probagator whose squared momentum is C KSCQ, while JCUT is the label of the loopcut propagator. C Consider a loop to the left of the cut. C For propagators with JPROBED < JCUT + 1 ,we do not have a singularity C protected by an i\epsilon, so we do not test. For propagators with C JPROBED > JCUT + 1, we are deforming around a singularity that is C protected with an i\epsilon, and we should check that we deformed C in the right direction. C For loops to the right of the cut, the situation is the opposite. C C 14 February 1996 C 16 March 1996 C 24 October 1996 (Checks also for NINLOOP.eq.2, 21 September 1997) C 2 January 1998 Simplify, look at all cases. C 14 March 1998 C 4 August 1998 Check all cases for 3 and 4 point subgraphs. C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT C REAL*8 XXIMAG,XXREAL C-- C IF ( LEFTLOOP ) THEN C IF (JPROBED.GT.(JCUT+1)) THEN IF ( - XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign for deformation with an L loop?' WRITE(NOUT,*) KCSQ STOP ENDIF ELSEIF (JPROBED.EQ.(JCUT+1)) THEN C C Check the special cases in which there is a three particle loop C that contains one of the current vertices. Then there is one C propagator coming out of the loop that is not cut, but connects C to two cut propagators. When the loopcut propagator is the one just C before this one, we need to check the deformation sign. Since this C is just a check, we find the cases to look at by hand instead of C writing code to do the logic. C IF (GRAPHNUMBER.EQ.5) THEN IF (LOOPCUT(5)) THEN IF ( - XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 5, L loop?' WRITE(NOUT,*) KCSQ STOP ENDIF ENDIF ELSE IF (GRAPHNUMBER.EQ.8) THEN IF (LOOPCUT(1).AND.CUT(3)) THEN IF ( - XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 8, L loop?' WRITE(NOUT,*) KCSQ STOP ENDIF ELSEIF (LOOPCUT(5).AND.CUT(4)) THEN IF ( - XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 8, L loop?' WRITE(NOUT,*) KCSQ STOP ENDIF ENDIF ENDIF C Close IF (JPROBED.GT.(JCUT+1)) ... ELSEIF (JPROBED.EQ.(JCUT+1)) ... ENDIF C ELSEIF ( RIGHTLOOP ) THEN C IF (JPROBED.LT.(JCUT-1)) THEN IF ( XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign for deformation with an R loop?' WRITE(NOUT,*) KCSQ STOP ENDIF C C Special case checks as for left loop. C ELSE IF (JPROBED.EQ.(JCUT-1)) THEN C IF (GRAPHNUMBER.EQ.6) THEN IF (LOOPCUT(4)) THEN IF ( XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 6, R loop?' WRITE(NOUT,*) KCSQ STOP ENDIF ENDIF ELSE IF (GRAPHNUMBER.EQ.8) THEN IF (LOOPCUT(8).AND.CUT(1)) THEN IF ( XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 8, R loop?' WRITE(NOUT,*) KCSQ STOP ENDIF ELSE IF (LOOPCUT(4).AND.CUT(2)) THEN IF ( XXIMAG(KCSQ) .GT. ABS(XXREAL(KCSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 8, R loop?' WRITE(NOUT,*) KCSQ STOP ENDIF ENDIF ENDIF C Close IF (JPROBED.LT.(JCUT-1)) ... ELSEIF (JPROBED.EQ.(JCUT-1)) ... ENDIF C ELSE C WRITE(NOUT,*) 'If there is a loop, it must be L or R.' STOP C Close IF (LEFTLOOP) ... ELSE IF (RIGHTLOOP) ... ELSE ENDIF RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE CHECKDEFORM2(QBARSQ,CUT,P1,P2) C INTEGER SIZE PARAMETER (SIZE = 3) C In: COMPLEX*16 QBARSQ LOGICAL CUT(3*SIZE-1) INTEGER P1,P2 C Out: None C C This subroutine just checks that the contour deformation was in C the right direction in the case of a self-energy virtual correction C on a quark propagator with q^2 >0. Here QBARSQ is the squared C four-momentum in the dispersive integral. We want to check that the C sign of the imaginary part of QBARSQ is negative for a virtual loop C to the left of the final state cut and positive for a virtual loop C to the right of the final state cut. We use an ad-hoc method rather C than programming the logic to distinguish the cases since this is C just a check. C C 4 August 1998 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT C REAL*8 XXIMAG,XXREAL C-- C IF ( CUT(2).AND.CUT(3).AND.CUT(5) ) THEN C C This is graph 5, with a virtual self energy to the right of the cut. C IF ( - XXIMAG(QBARSQ) .GT. ABS(XXREAL(QBARSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 5, R loop?' WRITE(NOUT,*) QBARSQ STOP ENDIF C ELSE IF ( CUT(1).AND.CUT(6).AND.CUT(7) ) THEN C C This is graph 2, with a virtual self energy to the left of the cut. C IF ( XXIMAG(QBARSQ) .GT. ABS(XXREAL(QBARSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 2, L loop?' WRITE(NOUT,*) QBARSQ STOP ENDIF C ELSE IF ( CUT(1).AND.CUT(4).AND.CUT(5) ) THEN IF( (P1.EQ.2).OR.(P2.EQ.2) ) THEN C C This is graph 6, with a virtual self-energy to the left of the cut. C IF ( XXIMAG(QBARSQ) .GT. ABS(XXREAL(QBARSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 6, L loop?' WRITE(NOUT,*) QBARSQ STOP ENDIF C ELSE IF( (P1.EQ.3).OR.(P2.EQ.3) ) THEN C C This is graph 2, with a virtual self energy to the right of the cut. C IF ( - XXIMAG(QBARSQ) .GT. ABS(XXREAL(QBARSQ)) ) THEN WRITE(NOUT,*) 'Wrong sign deformation in graph 2, R loop?' WRITE(NOUT,*) QBARSQ STOP ENDIF C ELSE WRITE(NOUT,*) 'Snafu 2 in CHECKDEFORM2' STOP ENDIF C ELSE WRITE(NOUT,*) 'Snafu 1 in CHECKDEFORM2' STOP ENDIF RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION SMEAR(RTS) C REAL*8 RTS C C A smearing function that may do a good job of optimizing C the integration accuracy. It satisfies C C Int_0^\infty dE SMEAR(E) = 1 C C We take C C SMEAR(E) = (N-1)!/[M! (N-M-2)!] (A E_0 )**(N-M-1) C * E**M / [E + A * E_0]**N C C where E_0 = ENERGYSCALE. C C REAL*8 ENERGYSCALE COMMON /MSCALE/ ENERGYSCALE C REAL*8 SMEARFCTR INTEGER LOWPWR,HIGHPWR COMMON /SMEARPARMS/ SMEARFCTR,LOWPWR,HIGHPWR C REAL*8 FACTORIAL C SMEAR = FACTORIAL(HIGHPWR-1) > /(FACTORIAL(LOWPWR) * FACTORIAL(HIGHPWR-LOWPWR-2)) SMEAR = SMEAR * (SMEARFCTR * ENERGYSCALE)**(HIGHPWR-LOWPWR-1) SMEAR = SMEAR * RTS**LOWPWR > /( RTS + SMEARFCTR*ENERGYSCALE )**HIGHPWR RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C Numerator Function C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION NUMERATOR(GRAPHNUMBER,KC,CUT) C C In: INTEGER SIZE PARAMETER (SIZE = 3) INTEGER GRAPHNUMBER COMPLEX*16 KC(0:3*SIZE-1,0:3) LOGICAL CUT(3*SIZE-1) C C Numerator function for graph GRAPHNUMBER, with complex momenta KC, C including the color factor. C Early version: 17 July 1994 C This version written by Mathematica code of 4 September 1998 on C 4 Sep 1998 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C REAL*8 METRIC(0:3) DATA METRIC /+1.0D0,-1.0D0,-1.0D0,-1.0D0/ C INTEGER MU,NU C COMPLEX*16 K12,K12G5678,K13,K13G5678 COMPLEX*16 K14,K14G5678,K15,K16 COMPLEX*16 K17,K18,K1Q2854,K1Q3876 COMPLEX*16 K1Q4687,K1Q5687,K1Q6487,K1Q6587 COMPLEX*16 K1Q8254,K1Q8376,K22,K23 COMPLEX*16 K24,K24G5678,K25,K26 COMPLEX*16 K27,K28,K2Q2854,K2Q3876 COMPLEX*16 K2Q5687,K2Q6587,K2Q8254,K2Q8376 COMPLEX*16 K34,K34G5678,K35,K36 COMPLEX*16 K37,K38,K3Q2876,K3Q8276 COMPLEX*16 K45,K46,K47,K48 COMPLEX*16 K56,K57,K58,K67 COMPLEX*16 K68,K78,Q2487Q3165,Q2854Q8376 COMPLEX*16 Q3876Q8254,TRG5678 COMPLEX*16 Q2487(0:3),Q2854(0:3),Q2876(0:3),Q3165(0:3) COMPLEX*16 Q3876(0:3),Q4687(0:3),Q5687(0:3),Q6487(0:3) COMPLEX*16 Q6587(0:3),Q8254(0:3),Q8276(0:3),Q8376(0:3) COMPLEX*16 G5678(0:3,0:3) C NUMERATOR = (0.0D0,0.0D0) C IF (GRAPHNUMBER.EQ.1) THEN C CALL QPROP(CUT,KC,5,6,8,7,-1,-1,Q5687) CALL QPROP(CUT,KC,6,5,8,7,-1,1,Q6587) CALL GPROP(CUT,KC,5,6,7,8,1,1,G5678) K12 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K22 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K1Q5687 = (0.0D0,0.0D0) K1Q6587 = (0.0D0,0.0D0) K2Q5687 = (0.0D0,0.0D0) K2Q6587 = (0.0D0,0.0D0) TRG5678 = (0.0D0,0.0D0) K14G5678 = (0.0D0,0.0D0) K24G5678 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K22 = K22 + KC(2,MU)*KC(2,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K1Q5687 = K1Q5687 + KC(1,MU)*Q5687(MU)*METRIC(MU) K1Q6587 = K1Q6587 + KC(1,MU)*Q6587(MU)*METRIC(MU) K2Q5687 = K2Q5687 + KC(2,MU)*Q5687(MU)*METRIC(MU) K2Q6587 = K2Q6587 + KC(2,MU)*Q6587(MU)*METRIC(MU) TRG5678 = TRG5678 + G5678(MU,MU)*METRIC(MU) ENDDO DO MU = 0,3 DO NU = 0,3 K14G5678 = K14G5678 > + KC(1,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K24G5678 = K24G5678 > + KC(2,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) ENDDO ENDDO NUMERATOR = 8*(-1 + NC)*(1 + NC)*(-2*K14G5678*K22 + K1Q5687*K22 > - K1Q6587*K22 + 4*K12*K24G5678 - 2*K12*K2Q5687 + 2*K12*K2Q6587 > + K14*K22*TRG5678 - 2*K12*K24*TRG5678) C ELSE IF (GRAPHNUMBER.EQ.2) THEN C CALL QPROP(CUT,KC,2,8,5,4,-1,-1,Q2854) CALL QPROP(CUT,KC,3,8,7,6,-1,-1,Q3876) CALL QPROP(CUT,KC,8,2,5,4,1,1,Q8254) CALL QPROP(CUT,KC,8,3,7,6,-1,1,Q8376) K12 = (0.0D0,0.0D0) K1Q2854 = (0.0D0,0.0D0) K1Q3876 = (0.0D0,0.0D0) K1Q8254 = (0.0D0,0.0D0) K1Q8376 = (0.0D0,0.0D0) K2Q2854 = (0.0D0,0.0D0) K2Q3876 = (0.0D0,0.0D0) K2Q8254 = (0.0D0,0.0D0) K2Q8376 = (0.0D0,0.0D0) Q2854Q8376 = (0.0D0,0.0D0) Q3876Q8254 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K1Q2854 = K1Q2854 + KC(1,MU)*Q2854(MU)*METRIC(MU) K1Q3876 = K1Q3876 + KC(1,MU)*Q3876(MU)*METRIC(MU) K1Q8254 = K1Q8254 + KC(1,MU)*Q8254(MU)*METRIC(MU) K1Q8376 = K1Q8376 + KC(1,MU)*Q8376(MU)*METRIC(MU) K2Q2854 = K2Q2854 + KC(2,MU)*Q2854(MU)*METRIC(MU) K2Q3876 = K2Q3876 + KC(2,MU)*Q3876(MU)*METRIC(MU) K2Q8254 = K2Q8254 + KC(2,MU)*Q8254(MU)*METRIC(MU) K2Q8376 = K2Q8376 + KC(2,MU)*Q8376(MU)*METRIC(MU) Q2854Q8376 = Q2854Q8376 + Q2854(MU)*Q8376(MU)*METRIC(MU) Q3876Q8254 = Q3876Q8254 + Q3876(MU)*Q8254(MU)*METRIC(MU) ENDDO NUMERATOR = 8*NC*(K1Q8376*K2Q2854 + K1Q8254*K2Q3876 > + K1Q3876*K2Q8254 + K1Q2854*K2Q8376 - K12*Q2854Q8376 > - K12*Q3876Q8254) C ELSE IF (GRAPHNUMBER.EQ.3) THEN C K12 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K15 = (0.0D0,0.0D0) K16 = (0.0D0,0.0D0) K17 = (0.0D0,0.0D0) K18 = (0.0D0,0.0D0) K22 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K25 = (0.0D0,0.0D0) K26 = (0.0D0,0.0D0) K27 = (0.0D0,0.0D0) K28 = (0.0D0,0.0D0) K45 = (0.0D0,0.0D0) K56 = (0.0D0,0.0D0) K57 = (0.0D0,0.0D0) K58 = (0.0D0,0.0D0) K67 = (0.0D0,0.0D0) K78 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K15 = K15 + KC(1,MU)*KC(5,MU)*METRIC(MU) K16 = K16 + KC(1,MU)*KC(6,MU)*METRIC(MU) K17 = K17 + KC(1,MU)*KC(7,MU)*METRIC(MU) K18 = K18 + KC(1,MU)*KC(8,MU)*METRIC(MU) K22 = K22 + KC(2,MU)*KC(2,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K25 = K25 + KC(2,MU)*KC(5,MU)*METRIC(MU) K26 = K26 + KC(2,MU)*KC(6,MU)*METRIC(MU) K27 = K27 + KC(2,MU)*KC(7,MU)*METRIC(MU) K28 = K28 + KC(2,MU)*KC(8,MU)*METRIC(MU) K45 = K45 + KC(4,MU)*KC(5,MU)*METRIC(MU) K56 = K56 + KC(5,MU)*KC(6,MU)*METRIC(MU) K57 = K57 + KC(5,MU)*KC(7,MU)*METRIC(MU) K58 = K58 + KC(5,MU)*KC(8,MU)*METRIC(MU) K67 = K67 + KC(6,MU)*KC(7,MU)*METRIC(MU) K78 = K78 + KC(7,MU)*KC(8,MU)*METRIC(MU) ENDDO NUMERATOR = (16*(-1 + NC)*(1 + NC)*(-2*K18*K22*K56 > + 4*K12*K28*K56 + K17*K22*K45*NC**2 - 2*K12*K27*K45*NC**2 > - K14*K22*K57*NC**2 + K16*K22*K57*NC**2 + 2*K12*K24*K57*NC**2 > - 2*K12*K26*K57*NC**2 + K17*K22*K58*NC**2 - 2*K12*K27*K58*NC**2 > - K15*K22*K67*NC**2 + 2*K12*K25*K67*NC**2 - K15*K22*K78*NC**2 > + 2*K12*K25*K78*NC**2))/NC C ELSE IF (GRAPHNUMBER.EQ.4) THEN C CALL GPROP(CUT,KC,5,6,7,8,1,1,G5678) K12 = (0.0D0,0.0D0) K13 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K23 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K34 = (0.0D0,0.0D0) TRG5678 = (0.0D0,0.0D0) K12G5678 = (0.0D0,0.0D0) K13G5678 = (0.0D0,0.0D0) K24G5678 = (0.0D0,0.0D0) K34G5678 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K13 = K13 + KC(1,MU)*KC(3,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K23 = K23 + KC(2,MU)*KC(3,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K34 = K34 + KC(3,MU)*KC(4,MU)*METRIC(MU) TRG5678 = TRG5678 + G5678(MU,MU)*METRIC(MU) ENDDO DO MU = 0,3 DO NU = 0,3 K12G5678 = K12G5678 > + KC(1,MU)*KC(2,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K13G5678 = K13G5678 > + KC(1,MU)*KC(3,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K24G5678 = K24G5678 > + KC(2,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K34G5678 = K34G5678 > + KC(3,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) ENDDO ENDDO NUMERATOR = -4*(-1 + NC)*(1 + NC)*(-2*K13G5678*K24 > - 2*K13*K24G5678 + 2*K12G5678*K34 + 2*K12*K34G5678 > + K14*K23*TRG5678 + K13*K24*TRG5678 - K12*K34*TRG5678) C ELSE IF (GRAPHNUMBER.EQ.5) THEN C CALL QPROP(CUT,KC,4,6,8,7,-1,-1,Q4687) CALL QPROP(CUT,KC,6,4,8,7,-1,1,Q6487) K23 = (0.0D0,0.0D0) K1Q4687 = (0.0D0,0.0D0) K1Q6487 = (0.0D0,0.0D0) DO MU = 0,3 K23 = K23 + KC(2,MU)*KC(3,MU)*METRIC(MU) K1Q4687 = K1Q4687 + KC(1,MU)*Q4687(MU)*METRIC(MU) K1Q6487 = K1Q6487 + KC(1,MU)*Q6487(MU)*METRIC(MU) ENDDO NUMERATOR = -16*(K1Q4687 - K1Q6487)*K23*(-1 + NC)*(1 + NC) C ELSE IF (GRAPHNUMBER.EQ.6) THEN C CALL QPROP(CUT,KC,2,8,7,6,-1,-1,Q2876) CALL QPROP(CUT,KC,8,2,7,6,-1,1,Q8276) K14 = (0.0D0,0.0D0) K3Q2876 = (0.0D0,0.0D0) K3Q8276 = (0.0D0,0.0D0) DO MU = 0,3 K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K3Q2876 = K3Q2876 + KC(3,MU)*Q2876(MU)*METRIC(MU) K3Q8276 = K3Q8276 + KC(3,MU)*Q8276(MU)*METRIC(MU) ENDDO NUMERATOR = -16*K14*(K3Q2876 - K3Q8276)*(-1 + NC)*(1 + NC) C ELSE IF (GRAPHNUMBER.EQ.7) THEN C K12 = (0.0D0,0.0D0) K13 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K15 = (0.0D0,0.0D0) K16 = (0.0D0,0.0D0) K17 = (0.0D0,0.0D0) K18 = (0.0D0,0.0D0) K23 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K25 = (0.0D0,0.0D0) K26 = (0.0D0,0.0D0) K27 = (0.0D0,0.0D0) K28 = (0.0D0,0.0D0) K34 = (0.0D0,0.0D0) K35 = (0.0D0,0.0D0) K36 = (0.0D0,0.0D0) K37 = (0.0D0,0.0D0) K38 = (0.0D0,0.0D0) K45 = (0.0D0,0.0D0) K46 = (0.0D0,0.0D0) K47 = (0.0D0,0.0D0) K48 = (0.0D0,0.0D0) K67 = (0.0D0,0.0D0) K68 = (0.0D0,0.0D0) K78 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K13 = K13 + KC(1,MU)*KC(3,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K15 = K15 + KC(1,MU)*KC(5,MU)*METRIC(MU) K16 = K16 + KC(1,MU)*KC(6,MU)*METRIC(MU) K17 = K17 + KC(1,MU)*KC(7,MU)*METRIC(MU) K18 = K18 + KC(1,MU)*KC(8,MU)*METRIC(MU) K23 = K23 + KC(2,MU)*KC(3,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K25 = K25 + KC(2,MU)*KC(5,MU)*METRIC(MU) K26 = K26 + KC(2,MU)*KC(6,MU)*METRIC(MU) K27 = K27 + KC(2,MU)*KC(7,MU)*METRIC(MU) K28 = K28 + KC(2,MU)*KC(8,MU)*METRIC(MU) K34 = K34 + KC(3,MU)*KC(4,MU)*METRIC(MU) K35 = K35 + KC(3,MU)*KC(5,MU)*METRIC(MU) K36 = K36 + KC(3,MU)*KC(6,MU)*METRIC(MU) K37 = K37 + KC(3,MU)*KC(7,MU)*METRIC(MU) K38 = K38 + KC(3,MU)*KC(8,MU)*METRIC(MU) K45 = K45 + KC(4,MU)*KC(5,MU)*METRIC(MU) K46 = K46 + KC(4,MU)*KC(6,MU)*METRIC(MU) K47 = K47 + KC(4,MU)*KC(7,MU)*METRIC(MU) K48 = K48 + KC(4,MU)*KC(8,MU)*METRIC(MU) K67 = K67 + KC(6,MU)*KC(7,MU)*METRIC(MU) K68 = K68 + KC(6,MU)*KC(8,MU)*METRIC(MU) K78 = K78 + KC(7,MU)*KC(8,MU)*METRIC(MU) ENDDO NUMERATOR = (8*(-1 + NC)*(1 + NC)*(2*K18*K27*K34 + 2*K17*K28*K34 > - 2*K18*K24*K37 - 2*K14*K28*K37 - 2*K17*K24*K38 + 2*K14*K27*K38 > - 2*K18*K23*K47 - 2*K13*K28*K47 + 2*K12*K38*K47 + 2*K17*K23*K48 > - 2*K13*K27*K48 + 2*K12*K37*K48 - 2*K14*K23*K78 + 2*K13*K24*K78 > - 2*K12*K34*K78 + K17*K26*K34*NC**2 - K18*K26*K34*NC**2 > + K16*K27*K34*NC**2 - K16*K28*K34*NC**2 - 2*K14*K26*K35*NC**2 > - K17*K24*K36*NC**2 + K18*K24*K36*NC**2 + 2*K14*K25*K36*NC**2 > - K14*K27*K36*NC**2 - K14*K28*K36*NC**2 - K16*K24*K37*NC**2 > + K14*K26*K37*NC**2 + K16*K24*K38*NC**2 + K14*K26*K38*NC**2 > - 2*K16*K23*K45*NC**2 + 2*K15*K23*K46*NC**2 - K17*K23*K46*NC**2 > - K18*K23*K46*NC**2 - K13*K27*K46*NC**2 + K13*K28*K46*NC**2 > + K12*K37*K46*NC**2 - K12*K38*K46*NC**2 + K16*K23*K47*NC**2 > - K13*K26*K47*NC**2 + K12*K36*K47*NC**2 + K16*K23*K48*NC**2 > + K13*K26*K48*NC**2 - K12*K36*K48*NC**2 + 3*K14*K23*K67*NC**2 > + K13*K24*K67*NC**2 - K12*K34*K67*NC**2 - 3*K14*K23*K68*NC**2 > - K13*K24*K68*NC**2 + K12*K34*K68*NC**2))/NC C ELSE IF (GRAPHNUMBER.EQ.8) THEN C K12 = (0.0D0,0.0D0) K13 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K16 = (0.0D0,0.0D0) K17 = (0.0D0,0.0D0) K23 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K26 = (0.0D0,0.0D0) K27 = (0.0D0,0.0D0) K34 = (0.0D0,0.0D0) K36 = (0.0D0,0.0D0) K37 = (0.0D0,0.0D0) K46 = (0.0D0,0.0D0) K47 = (0.0D0,0.0D0) K67 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K13 = K13 + KC(1,MU)*KC(3,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K16 = K16 + KC(1,MU)*KC(6,MU)*METRIC(MU) K17 = K17 + KC(1,MU)*KC(7,MU)*METRIC(MU) K23 = K23 + KC(2,MU)*KC(3,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K26 = K26 + KC(2,MU)*KC(6,MU)*METRIC(MU) K27 = K27 + KC(2,MU)*KC(7,MU)*METRIC(MU) K34 = K34 + KC(3,MU)*KC(4,MU)*METRIC(MU) K36 = K36 + KC(3,MU)*KC(6,MU)*METRIC(MU) K37 = K37 + KC(3,MU)*KC(7,MU)*METRIC(MU) K46 = K46 + KC(4,MU)*KC(6,MU)*METRIC(MU) K47 = K47 + KC(4,MU)*KC(7,MU)*METRIC(MU) K67 = K67 + KC(6,MU)*KC(7,MU)*METRIC(MU) ENDDO NUMERATOR = (8*(K17*K26*K34 + K16*K27*K34 - K17*K24*K36 > - K14*K27*K36 - K16*K24*K37 + K14*K26*K37 + K17*K23*K46 > - K13*K27*K46 + K12*K37*K46 - K16*K23*K47 - K13*K26*K47 > + K12*K36*K47 + K14*K23*K67 + K13*K24*K67 - K12*K34*K67)*(-1 > + NC)**2*(1 + NC)**2)/NC C ELSE IF (GRAPHNUMBER.EQ.9) THEN C CALL QPROP(CUT,KC,2,4,8,7,-1,-1,Q2487) CALL QPROP(CUT,KC,3,1,6,5,-1,1,Q3165) Q2487Q3165 = (0.0D0,0.0D0) DO MU = 0,3 Q2487Q3165 = Q2487Q3165 + Q2487(MU)*Q3165(MU)*METRIC(MU) ENDDO NUMERATOR = -8*NC*Q2487Q3165 C ELSE IF (GRAPHNUMBER.EQ.10) THEN C K13 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K67 = (0.0D0,0.0D0) DO MU = 0,3 K13 = K13 + KC(1,MU)*KC(3,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K67 = K67 + KC(6,MU)*KC(7,MU)*METRIC(MU) ENDDO NUMERATOR = (-32*K13*K24*K67*(-1 + NC)*(1 + NC))/NC C ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C Numerator Function for Renormalization C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RNUMERATOR(GRAPHNUMBER,KC,MUMSBAR,CUT) C C In: INTEGER SIZE PARAMETER (SIZE = 3) INTEGER GRAPHNUMBER COMPLEX*16 KC(0:3*SIZE-1,0:3) REAL*8 MUMSBAR LOGICAL CUT(3*SIZE-1) C C Numerator function for renormalization of graph GRAPHNUMBER, C with complex momenta KC, including the color factor. C C For input here, we want CUT(P) = .TRUE. for the cut propagators C but not for a propagator with a loopcut. C C First version 28 August 1997. C This version written by Mathematica code of 4 September 1998 on C 4 Sep 1998 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C REAL*8 METRIC(0:3) DATA METRIC /+1.0D0,-1.0D0,-1.0D0,-1.0D0/ C INTEGER MU,NU C COMPLEX*16 K12,K12G5678,K13,K13G5678 COMPLEX*16 K14,K14G5678,K15,K16 COMPLEX*16 K17,K1L,K1Q2854,K1Q2876 COMPLEX*16 K1Q3876,K1Q4687,K1Q5687,K1Q6487 COMPLEX*16 K1Q6587,K1Q8254,K1Q8276,K1Q8376 COMPLEX*16 K22,K23,K24,K24G5678 COMPLEX*16 K25,K26,K27,K2L COMPLEX*16 K2Q2854,K2Q3876,K2Q5687,K2Q6587 COMPLEX*16 K2Q8254,K2Q8376,K34,K34G5678 COMPLEX*16 K36,K37,K3L,K3Q2876 COMPLEX*16 K3Q4687,K3Q6487,K3Q8276,K46 COMPLEX*16 K47,K4L,K5L,K67 COMPLEX*16 K6L,K7L,KLL,KLQ2876 COMPLEX*16 KLQ4687,KLQ6487,KLQ8276,NK1 COMPLEX*16 NK2,NK3,NK4,NK5 COMPLEX*16 NK6,NK7,NQ2876,NQ4687 COMPLEX*16 NQ6487,NQ8276,Q2487Q3165,Q2854Q8376 COMPLEX*16 Q3165Q2487,Q3876Q8254,Q8254Q3876,Q8376Q2854 COMPLEX*16 RQQG3AG,RQQG3AGV,RQQG3AQ,RQQG3AQV COMPLEX*16 RQQG3BG,RQQG3BGV,RQQG3BQ,RQQG3BQV COMPLEX*16 RQQG3CG,RQQG3CGV,RQQG3CQ,RQQG3CQV COMPLEX*16 RQQP3AQ,RQQP3AQV,RQQP3BQ,RQQP3BQV COMPLEX*16 RQQP3CQ,RQQP3CQV,TRG5678 COMPLEX*16 KLOOP(0:3),Q2487(0:3),Q2854(0:3),Q2876(0:3) COMPLEX*16 Q3165(0:3),Q3876(0:3),Q4687(0:3),Q5687(0:3) COMPLEX*16 Q6487(0:3),Q6587(0:3),Q8254(0:3),Q8276(0:3) COMPLEX*16 Q8376(0:3) COMPLEX*16 G5678(0:3,0:3) C RNUMERATOR = (0.0D0,0.0D0) C IF (GRAPHNUMBER.EQ.1) THEN C C ***** Divergent one loop subgraph is {7, 8}. ***** C IF (CUT(7).OR.CUT(8)) THEN C C The divergent subgraph is cut. C RNUMERATOR = (0.00D0,0.00D0) C ELSE C C The divergent subgraph {7, 8} is virtual. C CALL QPROPR(CUT,MUMSBAR,KC,5,6,8,7,-1,-1,Q5687) CALL QPROPR(CUT,MUMSBAR,KC,6,5,8,7,-1,1,Q6587) CALL GPROPR(CUT,MUMSBAR,KC,5,6,7,8,1,1,G5678) K12 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K22 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K1Q5687 = (0.0D0,0.0D0) K1Q6587 = (0.0D0,0.0D0) K2Q5687 = (0.0D0,0.0D0) K2Q6587 = (0.0D0,0.0D0) TRG5678 = (0.0D0,0.0D0) K14G5678 = (0.0D0,0.0D0) K24G5678 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K22 = K22 + KC(2,MU)*KC(2,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K1Q5687 = K1Q5687 + KC(1,MU)*Q5687(MU)*METRIC(MU) K1Q6587 = K1Q6587 + KC(1,MU)*Q6587(MU)*METRIC(MU) K2Q5687 = K2Q5687 + KC(2,MU)*Q5687(MU)*METRIC(MU) K2Q6587 = K2Q6587 + KC(2,MU)*Q6587(MU)*METRIC(MU) TRG5678 = TRG5678 + G5678(MU,MU)*METRIC(MU) ENDDO DO MU = 0,3 DO NU = 0,3 K14G5678 = K14G5678 > + KC(1,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K24G5678 = K24G5678 > + KC(2,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) ENDDO ENDDO RNUMERATOR = 8*(-1 + NC)*(1 + NC)*(-2*K14G5678*K22 + K1Q5687*K22 > - K1Q6587*K22 + 4*K12*K24G5678 - 2*K12*K2Q5687 + 2*K12*K2Q6587 > + K14*K22*TRG5678 - 2*K12*K24*TRG5678) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.2) THEN C C ***** Divergent one loop subgraphs are {{4, 5}, {6, 7}}. ***** C IF (CUT(4).OR.CUT(5)) THEN C C The divergent subgraph {6, 7} is virtual. C CALL QPROP(CUT,KC,2,8,5,4,-1,-1,Q2854) CALL QPROP(CUT,KC,8,2,5,4,1,1,Q8254) CALL QPROPR(CUT,MUMSBAR,KC,3,8,7,6,-1,-1,Q3876) CALL QPROPR(CUT,MUMSBAR,KC,8,3,7,6,-1,1,Q8376) K12 = (0.0D0,0.0D0) K1Q2854 = (0.0D0,0.0D0) K1Q3876 = (0.0D0,0.0D0) K1Q8254 = (0.0D0,0.0D0) K1Q8376 = (0.0D0,0.0D0) K2Q2854 = (0.0D0,0.0D0) K2Q3876 = (0.0D0,0.0D0) K2Q8254 = (0.0D0,0.0D0) K2Q8376 = (0.0D0,0.0D0) Q2854Q8376 = (0.0D0,0.0D0) Q8254Q3876 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K1Q2854 = K1Q2854 + KC(1,MU)*Q2854(MU)*METRIC(MU) K1Q3876 = K1Q3876 + KC(1,MU)*Q3876(MU)*METRIC(MU) K1Q8254 = K1Q8254 + KC(1,MU)*Q8254(MU)*METRIC(MU) K1Q8376 = K1Q8376 + KC(1,MU)*Q8376(MU)*METRIC(MU) K2Q2854 = K2Q2854 + KC(2,MU)*Q2854(MU)*METRIC(MU) K2Q3876 = K2Q3876 + KC(2,MU)*Q3876(MU)*METRIC(MU) K2Q8254 = K2Q8254 + KC(2,MU)*Q8254(MU)*METRIC(MU) K2Q8376 = K2Q8376 + KC(2,MU)*Q8376(MU)*METRIC(MU) Q2854Q8376 = Q2854Q8376 + Q2854(MU)*Q8376(MU)*METRIC(MU) Q8254Q3876 = Q8254Q3876 + Q8254(MU)*Q3876(MU)*METRIC(MU) ENDDO RNUMERATOR = 8*NC*(K1Q8376*K2Q2854 + K1Q8254*K2Q3876 > + K1Q3876*K2Q8254 + K1Q2854*K2Q8376 - K12*Q2854Q8376 > - K12*Q8254Q3876) C ELSE C C The divergent subgraph {4, 5} is virtual. C CALL QPROP(CUT,KC,3,8,7,6,-1,-1,Q3876) CALL QPROP(CUT,KC,8,3,7,6,-1,1,Q8376) CALL QPROPR(CUT,MUMSBAR,KC,2,8,5,4,-1,-1,Q2854) CALL QPROPR(CUT,MUMSBAR,KC,8,2,5,4,1,1,Q8254) K12 = (0.0D0,0.0D0) K1Q2854 = (0.0D0,0.0D0) K1Q3876 = (0.0D0,0.0D0) K1Q8254 = (0.0D0,0.0D0) K1Q8376 = (0.0D0,0.0D0) K2Q2854 = (0.0D0,0.0D0) K2Q3876 = (0.0D0,0.0D0) K2Q8254 = (0.0D0,0.0D0) K2Q8376 = (0.0D0,0.0D0) Q3876Q8254 = (0.0D0,0.0D0) Q8376Q2854 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K1Q2854 = K1Q2854 + KC(1,MU)*Q2854(MU)*METRIC(MU) K1Q3876 = K1Q3876 + KC(1,MU)*Q3876(MU)*METRIC(MU) K1Q8254 = K1Q8254 + KC(1,MU)*Q8254(MU)*METRIC(MU) K1Q8376 = K1Q8376 + KC(1,MU)*Q8376(MU)*METRIC(MU) K2Q2854 = K2Q2854 + KC(2,MU)*Q2854(MU)*METRIC(MU) K2Q3876 = K2Q3876 + KC(2,MU)*Q3876(MU)*METRIC(MU) K2Q8254 = K2Q8254 + KC(2,MU)*Q8254(MU)*METRIC(MU) K2Q8376 = K2Q8376 + KC(2,MU)*Q8376(MU)*METRIC(MU) Q3876Q8254 = Q3876Q8254 + Q3876(MU)*Q8254(MU)*METRIC(MU) Q8376Q2854 = Q8376Q2854 + Q8376(MU)*Q2854(MU)*METRIC(MU) ENDDO RNUMERATOR = 8*NC*(K1Q8376*K2Q2854 + K1Q8254*K2Q3876 > + K1Q3876*K2Q8254 + K1Q2854*K2Q8376 - K12*Q3876Q8254 > - K12*Q8376Q2854) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.3) THEN C C ***** Divergent one loop subgraphs are {{4, 5, 8}, {6, 7, 8}}. ***** C IF (CUT(4).OR.CUT(5)) THEN C C The divergent subgraph {6, 7, 8} is virtual. C NK1 = KC(1,0) NK2 = KC(2,0) NK5 = KC(5,0) K12 = (0.0D0,0.0D0) K15 = (0.0D0,0.0D0) K1L = (0.0D0,0.0D0) K22 = (0.0D0,0.0D0) K25 = (0.0D0,0.0D0) K2L = (0.0D0,0.0D0) K5L = (0.0D0,0.0D0) KLL = (0.0D0,0.0D0) KLOOP(0) = (0.0D0,0.0D0) DO MU = 1,3 KLOOP(MU) = (-KC(6,MU) + KC(7,MU) - KC(8,MU))/3.0D0 ENDDO DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K15 = K15 + KC(1,MU)*KC(5,MU)*METRIC(MU) K1L = K1L + KC(1,MU)*KLOOP(MU)*METRIC(MU) K22 = K22 + KC(2,MU)*KC(2,MU)*METRIC(MU) K25 = K25 + KC(2,MU)*KC(5,MU)*METRIC(MU) K2L = K2L + KC(2,MU)*KLOOP(MU)*METRIC(MU) K5L = K5L + KC(5,MU)*KLOOP(MU)*METRIC(MU) KLL = KLL + KLOOP(MU)*KLOOP(MU)*METRIC(MU) ENDDO RQQG3AGV = RQQG3AG(KLOOP, MUMSBAR) RQQG3AQV = RQQG3AQ(KLOOP, MUMSBAR) RQQG3BGV = RQQG3BG(KLOOP, MUMSBAR) RQQG3BQV = RQQG3BQ(KLOOP, MUMSBAR) RQQG3CGV = RQQG3CG(KLOOP, MUMSBAR) RQQG3CQV = RQQG3CQ(KLOOP, MUMSBAR) RNUMERATOR = -8*(-1 + NC)*(1 + NC)*(2*K15*K22*RQQG3AGV > - 4*K12*K25*RQQG3AGV + 2*K15*K22*RQQG3AQV - 4*K12*K25*RQQG3AQV > + K15*K22*RQQG3BGV - 2*K12*K25*RQQG3BGV - 2*K22*NK1*NK5*RQQG3BGV > + 4*K12*NK2*NK5*RQQG3BGV + K15*K22*RQQG3BQV - 2*K12*K25*RQQG3BQV > - 2*K22*NK1*NK5*RQQG3BQV + 4*K12*NK2*NK5*RQQG3BQV > - 2*K1L*K22*K5L*RQQG3CGV + 4*K12*K2L*K5L*RQQG3CGV > + K15*K22*KLL*RQQG3CGV - 2*K12*K25*KLL*RQQG3CGV > - 2*K1L*K22*K5L*RQQG3CQV + 4*K12*K2L*K5L*RQQG3CQV > + K15*K22*KLL*RQQG3CQV - 2*K12*K25*KLL*RQQG3CQV) C ELSE C C The divergent subgraph {4, 5, 8} is virtual. C NK1 = KC(1,0) NK2 = KC(2,0) NK6 = KC(6,0) NK7 = KC(7,0) K12 = (0.0D0,0.0D0) K16 = (0.0D0,0.0D0) K17 = (0.0D0,0.0D0) K1L = (0.0D0,0.0D0) K22 = (0.0D0,0.0D0) K26 = (0.0D0,0.0D0) K27 = (0.0D0,0.0D0) K2L = (0.0D0,0.0D0) K6L = (0.0D0,0.0D0) K7L = (0.0D0,0.0D0) KLL = (0.0D0,0.0D0) KLOOP(0) = (0.0D0,0.0D0) DO MU = 1,3 KLOOP(MU) = (-KC(4,MU) + KC(5,MU) - KC(8,MU))/3.0D0 ENDDO DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K16 = K16 + KC(1,MU)*KC(6,MU)*METRIC(MU) K17 = K17 + KC(1,MU)*KC(7,MU)*METRIC(MU) K1L = K1L + KC(1,MU)*KLOOP(MU)*METRIC(MU) K22 = K22 + KC(2,MU)*KC(2,MU)*METRIC(MU) K26 = K26 + KC(2,MU)*KC(6,MU)*METRIC(MU) K27 = K27 + KC(2,MU)*KC(7,MU)*METRIC(MU) K2L = K2L + KC(2,MU)*KLOOP(MU)*METRIC(MU) K6L = K6L + KC(6,MU)*KLOOP(MU)*METRIC(MU) K7L = K7L + KC(7,MU)*KLOOP(MU)*METRIC(MU) KLL = KLL + KLOOP(MU)*KLOOP(MU)*METRIC(MU) ENDDO RQQG3AGV = RQQG3AG(KLOOP, MUMSBAR) RQQG3AQV = RQQG3AQ(KLOOP, MUMSBAR) RQQG3BGV = RQQG3BG(KLOOP, MUMSBAR) RQQG3BQV = RQQG3BQ(KLOOP, MUMSBAR) RQQG3CGV = RQQG3CG(KLOOP, MUMSBAR) RQQG3CQV = RQQG3CQ(KLOOP, MUMSBAR) RNUMERATOR = 8*(-1 + NC)*(1 + NC)*(2*K17*K22*RQQG3AGV > - 4*K12*K27*RQQG3AGV + 2*K16*K22*RQQG3AQV - 4*K12*K26*RQQG3AQV > + K17*K22*RQQG3BGV - 2*K12*K27*RQQG3BGV - 2*K22*NK1*NK7*RQQG3BGV > + 4*K12*NK2*NK7*RQQG3BGV + K16*K22*RQQG3BQV - 2*K12*K26*RQQG3BQV > - 2*K22*NK1*NK6*RQQG3BQV + 4*K12*NK2*NK6*RQQG3BQV > - 2*K1L*K22*K7L*RQQG3CGV + 4*K12*K2L*K7L*RQQG3CGV > + K17*K22*KLL*RQQG3CGV - 2*K12*K27*KLL*RQQG3CGV > - 2*K1L*K22*K6L*RQQG3CQV + 4*K12*K2L*K6L*RQQG3CQV > + K16*K22*KLL*RQQG3CQV - 2*K12*K26*KLL*RQQG3CQV) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.4) THEN C C ***** Divergent one loop subgraph is {7, 8}. ***** C IF (CUT(7).OR.CUT(8)) THEN C C The divergent subgraph is cut. C RNUMERATOR = (0.00D0,0.00D0) C ELSE C C The divergent subgraph {7, 8} is virtual. C CALL GPROPR(CUT,MUMSBAR,KC,5,6,7,8,1,1,G5678) K12 = (0.0D0,0.0D0) K13 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K23 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K34 = (0.0D0,0.0D0) TRG5678 = (0.0D0,0.0D0) K12G5678 = (0.0D0,0.0D0) K13G5678 = (0.0D0,0.0D0) K24G5678 = (0.0D0,0.0D0) K34G5678 = (0.0D0,0.0D0) DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K13 = K13 + KC(1,MU)*KC(3,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K23 = K23 + KC(2,MU)*KC(3,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K34 = K34 + KC(3,MU)*KC(4,MU)*METRIC(MU) TRG5678 = TRG5678 + G5678(MU,MU)*METRIC(MU) ENDDO DO MU = 0,3 DO NU = 0,3 K12G5678 = K12G5678 > + KC(1,MU)*KC(2,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K13G5678 = K13G5678 > + KC(1,MU)*KC(3,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K24G5678 = K24G5678 > + KC(2,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) K34G5678 = K34G5678 > + KC(3,MU)*KC(4,NU)*G5678(MU,NU)*METRIC(MU)*METRIC(NU) ENDDO ENDDO RNUMERATOR = -4*(-1 + NC)*(1 + NC)*(-2*K13G5678*K24 > - 2*K13*K24G5678 + 2*K12G5678*K34 + 2*K12*K34G5678 > + K14*K23*TRG5678 + K13*K24*TRG5678 - K12*K34*TRG5678) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.5) THEN C C ***** Divergent one loop subgraphs are {{7, 8}, {1, 2, 5}}. ***** C IF (CUT(7).OR.CUT(8)) THEN C C The divergent subgraph {1, 2, 5} is virtual. C CALL QPROP(CUT,KC,4,6,8,7,-1,-1,Q4687) CALL QPROP(CUT,KC,6,4,8,7,-1,1,Q6487) NK3 = KC(3,0) NQ4687 = Q4687(0) NQ6487 = Q6487(0) K3L = (0.0D0,0.0D0) KLL = (0.0D0,0.0D0) K3Q4687 = (0.0D0,0.0D0) K3Q6487 = (0.0D0,0.0D0) KLQ4687 = (0.0D0,0.0D0) KLQ6487 = (0.0D0,0.0D0) KLOOP(0) = (0.0D0,0.0D0) DO MU = 1,3 KLOOP(MU) = (-KC(1,MU) + KC(2,MU) - KC(5,MU))/3.0D0 ENDDO DO MU = 0,3 K3L = K3L + KC(3,MU)*KLOOP(MU)*METRIC(MU) KLL = KLL + KLOOP(MU)*KLOOP(MU)*METRIC(MU) K3Q4687 = K3Q4687 + KC(3,MU)*Q4687(MU)*METRIC(MU) K3Q6487 = K3Q6487 + KC(3,MU)*Q6487(MU)*METRIC(MU) KLQ4687 = KLQ4687 + KLOOP(MU)*Q4687(MU)*METRIC(MU) KLQ6487 = KLQ6487 + KLOOP(MU)*Q6487(MU)*METRIC(MU) ENDDO RQQP3AQV = RQQP3AQ(KLOOP, MUMSBAR) RQQP3BQV = RQQP3BQ(KLOOP, MUMSBAR) RQQP3CQV = RQQP3CQ(KLOOP, MUMSBAR) RNUMERATOR = -4*NC*(2*K3Q4687*RQQP3AQV - 2*K3Q6487*RQQP3AQV > + K3Q4687*RQQP3BQV - K3Q6487*RQQP3BQV - 2*NK3*NQ4687*RQQP3BQV > + 2*NK3*NQ6487*RQQP3BQV + K3Q4687*KLL*RQQP3CQV > - K3Q6487*KLL*RQQP3CQV - 2*K3L*KLQ4687*RQQP3CQV > + 2*K3L*KLQ6487*RQQP3CQV) C ELSE C C The divergent subgraph {7, 8} is virtual. C CALL QPROPR(CUT,MUMSBAR,KC,4,6,8,7,-1,-1,Q4687) CALL QPROPR(CUT,MUMSBAR,KC,6,4,8,7,-1,1,Q6487) K23 = (0.0D0,0.0D0) K1Q4687 = (0.0D0,0.0D0) K1Q6487 = (0.0D0,0.0D0) DO MU = 0,3 K23 = K23 + KC(2,MU)*KC(3,MU)*METRIC(MU) K1Q4687 = K1Q4687 + KC(1,MU)*Q4687(MU)*METRIC(MU) K1Q6487 = K1Q6487 + KC(1,MU)*Q6487(MU)*METRIC(MU) ENDDO RNUMERATOR = -16*(K1Q4687 - K1Q6487)*K23*(-1 + NC)*(1 + NC) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.6) THEN C C ***** Divergent one loop subgraphs are {{6, 7}, {3, 4, 5}}. ***** C IF (CUT(6).OR.CUT(7)) THEN C C The divergent subgraph {3, 4, 5} is virtual. C CALL QPROP(CUT,KC,2,8,7,6,-1,-1,Q2876) CALL QPROP(CUT,KC,8,2,7,6,-1,1,Q8276) NK1 = KC(1,0) NQ2876 = Q2876(0) NQ8276 = Q8276(0) K1L = (0.0D0,0.0D0) KLL = (0.0D0,0.0D0) K1Q2876 = (0.0D0,0.0D0) K1Q8276 = (0.0D0,0.0D0) KLQ2876 = (0.0D0,0.0D0) KLQ8276 = (0.0D0,0.0D0) KLOOP(0) = (0.0D0,0.0D0) DO MU = 1,3 KLOOP(MU) = (-KC(3,MU) + KC(4,MU) - KC(5,MU))/3.0D0 ENDDO DO MU = 0,3 K1L = K1L + KC(1,MU)*KLOOP(MU)*METRIC(MU) KLL = KLL + KLOOP(MU)*KLOOP(MU)*METRIC(MU) K1Q2876 = K1Q2876 + KC(1,MU)*Q2876(MU)*METRIC(MU) K1Q8276 = K1Q8276 + KC(1,MU)*Q8276(MU)*METRIC(MU) KLQ2876 = KLQ2876 + KLOOP(MU)*Q2876(MU)*METRIC(MU) KLQ8276 = KLQ8276 + KLOOP(MU)*Q8276(MU)*METRIC(MU) ENDDO RQQP3AQV = RQQP3AQ(KLOOP, MUMSBAR) RQQP3BQV = RQQP3BQ(KLOOP, MUMSBAR) RQQP3CQV = RQQP3CQ(KLOOP, MUMSBAR) RNUMERATOR = -4*NC*(2*K1Q2876*RQQP3AQV - 2*K1Q8276*RQQP3AQV > + K1Q2876*RQQP3BQV - K1Q8276*RQQP3BQV - 2*NK1*NQ2876*RQQP3BQV > + 2*NK1*NQ8276*RQQP3BQV + K1Q2876*KLL*RQQP3CQV > - K1Q8276*KLL*RQQP3CQV - 2*K1L*KLQ2876*RQQP3CQV > + 2*K1L*KLQ8276*RQQP3CQV) C ELSE C C The divergent subgraph {6, 7} is virtual. C CALL QPROPR(CUT,MUMSBAR,KC,2,8,7,6,-1,-1,Q2876) CALL QPROPR(CUT,MUMSBAR,KC,8,2,7,6,-1,1,Q8276) K14 = (0.0D0,0.0D0) K3Q2876 = (0.0D0,0.0D0) K3Q8276 = (0.0D0,0.0D0) DO MU = 0,3 K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K3Q2876 = K3Q2876 + KC(3,MU)*Q2876(MU)*METRIC(MU) K3Q8276 = K3Q8276 + KC(3,MU)*Q8276(MU)*METRIC(MU) ENDDO RNUMERATOR = -16*K14*(K3Q2876 - K3Q8276)*(-1 + NC)*(1 + NC) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.7) THEN C C ***** Divergent one loop subgraph is {6, 7, 8}. ***** C IF (CUT(6).OR.CUT(7)) THEN C C The divergent subgraph is cut. C RNUMERATOR = (0.00D0,0.00D0) C ELSE C C The divergent subgraph {6, 7, 8} is virtual. C NK1 = KC(1,0) NK2 = KC(2,0) NK3 = KC(3,0) NK4 = KC(4,0) K12 = (0.0D0,0.0D0) K13 = (0.0D0,0.0D0) K14 = (0.0D0,0.0D0) K1L = (0.0D0,0.0D0) K23 = (0.0D0,0.0D0) K24 = (0.0D0,0.0D0) K2L = (0.0D0,0.0D0) K34 = (0.0D0,0.0D0) K3L = (0.0D0,0.0D0) K4L = (0.0D0,0.0D0) KLL = (0.0D0,0.0D0) KLOOP(0) = (0.0D0,0.0D0) DO MU = 1,3 KLOOP(MU) = (-KC(6,MU) + KC(7,MU) - KC(8,MU))/3.0D0 ENDDO DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K13 = K13 + KC(1,MU)*KC(3,MU)*METRIC(MU) K14 = K14 + KC(1,MU)*KC(4,MU)*METRIC(MU) K1L = K1L + KC(1,MU)*KLOOP(MU)*METRIC(MU) K23 = K23 + KC(2,MU)*KC(3,MU)*METRIC(MU) K24 = K24 + KC(2,MU)*KC(4,MU)*METRIC(MU) K2L = K2L + KC(2,MU)*KLOOP(MU)*METRIC(MU) K34 = K34 + KC(3,MU)*KC(4,MU)*METRIC(MU) K3L = K3L + KC(3,MU)*KLOOP(MU)*METRIC(MU) K4L = K4L + KC(4,MU)*KLOOP(MU)*METRIC(MU) KLL = KLL + KLOOP(MU)*KLOOP(MU)*METRIC(MU) ENDDO RQQG3AGV = RQQG3AG(KLOOP, MUMSBAR) RQQG3AQV = RQQG3AQ(KLOOP, MUMSBAR) RQQG3BGV = RQQG3BG(KLOOP, MUMSBAR) RQQG3BQV = RQQG3BQ(KLOOP, MUMSBAR) RQQG3CGV = RQQG3CG(KLOOP, MUMSBAR) RQQG3CQV = RQQG3CQ(KLOOP, MUMSBAR) RNUMERATOR = 8*(-1 + NC)*(1 + NC)*(4*K14*K23*RQQG3AGV > + 4*K14*K23*RQQG3AQV + K14*K23*RQQG3BGV + K13*K24*RQQG3BGV > - K12*K34*RQQG3BGV + 2*K34*NK1*NK2*RQQG3BGV > - 2*K24*NK1*NK3*RQQG3BGV - 2*K13*NK2*NK4*RQQG3BGV > + 2*K12*NK3*NK4*RQQG3BGV + K14*K23*RQQG3BQV + K13*K24*RQQG3BQV > - K12*K34*RQQG3BQV + 2*K34*NK1*NK2*RQQG3BQV > - 2*K24*NK1*NK3*RQQG3BQV - 2*K13*NK2*NK4*RQQG3BQV > + 2*K12*NK3*NK4*RQQG3BQV + 2*K1L*K2L*K34*RQQG3CGV > - 2*K1L*K24*K3L*RQQG3CGV - 2*K13*K2L*K4L*RQQG3CGV > + 2*K12*K3L*K4L*RQQG3CGV + K14*K23*KLL*RQQG3CGV > + K13*K24*KLL*RQQG3CGV - K12*K34*KLL*RQQG3CGV > + 2*K1L*K2L*K34*RQQG3CQV - 2*K1L*K24*K3L*RQQG3CQV > - 2*K13*K2L*K4L*RQQG3CQV + 2*K12*K3L*K4L*RQQG3CQV > + K14*K23*KLL*RQQG3CQV + K13*K24*KLL*RQQG3CQV > - K12*K34*KLL*RQQG3CQV) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.8) THEN C C ***** Divergent one loop subgraphs are {{1, 2, 5}, {3, 4, 8}}. ***** C IF (CUT(1).OR.CUT(2)) THEN C C The divergent subgraph {3, 4, 8} is virtual. C NK1 = KC(1,0) NK2 = KC(2,0) NK6 = KC(6,0) NK7 = KC(7,0) K12 = (0.0D0,0.0D0) K16 = (0.0D0,0.0D0) K17 = (0.0D0,0.0D0) K1L = (0.0D0,0.0D0) K26 = (0.0D0,0.0D0) K27 = (0.0D0,0.0D0) K2L = (0.0D0,0.0D0) K67 = (0.0D0,0.0D0) K6L = (0.0D0,0.0D0) K7L = (0.0D0,0.0D0) KLL = (0.0D0,0.0D0) KLOOP(0) = (0.0D0,0.0D0) DO MU = 1,3 KLOOP(MU) = (-KC(3,MU) + KC(4,MU) - KC(8,MU))/3.0D0 ENDDO DO MU = 0,3 K12 = K12 + KC(1,MU)*KC(2,MU)*METRIC(MU) K16 = K16 + KC(1,MU)*KC(6,MU)*METRIC(MU) K17 = K17 + KC(1,MU)*KC(7,MU)*METRIC(MU) K1L = K1L + KC(1,MU)*KLOOP(MU)*METRIC(MU) K26 = K26 + KC(2,MU)*KC(6,MU)*METRIC(MU) K27 = K27 + KC(2,MU)*KC(7,MU)*METRIC(MU) K2L = K2L + KC(2,MU)*KLOOP(MU)*METRIC(MU) K67 = K67 + KC(6,MU)*KC(7,MU)*METRIC(MU) K6L = K6L + KC(6,MU)*KLOOP(MU)*METRIC(MU) K7L = K7L + KC(7,MU)*KLOOP(MU)*METRIC(MU) KLL = KLL + KLOOP(MU)*KLOOP(MU)*METRIC(MU) ENDDO RQQP3AQV = RQQP3AQ(KLOOP, MUMSBAR) RQQP3BQV = RQQP3BQ(KLOOP, MUMSBAR) RQQP3CQV = RQQP3CQ(KLOOP, MUMSBAR) RNUMERATOR = 4*(-1 + NC)*(1 + NC)*(4*K17*K26*RQQP3AQV > + K17*K26*RQQP3BQV - K16*K27*RQQP3BQV + K12*K67*RQQP3BQV > - 2*K67*NK1*NK2*RQQP3BQV + 2*K27*NK1*NK6*RQQP3BQV > + 2*K16*NK2*NK7*RQQP3BQV - 2*K12*NK6*NK7*RQQP3BQV > - 2*K1L*K2L*K67*RQQP3CQV + 2*K1L*K27*K6L*RQQP3CQV > + 2*K16*K2L*K7L*RQQP3CQV - 2*K12*K6L*K7L*RQQP3CQV > + K17*K26*KLL*RQQP3CQV - K16*K27*KLL*RQQP3CQV > + K12*K67*KLL*RQQP3CQV) C ELSE C C The divergent subgraph {1, 2, 5} is virtual. C NK3 = KC(3,0) NK4 = KC(4,0) NK6 = KC(6,0) NK7 = KC(7,0) K34 = (0.0D0,0.0D0) K36 = (0.0D0,0.0D0) K37 = (0.0D0,0.0D0) K3L = (0.0D0,0.0D0) K46 = (0.0D0,0.0D0) K47 = (0.0D0,0.0D0) K4L = (0.0D0,0.0D0) K67 = (0.0D0,0.0D0) K6L = (0.0D0,0.0D0) K7L = (0.0D0,0.0D0) KLL = (0.0D0,0.0D0) KLOOP(0) = (0.0D0,0.0D0) DO MU = 1,3 KLOOP(MU) = (-KC(1,MU) + KC(2,MU) - KC(5,MU))/3.0D0 ENDDO DO MU = 0,3 K34 = K34 + KC(3,MU)*KC(4,MU)*METRIC(MU) K36 = K36 + KC(3,MU)*KC(6,MU)*METRIC(MU) K37 = K37 + KC(3,MU)*KC(7,MU)*METRIC(MU) K3L = K3L + KC(3,MU)*KLOOP(MU)*METRIC(MU) K46 = K46 + KC(4,MU)*KC(6,MU)*METRIC(MU) K47 = K47 + KC(4,MU)*KC(7,MU)*METRIC(MU) K4L = K4L + KC(4,MU)*KLOOP(MU)*METRIC(MU) K67 = K67 + KC(6,MU)*KC(7,MU)*METRIC(MU) K6L = K6L + KC(6,MU)*KLOOP(MU)*METRIC(MU) K7L = K7L + KC(7,MU)*KLOOP(MU)*METRIC(MU) KLL = KLL + KLOOP(MU)*KLOOP(MU)*METRIC(MU) ENDDO RQQP3AQV = RQQP3AQ(KLOOP, MUMSBAR) RQQP3BQV = RQQP3BQ(KLOOP, MUMSBAR) RQQP3CQV = RQQP3CQ(KLOOP, MUMSBAR) RNUMERATOR = 4*(-1 + NC)*(1 + NC)*(4*K37*K46*RQQP3AQV > + K37*K46*RQQP3BQV - K36*K47*RQQP3BQV + K34*K67*RQQP3BQV > - 2*K67*NK3*NK4*RQQP3BQV + 2*K47*NK3*NK6*RQQP3BQV > + 2*K36*NK4*NK7*RQQP3BQV - 2*K34*NK6*NK7*RQQP3BQV > - 2*K3L*K4L*K67*RQQP3CQV + 2*K3L*K47*K6L*RQQP3CQV > + 2*K36*K4L*K7L*RQQP3CQV - 2*K34*K6L*K7L*RQQP3CQV > + K37*K46*KLL*RQQP3CQV - K36*K47*KLL*RQQP3CQV > + K34*K67*KLL*RQQP3CQV) C ENDIF C ELSE IF (GRAPHNUMBER.EQ.9) THEN C C ***** Divergent one loop subgraphs are {{5, 6}, {7, 8}}. ***** C IF (CUT(5).OR.CUT(6)) THEN C C The divergent subgraph {7, 8} is virtual. C CALL QPROP(CUT,KC,3,1,6,5,-1,1,Q3165) CALL QPROPR(CUT,MUMSBAR,KC,2,4,8,7,-1,-1,Q2487) Q3165Q2487 = (0.0D0,0.0D0) DO MU = 0,3 Q3165Q2487 = Q3165Q2487 + Q3165(MU)*Q2487(MU)*METRIC(MU) ENDDO RNUMERATOR = -8*NC*Q3165Q2487 C ELSE C C The divergent subgraph {5, 6} is virtual. C CALL QPROP(CUT,KC,2,4,8,7,-1,-1,Q2487) CALL QPROPR(CUT,MUMSBAR,KC,3,1,6,5,-1,1,Q3165) Q2487Q3165 = (0.0D0,0.0D0) DO MU = 0,3 Q2487Q3165 = Q2487Q3165 + Q2487(MU)*Q3165(MU)*METRIC(MU) ENDDO RNUMERATOR = -8*NC*Q2487Q3165 C ENDIF C ELSE IF (GRAPHNUMBER.EQ.10) THEN C C ***** There is no divergent one loop subgraph. ***** C RNUMERATOR = (0.00D0,0.00D0) C C ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C Propagator Functions C C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE QPROP(CUT,KC,P1,P2,PQ,PG,S1,SQ,RESULT) C In: INTEGER SIZE PARAMETER (SIZE = 3) LOGICAL CUT(3*SIZE-1) COMPLEX*16 KC(0:3*SIZE-1,0:3) INTEGER P1,P2,PQ,PG,S1,SQ C Out: COMPLEX*16 RESULT(0:3) C C A quark self-energy subgraph, together with its adjoining propagator C factors, becomes C C sequark[p1,p2,pq,pg,s1,sq]_mu * gamma^mu C C where sequark is a four-vector, the RESULT here, and C p1 is the label of the outgoing quark line, C s1 = +1 if this line carries momentum in the quark direction, C s1 = -1 otherwise; C p2 is the label of the ingoing quark line; C pq is the label of the internal quark line, C sq = +1 if this line carries momentum in the quark direstion, C sq = -1 otherwise; C pg is the label of the internal gluon line. C C This subroutine calculates the integrand for sequark^mu. There C are three cases: C C 1) self-energy subgraph is cut. C 2) self-energy subgraph is virtual, an adjoining propagator is cut. C 3) self-energy subgraph is virtual, adnoining propagators not cut. C C We let Q^mu be the incoming quark 3-momentum, C KQ^mu be the 3-momentum carried by the internal quark line, C KG^mu be the 3-momentum carried by the internal gluon line. C C 4 January 1998 C 4 August 1998 Add call to CHECKDEFORM2 C 4 September 1998 Add color factors C 21 September 1998 Fix color factors C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 Q(1:3),KQ(1:3),KG(1:3) COMPLEX*16 QSQ,KQSQ,KGSQ,ABSQ,ABSKQ,ABSKG COMPLEX*16 EQ,EBAR,QBARSQ,Q4SQ COMPLEX*16 PREFACTOR,SPECIAL COMPLEX*16 COMPLEXSQRT INTEGER MU C EQ = S1*KC(P1,0) QSQ = (0.0D0,0.0D0) KQSQ = (0.0D0,0.0D0) KGSQ = (0.0D0,0.0D0) DO MU = 1,3 Q(MU) = S1*KC(P1,MU) KQ(MU) = SQ*KC(PQ,MU) KG(MU) = Q(MU) - KQ(MU) QSQ = QSQ + Q(MU)**2 KQSQ = KQSQ + KQ(MU)**2 KGSQ = KGSQ + KG(MU)**2 ENDDO C ABSQ = COMPLEXSQRT(QSQ) ABSKQ = COMPLEXSQRT(KQSQ) ABSKG = COMPLEXSQRT(KGSQ) C EBAR = ABSKG + ABSKQ QBARSQ = EBAR**2 - QSQ PREFACTOR = (NC**2 - 1)/(4*NC*ABSKG*ABSKQ*ABSQ*QBARSQ) C C Now we need the proper special factor for each special case. C The first case is that of a cut self-energy. C In the second case, the self-energy subgraph is virtual, but C one of the adjoing propagators is cut. There are two such cuts C possible. For each of them, we have a factor 1/2 (and a minus sign). C In the third case, the self-energy subgraph and both of the C adjoinging propagaators is virtual. Here Q4SQ is the square of the C four-vector q^mu. C IF (CUT(PQ).AND.CUT(PG)) THEN SPECIAL = ABSQ/EBAR ELSE IF (CUT(P1).OR.CUT(P2)) THEN SPECIAL = (-0.5D0,0.0D0) ELSE CALL CHECKDEFORM2(QBARSQ,CUT,P1,P2) Q4SQ = EQ**2 - QSQ SPECIAL = - 2.0D0*ABSQ*QBARSQ/(Q4SQ*(QBARSQ-Q4SQ)) ENDIF C C Finally, we assemble the result. C RESULT(0) = SPECIAL*PREFACTOR*EQ*ABSKG DO MU = 1,3 RESULT(MU) = SPECIAL*PREFACTOR*EBAR*KG(MU) ENDDO C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE QPROPR(CUT,MUMSBAR,KC,P1,P2,PQ,PG,S1,SQ,RESULT) C In: INTEGER SIZE PARAMETER (SIZE = 3) LOGICAL CUT(3*SIZE-1) REAL*8 MUMSBAR COMPLEX*16 KC(0:3*SIZE-1,0:3) INTEGER P1,P2,PQ,PG,S1,SQ C Out: COMPLEX*16 RESULT(0:3) C C A quark self-energy subgraph, together with its adjoining propagator C factors, becomes C C sequark[p1,p2,pq,pg,s1,sq]_mu * gamma^mu C C where sequark is a four-vector, the RESULT here, and C p1 is the label of the outgoing quark line, C s1 = +1 if this line carries momentum in the quark direction, C s1 = -1 otherwise; C p2 is the label of the ingoing quark line; C pq is the label of the internal quark line, C sq = +1 if this line carries momentum in the quark direstion, C sq = -1 otherwise; C pg is the label of the internal gluon line. C C This subroutine calculates the integrand for the renormalization C counterterm sequark^mu. There are three cases: C C 1) self-energy subgraph is cut. (QPROPR not called in this case) C 2) self-energy subgraph is virtual, an adjoining propagator is cut. C 3) self-energy subgraph is virtual, adjoining propagators not cut. C C We let Q^mu be the incoming quark 3-momentum, C KQ^mu be the 3-momentum carried by the internal quark line, C KG^mu be the 3-momentum carried by the internal gluon line. C C 10 January 1998 C 4 September 1998 Add color factors C 28 October 1998 Add another subtraction term. C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 Q(1:3),KQ(1:3),ELL(1:3) COMPLEX*16 QSQ,ELLSQ,ABSQ COMPLEX*16 EQ,Q4SQ,DOTQL COMPLEX*16 DENOM,PREFACTOR,SPECIAL COMPLEX*16 FACTOREQ,FACTORQ,FACTORL COMPLEX*16 COMPLEXSQRT INTEGER MU C EQ = S1*KC(P1,0) QSQ = (0.0D0,0.0D0) ELLSQ = (0.0D0,0.0D0) DOTQL = (0.0D0,0.0D0) DO MU = 1,3 Q(MU) = S1*KC(P1,MU) KQ(MU) = SQ*KC(PQ,MU) ELL(MU) = 0.5D0*Q(MU) - KQ(MU) QSQ = QSQ + Q(MU)**2 ELLSQ = ELLSQ + ELL(MU)**2 DOTQL = DOTQL + Q(MU)*ELL(MU) ENDDO C ABSQ = COMPLEXSQRT(QSQ) C DENOM = ELLSQ + MUMSBAR**2 PREFACTOR = (NC**2 - 1)/(16*NC*ABSQ*COMPLEXSQRT(DENOM)**3) FACTOREQ = (DOTQL + 3*MUMSBAR**2)/(2*DENOM) FACTORQ = 3*MUMSBAR**2/(2*DENOM) FACTORL = (2*EQ**2 - 3*QSQ + 12*MUMSBAR**2)*ELLSQ + 5*DOTQL**2 FACTORL = FACTORL /(4*DENOM**2) C C Now we need the proper special factor for each special case. C The first case is that of a cut self-energy. C In the second case, the self-energy subgraph is virtual, but C one of the adjoing propagators is cut. There are two such cuts C possible. For each of them, we have a factor 1/2 (and a minus sign). C In the third case, the self-energy subgraph and both of the C adjoinging propagaators is virtual. Here Q4SQ is the square of the C four-vector q^mu. C IF (CUT(PQ).AND.CUT(PG)) THEN WRITE(NOUT,*)'Something is rotten in QPROPR.' STOP ELSE IF (CUT(P1).OR.CUT(P2)) THEN SPECIAL = (-0.5D0,0.0D0) ELSE Q4SQ = EQ**2 - QSQ SPECIAL = - 2.0D0*ABSQ/Q4SQ ENDIF C C Finally, we assemble the result. C RESULT(0) = SPECIAL*PREFACTOR*(EQ + FACTOREQ*EQ) DO MU = 1,3 RESULT(MU) = SPECIAL*PREFACTOR* > (2*ELL(MU) + Q(MU) + FACTORQ*Q(MU) + FACTORL*ELL(MU)) ENDDO C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE GPROP(CUT,KC,P1,P2,PG1,PG2,S1,SG1,RESULT) C In: INTEGER SIZE PARAMETER (SIZE = 3) LOGICAL CUT(3*SIZE-1) COMPLEX*16 KC(0:3*SIZE-1,0:3) INTEGER P1,P2,PG1,PG2,S1,SG1 C Out: COMPLEX*16 RESULT(0:3,0:3) C C A gluon self-energy subgraph, together with its adjoining propagator C factors, becomes C C seglue[p1,p2,pg1,pg2,s1,sg1]^{mu,nu} C C where seglue is a second rank tensor, the RESULT here, with C seglue^{0i} = seglue{i,0} = seglue(0,0) = 0 for a virtual C graph at q^2 = 0, but not for the cut graph. Here C p1 is the label of the incoming gluon line, C s1 = +1 if this line carries momentum into the graph, C s1 = -1 otherwise; C p2 is the label of the outgoing line; C pg1 is the label of one of the the internal parton lines, C sg1 = +1 if this line carries momentum away from the vertex C at which p1 joins the subgraph, C sg1 = -1 otherwise; C pg2 is the label of the other internal parton line. C C C 13 January 1998 C 7 April 1998 C 4 September 1998 Add color factors C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C REAL*8 DELTA COMPLEX*16 QSQ,Q4SQ,ELLSQ,ELL4SQ,KG1SQ,KG2SQ,DOTQELL COMPLEX*16 ABSQ,ABSKG1,ABSKG2 COMPLEX*16 KG1(3),KG2(3),Q(3),Q0,ELL(3),ELL0 COMPLEX*16 ELLP(0:3),T(0:3,0:3) COMPLEX*16 PREFACTOR,SPECIAL,TEMP COMPLEX*16 COMPLEXSQRT INTEGER CUTSIGN REAL*8 XXREAL INTEGER MU,NU C C We define CUTSIGN to be +1 if the self-energy subgraph crosses C the final state cut with propagator P1 to the left and P2 to C the right, -1 if it is the other way around. This allows us C to choose the signs so that q^mu is flowing toward the cut. C IF (S1*XXREAL(KC(P1,0)).GT.0.0D0) THEN CUTSIGN = 1 ELSE CUTSIGN = -1 ENDIF C QSQ = (0.0D0,0.0D0) ELLSQ = (0.0D0,0.0D0) KG1SQ = (0.0D0,0.0D0) KG2SQ = (0.0D0,0.0D0) DOTQELL = (0.0D0,0.0D0) DO MU = 1,3 Q(MU) = CUTSIGN*S1*KC(P1,MU) KG1(MU) = CUTSIGN*SG1*KC(PG1,MU) KG2(MU) = Q(MU) - KG1(MU) ELL(MU) = KG1(MU) - 0.5D0*Q(MU) QSQ = QSQ + Q(MU)**2 ELLSQ = ELLSQ + ELL(MU)**2 DOTQELL = DOTQELL + Q(MU)*ELL(MU) KG1SQ = KG1SQ + KG1(MU)**2 KG2SQ = KG2SQ + KG2(MU)**2 ENDDO ABSQ = COMPLEXSQRT(QSQ) ABSKG1 = COMPLEXSQRT(KG1SQ) ABSKG2 = COMPLEXSQRT(KG2SQ) C Q0 = ABSKG1 + ABSKG2 ELL0 = 0.5D0*(ABSKG1 - ABSKG2) Q4SQ = Q0**2 - QSQ ELL4SQ = ELL0**2 - ELLSQ C ELLP(0) = ELL0 - Q0*DOTQELL/QSQ T(0,0) = - Q4SQ/QSQ DO MU = 1,3 ELLP(MU) = ELL(MU) - Q(MU)*DOTQELL/QSQ T(0,MU) = 0.0D0 T(MU,0) = 0.0D0 DO NU = 1,3 T(MU,NU) = - DELTA(MU,NU) + Q(MU)*Q(NU)/QSQ ENDDO ENDDO C C PREFACTOR = 1.0D0/(4.0D0*ABSKG1*ABSKG2*Q4SQ**2) C C Now we need the proper special factor for each special case. C The first case is that of a cut self-energy. C In the second case, the self-energy subgraph is virtual, but C one of the adjoing propagators is cut. There are two such cuts C possible. For each of them, we have a factor 1/2 (and a minus sign). C In the third case, the self-energy subgraph and both of the C adjoining propagators are virtual. This does not happen for gluon C self energies in three loop graphs. Here Q4SQ is the square of the C four-vector q^mu. C IF (CUT(PG1).AND.CUT(PG2)) THEN SPECIAL = (1.0D0,0.0D0) TEMP = (NC + 2*NF)*ELL4SQ + (2.25D0*NC - 0.5D0*NF)*Q4SQ DO MU = 0,3 DO NU = 0,3 RESULT(MU,NU) = SPECIAL*PREFACTOR* > ( 4*(NC - NF)*ELLP(MU)*ELLP(NU) + TEMP*T(MU,NU)) ENDDO ENDDO ELSE IF (CUT(P1).OR.CUT(P2)) THEN SPECIAL = - 0.5D0*Q0/ABSQ TEMP = (NC + 2*NF)*ELL4SQ + (2.25D0*NC - 0.5D0*NF)*Q4SQ RESULT(0,0) = (0.0D0,0.0D0) DO MU = 1,3 RESULT(0,MU) = (0.0D0,0.0D0) RESULT(MU,0) = (0.0D0,0.0D0) DO NU = 1,3 RESULT(MU,NU) = SPECIAL*PREFACTOR* > ( 4*(NC - NF)*ELLP(MU)*ELLP(NU) + TEMP*T(MU,NU)) ENDDO ENDDO ELSE WRITE(NOUT,*)'Something is rotten in GPROP.' STOP ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE GPROPR(CUT,MUMSBAR,KC,P1,P2,PG1,PG2,S1,SG1,RESULT) C In: INTEGER SIZE PARAMETER (SIZE = 3) LOGICAL CUT(3*SIZE-1) REAL*8 MUMSBAR COMPLEX*16 KC(0:3*SIZE-1,0:3) INTEGER P1,P2,PG1,PG2,S1,SG1 C Out: COMPLEX*16 RESULT(0:3,0:3) C C A gluon self-energy subgraph, together with its adjoining propagator C factors, becomes C C seglue[p1,p2,pg1,pg2,s1,sg1]^{mu,nu} C C where seglue is a second rank tensor, the RESULT here, with C seglue^{0i} = seglue{i,0} = seglue(0,0) = 0. Here C p1 is the label of the incoming gluon line, C s1 = +1 if this line carries momentum into the graph, C s1 = -1 otherwise; C p2 is the label of the outgoing line; C pg1 is the label of one of the the internal gluon lines, C sg1 = +1 if this line carries away from the vertex C at which p1 joins the subgraph, C sg1 = -1 otherwise; C pg2 is the label of the other internal gluon line. C C This subroutine calculates the integrand for the renormalization C counterterm for sequark^mu. There are two cases C C 1) self-energy subgraph is cut. (Then GPROPR is not called.) C 2) self-energy subgraph is virtual, an adjoining propagator is cut. C C We let Q^mu be the incoming quark 3-momentum, C KG1^mu be the 3-momentum carried by the first internal line, C KG2^mu be the 3-momentum carried by the other internal line. C C 13 January 1998 C 7 April 1998 C 4 September 1998 Add color factors C 21 September 1998 Fix color factors C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C REAL*8 DELTA COMPLEX*16 QSQ,ELLSQ,DOTQELL,ABSQ COMPLEX*16 Q(3),KG1(3),ELL(3),ELLP(3) COMPLEX*16 T(3,3) COMPLEX*16 PREFACTOR,SPECIAL,TEMP COMPLEX*16 COMPLEXSQRT INTEGER CUTSIGN REAL*8 XXREAL INTEGER MU,NU C C We define CUTSIGN to be +1 if the self-energy subgraph crosses C the final state cut with propagator P1 to the left and P2 to C the right, -1 if it is the other way around. This allows us C to choose the signs so that q^mu is flowing toward the cut. C IF (S1*XXREAL(KC(P1,0)).GT.0.0D0) THEN CUTSIGN = 1 ELSE CUTSIGN = -1 ENDIF C QSQ = (0.0D0,0.0D0) ELLSQ = (0.0D0,0.0D0) DOTQELL = (0.0D0,0.0D0) DO MU = 1,3 Q(MU) = S1*KC(P1,MU) KG1(MU) = SG1*KC(PG1,MU) ELL(MU) = KG1(MU) - 0.5D0*Q(MU) QSQ = QSQ + Q(MU)**2 ELLSQ = ELLSQ + ELL(MU)**2 DOTQELL = DOTQELL + Q(MU)*ELL(MU) ENDDO DO MU = 1,3 ELLP(MU) = ELL(MU) - Q(MU)*DOTQELL/QSQ DO NU = 1,3 T(MU,NU) = - DELTA(MU,NU) + Q(MU)*Q(NU)/QSQ ENDDO ENDDO C ABSQ = COMPLEXSQRT(QSQ) C PREFACTOR = 1.0D0/(32.0D0*ABSQ*COMPLEXSQRT(ELLSQ + MUMSBAR**2)**5) C C Now we need the proper special factor for each special case. C The first case is that of a cut self-energy, which should not C happen since this is the subroutine for renormalization. C In the second case, the self-energy subgraph is virtual, but C one of the adjoing propagators is cut. There are two such cuts C possible. For each of them, we have a factor 1/2 (and a minus sign). C In the third case, the self-energy subgraph and both of the C adjoining propagators is virtual. Here Q4SQ is the square of the C four-vector q^mu. C IF (CUT(PG1).AND.CUT(PG2)) THEN WRITE(NOUT,*)'Something is rotten in GPROPR.' STOP ELSE IF (CUT(P1).OR.CUT(P2)) THEN SPECIAL = (-0.5D0,0.0D0) ELSE WRITE(NOUT,*)'Something else is rotten in GPROPR.' STOP ENDIF C C Finally, we assemble the result. C TEMP = (8*NC - 4*NF)*ELLSQ + (6*NC - 4*NF)*MUMSBAR**2 RESULT(0,0) = (0.0D0,0.0D0) DO MU = 1,3 RESULT(0,MU) = (0.0D0,0.0D0) RESULT(MU,0) = (0.0D0,0.0D0) DO NU = 1,3 RESULT(MU,NU) = SPECIAL*PREFACTOR* > ( 4*(NC - NF)*ELLP(MU)*ELLP(NU) + TEMP*T(MU,NU) ) ENDDO ENDDO C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C Renormalization functions for three point subgraphs. C 15 September 1997 C 4 September 1998 Add color factors C 21 September 1998 Fix color factors C Notation: C h[index, internal lines, external lines] C C C h[a,{gluon,gluon,quark},{gluon,qbar,quark}] = RQQG3AG(KLOOP,MUMSBAR), C h[a,{gluon,quark,quark},{gluon,qbar,quark}] = RQQG3AQ(KLOOP,MUMSBAR), C h[b,{gluon,gluon,quark},{gluon,qbar,quark}] = RQQG3BG(KLOOP,MUMSBAR), C h[b,{gluon,quark,quark},{gluon,qbar,quark}] = RQQG3BQ(KLOOP,MUMSBAR), C h[c,{gluon,gluon,quark},{gluon,qbar,quark}] = RQQG3CG(KLOOP,MUMSBAR), C h[c,{gluon,quark,quark},{gluon,qbar,quark}] = RQQG3CQ(KLOOP,MUMSBAR), C h[a,{gluon,quark,quark}, {phot,qbar,quark}] = RQQP3AQ(KLOOP,MUMSBAR), C h[b,{gluon,quark,quark}, {phot,qbar,quark}] = RQQP3BQ(KLOOP,MUMSBAR), C h[c,{gluon,quark,quark}, {phot,qbar,quark}] = RQQP3CQ(KLOOP,MUMSBAR) C C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQP3AQ(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gamma vertex, coefficient of gamma[mu] C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQP3AQ = 1 /(2*OMEGA**3) + 3*MUMSBAR**2 /(8*OMEGA**5) RQQP3AQ = RQQP3AQ * (NC**2 -1)/(2*NC) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQP3BQ(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gamma vertex, coefficient of N[mu] N.gamma C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQP3BQ = - (NC**2 - 1) /(8*NC*OMEGA**3) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQP3CQ(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gamma vertex, coefficient of L3[mu] L3.gamma C The loop momentum here is just the space part, but the dot products C are the 4-D Minkowski dot products. C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQP3CQ = 3*(NC**2 - 1) /(8*NC*OMEGA**5) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQG3AG(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gluon vertex, gluon connects to gluon line, C coefficient of gamma[mu] C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQG3AG = NC /(4*OMEGA**3) + 3*NC*MUMSBAR**2 /(16*OMEGA**5) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQG3BG(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gluon vertex, gluon connects to gluon line, C coefficient of N[mu] N.gamma C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQG3BG = NC /(8*OMEGA**3) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQG3CG(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gluon vertex, gluon connects to gluon line, C coefficient of L3[mu] L3.gamma C The loop momentum here is just the space part, but the dot products C are the 4-D Minkowski dot products. C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQG3CG = - 3*NC /(8*OMEGA**5) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQG3AQ(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gluon vertex, gluon connects to quark line, C coefficient of gamma[mu] C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQG3AQ = - 1 /(4*NC*OMEGA**3) - 3*MUMSBAR**2 /(16*NC*OMEGA**5) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQG3BQ(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gluon vertex, gluon connects to quark line, C coefficient of N[mu] N.gamma C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQG3BQ = 1 /(8*NC*OMEGA**3) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION RQQG3CQ(KLOOP,MUMSBAR) C C In: COMPLEX*16 KLOOP(0:3) REAL*8 MUMSBAR C C q-qbar-gluon vertex, gluon connects to quark line, C coefficient of L3[mu] L3.gamma C The loop momentum here is just the space part, but the dot products C are the 4-D Minkowski dot products. C REAL*8 NC,NF COMMON /COLORFACTORS/ NC,NF C COMPLEX*16 LSQ,OMEGA COMPLEX*16 COMPLEXSQRT INTEGER MU C LSQ = 0.0D0 DO MU = 1,3 LSQ = LSQ + KLOOP(MU)**2 ENDDO OMEGA = COMPLEXSQRT(LSQ + MUMSBAR**2) C RQQG3CQ = - 3 /(8*NC*OMEGA**5) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C C Miscellaneous Functions C C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION XXREAL(Z) C COMPLEX*16 Z C XXREAL = DBLE(Z) RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION XXIMAG(Z) C COMPLEX*16 Z C COMPLEX*16 ZZ C ZZ = (0.0D0,-1.0D0) * Z XXIMAG = DBLE(ZZ) RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C COMPLEX*16 FUNCTION COMPLEXSQRT(Z) COMPLEX*16 Z C C Square root for complex numbers with error checking. C 1 February 1998 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT REAL*8 XXREAL,XXIMAG C IF (XXREAL(Z).LT.0.0D0) THEN IF(ABS(XXIMAG(Z)).LT. (1.0D-1 * ABS(XXREAL(Z))) ) THEN WRITE(NOUT,*)'Too near cut of Sqrt(Z) in COMPLEXSQRT(Z)' WRITE(NOUT,*) Z STOP ENDIF ENDIF COMPLEXSQRT = SQRT(Z) RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION FACTORIAL(N) C INTEGER N C INTEGER J C FACTORIAL = 1.0D0 DO J = 1,N FACTORIAL = FACTORIAL * J ENDDO RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION SINHINV(Z) C REAL*8 Z REAL*8 Z2,Z3,Z5,Z7 C C Evaluate arcsinh(z) = log( z + Sqrt(1+z^2) ). C For small z, we use the first four terms in the power C series expansion of this function so that we do not C lose precision in evaluating log(1 + z + z^2/2 + ...). C At z = 0.03, the series evaluation is accurate to a C fractional error of 2E-14. At this same point, the C precision of the Log form of the function should be C about that for representing 1.03 - 1.00, or 14 digits C minus 1.5 digits, or 3E-12. C IF (Z.LT.3.0D-2) THEN Z2 = Z**2 Z3 = Z * Z2 Z5 = Z3 * Z2 Z7 = Z5 * Z2 SINHINV = Z - 1.66666666667D-1 * Z3 > + 7.5D-2 * Z5 - 4.46428571429D-2 * Z7 ELSE SINHINV = LOG(Z + SQRT(1.0D0 + Z**2)) ENDIF RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** C*********************************************************************** C SUBROUTINE DIAGNOSTIC(K,BADGRAPHNUMBER) INTEGER SIZE PARAMETER (SIZE = 3) C In: REAL*8 K(0:3*SIZE-1,0:3) INTEGER BADGRAPHNUMBER C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX 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 MAP variables: INTEGER NMAPS INTEGER QS(256,0:SIZE) INTEGER A1S(256,2:SIZE),A2S(256,2:SIZE) CHARACTER*6 MAPTYPES(256) C Calculate variables: LOGICAL REPORT,DETAILS COMMON /CALCULOOK/ REPORT,DETAILS COMPLEX*16 VALUE COMPLEX*16 VALUECHK REAL*8 MAXPART C REAL*8 ABSK(0:3*SIZE-1) INTEGER P,MU,V REAL*8 COS12,COS23,COS31,SIN12,SIN23,SIN31 REAL*8 BADNESS C C This finds data for the kind of graph with the worst value C Latest revision: 4 January 1999 C ------- C C First we have to run NEWGRAPH through all the graphs so that it C initializes itself. C GRAPHFOUND = .TRUE. DO WHILE (GRAPHFOUND) CALL NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND) ENDDO C C Now we find the graph we wanted. C GRAPHFOUND = .TRUE. GRAPHNUMBER = 0 DO WHILE (GRAPHFOUND) CALL NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND) IF (GRAPHFOUND) THEN GRAPHNUMBER = GRAPHNUMBER + 1 ENDIF IF (GRAPHFOUND.EQ..FALSE.) THEN WRITE(NOUT,*)'Oops, snafu in DIAGNOSTIC' ENDIF IF (GRAPHNUMBER.EQ.BADGRAPHNUMBER) THEN GRAPHFOUND = .FALSE. ENDIF ENDDO 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) C WRITE(NOUT,*)' ' WRITE(NOUT,*)'Analysis by subroutine DIAGNOSTIC' WRITE(NOUT,*)' ' WRITE(NOUT,702)GRAPHNUMBER 702 FORMAT('Graph number ',I3) WRITE(NOUT,*)'Point:' DO P=1,8 WRITE(NOUT,703) P,K(P,1),K(P,2),K(P,3) 703 FORMAT('P =',I2,' K = ',3(1P G12.3)) ENDDO C WRITE(NOUT,*)' ' WRITE(NOUT,*)'Softness:' DO P = 1,NPROPS ABSK(P) = 0.0D0 DO MU = 1,3 ABSK(P) = ABSK(P) + K(P,MU)**2 ENDDO ABSK(P) = SQRT(ABSK(P)) WRITE(NOUT,704) P,ABSK(P) 704 FORMAT('P =',I2,' |K| = ',1P G12.3) ENDDO C WRITE(NOUT,*)' ' WRITE(NOUT,*)'Collinearity:' DO V = 3,NVERTS COS12 = 0.0D0 COS23 = 0.0D0 COS31 = 0.0D0 DO MU = 1,3 COS12 = COS12 + K(PROP(V,1),MU) * K(PROP(V,2),MU) COS23 = COS23 + K(PROP(V,2),MU) * K(PROP(V,3),MU) COS31 = COS31 + K(PROP(V,3),MU) * K(PROP(V,1),MU) ENDDO COS12 = COS12 /ABSK(PROP(V,1))/ABSK(PROP(V,2)) COS23 = COS23 /ABSK(PROP(V,2))/ABSK(PROP(V,3)) COS31 = COS31 /ABSK(PROP(V,3))/ABSK(PROP(V,1)) SIN12 = SQRT(1.0D0 - COS12**2) SIN23 = SQRT(1.0D0 - COS23**2) SIN31 = SQRT(1.0D0 - COS31**2) WRITE(NOUT,705)V,PROP(V,1),PROP(V,2),PROP(V,3), > SIN12,SIN23,SIN31 705 FORMAT('V =',I2,' Ps =',3I2,' sines =', 3F10.5) ENDDO C WRITE(NOUT,*)' ' CALL CHECKPOINT(K,ABSK,PROP,BADNESS) WRITE(NOUT,706)BADNESS 706 FORMAT('Badness of this point is',1P G10.2) C WRITE(NOUT,*)' ' WRITE(NOUT,*)'CALCULATE finds the folowing:' REPORT = .TRUE. C CALL CALCULATE(VRTX,SELFPROP,GRAPHNUMBER,K,ABSK, > QS,A1S,A2S,MAPTYPES,NMAPS,VALUE,MAXPART,VALUECHK) C WRITE(NOUT,707)VALUE, ABS(VALUE),MAXPART 707 FORMAT('VALUE =',2(1P G12.4),' ABS(VALUE) = ',1P G12.4,/, > 'BIGGEST contribution was ',1P G12.4) RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** C Random Number Generator C*********************************************************************** C SUBROUTINE RANDOMINIT(IRAN) INTEGER IRAN C C Code from CERN, 1991. C Modified 19 Dec 1995 to replace .AND.,.OR. and .XOR. by C IAND(,), IOR(,), and IXOR(,) respectively. C C Initialize the shift-register random number generator from a random C seed IRAN, 0 <= IRAN < 259200, with the help of a portable "quick C and dirty" generator C INTEGER LBIT REAL*8 FAC PARAMETER(LBIT=31,FAC= 2.0D0**(-LBIT)) C Here LBIT is the number of bits used in the shift register generator INTEGER J,IR(250) REAL*8 RR(250) COMMON/RANDO/ RR,J,IR SAVE /RANDO/ INTEGER I,LB,JRAN,IFAC1,ISK,IDI,I1S C Configuration of the shift register generator INTEGER IM,IA,IC DATA IM,IA,IC /259200,7141,54773/ C C IF(IRAN.LT.0.OR.IRAN.GE.IM)STOP > 'RINI: IRAN OUT OF RANGE' JRAN=IRAN C Warm up the auxiliary generator a little DO I=1,10 JRAN=MOD(JRAN*IA+IC,IM) ENDDO IFAC1=((2**(LBIT-1)-1)*2+1)/IM DO I=2,250 JRAN=MOD(JRAN*IA+IC,IM) IR(I)=IFAC1*JRAN ENDDO C Guarantee LBIT linearly independent (over the field of 2 el.) C elements in IR(I): IR(1)=1 I=1 ISK=250/LBIT IDI=1 I1S=1 DO LB=1,LBIT-1 I=I+ISK IDI=IDI*2 I1S=I1S+IDI IR(I)=IR(I).OR.IDI IR(I)=IR(I).AND.I1S CM IR(I)=IOR(IR(I),IDI) CM IR(I)=IAND(IR(I),I1S) ENDDO CALL NEWRAN END C C*********************************************************************** C SUBROUTINE NEWRAN C C Code from CERN, 1991. C Modified 19 Dec 1995 to replace .AND.,.OR. and .XOR. by C IAND(,), IOR(,), and IXOR(,) respectively. C C Fills IR(I),RR(I) with 250 new random numbers, resets J=0. C Increment J before use! C INTEGER LBIT REAL*8 FAC PARAMETER(LBIT=31,FAC= 2.0D0**(-LBIT)) INTEGER J,IR(250) REAL*8 RR(250) COMMON/RANDO/ RR,J,IR SAVE /RANDO/ C INTEGER N C DO N=1,103 IR(N)=IR(N+147).XOR.IR(N) CM IR(N)=IXOR(IR(N+147),IR(N)) RR(N)=FAC*(DBLE(IR(N)) + 0.5D0) ENDDO C DO N=104,250 IR(N)=IR(N-103).XOR.IR(N) CM IR(N)=IXOR(IR(N-103),IR(N)) RR(N)=FAC*(DBLE(IR(N)) + 0.5D0) ENDDO C J=0 END C C*********************************************************************** C REAL*8 FUNCTION RANDOM(DUMMY) INTEGER DUMMY C C Code from CERN, 1991. C C Random number between 2**(-32) and 1 - 2**(-32): C INTEGER J,IR(250) REAL*8 RR(250) COMMON/RANDO/ RR,J,IR SAVE /RANDO/ C IF(J.GE.250)CALL NEWRAN J=J+1 RANDOM = RR(J) END C C*********************************************************************** C END OF LIBRARY ROUTINES FOR E+E- CALCULATION C***********************************************************************