C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C C --------------------- C beowulfsubs.f Version 0.7 C --------------------- C C23456789012345678901234567890123456789012345678901234567890123456789012 SUBROUTINE VERSION INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT WRITE(NOUT,*)'Beowulf 0.7 subroutines 9 April 1998' RETURN END C23456789012345678901234567890123456789012345678901234567890123456789012 C C Subroutines for numerical integration of jet cross sections in C electron-positron annihilation. C C D. E. Soper C C First code: 29 November 1992 C Latest 1993 revision: 27 September 1993 C Revision: 2 August 1994 C Multiple revisions: 1996 C Latest revision: see Version subroutine above 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(ABSK,Q,MAPTYPE,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 Results: C For each point, subroutine CALCULATE returns C VALUE = value of the integrand at that point (including jacobian); C BIGGEST = largest of the absolute values of the C contributions to VALUE. C The results VALUE are summed in C INTEGRAL for those points with BIGGEST < CANCELLIMIT * ABS(VALUE). C EXTRA for those points with C CANCELLIMIT * ABS(VALUE) < BIGGEST < 10 * CANCELLIMIT * ABS(VALUE). C The real and imaginary parts of INTEGRAL are summed in SUMR,SUMI. C The real and imaginary parts of EXTRA are summed in SUMXR,SUMXI. C The errors for SUMR and SUMI are summed in SQRSUMR,SQRSUMI and C reported in ERROR and ERRORIMAG. C C23456789012345678901234567890123456789012345678901234567890123456789012 C Subroutines C C Subroutine RENO(NRENO,SUMR,ERROR,SUMI,ERRORIMAG,SUMXR,SUMXI): C subroutine NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND) C subroutine NEWCHOICE(C,COUNT,V,FAIL) C subroutine NEXTCHOICE(C,COUNT,V,FAIL) C subroutine CHECKSTATUS(C,COUNT) C function ONEPI(C) C subroutine CHECK(C,NPERMS,OK) C subroutine CHECKOUT(C,CIN,NPERMS,OK) C subroutine EXCHANGE(V1,V2,C) C subroutine FINDA(VRTX,Q,NQ,A,QOK) C function PROPSIGN(VRTX,P,V) C subroutine FINDQS(VRTX,PROP,NMAPS,QS,A1S,A2S,MAPTYPES) C function DEPENDENT(VRTX,Q,NQ) C subroutine SOFTHELP(Q1,VRTX,PROP,P1,P2,SIGN1,SIGN2,YES) C uses NEWCUT(...,CUTINDEX,...,LOOPINDEX,...,CUTFOUND) C subroutine NEWPOINT(Q,A,VRTX,PROP,MAPTYPE,A1,A2,K,ABSK) C subroutine AXES(EA,EB,EC) C function RANDOM(1) C subroutine CHECKPOINT(K,ABSK,PROP,BADNESS) C subroutine CALCULATE(VRTX,SELFPROP,GRAPHNUMBER,K,ABSK, C QS,A1S,A2S,MAPTYPES,NMAPS,VALUE,BIGGEST,CHECKVALUE) C subroutine REFLECT(VRTX,ABSK,K,NEWABSK,NEWK,NREFLECT) C subroutine MOSTCOLLINEAR(VRTX,ABSK,K,VSTAR,P1,P2,P3) C function DENSITY(K,ABSK,QS,A1S,A2S,MAPTYPES,NMAPS) C subroutine NEWCUT(VRTX,NEWCUTINIT,NCUT,ISIGN, C CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP, C NINLOOP,LOOPINDEX,LOOPSIGN,CUTFOUND) C subroutine DEFORM(VRTX,LOOPINDEX,RTS,LEFTLOOP, C RIGHTLOOP,NINLOOP,KINLOOP,NEWKINLOOP,JACDEFORM) C subroutine CHECKDEFORM(P,KCSQ, C NINLOOP,LEFTLOOP,RIGHTLOOP,LOOPINDEX,JCUT) C function NUMERATOR(GRAPHNUMBER,KC,CUT) C subroutine QPROP(CUT,KC,P1,P2,PQ,PG,S1,SQ,RESULT) C subroutine PROP(CUT,KC,P1,P2,PG1,PG2,S1,SG1,RESULT) C function RNUMERATOR(GRAPHNUMBER,KC,MUMSBAR,CUT) C subroutine QPROPR(CUT,MUMSBAR, C KC,P1,P2,PQ,PG,S1,SQ,RESULT) C subroutine GPROPR(CUT,MUMSBAR, C KC,P1,P2,PG1,PG2,S1,SG1,RESULT) C functions RQQG3AG(KLOOP,MUMSBAR) ... C function CALS0(KCUT,NCUT) C function MOCKDELTA(Z,Z0) C function FEXP(Y) C function SMEAR(RTS) C function FACTORIAL(N) C C function XXREAL(Z) C function XXIMAG(Z) C - plus random number routines C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE RENO( > SUMR,ERRORR,SUMI,ERRORI, > SUMBISR,ERRORBISR,SUMBISI,ERRORBISI, > SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI, > INCLUDED,EXTRACOUNT,OMITTED, > NVALUE,VALUEMAX,KBAD,BADGRAPHNUMBER) C INTEGER SIZE PARAMETER (SIZE = 3) C Out: REAL*8 SUMR,ERRORR,SUMI,ERRORI REAL*8 SUMBISR,ERRORBISR,SUMBISI,ERRORBISI REAL*8 SUMCHKR,ERRORCHKR,SUMCHKI,ERRORCHKI INTEGER INCLUDED,EXTRACOUNT,OMITTED INTEGER NVALUE(-9:6) REAL*8 VALUEMAX REAL*8 KBAD(0:3*SIZE-1,0:3) INTEGER BADGRAPHNUMBER C C Computes the cross section integral by Monte Carlo integration. C C Latest revision 5 March 1998 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 C Graphs to include LOGICAL USEGRAPH(10) COMMON /WHICHGRAPHS/ USEGRAPH 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 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 variables: INTEGER GROUPSIZE,NRENO COMMON /MONTECARLO/GROUPSIZE,NRENO INTEGER GROUP,POINT C Reno results variables: REAL*8 SQRSUMR,SQRSUMBISR,SQRSUMCHKR REAL*8 SQRSUMI,SQRSUMBISI,SQRSUMCHKI COMPLEX*16 INTEGRAL,INTEGRALBIS,INTEGRALCHK C Calculate variables: COMPLEX*16 VALUE,VALUECHK REAL*8 MAXPART REAL*8 LOGVALUE LOGICAL REPORT,DETAILS COMMON /CALCULOOK/ REPORT,DETAILS 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 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 SUMBISR = 0.0D0 SUMBISI = 0.0D0 SUMCHKR = 0.0D0 SUMCHKI = 0.0D0 C SQRSUMR = 0.0D0 SQRSUMI = 0.0D0 SQRSUMBISR = 0.0D0 SQRSUMBISI = 0.0D0 SQRSUMCHKR = 0.0D0 SQRSUMCHKI = 0.0D0 C DO L = -9,6 NVALUE(L) = 0 ENDDO VALUEMAX = 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,0.0D0) INTEGRALCHK = (0.0D0,0.0D0) C C Loop over groups of points. C DO GROUP = 1,NRENO C C Call Hrothgar to tell him to that we are starting a new group. C CALL HROTHGAR(1,KCUT0,1.0D0,'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) 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 C C Call Hrothgar to tell him that we are starting a new point. C CALL HROTHGAR(1,KCUT0,1.0D0,'STARTPOINT') C C Get a new point. Check on its badness. If it is too bad, we omit C the point after notifying Hrothgar. C CALL NEWPOINT(Q,A,VRTX,PROP,MAPTYPE,A1,A2,K,ABSK) CALL CHECKPOINT(K,ABSK,PROP,BADNESS) IF (BADNESS.GT.BADNESSLIMIT) THEN OMITTED = OMITTED + 1 CALL HROTHGAR(1,KCUT0,1.0D0,'BADPOINT ') ELSE C C We know that the badness is not too bad. C Calculate the corresponding contribution. 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 CALL CALCULATE(VRTX,SELFPROP,GRAPHNUMBER,K,ABSK, > QS,A1S,A2S,MAPTYPES,NMAPS,VALUE,MAXPART,VALUECHK) C C Add contribution from this point to integral. C We count the point if Maxvalue/|Value| < Cancellimit. C IF ( MAXPART.LT.CANCELLIMIT*ABS(XXREAL(VALUE)) ) THEN INTEGRAL = INTEGRAL + VALUE INTEGRALCHK = INTEGRALCHK + VALUECHK INCLUDED = INCLUDED + 1 LOGVALUE = LOG10(ABS(XXREAL(VALUE))) DO L = -9,6 IF((LOGVALUE.GE.L).AND.(LOGVALUE.LT.(L+1))) THEN NVALUE(L) = NVALUE(L) + 1 ENDIF ENDDO IF (ABS(XXREAL(VALUE)).GT.VALUEMAX) THEN VALUEMAX = ABS(XXREAL(VALUE)) DO P = 1,NPROPS DO MU = 1,3 KBAD(P,MU) = K(P,MU) ENDDO ENDDO BADGRAPHNUMBER = GRAPHNUMBER ENDIF ELSE IF (MAXPART.LT.1.0D2*CANCELLIMIT*ABS(XXREAL(VALUE))) THEN INTEGRALBIS = INTEGRALBIS + VALUE INTEGRALCHK = INTEGRALCHK + VALUECHK EXTRACOUNT = EXTRACOUNT + 1 CALL HROTHGAR(1,KCUT0,1.0D0,'XTRAPOINT ') ELSE OMITTED = OMITTED + 1 CALL HROTHGAR(1,KCUT0,1.0D0,'BADPOINT ') ENDIF C C End IF (BADNESS.GT.BADNESSLIMIT) ... ELSE: C ENDIF C C End of loop over POINT. C CALL HROTHGAR(1,KCUT0,1.0D0,'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,'GROUPDONE ') C C Add results from this group to the SUM variables. C INTEGRAL = INTEGRAL/GROUPSIZE INTEGRALBIS = INTEGRALBIS/GROUPSIZE INTEGRALCHK = INTEGRALCHK/GROUPSIZE C SUMR = SUMR + XXREAL(INTEGRAL) SUMI = SUMI + XXIMAG(INTEGRAL) SUMBISR = SUMBISR + XXREAL(INTEGRALBIS) SUMBISI = SUMBISI + XXIMAG(INTEGRALBIS) SUMCHKR = SUMCHKR + XXREAL(INTEGRALCHK) SUMCHKI = SUMCHKI + XXIMAG(INTEGRALCHK) C SQRSUMR = SQRSUMR + XXREAL(INTEGRAL)**2 SQRSUMI = SQRSUMI + XXIMAG(INTEGRAL)**2 SQRSUMBISR = SQRSUMBISR + XXREAL(INTEGRALBIS)**2 SQRSUMBISI = SQRSUMBISI + XXIMAG(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 GROUP = 1,NRENO C ENDDO C C Calculate the SUM results. C SUMR = SUMR/NRENO SUMI = SUMI/NRENO SUMBISR = SUMBISR/NRENO SUMBISI = SUMBISI/NRENO SUMCHKR = SUMCHKR/NRENO SUMCHKI = SUMCHKI/NRENO C SQRSUMR = SQRSUMR/NRENO SQRSUMI = SQRSUMI/NRENO SQRSUMBISR = SQRSUMBISR/NRENO SQRSUMBISI = SQRSUMBISI/NRENO SQRSUMCHKR = SQRSUMCHKR/NRENO SQRSUMCHKI = SQRSUMCHKI/NRENO C ERRORR = SQRT((SQRSUMR - SUMR**2)/(NRENO - 1)) ERRORI = SQRT((SQRSUMI - SUMI**2)/(NRENO - 1)) ERRORBISR = SQRT((SQRSUMBISR - SUMBISR**2)/(NRENO - 1)) ERRORBISI = SQRT((SQRSUMBISI - SUMBISI**2)/(NRENO - 1)) ERRORCHKR = SQRT((SQRSUMCHKR - SUMCHKR**2)/(NRENO-1)) ERRORCHKI = SQRT((SQRSUMCHKI - SUMCHKI**2)/(NRENO-1)) C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C Subroutines associated with NEWGRAPH C C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE NEWGRAPH(VRTX,PROP,SELFPROP,GRAPHFOUND) c implicit none 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 C CALL CHECKSTATUS(C,COUNT) 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 implicit none 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 implicit none 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 SUBROUTINE CHECKSTATUS(C,COUNT) c implicit none INTEGER SIZE PARAMETER (SIZE = 3) C INTEGER C(2*SIZE,3),COUNT(2*SIZE) C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX C INTEGER CHECKSUM INTEGER NCONN(2*SIZE,2*SIZE) INTEGER V,V1,V2,I SAVE C DO V1 = 1,NVERTS DO V2 = 1,NVERTS NCONN(V1,V2) = 0 ENDDO ENDDO C DO V1 = 1,NVERTS DO I = 1,3 V2 = C(V1,I) IF (V2.GT.0) THEN NCONN(V1,V2) = NCONN(V1,V2) + 1 ENDIF ENDDO ENDDO C DO V = 1,NVERTS IF(NCONN(V,V).NE.0) THEN WRITE(NOUT,*)'NCONN(V,V).NE.0' STOP ENDIF ENDDO C DO V1 = 1,NVERTS CHECKSUM = 0 IF ((V1.EQ.1).OR.(V1.EQ.2)) THEN CHECKSUM = 1 ENDIF DO V2 = 1,NVERTS IF (NCONN(V1,V2).NE.NCONN(V2,V1)) THEN WRITE(NOUT,*)'NCONN(V1,V2).NE.NCONN(V2,V1)' STOP ENDIF CHECKSUM = CHECKSUM + NCONN(V1,V2) ENDDO IF(COUNT(V1).NE.CHECKSUM) THEN WRITE(NOUT,*)'COUNT(V1).NE.CHECKSUM' STOP ENDIF ENDDO C DO V = 1,NVERTS DO I = 1,2 IF ((C(V,I).GT.C(V,I+1)).AND.(C(V,I+1).NE.0)) THEN WRITE(NOUT,*)'C(V,I).GT.C(V,I+1)' STOP ENDIF ENDDO ENDDO C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C LOGICAL FUNCTION ONEPI(CIN) c implicit none 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 implicit none 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 implicit none 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 implicit none 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 implicit none 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 implicit none 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 implicit none 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 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 uniformly spread in angle. 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 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. Actually, some C of the tests are performed twice, so the subroutine could run faster C if the results of these tests were recorded to be used later. 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 uniformly spread in angle. 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) The IRPART case gives a soft map without special concentration C of points in a plane. We want this when there is no C cut for which propagator Q(1) is part of a virtual loop C and connects to cut lines. In this case subrountine SOFTHELP C returns YES = .FALSE. C 3) The normal IRFULL case requires a concentration of points near C 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 4) 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 CALL SOFTHELP(Q1,VRTX,PROP,P1,P2,SIGN1,SIGN2,YES) IF (YES) THEN 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 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 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 NMAPS = NMAPS + 1 MAPTYPES(NMAPS) = 'UVHARD' DO L = 2,NLOOPS A1S(NMAPS,L) = 0 A2S(NMAPS,L) = 0 ENDDO QS(NMAPS,0) = 0 QS(NMAPS,NLOOPS) = LOOPINDEX(1) DO L = 1,NLOOPS-1 QS(NMAPS,L) = CUTINDEX(L) 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 implicit none 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 implicit none 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) c implicit none 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) 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 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,CEPS COMMON /POINTS/ AS,BS,AM,BM,AH,BH,AUV,BUV,CEPS C REAL*8 PI PARAMETER (PI = 3.141592653589793239D0) C REAL*8 MAGK,EPSILON,GAMMA,Z REAL*8 TEMP1,TEMP2 REAL*8 DOT12,COS12 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 X,RANDOM REAL*8 TEMP,KSQ,PREVIOUSK INTEGER L,P,MU 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 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 ELL(NLOOPS,1) = MAGK * SINTHETA * COS(PHI) ELL(NLOOPS,2) = MAGK * SINTHETA * SIN(PHI) ELL(NLOOPS,3) = MAGK * COSTHETA 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 uniformly spread in angle. 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 For the other k_{Q(i)} except K_{Q(1)}, we use a nested "soft" map C with AS = 1, BS = 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**(-AS/BS) - 1.0D0 )**(1.0D0/AS) 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.'IRFULL') 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. C If NLOOPS is 1 or 2, the DO loop will not be C executed. For NLOOPS = 3, this evaluates ELL(2,mu). C PREVIOUSK = 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) 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 C AS, BS map, 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 DOT12 = 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 DOT12 = DOT12 + TEMP1*TEMP2 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) COS12 = DOT12/(K1ABS*K2ABS) DO MU = 1,3 K1HAT(MU) = K1(MU)/K1ABS K2HAT(MU) = K2(MU)/K2ABS ENDDO C NORMA = SQRT(2.0D0*(1.0D0 - COS12)) IF ((NORMA.LT.1.0D-30)) THEN WRITE(NOUT,*)'NORMA too small in NEWPOINT' STOP 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 Here we use arcsinh(z) = log( z + Sqrt(1+z^2) ). C X = RANDOM(1) MAGK = PREVIOUSK * ( X**(-AS/BS) - 1.0D0 )**(1.0D0/AS) EPSILON = CEPS * MAGK**2/PREVIOUSK**2 Z = 1.0D0/EPSILON GAMMA = LOG(Z + SQRT(1.0D0 + Z**2)) 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 IF (ABS(SINH(GAMMA)/Z-1.0D0).GT.1.0D-10)THEN IF(ABS(Z).GT.0.0001D0) THEN WRITE(NOUT,*)'Something is wrong in NEWPOINT',SINH(GAMMA),Z STOP ENDIF ENDIF 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 Now we generate the maps for nearly UV divergent loops. They 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 This completes the IF (MAPTYPE.NE.'UVHARD') THEN ... C ELSE C C The AH = 3, BH = 3 (normally) map for loop momenta 1,..., NLOOPS -1 C DO L = 1,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 ELL(L,1) = MAGK * SINTHETA * COS(PHI) ELL(L,2) = MAGK * SINTHETA * SIN(PHI) ELL(L,3) = MAGK * COSTHETA ENDDO 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 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 implicit none 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 + Kmid - Kmax )/Kmax C C is the 1/badness^2 for that vertex. The badness for the point is the C product of the badness values of all the vertices. C C Revised 27 December 1996 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX C REAL*8 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 SMALLNESS = SMALLNESS * (KMIN + KMID - KMAX)/KMAX ENDDO IF (SMALLNESS.LT.1.0D-30) THEN BADNESS = 1.0D30 ELSE BADNESS = SQRT(1.0D0/SMALLNESS) ENDIF RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE AXES(EA,EB,EC) c implicit none 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 c implicit none 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 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 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 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 LOGICAL NEWCUTINIT,CUTFOUND 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 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---------------------------------- 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 Also get a reflected point (to further kill collinear singularities). C If there is a reflection, REFLECT returns NREFLECT = 2. Else it is 1. C CALL REFLECT(VRTX,ABSK,K,ABSKREFLECTED,KREFLECTED,NREFLECT) C C Calculate jacobian. C RHO = DENSITY(K,ABSK,QS,A1S,A2S,MAPTYPES,NMAPS) IF (NREFLECT.EQ.2) THEN RHO = RHO + > DENSITY(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 Get a new cut. We set NEWCUTINIT to .TRUE. to tell NEWCUT to C initialize itself. Then NEWCUT will set this variable to .FALSE. C so that it remembers its state on successive calls in this DO loop. C CUTFOUND = .TRUE. NEWCUTINIT = .TRUE. DO WHILE (CUTFOUND) CALL NEWCUT(VRTX,NEWCUTINIT,NCUT,ISIGN, > CUTINDEX,CUTSIGN,LEFTLOOP,RIGHTLOOP, > NINLOOP,LOOPINDEX,LOOPSIGN,CUTFOUND) IF (CUTFOUND) THEN 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(KCUT,NCUT) PREFACTOR = 1.0D0 / (3.0D0 * 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) 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 C Just to be sure, let's check the direction of deformation also for C the uncut propagator in a virtual self-energy loop. C C IF ( SELFPROP(P).AND.INLOOP(P).AND.(.NOT.LOOPCUT(P)) ) THEN C KCSQ = 0.0D0 C DO MU = 0,3 C KCSQ = KCSQ + METRIC(MU) * KC(P,MU)**2 C ENDDO C CALL CHECKDEFORM(KCSQ,LEFTLOOP,RIGHTLOOP) C ENDIF C 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,'NEWRESULT ') C C End of loop DO WHILE (CUTFOUND)/ IF (CUTFOUND). C ENDIF 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 c implicit none 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]^(-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 INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT 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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(6,MU)*KC(6,MU) TEMP3 = TEMP3 + KC(7,MU)*KC(7,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU) TEMP3 = TEMP3 + KC(8,MU)*KC(8,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(8,MU)*KC(8,MU) TEMP3 = TEMP3 + KC(5,MU)*KC(5,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(3,MU)*KC(3,MU) TEMP3 = TEMP3 + KC(7,MU)*KC(7,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(3,MU)*KC(3,MU) TEMP2 = TEMP2 + KC(7,MU)*KC(7,MU) TEMP3 = TEMP3 + KC(8,MU)*KC(8,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU) TEMP3 = TEMP3 + KC(7,MU)*KC(7,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(5,MU)*KC(5,MU) TEMP3 = TEMP3 + KC(8,MU)*KC(8,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(6,MU)*KC(6,MU) TEMP3 = TEMP3 + KC(8,MU)*KC(8,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(6,MU)*KC(6,MU) TEMP3 = TEMP3 + KC(4,MU)*KC(4,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(1,MU)*KC(1,MU) TEMP2 = TEMP2 + KC(4,MU)*KC(4,MU) TEMP3 = TEMP3 + KC(8,MU)*KC(8,MU) ENDDO ELSE RETURN ENDIF C ELSE WRITE(NOUT,*)'Problem with graph number in CHECKCALC' STOP ENDIF C CHECK = (TEMP1 + MM**2)**2 CHECK = CHECK * (TEMP2 + MM**2)**2 CHECK = CHECK * (TEMP3 + MM**2)**2 CHECK = (MM/PI**2)**3 /CHECK CHECK = JACDEFORM * JACNEWPOINT * CHECK C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE REFLECT(VRTX,ABSK,K,NEWABSK,NEWK,NREFLECT) c implicit none 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 transformation K --> NEWK is 'valid' if the most collinear C vertex NEWVSTAR constructed from NEWK is the same as the original most C collinear vertex VSTAR. 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 INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX C INTEGER P,L,MU INTEGER VSTAR,P1,P2,P3 INTEGER NEWVSTAR,NEWP1,NEWP2,NEWP3 INTEGER Q(0:SIZE) REAL*8 ABSKBEST INTEGER PBEST INTEGER A(0:3*SIZE-1,0:SIZE) REAL*8 TEMP,DOT12 REAL*8 FACTOR REAL*8 KSQUARE LOGICAL DEPENDENT,LOOKMORE,QOK C C NREFLECT set to 1. C NREFLECT = 1 C IF (NLOOPS.LT.2) RETURN C CALL MOSTCOLLINEAR(VRTX,ABSK,K,VSTAR,P1,P2,P3) C C Now we have the most collinear vertex VSTAR, attached to C propagators P1,P2,P3. We set Q(1) = P1 and Q(2) = P2, C then 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 Q(0) = 0 Q(1) = P1 Q(2) = P3 DO L = 3,NLOOPS C C First, we need to find at least one possible choice. C LOOKMORE = .TRUE. P=0 DO WHILE (LOOKMORE) P = P+1 IF(P.GT.NPROPS) THEN WRITE(NOUT,*)'Failed to find independent props in REFLECT' 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. C LOOKMORE = .TRUE. DO WHILE (LOOKMORE) P = P+1 IF(P.GT.NPROPS) THEN LOOKMORE = .FALSE. ELSE IF (ABSK(P).LT.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 Close loop over L = 3,NLOOPS. C Q(L) = PBEST ENDDO C C Now we have (Q(1),...,Q(NLOOPS)) and we can constructd NEWK. 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 we change K(Q(2),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 most collinear C vertex NEWVSTAR constructed from NEWL is the same as the original most C collinear vertex VSTAR. C CALL MOSTCOLLINEAR(VRTX,NEWABSK,NEWK,NEWVSTAR,NEWP1,NEWP2,NEWP3) C IF (NEWVSTAR.EQ.VSTAR) THEN NREFLECT = 2 ELSE NREFLECT = 1 ENDIF C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE MOSTCOLLINEAR(VRTX,ABSK,K,VSTAR,P1STAR,P2STAR,P3STAR) c implicit none 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 VSTAR,P1STAR,P2STAR,P3STAR C C 18 July 1994 C 19 Dec 1995 C INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX C INTEGER V INTEGER P,P1,P2,P3 INTEGER MU LOGICAL NONESOFAR INTEGER MATCHES INTEGER PROPSIGN REAL*8 COS23,COSMAX,MAXSOFAR,MINSOFAR 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 ABSOLDK(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 Calculate cos( theta_{23} ). This quantity is nearly 1 for C a highly collinear splitting. Compare to the previous C highest value of COS23, corresponding to vertex VSTAR. 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 IF (COS23.GT.COSMAX) THEN C C Make sure that vertex V is not equivalent to VSTAR. C MATCHES = 0 IF( (VRTX(P1,1).EQ.VSTAR).OR.(VRTX(P1,2).EQ.VSTAR) ) THEN MATCHES = MATCHES + 1 ENDIF IF( (VRTX(P2,1).EQ.VSTAR).OR.(VRTX(P2,2).EQ.VSTAR) ) THEN MATCHES = MATCHES + 1 ENDIF IF( (VRTX(P3,1).EQ.VSTAR).OR.(VRTX(P3,2).EQ.VSTAR) ) THEN MATCHES = MATCHES + 1 ENDIF IF (MATCHES.LT.2) THEN C C Get to here if COS23.GT.COSMAX and V is not equivalent to VSTAR. C VSTAR = V P1STAR = P1 P2STAR = P2 P3STAR = P3 COSMAX = COS23 C C Close IF (COS23.GT.COSMAX) and IF (MATCHES.LT.2). C ENDIF ENDIF C C End DO V = 3,NVERTS. C ENDDO C RETURN END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C REAL*8 FUNCTION DENSITY(K,ABSK,QS,A1S,A2S,MAPTYPES,NMAPS) c implicit none INTEGER SIZE PARAMETER (SIZE = 3) C In: 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 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,CEPS COMMON /POINTS/ AS,BS,AM,BM,AH,BH,AUV,BUV,CEPS C INTEGER Q(0:SIZE),A1(2:SIZE),A2(2:SIZE) CHARACTER*6 MAPTYPE REAL*8 MAGK,K0,PREVIOUSK REAL*8 K1SQ,K2SQ,K1ABS,K2ABS,DOT12,COS12,TEMP1,TEMP2,NORMA REAL*8 K1(3),K2(3),K1HAT(3),K2HAT(3),EA(3) REAL*8 EPSILON,Z,COSTHETA,GAMMA INTEGER MAPNUMBER,L,MU REAL*8 J REAL*8 KRATIO REAL*8 FOURPI PARAMETER (FOURPI = 12.5663706143591730D0) C IF (NLOOPS.LT.1) THEN WRITE(NOUT,*) 'NLOOPS less than 1 in DENSITY' STOP ENDIF 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 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) 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 Q(1) soft singularity. C IRPART for infrared map with Q(1) points uniformly spread in angle. 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 "soft" 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. 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 /BS * K0**3 * KRATIO**(3.0D0 - AS) > * ( 1.0D0 + KRATIO**AS )**(1.0D0 + BS/AS) PREVIOUSK = MAGK ENDDO C ELSE IF (MAPTYPE.EQ.'IRFULL') 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. C If NLOOPS is 1 or 2, the DO loop will not be C executed. For NLOOPS = 3, this evaluates the jacobian for k_{Q(2)}. C PREVIOUSK = ENERGYSCALE DO L = NLOOPS-1,2,-1 MAGK = ABSK(Q(L)) K0 = PREVIOUSK KRATIO = MAGK/K0 J = J * FOURPI /BS * K0**3 * KRATIO**(3.0D0 - AS) > * ( 1.0D0 + KRATIO**AS )**(1.0D0 + BS/AS) PREVIOUSK = MAGK ENDDO C C AS, BS map, 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 DOT12 = 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 DOT12 = DOT12 + TEMP1*TEMP2 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) COS12 = DOT12/(K1ABS*K2ABS) DO MU = 1,3 K1HAT(MU) = K1(MU)/K1ABS K2HAT(MU) = K2(MU)/K2ABS ENDDO C NORMA = SQRT(2.0D0*(1.0D0 - COS12)) IF ((NORMA.LT.1.0D-30)) 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 Here we use arcsinh(z) = log( z + Sqrt(1+z^2) ). C EPSILON = CEPS * (MAGK/K0)**2 COSTHETA = 0.0D0 DO MU = 1,3 COSTHETA = COSTHETA + K(Q(1),MU) * EA(MU) ENDDO COSTHETA = COSTHETA/MAGK Z = 1.0D0/EPSILON GAMMA = LOG(Z + SQRT(1.0D0 + Z**2)) 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 NEWPOINT' STOP ENDIF 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 This completes the IF (MAPTYPE.NE.'UVHARD') THEN ... C ELSE C C The AH = 3, BH = 3 (normally) map for loop momenta 1,..., NLOOPS -1 C K0 = ENERGYSCALE DO L = 1,NLOOPS-1 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) ENDDO C C The AUV = 3, BUV = 1 map for loop momentum NLOOPS. C 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 End of IF (MAPTYPE.NE.'UVHARD') THEN ... ELSE ... C ENDIF C C End of DO MAPNUMBER = 1,NMAPS C DENSITY = DENSITY + 1.0D0/J ENDDO RETURN C END C C23456789012345678901234567890123456789012345678901234567890123456789012 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE NEWCUT(XVRTX,NEWCUTINIT,XNCUT,XISIGN, > XCUTINDEX,XCUTSIGN,XLEFTLOOP,XRIGHTLOOP, > XNINLOOP,XLOOPINDEX,XLOOPSIGN,CUTFOUND) C c implicit none 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 C23456789012345678901234567890123456789012345678901234567890123456789012 C SUBROUTINE DEFORM(VRTX,LOOPINDEX,RTS,LEFTLOOP,RIGHTLOOP, > NINLOOP,KINLOOP,NEWKINLOOP,JACDEFORM) C C implicit none 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 INTEGER NIN,NOUT COMMON /IOUNIT/ NIN,NOUT INTEGER NLOOPS,NPROPS,NVERTS,CUTMAX COMMON /SIZES/ NLOOPS,NPROPS,NVERTS,CUTMAX REAL*8 DFRM1,DFRM2,DFRM3,DFRM4 COMMON /DEFORMSCALES/DFRM1,DFRM2,DFRM3,DFRM4 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