gchiurato
/
EVENT2
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
EVENT2/event2_03.f

3627 lines
136 KiB
Fortran
Raw Blame History

PROGRAM MAIN
IMPLICIT NONE
INTEGER NF,NEV
DOUBLE PRECISION EM
EXTERNAL DEMO
NF=5
EM=91.1876
C NEV=1 000 000 000
C NEV=5 000 000
NEV= 500 000
CALL DEMOIN
CALL EVENT2(NEV,EM,NF,DEMO)
CALL DEMOUT(NEV)
END PROGRAM
C-----------------------------------------------------------------------
SUBROUTINE EVENT2(NEV,EM,NFL,USER)
IMPLICIT NONE
C---CALCULATE E+E- EVENT FEATURES TO NEXT-TO-LEADING ORDER
C ACCORDING TO THE METHOD OF CATANI AND SEYMOUR PLB378 (1996) 287
C
C VERSION 0.0 - 18th January 1996
C VERSION 0.1 - 26th February 1996 - bug fix in event orientation
C VERSION 0.2 - 15th March 1996 - removed need for Lorentz boosts
C VERSION 0.3 - 10th November 1997 - improved numerical convergence
C and safety against arithmetic errors
C - added CUTUP variable
C
C - NEV IS THE NUMBER OF EVENTS TO GENERATE
C - EM IS THE TOTAL CENTRE-OF-MASS ENERGY
C - NFL IS THE NUMBER OF FLAVOURS
C - USER IS THE NAME OF ROUTINE TO ANALYSE THE EVENTS
C IT TAKES FIVE ARGUMENTS:
C SUBROUTINE USER(N,NA,ITYPE,P,WEIGHT)
C N = THE NUMBER OF PARTONS
C NA = THE ORDER IN ALPHA
C ITYPE = THE TYPE: 0=TREE LEVEL, 1=SUBTRACTION, 2=FINITE VIRTUAL
C P(4,7) = THEIR MOMENTA (P(1-4,I) IS 4-MOMENTUM OF PARTICLE I,
C IN ORDER Q,QBAR, Q,QBAR,G, Q,QBAR,G,G OR Q',Q,Q'BAR,QBAR
C P(1-4,5) IS TOTAL 4-MOM OF THE N PARTICLES,
C P(1-4,6/7) IS THE 4-MOM OF THE ELECTRON/POSITRON)
C WEIGHT = THE WEIGHT IN THE INTEGRAL
C IT IS ALSO CALLED AT THE END OF EACH FULL EVENT WITH N=0
C A SIMPLE EXAMPLE (CALLED DEMO) IS GIVEN BELOW
C
C IT IS IMPORTANT TO THINK ABOUT THE EFFICIENCY OF THE USER ROUTINE,
C AS IT IS CALLED 17 TIMES PER EVENT
C (THIS IS 1 4-PARTON CALL, 12 3-PARTON CALLS AND 4 2-PARTON CALLS)
C
INTEGER NEV,NFL,NPERM3,NPERM4,NU,ND,I,J,K
LOGICAL KEEP
PARAMETER (NPERM3=2,NPERM4=10)
DOUBLE PRECISION EM,WTWO,WTHR,WFOR,JTWO,JTHR,JFOR,JTMP,
$ MTWO,MTHR,MFOR,VTWO,VTHR,STHR(NPERM3),SFOR(NPERM4),WEIGHT,
$ ZERO,ONE,P(4,7),Q(4,7,NPERM4),NRM,QD,QU,AEM,SUBS(6,NPERM4)
PARAMETER (ZERO=0,ONE=1)
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
INTEGER NPOW1,NPOW2
DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ
COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
EXTERNAL USER
C---PRINT OPENING MESSAGE
WRITE (*,'(/2A)') ' This is EVENT2, a program for calculating',
$ ' jet quantities'
WRITE (*,'(2A)') ' in e+e- annihilation to next-to-leading',
$ ' order in alpha_s'
WRITE (*,'(A)') ' If you use this program, please reference:'
WRITE (*,'(2A)') ' S.Catani & M.H.Seymour,',
$ ' Phys. Lett. B378 (1996) 287.'
WRITE (*,'(/A)') ' Written by Mike Seymour, January 1996'
WRITE (*,'(A/)') ' Version 0.3, November 1997'
WRITE (*,'(A,I10)')' NEV=',NEV
C---INITIALIZE CONSTANTS
PI=ATAN(ONE)*4
PISQ=PI**2
HF=0.5
CA=3
CF=(CA**2-1)/(2*CA)
TR=HF
NF=NFL
ONF=1
IF (NF.NE.0) ONF=ONE/NF
ND=INT((NF+1)/2)
NU=INT(NF/2)
QD=-ONE/3
QU=QD+1
CQ=ND*QD**2+NU*QU**2
AEM=ONE/128
FPAL2=(4*PI*AEM)**2
CUTOFF=1D-8
CUTUP=1
C---METYPE IS 0 FOR ERT, 1 FOR LEIDEN AND 2 FOR GIELE-GLOVER
C (CROSS-SECTION NORMALIZATION IS DIFFERENT IN THE DIFFERENT CASES.
C FOR ERT, 2-PARTON CROSS-SECTION IS 1, FOR THE OTHERS IT IS THE
C BORN CROSS-SECTION FOR E+E- --> QQBAR IN NANOBARNS.)
METYPE=1
IF (METYPE.EQ.0) THEN
NRM=16*PI*EM**2
ELSE
NRM=389385
ENDIF
C---PARAMETERS RELATED TO IMPORTANCE SAMPLING
NPOW1=2
NPOW2=2
XPOW1=1-ONE/NPOW1
XPOW2=1-ONE/NPOW2
COSP=1
COSM=1
COSZ=0
C---START MAIN LOOP
DO I=1,NEV
IF (MOD(I,100 000).EQ.1) CALL RANGEN(-I,SFOR)
C---GENERATE A TWO-PARTON STATE
CALL GENTWO(EM,P,WTWO,*1000)
C---CALCULATE THE JACOBIAN FACTOR
CALL JACTWO(P,JTWO)
C---EVALUATE THE TWO-PARTON TREE-LEVEL MATRIX ELEMENT
CALL MATTWO(P,MTWO)
C---GIVE IT TO THE USER
WEIGHT=NRM*MTWO*WTWO/JTWO/NEV
CALL USER(2,0,0,P,WEIGHT)
IF (CUTUP.GE.1) THEN
C---EVALUATE THE TWO-PARTON ONE-LOOP MATRIX ELEMENT
CALL VIRTWO(P,VTWO)
C---GIVE IT TO THE USER
WEIGHT=NRM*VTWO*WTWO/JTWO/NEV
CALL USER(2,1,2,P,WEIGHT)
ENDIF
C---GENERATE A THREE-PARTON STATE
KEEP=.FALSE.
CALL GENTHR(P,WTHR,KEEP,*1000)
C---CALCULATE THE JACOBIAN FACTOR AND SUBTRACTION CONFIGURATIONS
JTHR=0
DO J=1,NPERM3
CALL SUBTHR(J,P,Q(1,1,J),STHR(J),JTMP,*1000)
JTHR=JTHR+JTMP
ENDDO
IF (KEEP) THEN
C---EVALUATE THE THREE-PARTON TREE-LEVEL MATRIX ELEMENT
CALL MATTHR(P,MTHR)
C---GIVE IT TO THE USER
WEIGHT=NRM*MTHR*WTWO*WTHR/JTHR/NEV
CALL USER(3,1,0,P,WEIGHT)
C---GIVE THE SUBTRACTION CONFIGURATIONS TO THE USER
DO J=1,NPERM3
WEIGHT=-NRM*STHR(J)*WTWO*WTHR/JTHR/NEV
CALL USER(2,1,1,Q(1,1,J),WEIGHT)
ENDDO
C---EVALUATE THE THREE-PARTON ONE-LOOP MATRIX ELEMENT
CALL VIRTHR(P,VTHR)
C---GIVE IT TO THE USER
WEIGHT=NRM*VTHR*WTWO*WTHR/JTHR/NEV
CALL USER(3,2,2,P,WEIGHT)
ENDIF
C---GENERATE A FOUR-PARTON STATE
CALL GENFOR(P,WFOR,KEEP,*1000)
C---CALCULATE THE JACOBIAN FACTOR AND SUBTRACTION CONFIGURATIONS
JFOR=0
DO J=1,NPERM4
CALL SUBFOR(J,P,Q(1,1,J),SFOR(J),JTMP,KEEP,*1000)
JFOR=JFOR+JTMP
SUBS(1,J)=CFSUB
SUBS(2,J)=CASUB
SUBS(3,J)=TFSUB
SUBS(4,J)=GGSUB
SUBS(5,J)=QQSUB
SUBS(6,J)=QPSUB
ENDDO
IF (KEEP) THEN
C---EVALUATE THE FOUR-PARTON TREE-LEVEL MATRIX ELEMENT
CALL MATFOR(P,MFOR)
C---GIVE IT TO THE USER
WEIGHT=NRM*MFOR*WTWO*WTHR*WFOR/JFOR/NEV
CALL USER(4,2,0,P,WEIGHT)
C---GIVE THE SUBTRACTION CONFIGURATIONS TO THE USER
DO J=1,NPERM4
CFSUB=SUBS(1,J)
CASUB=SUBS(2,J)
TFSUB=SUBS(3,J)
GGSUB=SUBS(4,J)
QQSUB=SUBS(5,J)
QPSUB=SUBS(6,J)
WEIGHT=-NRM*SFOR(J)*WTWO*WTHR*WFOR/JFOR/NEV
CALL USER(3,2,1,Q(1,1,J),WEIGHT)
ENDDO
ENDIF
C---TELL THE USER THAT THE EVENT IS COMPLETE
1000 CALL USER(0,0,0,P,ZERO)
ENDDO
END
C-----------------------------------------------------------------------
SUBROUTINE READCARD
END
C-----------------------------------------------------------------------
SUBROUTINE DEMOIN
IMPLICIT NONE
INTEGER I
DOUBLE PRECISION CSUM,CSQR,BSUM,BSQR
COMMON /DEMCOM/CSUM(8),CSQR(8),BSUM(2,6),BSQR(2,6)
DO I=1,8
CSUM(I)=0
CSQR(I)=0
IF (I.LE.6) THEN
BSUM(1,I)=0
BSUM(2,I)=0
BSQR(1,I)=0
BSQR(2,I)=0
ENDIF
ENDDO
CALL GBOOK1( 1,'C-PARAMETER',50,0D0,1D0)
CALL GBOOK1(101,'C-PARAMETER',50,0D0,1D0)
CALL GBOOK1(201,'C-PARAMETER',50,0D0,1D0)
CALL GBOOK1( 11,'C-PARAMETER',37,0D0,0.74D0)
CALL GBOOK1(111,'C-PARAMETER',37,0D0,0.74D0)
CALL GBOOK1(211,'C-PARAMETER',37,0D0,0.74D0)
CALL GBOOK1( 2,'D-PARAMETER',50,0D0,1D0)
CALL GBOOK1(102,'D-PARAMETER',50,0D0,1D0)
CALL GBOOK1(202,'D-PARAMETER',50,0D0,1D0)
CALL GBOOK1( 12,'D-PARAMETER',50,0D0,1D0)
CALL GBOOK1(112,'D-PARAMETER',50,0D0,1D0)
CALL GBOOK1(212,'D-PARAMETER',50,0D0,1D0)
CALL GBOOK1( 3,'THRUST',50,0D0,0.5D0)
CALL GBOOK1(103,'THRUST',50,0D0,0.5D0)
CALL GBOOK1(203,'THRUST',50,0D0,0.5D0)
CALL GBOOK1( 13,'THRUST',50,0D0,0.5D0)
CALL GBOOK1(113,'THRUST',50,0D0,0.5D0)
CALL GBOOK1(213,'THRUST',50,0D0,0.5D0)
CALL GBOOK1( 4,'y3(JADE,P)',40,0D0,0.4D0)
CALL GBOOK1(104,'y3(JADE,P)',40,0D0,0.4D0)
CALL GBOOK1(204,'y3(JADE,P)',40,0D0,0.4D0)
CALL GBOOK1( 14,'y3(JADE,P)',40,0D0,0.4D0)
CALL GBOOK1(114,'y3(JADE,P)',40,0D0,0.4D0)
CALL GBOOK1(214,'y3(JADE,P)',40,0D0,0.4D0)
CALL GBOOK1( 5,'EEC',50,-1D0,1D0)
CALL GBOOK1(105,'EEC',50,-1D0,1D0)
CALL GBOOK1(205,'EEC',50,-1D0,1D0)
CALL GBOOK1( 15,'EEC',50,-1D0,1D0)
CALL GBOOK1(115,'EEC',50,-1D0,1D0)
CALL GBOOK1(215,'EEC',50,-1D0,1D0)
END
C-----------------------------------------------------------------------
SUBROUTINE DEMOUT(NEV)
IMPLICIT NONE
INTEGER I,J,NEV
DOUBLE PRECISION CSUM,CSQR,BSUM,BSQR,SF,SQR
COMMON /DEMCOM/CSUM(8),CSQR(8),BSUM(2,6),BSQR(2,6)
CHARACTER*10 NAME(8)
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
SQR(SF)=SIGN(SQRT(ABS(SF)),SF)
DATA NAME
$/' sigma_0',' sigma_1',' C-PARAM',' D-PARAM',
$ ' 1-THRUST','y3(JADE,p)','EEC(90deg)','EEC(total)'/
SF=CSUM(1)
IF (SF.NE.0) THEN
DO I=1,5
CALL GSCALE( I,1/SF)
CALL GSCALE(100+I,1/SF**2)
CALL GSCALE( 10+I,1/SF)
CALL GSCALE(110+I,1/SF**2)
ENDDO
DO I=1,7
CSUM(I+1)=CSUM(I+1)/SF
CSQR(I+1)=CSQR(I+1)/SF**2
ENDDO
DO I=1,6
BSUM(1,I)=BSUM(1,I)/SF
BSQR(1,I)=BSQR(1,I)/SF**2
BSUM(2,I)=BSUM(2,I)/SF
BSQR(2,I)=BSQR(2,I)/SF**2
ENDDO
ENDIF
DO I=1,5
CALL GOPERA( I,'*', I,200+I,1D0,1D0)
CALL GOPERA(100+I,'-',200+I,100+I,1D0,1D0/NEV)
CALL GOPERA(100+I,'S', 0,100+I,1D0,0D0)
CALL GOPERA( 10+I,'*', 10+I,210+I,1D0,1D0)
CALL GOPERA(110+I,'-',210+I,110+I,1D0,1D0/NEV)
CALL GOPERA(110+I,'S', 0,110+I,1D0,0D0)
CALL GTOPER( I,1,1,0,100+I)
CALL GTOPER(10+I,0,0,0,110+I)
ENDDO
WRITE (*,'(8(1X,2(A,F15.5)/))')
$ (NAME(I),CSUM(I),' +-',SQR(CSQR(I)-CSUM(I)**2/NEV),I=1,8)
WRITE (*,'(2(/A))') ' Moments of EEC function:',' m n'
WRITE (*,'(2I4,F15.5,A,F15.5)')
$ ((I-1,J-1,BSUM(I,J),' +-',SQR(BSQR(I,J)-BSUM(I,J)**2/NEV)
$ ,J=1,6),I=1,2)
END
C-----------------------------------------------------------------------
SUBROUTINE DEMO(N,NA,ITYPE,P,WEIGHT)
IMPLICIT NONE
INTEGER N,NA,ITYPE,I,J,K,L,ORD
DOUBLE PRECISION P(4,7),WEIGHT,DOT,C,D,OS,ORS,E(7),OE(7),COSANG,
$ COSLIM,BEEC,T,TAU,TL,TA(3),PT(3),PC(3,4),Y3,Y,Q(4,7),
$ CSUM(8),CSQR(8),CS(8),BSUM(2,6),BSQR(2,6),BS(2,6)
DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /DEMCOM/CSUM,CSQR,BSUM,BSQR
DATA CS,BS/20*0/
C---FILL THE HISTOGRAMS AT THE END OF EACH EVENT
C (THIS IS SO THAT PAIRS OF LARGE POSITIVE AND NEGATIVE WEIGHTS
C GET CANCELLED BEFORE CALCULATING THE VARIANCE)
IF (N.EQ.0) THEN
DO I=0,10,10
DO J=1,5
CALL GOPERA(200+I+J,'+', I+J, I+J,1D0,1D0)
CALL GOPERA(200+I+J,'*',200+I+J,200+I+J,1D0,1D0)
CALL GOPERA(200+I+J,'+',100+I+J,100+I+J,1D0,1D0)
CALL GRESET(200+I+J)
ENDDO
ENDDO
DO I=1,8
CSUM(I)=CSUM(I)+CS(I)
CSQR(I)=CSQR(I)+CS(I)**2
CS(I)=0
IF (I.LE.6) THEN
BSUM(1,I)=BSUM(1,I)+BS(1,I)
BSUM(2,I)=BSUM(2,I)+BS(2,I)
BSQR(1,I)=BSQR(1,I)+BS(1,I)**2
BSQR(2,I)=BSQR(2,I)+BS(2,I)**2
BS(1,I)=0
BS(2,I)=0
ENDIF
ENDDO
RETURN
ENDIF
C---CALCULATE THE ENERGY OF EACH PARTICLE
OS=1/DOT(P,5,5)
ORS=SQRT(OS)
DO I=1,7
IF (I.LE.N.OR.I.GE.5) THEN
E(I)=DOT(P,I,5)*ORS
OE(I)=1/E(I)
ENDIF
ENDDO
C---CALCULATE THE TOTAL CROSS-SECTION
IF (NA.EQ.0) CS(1)=CS(1)+WEIGHT
IF (NA.EQ.1) CS(2)=CS(2)+WEIGHT
IF (N.LT.3) RETURN
ORD=20-10*NA
c$$$C---ONLY INCLUDE ONE OF THE COLOUR PIECES IN THE ALPHA_S**2 TERMS
c$$$ IF (NA.EQ.2) WEIGHT=WEIGHT*CASUB
C---CALCULATE THE C-PARAMETER
C=3
DO I=2,N
DO J=1,I-1
C=C-3*DOT(P,I,J)**2*OS/(E(I)*E(J))
ENDDO
ENDDO
CALL GFILL1(201+ORD,C,C*WEIGHT*50)
IF (ORD.EQ.0) CS(3)=CS(3)+WEIGHT*C
C---CALCULATE THE D-PARAMETER
IF (N.EQ.4) THEN
D=27*(2*(DOT(P,1,2)*DOT(P,1,3)*DOT(P,2,4)*DOT(P,3,4)+
$ DOT(P,1,2)*DOT(P,1,4)*DOT(P,2,3)*DOT(P,3,4)+
$ DOT(P,1,3)*DOT(P,1,4)*DOT(P,2,3)*DOT(P,2,4))-(
$ DOT(P,1,2)**2*DOT(P,3,4)**2+DOT(P,1,3)**2*DOT(P,2,4)**2+
$ DOT(P,1,4)**2*DOT(P,2,3)**2))/(
$ (DOT(P,1,2)+DOT(P,1,3)+DOT(P,1,4))*
$ (DOT(P,1,2)+DOT(P,2,3)+DOT(P,2,4))*
$ (DOT(P,1,3)+DOT(P,2,3)+DOT(P,3,4))*
$ (DOT(P,1,4)+DOT(P,2,4)+DOT(P,3,4)))
ELSE
D=0
ENDIF
CALL GFILL1(202+ORD,D,D*WEIGHT*50)
IF (ORD.EQ.0) CS(4)=CS(4)+WEIGHT*D
C---CALCULATE THE THRUST
IF (N.EQ.4) THEN
T=0
DO K=2,N
DO J=1,K-1
TA(1)=P(2,J)*P(3,K)-P(3,J)*P(2,K)
TA(2)=P(3,J)*P(1,K)-P(1,J)*P(3,K)
TA(3)=P(1,J)*P(2,K)-P(2,J)*P(1,K)
PT(1)=0
PT(2)=0
PT(3)=0
DO I=1,N
IF (I.NE.J.AND.I.NE.K) THEN
IF (P(1,I)*TA(1)+P(2,I)*TA(2)+P(3,I)*TA(3).GE.0.) THEN
PT(1)=PT(1)+P(1,I)
PT(2)=PT(2)+P(2,I)
PT(3)=PT(3)+P(3,I)
ELSE
PT(1)=PT(1)-P(1,I)
PT(2)=PT(2)-P(2,I)
PT(3)=PT(3)-P(3,I)
ENDIF
ENDIF
ENDDO
DO I=1,3
PC(I,1)=PT(I)+P(I,J)+P(I,K)
PC(I,2)=PT(I)+P(I,J)-P(I,K)
PC(I,3)=PT(I)-P(I,J)+P(I,K)
PC(I,4)=PT(I)-P(I,J)-P(I,K)
ENDDO
DO I=1,4
TL=PC(1,I)**2+PC(2,I)**2+PC(3,I)**2
IF (TL.GT.T) T=TL
ENDDO
ENDDO
ENDDO
T=SQRT(T)*ORS
ELSE
T=MAX(E(1),E(2),E(3))*ORS*2
ENDIF
TAU = 1 - T
CALL GFILL1(203+ORD,TAU,TAU*WEIGHT*100)
IF (ORD.EQ.0) CS(5)=CS(5)+TAU*WEIGHT
C---CALCULATE THE Y3 VALUE (USES P SCHEME FOR NO PARTICULAR REASON)
IF (N.EQ.4) THEN
Y3=1D6
DO I=2,4
DO J=1,I-1
Y=2*DOT(P,I,J)*OS
IF (Y.LT.Y3) THEN
Y3=Y
K=I
L=J
ENDIF
ENDDO
ENDDO
DO I=1,3
Q(I,L)=P(I,L)+P(I,K)
ENDDO
Q(4,L)=SQRT(Q(1,L)**2+Q(2,L)**2+Q(3,L)**2)
IF (K.LT.4) THEN
DO I=1,4
Q(I,K)=P(I,4)
ENDDO
ENDIF
DO I=1,3
IF (I.NE.K.AND.I.NE.L) THEN
DO J=1,4
Q(J,I)=P(J,I)
ENDDO
ENDIF
ENDDO
Y3=1D6
DO I=2,3
DO J=1,I-1
Y=2*DOT(Q,I,J)*OS
IF (Y.LT.Y3) Y3=Y
ENDDO
ENDDO
ELSE
Y3=1-T
ENDIF
CALL GFILL1(204+ORD,Y3,Y3*WEIGHT*100)
IF (ORD.EQ.0) CS(6)=CS(6)+Y3*WEIGHT
C---CALCULATE THE ENERGY-ENERGY CORRELATION
DO I=2,N
DO J=1,I-1
COSANG=1-DOT(P,I,J)*OE(I)*OE(J)
IF (ABS(COSANG).GE.1) COSANG=SIGN(1-1D-12,COSANG)
BEEC=(1-COSANG**2)*2*WEIGHT*E(I)*E(J)*OS
CALL GFILL1(205+ORD,COSANG,BEEC*50/2)
COSLIM=0.1D0
IF (ORD.EQ.0) THEN
IF (ABS(COSANG).LT.COSLIM)
& CS(7)=CS(7)+BEEC/(2*COSLIM)
CS(8)=CS(8)+BEEC
DO K=1,2
DO L=1,6
BS(K,L)=BS(K,L)+
$ BEEC*SQRT(1-COSANG**2)**(L-1)*COSANG**(K-1)
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
END
C-----------------------------------------------------------------------
FUNCTION DOT(P,I,J)
IMPLICIT NONE
C---RETURN THE DOT PRODUCT OF P(*,I) AND P(*,J)
INTEGER I,J
DOUBLE PRECISION DOT,P(4,7)
DOT=P(4,I)*P(4,J)-P(3,I)*P(3,J)-P(2,I)*P(2,J)-P(1,I)*P(1,J)
END
C-----------------------------------------------------------------------
SUBROUTINE GENTWO(EM,P,W,*)
IMPLICIT NONE
C---GENERATE A TWO-PARTON CONFIGURATION
C THE WEIGHT GIVES THE TOTAL VOLUME OF PHASE-SPACE
DOUBLE PRECISION EM,P(4,7),W,R(5),E,C,S,A
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
INTEGER NPOW1,NPOW2
DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ
COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
C---GENERATE C=COS(THETA) ACCORDING TO COSZ+COSP*(1+C)**2+COSM*(1-C)**2
CALL RANGEN(5,R)
E=EM/2
IF (R(1)*(4*COSP+4*COSM+3*COSZ).LT.4*COSP) THEN
C=2*MAX(R(2),R(3),R(4))-1
ELSEIF (R(1)*(4*COSP+4*COSM+3*COSZ).LT.4*COSP+4*COSM) THEN
C=2*MIN(R(2),R(3),R(4))-1
ELSE
C=2*R(2)-1
ENDIF
S=1-C**2
IF (S.LT.0) CALL WARNER('GENTWO',1,*999)
S=SQRT(S)
A=2*PI*R(5)
P(1,1)= E*COS(A)*S
P(2,1)= E*SIN(A)*S
P(3,1)= E*C
P(4,1)= E
P(1,2)= -E*COS(A)*S
P(2,2)= -E*SIN(A)*S
P(3,2)= -E*C
P(4,2)= E
P(1,5)= 0
P(2,5)= 0
P(3,5)= 0
P(4,5)=2*E
P(1,6)= 0
P(2,6)= 0
P(3,6)= E
P(4,6)= E
P(1,7)= 0
P(2,7)= 0
P(3,7)= -E
P(4,7)= E
W=1/(16*PI*EM**2)
RETURN
999 RETURN 1
END
C-----------------------------------------------------------------------
SUBROUTINE GENTHR(P,W,KEEP,*)
IMPLICIT NONE
C---GENERATE A THREE-PARTON CONFIGURATION FROM A GIVEN TWO-PARTON ONE
C IF ANY DOT PRODUCTS ARE BELOW CUTOFF, EVENT IS REJECTED,
C BY JUMPING TO THE SPECIFIED LABEL
C IF ANY DOT PRODUCTS ARE BELOW CUTUP, EVENT IS MARKED BY KEEP=.TRUE.
INTEGER EMIT,OMIT
LOGICAL KEEP
DOUBLE PRECISION P(4,7),W,X(2),R(4),EMSQ,DOT,T
INTEGER I
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
INTEGER NPOW1,NPOW2
DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ
COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
C---GENERATE X1,X2 VALUES
CALL RANGEN(4,R)
X(1)=1-(R(1)*R(2))**NPOW1
IF (X(1).GT.1) CALL WARNER('GENTHR',1,*999)
X(2)=2-X(1)-(1-X(1))**R(3)
C---ENFORCE INVARIANT MASS CUTOFFS
IF (1-X(1).LT.CUTOFF.OR.1-X(2).LT.CUTOFF.OR.X(1)+X(2)-1.LT.CUTOFF)
$ RETURN 1
IF (1-X(1).LT.CUTUP.OR.1-X(2).LT.CUTUP.OR.X(1)+X(2)-1.LT.CUTUP)
$ KEEP=.TRUE.
C---CHOOSE WHICH ONE EMITTED
EMIT=INT(R(4)*2)+1
OMIT=3-EMIT
C---GENERATE THEIR MOMENTA
CALL GENDIP(P,EMIT,OMIT,3,X,*999)
C---NOW SWAP EMITTER AND 3 HALF THE TIME
IF (MOD(INT(4*R(4)),2).EQ.1) THEN
DO I=1,4
T=P(I,EMIT)
P(I,EMIT)=P(I,3)
P(I,3)=T
ENDDO
ENDIF
C---CALCULATE WEIGHT
EMSQ=DOT(P,5,5)
W=HF*EMSQ/(16*PISQ)
RETURN
999 RETURN 1
END
C-----------------------------------------------------------------------
SUBROUTINE GENFOR(P,W,KEEP,*)
IMPLICIT NONE
C---GENERATE A FOUR-BODY CONFIGURATION FROM A GIVEN THREE-BODY ONE
C IF ANY DOT PRODUCTS ARE BELOW CUTOFF, EVENT IS REJECTED,
C BY JUMPING TO THE SPECIFIED LABEL
C IF ANY DOT PRODUCTS ARE BELOW CUTUP, EVENT IS MARKED BY KEEP=.TRUE.
INTEGER I,J,EMIT,OMIT
LOGICAL KEEP
DOUBLE PRECISION P(4,7),W,X(2),R(6),T,EMSQ,DOT
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
INTEGER NPOW1,NPOW2
DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ
COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
C---GENERATE X1,X2 VALUES
CALL RANGEN(6,R)
X(1)=1-(R(1)*R(2))**NPOW2
IF (X(1).GT.1) CALL WARNER('GENFOR',1,*999)
X(2)=2-X(1)-(1-X(1))**R(3)
C---ENFORCE INVARIANT MASS CUTOFFS
W=0
IF (1-X(1).LT.CUTOFF.OR.1-X(2).LT.CUTOFF.OR.X(1)+X(2)-1.LT.CUTOFF)
$ RETURN 1
C---CHOOSE WHICH ONE TO SPLIT UNIFORMLY
EMIT=INT(R(4)*3)+1
C---AND WHICH OF THE OTHERS TO BALANCE ITS MOMENTUM WITH
OMIT=INT(R(5)*2)+1
IF (OMIT.EQ.EMIT) OMIT=3
C---GENERATE THEIR MOMENTA
CALL GENDIP(P,EMIT,OMIT,4,X,*999)
C---HALF THE TIME, SWAP PARTICLES 3 AND 4 SO THEY ARE SYMMETRIC
IF (R(6).GT.0.5) THEN
DO I=1,4
T=P(I,3)
P(I,3)=P(I,4)
P(I,4)=T
ENDDO
ENDIF
C---REJECT THE EVENT IF ANY PAIR MASSES ARE BELOW CUTOFF*Q**2
EMSQ=DOT(P,5,5)
DO I=1,3
DO J=I+1,4
IF (2*DOT(P,I,J).LT.CUTOFF*EMSQ) RETURN 1
IF (2*DOT(P,I,J).LT.CUTUP*EMSQ) KEEP=.TRUE.
ENDDO
ENDDO
C---CALCULATE WEIGHT (NOTE: THIS IS NORMALIZED TO THE TOTAL CMF MASS,
C RATHER THAN THE SUBSYSTEM MASS. THE EXTRA FACTOR GETS INSERTED IN
C THE JACOBIAN FACTOR CALCULATED IN SUBFOR)
W=HF*EMSQ/(16*PISQ)
RETURN
999 RETURN 1
END
C-----------------------------------------------------------------------
SUBROUTINE GENDIP(P,EMIT,OMIT,THRD,X,*)
IMPLICIT NONE
C---GENERATE A 2->3 DIPOLE EMISSION
DOUBLE PRECISION P(4,7),X(2),R(1),Q(4,2),PTDQ,C(4),D(4)
INTEGER EMIT,OMIT,THRD,I
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
CALL RANGEN(1,R)
PTDQ=(X(1)+X(2)-1)*(1-X(1))*(1-X(2))/X(1)**2
IF (PTDQ.LT.0) CALL WARNER('GENDIP',1,*999)
PTDQ=SQRT(PTDQ)
CALL GTPERP(PTDQ,P,EMIT,OMIT,6,C,D,*999)
DO I=1,4
Q(I,1)=X(1)*P(I,OMIT)
Q(I,2)=(1-(1-X(2))/X(1))*P(I,EMIT)
$ +(1-X(1))*(1-X(2))/X(1)*P(I,OMIT)
$ +COS(R(1)*2*PI)*C(I)+SIN(R(1)*2*PI)*D(I)
ENDDO
DO I=1,4
P(I,THRD)=P(I,EMIT)+P(I,OMIT)-Q(I,1)-Q(I,2)
P(I,OMIT)=Q(I,1)
P(I,EMIT)=Q(I,2)
ENDDO
RETURN
999 RETURN 1
END
C-----------------------------------------------------------------------
SUBROUTINE GTPERP(PTDQ,P,I,J,K,C,D,*)
IMPLICIT NONE
C---FIND THE VECTORS PERPENDICULAR TO P(I) AND P(J)
C C AND D ARE PURELY SPACE-LIKE VECTORS IN THE P(I)+P(J) CMF,
C WITH C IN THE SAME PLANE AS P(K) AND D PERPENDICULAR TO IT,
C BOTH HAVING LENGTH PTDQ*SQRT(2*DOT(P,I,J))
DOUBLE PRECISION PTDQ,P(4,7),C(4),D(4),PTF,DIJ,DIK,DJK,DOT,EPS4
INTEGER I,J,K,L
DIJ=DOT(P,I,J)
DIK=DOT(P,I,K)
DJK=DOT(P,J,K)
IF (DIK.LE.0.OR.DJK.LE.0) CALL WARNER('GTPERP',1,*999)
PTF=PTDQ/SQRT(DIK*DJK)
DO L=1,4
C(L)=PTF*(DIJ*P(L,K)-DJK*P(L,I)-DIK*P(L,J))
ENDDO
DO L=1,4
D(L)=EPS4(L,P(1,I),P(1,J),C)/DIJ
ENDDO
RETURN
999 RETURN 1
END
C-----------------------------------------------------------------------
FUNCTION EPS4(I,A,B,C)
IMPLICIT NONE
DOUBLE PRECISION EPS4,EPS3,A(4),B(4),C(4),AA(3),BB(3),CC(3)
INTEGER I,J,K,S(4)
DATA S/+1,-1,+1,+1/
J=1
DO K=1,3
IF (I.EQ.J) J=J+1
AA(K)=A(J)
BB(K)=B(J)
CC(K)=C(J)
J=J+1
ENDDO
EPS4=0
DO J=1,3
EPS4=EPS4+CC(J)*EPS3(J,AA,BB)
ENDDO
EPS4=S(I)*EPS4
END
C-----------------------------------------------------------------------
FUNCTION EPS3(I,A,B)
IMPLICIT NONE
DOUBLE PRECISION EPS3,A(3),B(3),AA(2),BB(2)
INTEGER I,J,K,S(3)
DATA S/+1,-1,+1/
J=1
DO K=1,2
IF (I.EQ.J) J=J+1
AA(K)=A(J)
BB(K)=B(J)
J=J+1
ENDDO
EPS3=S(I)*(AA(1)*BB(2)-AA(2)*BB(1))
END
C-----------------------------------------------------------------------
SUBROUTINE JACTWO(P,J)
IMPLICIT NONE
C---CALCULATE THE TWO-PARTON JACOBIAN FACTOR
DOUBLE PRECISION P(4,7),J,CP,CM,DOT
INTEGER NPOW1,NPOW2
DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ
COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
CP=DOT(P,1,7)*DOT(P,5,5)/(DOT(P,1,5)*DOT(P,7,5))
CM=DOT(P,1,6)*DOT(P,5,5)/(DOT(P,1,5)*DOT(P,6,5))
J=(COSZ+COSP*CP**2+COSM*CM**2)*3/(3*COSZ+4*COSP+4*COSM)
END
C-----------------------------------------------------------------------
SUBROUTINE MATTWO(P,M)
IMPLICIT NONE
C---EVALUATE THE TWO-PARTON MATRIX ELEMENT SQUARED FOR THE GIVEN
C CONFIGURATION. USING ERT NORMALIZATION, THIS IS JUST UNITY.
DOUBLE PRECISION P(4,7),M,EMSQ,DOT
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
IF (METYPE.EQ.0) THEN
M=1
ELSEIF (METYPE.EQ.1) THEN
EMSQ=DOT(P,5,5)
M=2*(-8*DOT(P,1,6)*DOT(P,2,6))+EMSQ/2*(4*EMSQ)
C---INCLUDE EXTERNAL FACTORS
M=M/4/EMSQ
M=M*FPAL2*CA*CQ*4/EMSQ
ELSEIF (METYPE.EQ.2) THEN
CALL CONVERT(P)
CALL SETUP
CALL SFMEE(2,0,M)
ELSE
STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN MATTWO!'
ENDIF
END
C-----------------------------------------------------------------------
SUBROUTINE MATTHR(P,M)
IMPLICIT NONE
C---EVALUATE THE THREE-PARTON MATRIX ELEMENT SQUARED FOR THE GIVEN
C CONFIGURATION. INCLUDES COLOUR FACTORS, BUT NOT ALPHA_S/2PI.
DOUBLE PRECISION P(4,7),M,EMSQ,E1,E2,S13,S23,DOT,K1,K2,K3
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
IF (METYPE.EQ.0) THEN
EMSQ=DOT(P,5,5)
S13=2*DOT(P,1,3)
S23=2*DOT(P,2,3)
E1=EMSQ-S23
E2=EMSQ-S13
M=CF*16*PISQ*(E1**2+E2**2)/(EMSQ*S13*S23)
ELSEIF (METYPE.EQ.1) THEN
EMSQ=DOT(P,5,5)
S13=2*DOT(P,1,3)
S23=2*DOT(P,2,3)
E1=EMSQ-S23
E2=EMSQ-S13
K1=2*DOT(P,1,6)
K2=2*DOT(P,2,6)
K3=2*DOT(P,3,6)
M=CF*16*PISQ/(S13*S23)*(
$ 2*(2*(K1**2*S23+K2**2*S13
$ -K1*K2*(E1+E2)-K2*K3*E2-K3*K1*E1))
$ +EMSQ/2*(4*(E1**2+E2**2)))
C---INCLUDE EXTERNAL FACTORS
M=M/4/EMSQ
M=M*FPAL2*CA*CQ*4/EMSQ
ELSEIF (METYPE.EQ.2) THEN
CALL CONVERT(P)
CALL SETUP
CALL SFMEE(3,0,M)
ELSE
STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN MATTHR!'
ENDIF
END
C-----------------------------------------------------------------------
SUBROUTINE CONTHR(P,V,VV,M)
IMPLICIT NONE
C---EVALUATE THE CONTRACTION OF THE THREE-PARTON MATRIX ELEMENT SQUARED
C FOR THE GIVEN CONFIGURATION WITH THE VECTOR V. IT IS ASSUMED THAT
C V IS PERPENDICULAR TO THE GLUON, V.P3.EQ.0, AND THAT V.V.NE.0 .
C THE NORMALIZATION IS SUCH THAT THE AVERAGE OVER V'S AZIMUTH AROUND
C P3 IS EQUAL TO HALF OF MATTHR
INTEGER I,J
DOUBLE PRECISION P(4,7),V(4),M,DOT,VDOT,EMSQ,OMSQ,VV,Y1,Y2,L1,L2,
$ DOTS,T1,T2,T3,T4,T
C---EERAD COMMON BLOCKS
DOUBLE PRECISION CZUR,CZUL,CZDR,CZDL,CZLR,CZLL,CZNR,CZNL,
$ CGUR,CGUL,CGDR,CGDL,CGLR,CGLL,CGNR,CGNL,
$ W,EM,SW2,RMZ,RGZ,
$ CGGVV,CGZVV,CZZVV,CVV,CGGAA,CGZAA,CZZAA,CAA
INTEGER NUTF,NDTF,ICON
COMMON /ZCUP/CZUR,CZUL,CZDR,CZDL,CZLR,CZLL,CZNR,CZNL
COMMON /GCUP/CGUR,CGUL,CGDR,CGDL,CGLR,CGLL,CGNR,CGNL
COMMON /GENERA/ W,EM,SW2,RMZ,RGZ,NUTF,NDTF,ICON
C---END EERAD COMMON BLOCKS
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
VDOT(I)=P(4,I)*V(4)-P(3,I)*V(3)-P(2,I)*V(2)-P(1,I)*V(1)
IF (METYPE.EQ.0) THEN
EMSQ=DOT(P,5,5)
OMSQ=1/EMSQ
Y1=2*DOT(P,2,3)
Y2=2*DOT(P,1,3)
T=CF*16*PISQ/(Y1*Y2)**2
M=-T*(2*(VDOT(1)*Y1-VDOT(2)*Y2)**2/VV
$ -HF*(Y1**2+Y2**2)*Y1*Y2*OMSQ)
ELSEIF (METYPE.EQ.1.OR.METYPE.EQ.2) THEN
EMSQ=DOT(P,5,5)
OMSQ=1/EMSQ
Y1=2*DOT(P,2,3)
Y2=2*DOT(P,1,3)
L1=2*DOT(P,1,7)
L2=2*DOT(P,2,7)
DOTS=1/(VV*(Y1+Y2)**2)
C---INCLUDE COUPLING STRENGTHS
IF (METYPE.EQ.1) THEN
CVV=CQ*FPAL2
CAA=0
ELSE
cggvv=nutf*0.25*(cgur**2+cgul**2)*(cglr**2+cgll**2)
$ +ndtf*0.25*(cgdr**2+cgdl**2)*(cglr**2+cgll**2)
cggaa=nutf*0.25*(cgur**2-cgul**2)*(cglr**2-cgll**2)
$ +ndtf*0.25*(cgdr**2-cgdl**2)*(cglr**2-cgll**2)
cgzvv=nutf*0.5*(cgur*czur+cgul*czul)*(cglr*czlr+cgll*czll)
$ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2)
$ +ndtf*0.5*(cgdr*czdr+cgdl*czdl)*(cglr*czlr+cgll*czll)
$ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2)
cgzaa=nutf*0.5*(cgur*czur-cgul*czul)*(cglr*czlr-cgll*czll)
$ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2)
$ +ndtf*0.5*(cgdr*czdr-cgdl*czdl)*(cglr*czlr-cgll*czll)
$ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2)
czzvv=nutf*0.25*(czur**2+czul**2)*(czlr**2+czll**2)
$ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2)
$ +ndtf*0.25*(czdr**2+czdl**2)*(czlr**2+czll**2)
$ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2)
czzaa=nutf*0.25*(czur**2-czul**2)*(czlr**2-czll**2)
$ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2)
$ +ndtf*0.25*(czdr**2-czdl**2)*(czlr**2-czll**2)
$ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2)
if (icon.eq.1) then
cvv=cggvv
caa=cggaa
elseif (icon.eq.2) then
cvv=cgzvv
caa=cgzaa
elseif (icon.eq.3) then
cvv=czzvv
caa=czzaa
else
cvv=cggvv+cgzvv+czzvv
caa=cggaa+cgzaa+czzaa
endif
cvv=cvv*em**4
caa=caa*em**4
ENDIF
C---THE VECTOR-VECTOR PART
T=CF*16*PISQ*CA*CVV/(Y1*Y2)**2
T1=T*4*(Y1**2*(L2**2+(EMSQ-L2)**2)+Y2**2*(L1**2+(EMSQ-L1)**2)
$ +2*Y1*Y2*(L1*EMSQ+L2*EMSQ-2*L1*L2))*OMSQ**2*DOTS
T2=T*4*(Y2*(EMSQ-2*L1)-Y1*(EMSQ-2*L2))*Y1*Y2*OMSQ**2*DOTS
T3=T*8*(Y1*Y2*OMSQ)**2*DOTS
T4=T*((EMSQ-Y1-L1-L2)**2+(EMSQ-Y2-L1-L2)**2)
$ *Y1*Y2*OMSQ
C---THE AXIAL-AXIAL PART
T=CF*16*PISQ*CA*CAA/(Y1*Y2)**2
T1=T1+T*4*(Y1*(EMSQ-2*L2)-Y2*(EMSQ-2*L1))*(Y1+Y2)*OMSQ*DOTS
T2=T2-T*4*(Y1+Y2)*Y1*Y2*OMSQ*DOTS
T3=T3
T4=T4-T*(Y1+Y2+2*(L1+L2-EMSQ))*(Y1-Y2)*Y1*Y2*OMSQ
M=-T1*(VDOT(1)*Y1-VDOT(2)*Y2)**2
$ -T2*2*(VDOT(1)*Y1-VDOT(2)*Y2)*
$ (VDOT(7)*(Y2+Y1-EMSQ+L1+L2)-(EMSQ-L1-L2)*VDOT(6))
$ -T3*(VDOT(7)*(Y2+Y1-EMSQ+L1+L2)-(EMSQ-L1-L2)*VDOT(6))**2
$ +T4
ELSE
STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN CONTHR!'
ENDIF
END
C-----------------------------------------------------------------------
SUBROUTINE MATFOR(P,M)
IMPLICIT NONE
C---EVALUATE THE FOUR-PARTON MATRIX ELEMENT SQUARED FOR THE GIVEN
C CONFIGURATION. INCLUDES COLOUR FACTORS, BUT NOT (ALPHA_S/2PI)**2.
DOUBLE PRECISION P(4,7),M,ERTA,ERTB,ERTC,ERTD,ERTE,
$ LEIA,LEIB,LEIC,LEID,LEIE,
$ A,B,C,DS,DD,E,EMSQ,DOT,T
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
C---EERAD COMMON BLOCKS
INTEGER IPS
COMMON /PERM/IPS
C---END EERAD COMMON BLOCKS
IF (METYPE.EQ.0) THEN
A=HF*(
$ ERTA(P,1,2,3,4)+ERTA(P,2,1,3,4)+ERTA(P,1,2,4,3)+ERTA(P,2,1,4,3))
B=HF*(
$ ERTB(P,1,2,3,4)+ERTB(P,2,1,3,4)+ERTB(P,1,2,4,3)+ERTB(P,2,1,4,3))
C=HF*(
$ ERTC(P,1,2,3,4)+ERTC(P,2,1,3,4)+ERTC(P,1,2,4,3)+ERTC(P,2,1,4,3))
DD=HF*(
$ ERTD(P,1,2,3,4)+ERTD(P,2,1,4,3)+ERTD(P,3,4,1,2)+ERTD(P,4,3,2,1))
DS=HF*HF*(
$ ERTD(P,2,1,3,4)+ERTD(P,1,2,4,3)+ERTD(P,3,4,2,1)+ERTD(P,4,3,1,2))
E=HF*HF*(
$ ERTE(P,1,2,3,4)+ERTE(P,2,1,4,3)+ERTE(P,3,4,1,2)+ERTE(P,4,3,2,1)+
$ ERTE(P,2,1,3,4)+ERTE(P,1,2,4,3)+ERTE(P,3,4,2,1)+ERTE(P,4,3,1,2))
CFSUB=A+B+E
CASUB=C-(B+E)/2
TFSUB=DD+(DS-HF*DD)*ONF
M=CF*CFSUB+CA*CASUB+NF*TR*TFSUB
CFSUB=CFSUB/(M*CF)
CASUB=CASUB/(M*CF)
TFSUB=TFSUB/(M*CF)
GGSUB=(CF*A+(CF-CA/2)*B+CA*C)/M
QQSUB=(TR*(DD*HF+DS)+(CF-CA/2)*E)/M
QPSUB=TR*(NF-1)*DD/M
C---INCLUDE EXTERNAL FACTORS
EMSQ=DOT(P,5,5)
M=M*256*PI**4*CF/EMSQ
ELSEIF (METYPE.EQ.1) THEN
EMSQ=DOT(P,5,5)
A=2*(LEIA(P,P(1,6),1,2,3,4)+LEIA(P,P(1,6),2,1,3,4)
$ +LEIA(P,P(1,6),1,2,4,3)+LEIA(P,P(1,6),2,1,4,3))+EMSQ/2*(
$ ERTA(P,1,2,3,4)+ERTA(P,2,1,3,4)+ERTA(P,1,2,4,3)+ERTA(P,2,1,4,3))
B=2*(LEIB(P,P(1,6),1,2,3,4)+LEIB(P,P(1,6),2,1,3,4)
$ +LEIB(P,P(1,6),1,2,4,3)+LEIB(P,P(1,6),2,1,4,3))+EMSQ/2*(
$ ERTB(P,1,2,3,4)+ERTB(P,2,1,3,4)+ERTB(P,1,2,4,3)+ERTB(P,2,1,4,3))
C=2*(LEIC(P,P(1,6),1,2,3,4)+LEIC(P,P(1,6),2,1,3,4)
$ +LEIC(P,P(1,6),1,2,4,3)+LEIC(P,P(1,6),2,1,4,3))+EMSQ/2*(
$ ERTC(P,1,2,3,4)+ERTC(P,2,1,3,4)+ERTC(P,1,2,4,3)+ERTC(P,2,1,4,3))
DD=2*(LEID(P,P(1,6),1,2,3,4)+LEID(P,P(1,6),2,1,4,3)
$ +LEID(P,P(1,6),3,4,1,2)+LEID(P,P(1,6),4,3,2,1))+EMSQ/2*(
$ ERTD(P,1,2,3,4)+ERTD(P,2,1,4,3)+ERTD(P,3,4,1,2)+ERTD(P,4,3,2,1))
DS=2*(LEID(P,P(1,6),2,1,3,4)+LEID(P,P(1,6),1,2,4,3)
$ +LEID(P,P(1,6),3,4,2,1)+LEID(P,P(1,6),4,3,1,2))+EMSQ/2*(
$ ERTD(P,2,1,3,4)+ERTD(P,1,2,4,3)+ERTD(P,3,4,2,1)+ERTD(P,4,3,1,2))
E=2*(LEIE(P,P(1,6),1,2,3,4)+LEIE(P,P(1,6),2,1,4,3)
$ +LEIE(P,P(1,6),3,4,1,2)+LEIE(P,P(1,6),4,3,2,1)
$ +LEIE(P,P(1,6),2,1,3,4)+LEIE(P,P(1,6),1,2,4,3)
$ +LEIE(P,P(1,6),3,4,2,1)+LEIE(P,P(1,6),4,3,1,2))+EMSQ/2*(
$ ERTE(P,1,2,3,4)+ERTE(P,2,1,4,3)+ERTE(P,3,4,1,2)+ERTE(P,4,3,2,1)+
$ ERTE(P,2,1,3,4)+ERTE(P,1,2,4,3)+ERTE(P,3,4,2,1)+ERTE(P,4,3,1,2))
A=HF*A
B=HF*B
C=HF*C
DD=HF*DD
DS=HF*HF*DS
E=HF*HF*E
CFSUB=A+B+E
CASUB=C-(B+E)/2
TFSUB=DD+(DS-HF*DD)*ONF
M=CF*CFSUB+CA*CASUB+NF*TR*TFSUB
CFSUB=CFSUB/(M*CF)
CASUB=CASUB/(M*CF)
TFSUB=TFSUB/(M*CF)
GGSUB=(CF*A+(CF-CA/2)*B+CA*C)/M
QQSUB=(TR*(DD*HF+DS)+(CF-CA/2)*E)/M
QPSUB=TR*(NF-1)*DD/M
C---INCLUDE EXTERNAL FACTORS
M=M*256*PI**4*CF/EMSQ
M=M*FPAL2*CA*CQ*4/EMSQ
ELSEIF (METYPE.EQ.2) THEN
CALL CONVERT(P)
CALL SETUP
IPS=1
CALL SFMEE(4,0,M)
IPS=2
CALL SFMEE(4,0,T)
M=M+T
ELSE
STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN MATFOR!'
ENDIF
END
C-----------------------------------------------------------------------
SUBROUTINE VIRTWO(P,V)
IMPLICIT NONE
C---CALCULATE THE TWO-PARTON MATRIX-ELEMENT AT NEXT-TO-LEADING ORDER
DOUBLE PRECISION P(4,7),V,M,SUB
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
C---CALCULATE THE LOWEST-ORDER MATRIX-ELEMENT
CALL MATTWO(P,M)
IF (METYPE.EQ.0.OR.METYPE.EQ.1) THEN
C---VIRTUAL CROSS-SECTION, WITH THE POLES REMOVED AND THE LIMIT TAKEN
V=M*CF*(PISQ-8)
ELSEIF (METYPE.EQ.2) THEN
CALL CONVERT(P)
CALL SETUP
CALL SFMEE(2,1,V)
ELSE
STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN VIRTWO!'
ENDIF
C---THE SUBTRACTION
SUB=CF*(10-PISQ)
C---THE TOTAL
V=V + M*SUB
END
C-----------------------------------------------------------------------
SUBROUTINE VIRTHR(P,V)
IMPLICIT NONE
C---CALCULATE THE THREE-PARTON MATRIX-ELEMENT AT NEXT-TO-LEADING ORDER
DOUBLE PRECISION P(4,7),V,DOT,EMSQ,Y12,Y13,Y23,L12,L13,L23,
$ M,ERTV,LEIV,CFTOT,CATOT,TFTOT,TOT
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
EMSQ=DOT(P,5,5)
Y12=2*DOT(P,1,2)/EMSQ
Y13=2*DOT(P,1,3)/EMSQ
Y23=2*DOT(P,2,3)/EMSQ
L12=LOG(Y12)
L13=LOG(Y13)
L23=LOG(Y23)
C---CALCULATE THE LOWEST-ORDER MATRIX-ELEMENT
CALL MATTHR(P,M)
IF (METYPE.EQ.0) THEN
C---THE FIRST ERT EXPRESSION, WITH THE POLES REMOVED AND THE LIMIT TAKEN
CFTOT=M*(HF*PISQ*2 - 8 - L12**2)
CATOT=M*(HF*PISQ + HF*(L12**2-L13**2-L23**2))
TFTOT=0
C---THE FINITE TERMS NOT PROPORTIONAL TO TREE-LEVEL RESULT
TOT=ERTV(P)
CFTOT=CFTOT+CF*CFSUB*TOT
CATOT=CATOT+CF*CASUB*TOT
TFTOT=TFTOT+CF*TFSUB*TOT
ELSEIF (METYPE.EQ.1) THEN
C---THE FIRST ERT EXPRESSION, WITH THE POLES REMOVED AND THE LIMIT TAKEN
CFTOT=M*(HF*PISQ*2 - 8 - L12**2)
CATOT=M*(HF*PISQ + HF*(L12**2-L13**2-L23**2))
TFTOT=0
C---THE FINITE TERMS NOT PROPORTIONAL TO TREE-LEVEL RESULT
TOT=(FPAL2*CA*CQ*4/EMSQ)*(2*LEIV(P,P(1,6)))
CFTOT=CFTOT+CF*CFSUB*TOT
CATOT=CATOT+CF*CASUB*TOT
TFTOT=TFTOT+CF*TFSUB*TOT
TOT=(FPAL2*CA*CQ*4/EMSQ)*(EMSQ/2*ERTV(P))
CFTOT=CFTOT+CF*CFSUB*TOT
CATOT=CATOT+CF*CASUB*TOT
TFTOT=TFTOT+CF*TFSUB*TOT
ELSEIF (METYPE.EQ.2) THEN
CALL CONVERT(P)
CALL SETUP
CALL SFMEE(3,1,V)
CFTOT=V/CF
CATOT=0
TFTOT=0
ELSE
STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN VIRTHR!'
ENDIF
C---THE SUBTRACTION
CFTOT=CFTOT+M*(-HF*PISQ*2 + 10 + L12**2 - 3D0/2*(2*L12))
CATOT=CATOT+M*(-HF*PISQ + 50D0/9 - HF*(L12**2-L13**2-L23**2)
$ - HF*(11D0/6)*(L13+L23) - 3D0/2*(HF*(L13+L23-2*L12)))
TFTOT=TFTOT+M*( - 16D0/9 - HF*(-2D0/3)*(L13+L23))
C---THE TOTAL
V=CF*CFTOT+CA*CATOT+TR*NF*TFTOT
CFSUB=CFTOT/(CF*V)
CASUB=CATOT/(CF*V)
TFSUB=TFTOT/(CF*V)
GGSUB=CF*(CF*CFSUB+CA*CASUB)
QQSUB=CF*TR*TFSUB
QPSUB=CF*TR*(NF-1)*TFSUB
END
C-----------------------------------------------------------------------
DOUBLE PRECISION FUNCTION ERTV(P)
IMPLICIT NONE
C---RETURN THE ERT F FUNCTION FOR VIRTUAL TERMS THAT ARE NOT
C TRIVIALLY PROPORTIONAL TO TREE-LEVEL
DOUBLE PRECISION P(4,7),DOT,R,X,Y,DILOG,EMSQ,
$ Y12,Y13,Y23,L12,L13,L23
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
R(X,Y)=LOG(X)*LOG(Y)-LOG(X)*LOG(1-X)-LOG(Y)*LOG(1-Y)+PISQ/6
$ -DILOG(X)-DILOG(Y)
EMSQ=DOT(P,5,5)
Y12=2*DOT(P,1,2)/EMSQ
Y13=2*DOT(P,1,3)/EMSQ
Y23=2*DOT(P,2,3)/EMSQ
L12=LOG(Y12)
L13=LOG(Y13)
L23=LOG(Y23)
CFSUB=Y12/(Y12+Y13)+Y12/(Y12+Y23)+(Y12+Y23)/Y13+(Y12+Y13)/Y23
$ +L13*(4*Y12**2+2*Y12*Y13+4*Y12*Y23+Y13*Y23)/(Y12+Y23)**2
$ +L23*(4*Y12**2+2*Y12*Y23+4*Y12*Y13+Y13*Y23)/(Y12+Y13)**2
$ -2*((Y12**2+(Y12+Y13)**2)/(Y13*Y23)*R(Y12,Y23)
$ +(Y12**2+(Y12+Y23)**2)/(Y13*Y23)*R(Y12,Y13)
$ +(Y13**2+Y23**2)/(Y13*Y23*(Y13+Y23))
$ -2*L12*(Y12**2/(Y13+Y23)**2+2*Y12/(Y13+Y23)))
CASUB=L13*Y13/(Y12+Y23)+L23*Y23/(Y12+Y13)
$ +((Y12**2+(Y12+Y13)**2)/(Y13*Y23)*R(Y12,Y23)
$ +(Y12**2+(Y12+Y23)**2)/(Y13*Y23)*R(Y12,Y13)
$ +(Y13**2+Y23**2)/(Y13*Y23*(Y13+Y23))
$ -2*L12*(Y12**2/(Y13+Y23)**2+2*Y12/(Y13+Y23)))
$ -R(Y13,Y23)*(Y13/Y23+Y23/Y13+2*Y12/(Y13*Y23))
TFSUB=0
ERTV=CF*(CF*CFSUB+CA*CASUB)
CFSUB=CFSUB/ERTV
CASUB=CASUB/ERTV
ERTV=16*PISQ/EMSQ*ERTV
END
C-----------------------------------------------------------------------
DOUBLE PRECISION FUNCTION LEIV(P,K)
IMPLICIT NONE
C---RETURN THE PART OF THE ONE-LOOP HADRONIC TENSOR NOT TRIVIALLY
C PROPORTIONAL TO TREE-LEVEL ONE CONTRACTED WITH A TENSOR -K(MU)K(NU)
DOUBLE PRECISION P(4,7),K(4),DOT,R,X,Y,Z,DILOG,EMSQ,Y12,Y13,Y23,
$ L12,L13,L23,A,B,C,D,P1K,P2K,P3K,KK,
$ CF1,CF2,CF3,CF4,CA1,CA2,CA3,CA4
PARAMETER (Z=0)
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
R(X,Y)=LOG(X)*LOG(Y)-LOG(X)*LOG(1-X)-LOG(Y)*LOG(1-Y)+PISQ/6
$ -DILOG(X)-DILOG(Y)
EMSQ=DOT(P,5,5)
Y12=2*DOT(P,1,2)/EMSQ
Y13=2*DOT(P,1,3)/EMSQ
Y23=2*DOT(P,2,3)/EMSQ
L12=LOG(Y12)
L13=LOG(Y13)
L23=LOG(Y23)
C---FIRST DECOMPOSE K AS A*P1+B*P2+C*P3+D*E_(MU,P1,P2,P3)/Q^2
P1K=P(4,1)*K(4)-P(3,1)*K(3)-P(2,1)*K(2)-P(1,1)*K(1)
P2K=P(4,2)*K(4)-P(3,2)*K(3)-P(2,2)*K(2)-P(1,2)*K(1)
P3K=P(4,3)*K(4)-P(3,3)*K(3)-P(2,3)*K(2)-P(1,3)*K(1)
KK=K(4)*K(4)-K(3)*K(3)-K(2)*K(2)-K(1)*K(1)
A=(P3K*Y12+P2K*Y13-P1K*Y23)/(Y12*Y13*EMSQ)
B=(P1K*Y23+P3K*Y12-P2K*Y13)/(Y12*Y23*EMSQ)
C=(P2K*Y13+P1K*Y23-P3K*Y12)/(Y23*Y13*EMSQ)
D=SQRT(MAX(Z,4*(KK/EMSQ-A*B*Y12-A*C*Y13-B*C*Y23)/(-Y12*Y13*Y23)))
C---THEN CALCULATE THE CONTRACTIONS WITH -P1(MU)P1(NU) ETC
CF1=R(Y13,Y12)*Y12-R(Y23,Y12)*Y12**2/Y13
$ -L12*Y12*Y13/(Y13+Y23)**2-L23*Y12+Y12/2*(Y23-Y13)/(Y23+Y13)
CA1=-R(Y13,Y12)*Y12/2+R(Y23,Y12)*Y12**2/Y13/2+R(Y23,Y13)*Y12/2
$ +L12*Y12*Y13/(Y13+Y23)**2/2-Y12*Y23/(Y23+Y13)/2
CF2=R(Y23,Y12)*Y12-R(Y13,Y12)*Y12**2/Y23
$ -L12*Y12*Y23/(Y13+Y23)**2-L13*Y12+Y12/2*(Y13-Y23)/(Y13+Y23)
CA2=-R(Y23,Y12)*Y12/2+R(Y13,Y12)*Y12**2/Y23/2+R(Y23,Y13)*Y12/2
$ +L12*Y12*Y23/(Y13+Y23)**2/2-Y12*Y13/(Y13+Y23)/2
CF3=R(Y13,Y12)*Y12+R(Y23,Y12)*Y12-2*L12*Y12
$ -L13*Y12*Y23/2/(Y12+Y23)**2*(1+Y12+Y23)
$ -L23*Y12*Y13/2/(Y12+Y13)**2*(1+Y12+Y13)
$ -Y12/2*(Y13/(Y12+Y13)+Y23/(Y12+Y23))
CA3=-R(Y13,Y12)*Y12/2-R(Y23,Y12)*Y12/2+R(Y23,Y13)*Y12
$ +L12*Y12+L13*Y12*Y13/(Y12+Y23)/2+L23*Y12*Y23/(Y12+Y13)/2
CF4=R(Y13,Y12)*Y12/4/Y23*(Y12-Y12**2-Y12*Y23+2*Y12**2*Y23
$ +2*Y12*Y23**2+Y23**3)
$ +R(Y23,Y12)*Y12/4/Y13*(Y12-Y12**2-Y12*Y13+2*Y12**2*Y13
$ +2*Y12*Y13**2+Y13**3)
$ +L12*Y12**2/4/(Y13+Y23)*(Y13+Y23-2*Y13*Y23)
$ -L13*Y12/8/(Y12+Y23)*Y13*(Y12+Y13)*(Y23-2*Y12)
$ -L23*Y12/8/(Y12+Y13)*Y23*(Y12+Y23)*(Y13-2*Y12)
$ +Y12/8*(Y13+Y23-2*Y13*Y23)
CA4=-R(Y13,Y12)*Y12/8/Y23*(Y12-Y12**2-Y12*Y23+2*Y12**2*Y23
$ +2*Y12*Y23**2+Y23**3)
$ -R(Y23,Y12)*Y12/8/Y13*(Y12-Y12**2-Y12*Y13+2*Y12**2*Y13
$ +2*Y12*Y13**2+Y13**3)
$ +R(Y13,Y23)*Y12/8*((1-Y13)**2+(1-Y23)**2)
$ -L12*Y12**2/8/(Y13+Y23)*(Y13+Y23-2*Y13*Y23)
$ -L13*Y12/8*Y13*(Y12+Y13)-L23*Y12/8*Y23*(Y12+Y23)
$ -Y12/8*(Y13+Y23-2*Y13*Y23)
C---AND COMBINE THEM WITH THE APPROPRIATE WEIGHTS
CFSUB=(A-B)*(A-C)*CF1+(B-C)*(B-A)*CF2+(C-A)*(C-B)*CF3+D*D*CF4
CASUB=(A-B)*(A-C)*CA1+(B-C)*(B-A)*CA2+(C-A)*(C-B)*CA3+D*D*CA4
TFSUB=0
LEIV=CF*(CF*CFSUB+CA*CASUB)
CFSUB=CFSUB/LEIV
CASUB=CASUB/LEIV
LEIV=16*PISQ*LEIV
END
C-----------------------------------------------------------------------
SUBROUTINE SUBTHR(PERM,P,Q,S,JAC,*)
IMPLICIT NONE
C---GENERATE A TWO-PARTON STATE FROM A THREE-PARTON STATE,
C CALCULATE THE JACOBIAN FACTOR FOR THE CORRESPONDING CHANNEL,
C AND (IF PERM.GT.0) THE APPROXIMATE MATRIX-ELEMENT.
INTEGER NPERM,PERM,IPERM,IJF,KF,I,J,K,M
DOUBLE PRECISION P(4,7),Q(4,7),S,JAC,DOT,EMSQ,DEN,Y,ZI,ZJ,OY,ZTI
PARAMETER (NPERM=2)
DIMENSION IPERM(3,NPERM),IJF(NPERM),KF(NPERM)
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
INTEGER NPOW1,NPOW2
DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ
COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
DATA IPERM/1,3,2, 2,3,1/
DATA IJF/1,2/
DATA KF/2,1/
IF (ABS(PERM).GT.NPERM.OR.PERM.EQ.0)
$ STOP 'PERM TOO BIG IN SUBTHR!'
C---FIND WHICH PARTONS TO TREAT
I=IPERM(1,ABS(PERM))
J=IPERM(2,ABS(PERM))
K=IPERM(3,ABS(PERM))
C---FIND THE KINEMATIC VARIABLES
EMSQ=2*(DOT(P,I,J)+DOT(P,J,K)+DOT(P,K,I))
IF (EMSQ.LE.0) CALL WARNER('SUBTHR',1,*999)
DEN=2/EMSQ
Y=DOT(P,I,J)*DEN
IF (Y.LE.0.OR.Y.GE.1) CALL WARNER('SUBTHR',2,*999)
ZI=DOT(P,I,K)*DEN
IF (ZI.LE.0.OR.ZI.GE.1) CALL WARNER('SUBTHR',3,*999)
ZJ=DOT(P,J,K)*DEN
IF (ZJ.LE.0.OR.ZJ.GE.1) CALL WARNER('SUBTHR',4,*999)
OY=1/(1-Y)
ZTI=ZI*OY
C---COPY INTO Q, REPLACING K BY KTILDE AND I BY IJTILDE
DO M=1,4
Q(M,IJF(ABS(PERM)))=P(M,I)+P(M,J)-Y*OY*P(M,K)
Q(M,KF(ABS(PERM)))=OY*P(M,K)
Q(M,5)=P(M,5)
Q(M,6)=P(M,6)
Q(M,7)=P(M,7)
ENDDO
C---CALCULATE THE CORRESPONDING JACOBIAN FACTOR
CALL JACTWO(Q,JAC)
JAC=JAC/(4*NPOW1**2*(1-ZI)*(1-ZJ)*Y**XPOW1)*(2-ZI-ZJ)
C---INCLUDE A PRIORI CHANNEL WEIGHTS
JAC=JAC/2
IF (PERM.LT.0) RETURN
C---CALCULATE WEIGHT FOR QUARK-GLUON SPLITTING FUNCTION
CALL MATTWO(Q,S)
S=CF*S*(2/(1-ZTI*(1-Y))-(1+ZTI))
C---INCLUDE OVERALL FACTOR
DEN=1/(Y*EMSQ)
S=S*16*PISQ*DEN
RETURN
999 RETURN 1
END
C-----------------------------------------------------------------------
SUBROUTINE SUBFOR(PERM,P,Q,S,JAC,KEEP,*)
IMPLICIT NONE
C---GENERATE A THREE-PARTON STATE FROM A FOUR-PARTON STATE,
C CALCULATE THE JACOBIAN FACTOR FOR THE CORRESPONDING CHANNEL,
C AND (IF PERM.GT.0) THE APPROXIMATE MATRIX-ELEMENT.
C IF ANY DOT PRODUCTS ARE BELOW CUTOFF, EVENT IS REJECTED,
C BY JUMPING TO THE SPECIFIED LABEL
C IF ANY DOT PRODUCTS ARE BELOW CUTUP, EVENT IS MARKED BY KEEP=.TRUE.
INTEGER NPERM,PERM,IPERM,IJF,KF,LF,NPERM3,I,J,K,L,M,N
LOGICAL KEEP
DOUBLE PRECISION P(4,7),Q(4,7),S,S1,S2,S3,JAC,DOT,EMSQ,QUSQ,T,
$ DEN,Y,ZI,ZJ,OY,ZTI,ZTJ,OX1,OX2,JTMP,STMP,QTMP(4,7),V(4),VV
PARAMETER (NPERM=10,NPERM3=2)
DIMENSION IPERM(4,NPERM),IJF(NPERM),KF(NPERM),LF(NPERM)
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
INTEGER NPOW1,NPOW2
DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ
COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
DATA IPERM/1,3,2,4, 1,3,4,2, 1,4,2,3, 1,4,3,2,
$ 2,4,1,3, 2,4,3,1, 2,3,1,4, 2,3,4,1,
$ 3,4,1,2, 3,4,2,1/
DATA IJF/1,1,1,1,2,2,2,2,3,3/
DATA KF/2,3,2,3,1,3,1,3,1,2/
DATA LF/3,2,3,2,3,1,3,1,2,1/
IF (ABS(PERM).GT.NPERM.OR.PERM.EQ.0)
$ STOP 'PERM TOO BIG IN SUBFOR!'
C---FIND WHICH PARTONS TO TREAT
I=IPERM(1,ABS(PERM))
J=IPERM(2,ABS(PERM))
K=IPERM(3,ABS(PERM))
L=IPERM(4,ABS(PERM))
C---FIND THE KINEMATIC VARIABLES
EMSQ=2*(DOT(P,I,J)+DOT(P,J,K)+DOT(P,K,I))
IF (EMSQ.LE.0) CALL WARNER('SUBFOR',1,*999)
DEN=2/EMSQ
Y=DOT(P,I,J)*DEN
IF (Y.LE.0.OR.Y.GE.1) CALL WARNER('SUBFOR',2,*999)
ZI=DOT(P,I,K)*DEN
IF (ZI.LE.0.OR.ZI.GE.1) CALL WARNER('SUBFOR',3,*999)
ZJ=DOT(P,J,K)*DEN
IF (ZJ.LE.0.OR.ZJ.GE.1) CALL WARNER('SUBFOR',4,*999)
OY=1/(1-Y)
ZTI=ZI*OY
ZTJ=ZJ*OY
C---COPY INTO Q, REPLACING K BY KTILDE AND I BY IJTILDE
DO M=1,4
Q(M,IJF(ABS(PERM)))=P(M,I)+P(M,J)-Y*OY*P(M,K)
Q(M,KF(ABS(PERM)))=OY*P(M,K)
Q(M,LF(ABS(PERM)))=P(M,L)
Q(M,5)=P(M,5)
Q(M,6)=P(M,6)
Q(M,7)=P(M,7)
ENDDO
C---REENFORCE THE MINIMUM CUTOFF ON ALL PAIR MASSES
QUSQ=DOT(Q,5,5)
OX1=2*DOT(Q,2,3)/QUSQ
OX2=2*DOT(Q,1,3)/QUSQ
IF (OX1.LT.CUTOFF.OR.OX2.LT.CUTOFF.OR.1-OX1-OX2.LT.CUTOFF)
$ RETURN 1
IF (OX1.LT.CUTUP.OR.OX2.LT.CUTUP.OR.1-OX1-OX2.LT.CUTUP)
$ KEEP=.TRUE.
C---CALCULATE THE CORRESPONDING JACOBIAN FACTOR
JAC=0
DO M=1,NPERM3
CALL SUBTHR(-M,Q,QTMP,STMP,JTMP,*999)
JAC=JAC+JTMP
ENDDO
IF (ABS(PERM).LE.8) THEN
JAC=JAC/(2*NPOW2**2*(1-ZI)*Y**XPOW2)
ELSE
JAC=JAC/(4*NPOW2**2*(1-ZI)*(1-ZJ)*Y**XPOW2)*(2-ZI-ZJ)
ENDIF
JAC=JAC*QUSQ/EMSQ
C---INCLUDE A PRIORI CHANNEL WEIGHTS
JAC=JAC/12
IF (ABS(PERM).GT.8) JAC=JAC*2
IF (PERM.LT.0) RETURN
C---CALCULATE WEIGHT FOR QUARK-GLUON SPLITTING FUNCTION
IF (PERM.LE.8) THEN
CALL MATTHR(Q,S1)
S1=S1*(2/(1-ZTI*(1-Y))-(1+ZTI))
IF (KF(PERM).EQ.3) THEN
S1=S1*HF*CA
ELSE
S1=S1*(CF-HF*CA)
ENDIF
C---THEN PERMUTE Q
IF (PERM.LE.4) THEN
DO M=1,4
T=Q(M,1)
Q(M,1)=Q(M,2)
Q(M,2)=Q(M,3)
Q(M,3)=T
ENDDO
ELSE
DO M=1,4
T=Q(M,3)
Q(M,3)=Q(M,2)
Q(M,2)=T
ENDDO
ENDIF
ELSE
S1=0
ENDIF
C---CALCULATE WEIGHT FOR GLUON DECAY
DEN=1/(Y*EMSQ)
CALL MATTHR(Q,S2)
DO M=1,4
V(M)=ZTI*P(M,I)-ZTJ*P(M,J)
ENDDO
VV=-ZTI*ZTJ*Y*EMSQ
CALL CONTHR(Q,V,VV,S3)
S3=4*ZTI*ZTJ*S3
IF (PERM.LE.8) THEN
S2=S2-S3
IF (MOD(PERM-1,4).LE.1) THEN
S2=S2*HF*TR*(NF-HF)*2
ELSE
S2=S2*HF*TR
ENDIF
ELSE
S2=S2*(2/(1-ZTI*(1-Y))+2/(1-ZTJ*(1-Y))-4)+S3
S2=S2*HF*CA
ENDIF
C---INCLUDE SYMMETRY FACTORS
S1=S1/2
IF (PERM.LE.8) THEN
S2=S2/4
ELSE
S2=S2/2
ENDIF
C---INCLUDE OVERALL FACTOR
S=(S1+S2)*16*PISQ*DEN
IF (PERM.LE.8) THEN
CASUB=S1/(S1+S2)
TFSUB=(1-CASUB)/(TR*CF)*ONF
IF (KF(PERM).EQ.3) THEN
CFSUB=0
CASUB=CASUB/(CA*CF)
ELSE
CFSUB=CASUB/((CF-CA/2)*CF)
CASUB=CASUB/((CA-2*CF)*CF)
ENDIF
ELSE
CASUB=1/(CA*CF)
CFSUB=0
TFSUB=0
ENDIF
RETURN
999 RETURN 1
END
C-----------------------------------------------------------------------
FUNCTION ERTA(P,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE ERT A FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION ERTA,P(4,7),DOT,
$ S12,S13,S14,S23,S24,S34,S134,S234,S
S12=2*DOT(P,I,J)
S13=2*DOT(P,I,K)
S14=2*DOT(P,I,L)
S23=2*DOT(P,J,K)
S24=2*DOT(P,J,L)
S34=2*DOT(P,K,L)
S134=S13+S14+S34
S234=S23+S24+S34
S=S12+S13+S14+S23+S24+S34
ERTA=(S12*S34**2-S13*S24*S34+S14*S23*S34+3*S12*S23*S34+
$ 3*S12*S14*S34+4*S12**2*S34-S13*S23*S24+2*S12*S23*S24-
$ S13*S14*S24-2*S12*S13*S24+2*S12**2*S24+S14*S23**2+
$ 2*S12*S23**2+S14**2*S23+4*S12*S14*S23+4*S12**2*S23+
$ 2*S12*S14**2+2*S12*S13*S14+4*S12**2*S14+2*S12**2*S13+
$ 2*S12**3)/(2*S13*S134*S234*S24)+
$ (S24*S34+S12*S34+S13*S24-S14*S23+S12*S13)/(S13*S134**2)+
$ 2*S23*(S-S13)/(S13*S134*S24)+
$ S34/(2*S13*S24)
END
C-----------------------------------------------------------------------
FUNCTION ERTB(P,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE ERT B FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION ERTB,P(4,7),DOT,
$ S12,S13,S14,S23,S24,S34,S123,S124,S134,S234,S
S12=2*DOT(P,I,J)
S13=2*DOT(P,I,K)
S14=2*DOT(P,I,L)
S23=2*DOT(P,J,K)
S24=2*DOT(P,J,L)
S34=2*DOT(P,K,L)
S123=S12+S13+S23
S124=S12+S14+S24
S134=S13+S14+S34
S234=S23+S24+S34
S=S12+S13+S14+S23+S24+S34
ERTB=(S12*S24*S34+S12*S14*S34-S13*S24**2+S13*S14*S24+
$ 2*S12*S14*S24)/(S13*S134*S23*S14)+
$ S12*(S+S34)*S124/(S134*S234*S14*S24)-
$ (2*S13*S24+S14**2+S13*S23+2*S12*S13)/(S13*S134*S14)+
$ S12*S123*S124/(2*S13*S14*S23*S24)
END
C-----------------------------------------------------------------------
FUNCTION ERTC(P,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE ERT C FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION ERTC,P(4,7),DOT,
$ S12,S13,S14,S23,S24,S34,S134,S234
S12=2*DOT(P,I,J)
S13=2*DOT(P,I,K)
S14=2*DOT(P,I,L)
S23=2*DOT(P,J,K)
S24=2*DOT(P,J,L)
S34=2*DOT(P,K,L)
S134=S13+S14+S34
S234=S23+S24+S34
ERTC=-(5*S12*S34**2+2*S12*S24*S34+2*S12*S23*S34+2*S12*S14*S34+
$ 2*S12*S13*S34+4*S12**2*S34-S13*S24**2+S14*S23*S24+
$ S13*S23*S24+S13*S14*S24-S12*S14*S24-S13**2*S24-
$ 3*S12*S13*S24-S14*S23**2-S14**2*S23+S13*S14*S23-
$ 3*S12*S14*S23-S12*S13*S23)/(4*S134*S234*S34**2)+
$ (3*S12*S34**2-3*S13*S24*S34+3*S12*S24*S34+3*S14*S23*S34-
$ S13*S24**2-S12*S23*S34+6*S12*S14*S34+2*S12*S13*S34-
$ 2*S12**2*S34+S14*S23*S24-3*S13*S23*S24-2*S13*S14*S24+
$ 4*S12*S14*S24+2*S12*S13*S24+3*S14*S23**2+2*S14**2*S23+
$ 2*S14**2*S12+2*S12**2*S14+6*S12*S14*S23-2*S12*S13**2-
$ 2*S12**2*S13)/(4*S13*S134*S234*S34)+
$ (2*S12*S34**2-2*S13*S24*S34+S12*S24*S34+4*S13*S23*S34+
$ 4*S12*S14*S34+2*S12*S13*S34+2*S12**2*S34-S13*S24**2+
$ 3*S14*S23*S24+4*S13*S23*S24-2*S13*S14*S24+4*S12*S14*S24+
$ 2*S12*S13*S24+2*S14*S23**2+4*S13*S23**2+2*S13*S14*S23+
$ 2*S12*S14*S23+4*S12*S13*S23+2*S12*S14**2+4*S12**2*S13+
$ 4*S12*S13*S14+2*S12**2*S14)/(4*S13*S134*S24*S34)
ERTC=ERTC-
$ (S12*S34**2-2*S14*S24*S34-2*S13*S24*S34-S14*S23*S34+
$ S13*S23*S34+S12*S14*S34+2*S12*S13*S34-2*S14**2*S24-
$ 4*S13*S14*S24-4*S13**2*S24-S14**2*S23-S13**2*S23+
$ S12*S13*S14-S12*S13**2)/(2*S13*S34*S134**2)+
$ (S12*S34**2-4*S14*S24*S34-2*S13*S24*S34-2*S14*S23*S34-
$ 4*S13*S23*S34-4*S12*S14*S34-4*S12*S13*S34-2*S13*S14*S24+
$ 2*S13**2*S24+2*S14**2*S23-2*S13*S14*S23-S12*S14**2-
$ 6*S12*S13*S14-S12*S13**2)/(4*S34**2*S134**2)
END
C-----------------------------------------------------------------------
FUNCTION ERTD(P,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE ERT D FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION ERTD,P(4,7),DOT,
$ S12,S13,S14,S23,S24,S34,S123,S134
S12=2*DOT(P,I,J)
S13=2*DOT(P,I,K)
S14=2*DOT(P,I,L)
S23=2*DOT(P,J,K)
S24=2*DOT(P,J,L)
S34=2*DOT(P,K,L)
S123=S12+S13+S23
S134=S13+S14+S34
ERTD=(S13*S23*S34+S12*S23*S34-S12**2*S34+S13*S23*S24+
$ 2*S12*S23*S24-S14*S23**2+S12*S13*S24+S12*S14*S23+
$ S12*S13*S14)/(S13**2*S123**2)-
$ (S12*S34**2-S13*S24*S34+S12*S24*S34-S14*S23*S34-
$ S12*S23*S34-S13*S24**2+S14*S23*S24-S13*S23*S24-
$ S13**2*S24+S14*S23**2)/(S13**2*S123*S134)
END
C-----------------------------------------------------------------------
FUNCTION ERTE(P,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE ERT E FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION ERTE,P(4,7),DOT,
$ S12,S13,S14,S23,S24,S34,S123,S124,S134,S234
S12=2*DOT(P,I,J)
S13=2*DOT(P,I,K)
S14=2*DOT(P,I,L)
S23=2*DOT(P,J,K)
S24=2*DOT(P,J,L)
S34=2*DOT(P,K,L)
S123=S12+S13+S23
S124=S12+S14+S24
S134=S13+S14+S34
S234=S23+S24+S34
ERTE=(S12*S23*S34-S12*S24*S34+S12*S14*S34+S12*S13*S34+S13*S24**2-
$ S14*S23*S24+S13*S23*S24+S13*S14*S24+S13**2*S24-S14*S23**2-
$ S14**2*S23-S13*S14*S23)/(S13*S23*S123*S134)-
$ S12*(S12*S34-S23*S24-S13*S24-S14*S23-
$ S14*S13)/(S13*S23*S123**2)-
$ (S14+S13)*(S24+S23)*S34/(S13*S23*S134*S234)
END
C-----------------------------------------------------------------------
FUNCTION LEIA(P,Q,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU)
C WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION LEIA,P(4,7),Q(4),DOT,
$ P12,P13,P14,P23,P24,P34,P15,P25,P35,Q2,
$ P1K,P2K,P3K,P4K,KK
P12=2*DOT(P,I,J)
P13=2*DOT(P,I,K)
P14=2*DOT(P,I,L)
P23=2*DOT(P,J,K)
P24=2*DOT(P,J,L)
P34=2*DOT(P,K,L)
P15=P23+P24+P34
P25=P13+P14+P34
P35=P12+P14+P24
Q2=P12+P13+P14+P23+P24+P34
P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1)
P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1)
P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1)
P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1)
KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1)
LEIA = 0
LEIA = LEIA + P1K*P2K * ( - 2*Q2*P12*P25 - Q2*P14*
+ P25 - Q2*P23*P25 + P12*P13*P25 + P12*P24*P25 + 2*P12*P25*
+ P34 - P13*P23*P25 - 2*P13*P24*P25 - P13*P25*P34 + 2*P14*P15*
+ P24 - 2*P14*P23*P25 - P14*P24*P25 - 2*P15*P23*P25 - 2*P15*P24
+ *P25 - 2*P15*P25*P34 - P24*P25*P34 )
LEIA = LEIA + P1K*P3K * ( - Q2*P12 + Q2*P24 + P12*
+ P13 + P14*P24 + 2*P15*P24 - P24*P34 )*P25
LEIA = LEIA + P1K*P4K * ( - Q2*P12 - Q2*P23 + P12*
+ P13 + P12*P34 - P14*P23 )*P25
LEIA = LEIA + P1K**2 * ( Q2*P23 + Q2*P24 - P13*P23
+ - P13*P24 - P24*P34 )*P25
LEIA = LEIA + P2K*P3K * ( - Q2*P12*P25 - Q2*P14*
+ P25 - 2*P12*P15*P25 + P12*P24*P25 + P12*P25*P34 + 2*P14*P15*
+ P24 - P14*P23*P25 - 4*P15*P23*P25 - 2*P15*P24*P25 - 2*P15*P25
+ *P34 )
LEIA = LEIA + P2K*P4K * ( - Q2*P12*P25 + Q2*P13*
+ P25 + P12*P24*P25 + 2*P13*P15*P25 + P13*P23*P25 - P13*P25*P34
+ + 2*P14*P15*P24 - 2*P15*P23*P25 - 2*P15*P24*P25 )
LEIA = LEIA + P2K**2 * ( Q2*P13 + Q2*P14 + 2*P13*P15
+ - P13*P24 - P13*P34 - P14*P24 + 2*P15*P34 )*P25
LEIA = LEIA + P3K*P4K * ( - 2*P12*P15 + P12*P34 - P13*
+ P24 - P14*P23 - 2*P15*P23 - P15*P25 )*P25
LEIA = LEIA + P3K**2 * ( P14 + 2*P15 )*P24*P25
LEIA = LEIA + P4K**2 * ( P13*P23*P25 )
c$$$ LEIA = LEIA + KK * ( 2*Q2*P12*P14*P25 + 2*Q2*P12*P23
c$$$ + *P25 + 2*Q2*P12**2*P25 - 2*Q2*P14*P15*P24 + 2*Q2*P14*
c$$$ + P23*P25 + P12*P13*P25*P34 + 4*P12*P15*P23*P25 + 2*P12*P15*P25
c$$$ + *P34 + P12*P24*P25*P34 - P13*P14*P23*P25 - 2*P13*P14*P24*P25
c$$$ + - 2*P13*P15*P24*P25 - 2*P13*P23*P24*P25 - P13*P24**2*P25 -
c$$$ + P13**2*P24*P25 + 2*P14*P15*P23*P25 - P14*P23*P24*P25 + 4*P15*
c$$$ + P23*P24*P25 + 4*P15*P23*P25*P34 + 4*P15*P23**2*P25 + 2*P15*
c$$$ + P24*P25*P35 + P15*P25**2*P34 )/4
LEIA = LEIA/(P15*P25**2*P13*P24)
END
C-----------------------------------------------------------------------
FUNCTION LEIB(P,Q,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU)
C WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION LEIB,P(4,7),Q(4),DOT,
$ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2,
$ P1K,P2K,P3K,P4K,KK
P12=2*DOT(P,I,J)
P13=2*DOT(P,I,K)
P14=2*DOT(P,I,L)
P23=2*DOT(P,J,K)
P24=2*DOT(P,J,L)
P34=2*DOT(P,K,L)
P15=P23+P24+P34
P25=P13+P14+P34
P35=P12+P14+P24
P45=P12+P13+P23
Q2=P12+P13+P14+P23+P24+P34
P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1)
P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1)
P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1)
P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1)
KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1)
LEIB = 0
LEIB = LEIB + P1K*P2K * ( - 4*Q2*P14*P24*P45 - 2*P12
+ *P15*P25*P34 - P12*P15*P25*P35 - P12*P15*P25*P45 + 8*P13*P15*
+ P23*P24 + 4*P13*P15*P24*P34 + 4*P13*P15*P24**2 + 4*P14*P15*
+ P23*P24 + 2*P14*P15*P24*P45 + 2*P14*P24*P25*P45 - 4*P15*P23*
+ P24*P25 + 4*P15*P23*P24*P34 )
LEIB = LEIB + P1K*P3K * ( - P12*P14*P15*P25 - 2*P12*
+ P14*P24*P45 + P12*P15*P24*P25 + 4*P12*P15*P24*P34 - P12**2*
+ P15*P25 - 4*P13*P15*P24**2 + 4*P14*P15*P23*P24 )
LEIB = LEIB + P1K*P4K * ( - P12*P13*P15*P25 - 2*P12*
+ P14*P24*P45 + P12*P15*P23*P25 - P12**2*P15*P25 + 4*P13*P15*
+ P23*P24 + 4*P14*P15*P23*P24 )
LEIB = LEIB + P1K**2 * ( 2*Q2*P14*P24*P45 - 2*P12*P14*
+ P24*P45 + P12*P15**2*P25 - 2*P14*P24*P25*P45 - 4*P15*P23*P24*
+ P34 )
LEIB = LEIB + P2K*P3K * ( 4*P12*P13*P15*P24 + P12*P14*
+ P15*P25 - 2*P12*P14*P24*P45 - P12*P15*P24*P25 - P12**2*P15*
+ P25 + 4*P13*P15*P23*P24 + 4*P13*P15*P24**2 )
LEIB = LEIB + P2K*P4K * ( - 4*P12*P13*P15*P24 + P12*
+ P13*P15*P25 - 8*P12*P14*P15*P24 - 2*P12*P14*P24*P45 - P12*P15
+ *P23*P25 - 4*P12*P15*P24*P34 - P12**2*P15*P25 - 4*P13*P14*P15
+ *P24 + 4*P13*P15*P23*P24 + 4*P13*P15*P24**2 - 4*P13**2*P15*
+ P24 )
LEIB = LEIB + P2K**2 * ( 2*Q2*P14*P24*P45 - 2*P12*P14*
+ P24*P45 + P12*P15*P25**2 - 4*P13*P14*P15*P24 - 4*P13*P15*P24*
+ P34 - 4*P13**2*P15*P24 - 2*P14*P15*P24*P45 )
LEIB = LEIB + P3K*P4K * ( - 2*P12*P25 - 4*P14*P24 )
+ *P12*P15
LEIB = LEIB + P3K**2 * ( - 4*P12*P14*P15*P24 )
c$$$ LEIB = LEIB + KK * ( 2*Q2*P12*P14*P24*P45 + Q2*
c$$$ + P12**2*P15*P25 - 2*Q2*P13*P15*P23*P24 - 2*Q2*P14*P15*P23*
c$$$ + P24 + 2*Q2*P15*P23*P24*P25 + P12*P13*P14*P15*P25 - 4*P12*
c$$$ + P13*P15*P23*P24 - 2*P12*P13*P15*P24*P34 + 2*P12*P14*P15*P24*
c$$$ + P34 + 4*P12*P14*P15*P24**2 + P12*P15*P23*P24*P25 - 2*P12*P15*
c$$$ + P23*P24*P34 + 2*P12*P15*P24**2*P34 - 2*P13*P14*P15*P23*P24 +
c$$$ + 2*P13*P14*P15*P24**2 - 2*P13*P15*P23*P24**2 - 2*P13*P15*
c$$$ + P24**3 + 2*P13**2*P15*P24**2 + 2*P14*P15*P23*P24**2 + 2*P14*
c$$$ + P15*P23**2*P24 - 2*P14**2*P15*P23*P24 - 2*P15**2*P23*P24*P25
c$$$ + )/2
LEIB = LEIB/(2*P13*P14*P15*P23*P24*P25)
END
C-----------------------------------------------------------------------
FUNCTION LEIC(P,Q,I,J,K,L)
IMPLICIT NONE
C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU)
C WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L
DOUBLE PRECISION LEIC,P(4,7),Q(4),DOT,
$ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2,
$ P1K,P2K,P3K,P4K,KK
P12=2*DOT(P,I,J)
P13=2*DOT(P,I,K)
P14=2*DOT(P,I,L)
P23=2*DOT(P,J,K)
P24=2*DOT(P,J,L)
P34=2*DOT(P,K,L)
P15=P23+P24+P34
P25=P13+P14+P34
P35=P12+P14+P24
P45=P12+P13+P23
Q2=P12+P13+P14+P23+P24+P34
P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1)
P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1)
P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1)
P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1)
KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1)
LEIC = 0
LEIC = LEIC + P1K*P2K * ( - 4*P12*P13*P15*P25*P34 + 6*
+ P12*P13*P24*P25*P34 - 2*P12*P14*P15*P25*P34 - 2*P12*P14*P24*
+ P25*P34 - 2*P12*P15*P25*P34**2 + 2*P12*P24*P25*P34**2 + 4*P13
+ *P14*P15*P24*P34 - 2*P13*P14*P15*P25*P34 - 3*P13*P14*P23*P24*
+ P25 - P13*P14*P24**2*P25 + P13*P14**2*P15*P24 - 2*P13*P15*P23
+ *P25*P34 + P13*P15*P24*P25*P34 - P13*P15*P24*P34**2 - P13*P15
+ *P25*P34**2 + P13*P23*P24*P25*P34 - 2*P13*P24**2*P25*P34 + 6*
+ P13**2*P14*P15*P24 - P13**2*P23*P24*P25 + P13**2*P24*P25*P34
+ - 3*P13**2*P24**2*P25 + P13**3*P15*P24 - P14*P15*P23*P25*P34
+ - 2*P14*P15*P24*P25*P34 - 2*P14*P15*P25*P34**2 - 3*P14*P23*
+ P24*P25*P34 + P14*P24*P25*P34**2 - 2*P14*P24**2*P25*P34 -
+ P14**2*P15*P25*P34 - P14**2*P24*P25*P34 + P15*P23*P25*P34**2
+ + P15*P24*P25*P34**2 + P15*P25*P34**3 + P24*P25*P34**3 -
+ P24**2*P25*P34**2 )
LEIC = LEIC + P1K*P3K * ( 3*P12*P13*P34 - P12*P14*P34
+ - P12*P34**2 + P13*P14*P24 + 2*P13*P15*P34 - P13**2*P24 -
+ P14*P15*P34 + 3*P14*P24*P34 + P15*P34**2 - 3*P24*P34**2 )
+ *P24*P25
LEIC = LEIC + P1K*P4K * ( - 2*P12*P13*P15*P34 + 3*P12*
+ P13*P24*P34 - P12*P14*P15*P34 - P12*P14*P24*P34 - 3*P12*P15*
+ P34**2 + 3*P12*P24*P34**2 - P13*P14*P23*P24 - P13*P15*P23*P34
+ - 2*P13*P15*P24*P34 - 2*P13*P23*P24*P34 + P13**2*P23*P24 + 2
+ *P14*P15*P23*P34 - 3*P14*P23*P24*P34 + P23*P24*P34**2 )*P25
LEIC = LEIC + P1K**2 * ( 2*P13*P15 - 3*P13*P23 - 3*P13*
+ P24 + P14*P15 + P14*P23 + P14*P24 + 3*P15*P34 + P23*P34 - 3*
+ P24*P34 )*P24*P25*P34
LEIC = LEIC + P2K*P3K * ( - 2*P12*P13*P15*P25*P34 +
+ P12*P13*P24*P25*P34 - P12*P14*P15*P25*P34 - 3*P12*P14*P24*P25
+ *P34 - P12*P15*P25*P34**2 + P12*P24*P25*P34**2 + 2*P13*P14*
+ P15*P24*P34 - P13*P14*P15*P25*P34 - P13*P14*P23*P24*P25 - 2*
+ P13*P14*P24*P25*P34 + P13*P14*P24**2*P25 - 2*P13*P14**2*P15*
+ P24 - 4*P13*P15*P23*P25*P34 - 2*P13*P15*P25*P34**2 + 2*P13**2
+ *P14*P15*P24 - 2*P14*P15*P23*P25*P34 - 3*P14*P15*P24*P25*P34
+ - 2*P14*P15*P25*P34**2 - 3*P14*P23*P24*P25*P34 - P14*P24**2*
+ P25*P34 - P14**2*P24*P25*P34 )
LEIC = LEIC + P2K*P4K * ( - P12*P13*P15*P25*P34 + 2*
+ P12*P13*P24*P25*P34 - 2*P12*P14*P15*P25*P34 - 2*P12*P14*P24*
+ P25*P34 - 2*P12*P24*P25*P34**2 + 4*P13*P14*P15*P24*P34 - P13*
+ P14*P15*P25*P34 + P13*P14*P24*P25*P34 - 2*P13*P15*P23*P25*P34
+ - 3*P13*P15*P24*P25*P34 - P13*P15*P25*P34**2 + 3*P13*P23*P24
+ *P25*P34 - 2*P13*P24*P25*P34**2 + P13*P24**2*P25*P34 + 2*
+ P13**2*P14*P15*P24 - 2*P13**2*P15*P24*P34 + 2*P13**2*P15*P25*
+ P34 + P13**2*P23*P24*P25 - P13**2*P24**2*P25 - 2*P13**3*P15*
+ P24 - 2*P14*P15*P23*P25*P34 - 4*P14*P15*P24*P25*P34 )
LEIC = LEIC + P2K**2 * ( 2*P13*P14*P15 + P13*P14*P24 +
+ P13*P15*P34 - P13*P24*P34 + 2*P13**2*P15 - P13**2*P24 + 2*P14
+ *P15*P34 + 2*P14*P24*P34 + 2*P14**2*P15 + 2*P14**2*P24 )
+ *P25*P34
LEIC = LEIC + P3K*P4K * ( - P12*P13*P15*P34 + 2*P12*
+ P13*P24*P34 - P12*P14*P15*P34 - P12*P15*P34**2 + 2*P12*P24*
+ P34**2 - 2*P13*P14*P23*P24 - 2*P13*P15*P23*P34 - 2*P13*P24**2
+ *P34 - 2*P13**2*P24**2 - 2*P14*P23*P24*P34 )*P25
LEIC = LEIC + P3K**2 * ( 2*P13*P14*P24 + 2*P13*P15*P34
+ + 2*P14*P24*P34 )*P24*P25
LEIC = LEIC + P4K**2 * ( 2*P12*P15*P34 + 2*P13*P23*P24
+ + 2*P23*P24*P34 )*P13*P25
c$$$ LEIC = LEIC + KK * ( - 4*P12*P13*P14*P15*P24*P34 + 4*
c$$$ + P12*P13*P14*P15*P25*P34 + 3*P12*P13*P14*P23*P24*P25 - 2*P12*
c$$$ + P13*P14*P24*P25*P34 + P12*P13*P14*P24**2*P25 - P12*P13*P14**2
c$$$ + *P15*P24 + 4*P12*P13*P15*P23*P25*P34 + P12*P13*P15*P24*P34**2
c$$$ + + 2*P12*P13*P15*P25*P34**2 - 2*P12*P13*P23*P24*P25*P34 - 2*
c$$$ + P12*P13*P24*P25*P34**2 - 6*P12*P13**2*P14*P15*P24 + P12*
c$$$ + P13**2*P23*P24*P25 - 4*P12*P13**2*P24*P25*P34 + 3*P12*P13**2*
c$$$ + P24**2*P25 - P12*P13**3*P15*P24 + 2*P12*P14*P15*P23*P25*P34
c$$$ + + 4*P12*P14*P15*P24*P25*P34 + 5*P12*P14*P15*P25*P34**2 + 6*
c$$$ + P12*P14*P23*P24*P25*P34 + P12*P14*P24*P25*P34**2 + 4*P12*P14*
c$$$ + P24**2*P25*P34 + 2*P12*P14**2*P15*P25*P34 + 2*P12*P14**2*P24*
c$$$ + P25*P34 - P12*P15*P24*P25*P34**2 + P12*P15*P25*P34**3 - P12*
c$$$ + P23*P24*P25*P34**2 + 3*P12*P24**2*P25*P34**2 + 4*P12**2*P13*
c$$$ + P15*P25*P34 - 6*P12**2*P13*P24*P25*P34 + 2*P12**2*P14*P15*P25
c$$$ + *P34 + 2*P12**2*P14*P24*P25*P34 + 2*P12**2*P15*P25*P34**2 - 2
c$$$ + *P12**2*P24*P25*P34**2 )/4
c$$$ LEIC = LEIC + KK * ( - 2*P13*P14*P15*P23*P24*P34 + 2*
c$$$ + P13*P14*P15*P23*P25*P34 + P13*P14*P15*P24*P25*P34 - 4*P13*P14
c$$$ + *P15*P24**2*P34 + 3*P13*P14*P23*P24*P25*P34 - P13*P14*P23*
c$$$ + P24**2*P25 + P13*P14*P23**2*P24*P25 - 3*P13*P14*P24**2*P25*
c$$$ + P34 + 2*P13*P14**2*P15*P23*P24 + P13*P14**2*P23*P24*P25 + 2*
c$$$ + P13*P15*P23*P24*P25*P34 + 4*P13*P15*P23*P25*P34**2 + 4*P13*
c$$$ + P15*P23**2*P25*P34 - P13*P15*P24*P25*P34**2 + 3*P13*P15*
c$$$ + P24**2*P25*P34 - 3*P13*P23*P24**2*P25*P34 + 2*P13*P24**2*P25*
c$$$ + P34**2 - P13*P24**3*P25*P34 - 2*P13**2*P14*P15*P23*P24 - 2*
c$$$ + P13**2*P14*P15*P24**2 - P13**2*P14*P23*P24*P25 - P13**2*P14*
c$$$ + P24**2*P25 - 2*P13**2*P15*P24*P25*P34 + 2*P13**2*P15*P24**2*
c$$$ + P34 - P13**2*P23*P24**2*P25 + P13**2*P24**2*P25*P34 + P13**2*
c$$$ + P24**3*P25 + 2*P13**3*P15*P24**2 + P13**3*P24**2*P25 + 5*P14*
c$$$ + P15*P23*P24*P25*P34 + P14*P15*P23*P25*P34**2 + 2*P14*P15*
c$$$ + P23**2*P25*P34 + 4*P14*P15*P24**2*P25*P34 + 2*P14*P23*P24*P25
c$$$ + *P34**2 + )/4
c$$$ LEIC = LEIC + KK * ( P14*P23*P24**2*P25*P34 + 3*P14*
c$$$ + P23**2*P24*P25*P34 - P14**2*P15*P23*P25*P34 + 3*P14**2*P23*
c$$$ + P24*P25*P34 )/4
LEIC = LEIC/(2*P15*P25**2*P13*P34**2*P24)
END
C-----------------------------------------------------------------------
FUNCTION LEID(P,Q,II,JJ,KKK,LL)
IMPLICIT NONE
C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU)
C WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L,II,JJ,KKK,LL
DOUBLE PRECISION LEID,P(4,7),Q(4),DOT,
$ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2,
$ P1K,P2K,P3K,P4K,KK
C---CONVERT FROM ERT NOTATION (Q'QQ'BARQBAR) TO LEIDEN (QQBARQ'BARQ')
I=JJ
J=LL
K=KKK
L=II
P12=2*DOT(P,I,J)
P13=2*DOT(P,I,K)
P14=2*DOT(P,I,L)
P23=2*DOT(P,J,K)
P24=2*DOT(P,J,L)
P34=2*DOT(P,K,L)
P15=P23+P24+P34
P25=P13+P14+P34
P35=P12+P14+P24
P45=P12+P13+P23
Q2=P12+P13+P14+P23+P24+P34
P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1)
P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1)
P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1)
P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1)
KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1)
LEID = 0
LEID = LEID + P1K*P2K * ( - 2*P12*P25*P34 - 4*P13*P14*
+ P15 - 2*P13*P15*P34 + 2*P13*P24*P25 - 2*P14*P15*P34 + 2*P14*
+ P23*P25 )/P25
LEID = LEID + P1K*P3K * ( - P12*P34 - 2*P14*P24 + P24*P25 )
LEID = LEID + P1K*P4K * ( - P12*P34 - 2*P13*P23 + P23*P25 )
LEID = LEID + P1K**2 * ( P23 + P24 )*P34
LEID = LEID + P2K*P3K * ( - P12*P25*P34 - 4*P13*P14*P15
+ - 2*P13*P15*P34 + 3*P14*P15*P25 - 2*P14*P15*P34 - 2*P14*P24*
+ P25 )/P25
LEID = LEID + P2K*P4K * ( - P12*P25*P34 - 4*P13*P14*P15
+ + 3*P13*P15*P25 - 2*P13*P15*P34 - 2*P13*P23*P25 - 2*P14*P15*
+ P34 )/P25
LEID = LEID + P2K**2 * ( P13 + P14 )*P34
LEID = LEID + P3K*P4K * ( - 2*P12*P34 + 2*P13*P24 + 2*
+ P14*P23 )
LEID = LEID + P3K**2 * ( - 2*P14*P24 )
LEID = LEID + P4K**2 * ( - 2*P13*P23 )
c$$$ LEID = LEID + KK * ( Q2*P12*P25*P34 - Q2*P13*P24*P25
c$$$ + - Q2*P14*P23*P25 + 2*P12*P13*P15*P25 - P12*P13*P24*P25 - 2
c$$$ + *P12*P13**2*P15 + 2*P12*P14*P15*P25 - P12*P14*P23*P25 - 2*P12
c$$$ + *P14**2*P15 + P12*P25*P34**2 + P12**2*P25*P34 + 2*P13*P14*P23
c$$$ + *P25 + 2*P13*P14*P24*P25 + 2*P13*P15*P23*P25 + 2*P13*P23*P24*
c$$$ + P25 - P13*P24*P25*P34 - 2*P13**2*P15*P23 - 2*P13**2*P15*P24
c$$$ + + 2*P14*P15*P24*P25 + 2*P14*P23*P24*P25 - P14*P23*P25*P34 -
c$$$ + 2*P14**2*P15*P23 - 2*P14**2*P15*P24 )/(4*P25)
LEID = LEID/(P15*P25*P34**2)
END
C-----------------------------------------------------------------------
FUNCTION LEIE(P,Q,II,JJ,KKK,LL)
IMPLICIT NONE
C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU)
C WITH P(*,I)=P1, P(*,J)=P2, ETC
INTEGER I,J,K,L,II,JJ,KKK,LL
DOUBLE PRECISION LEIE,P(4,7),Q(4),DOT,
$ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2,
$ P1K,P2K,P3K,P4K,KK
C---CONVERT FROM ERT NOTATION (Q'QQ'BARQBAR) TO LEIDEN (QQBARQ'BARQ')
I=JJ
J=LL
K=KKK
L=II
P12=2*DOT(P,I,J)
P13=2*DOT(P,I,K)
P14=2*DOT(P,I,L)
P23=2*DOT(P,J,K)
P24=2*DOT(P,J,L)
P34=2*DOT(P,K,L)
P15=P23+P24+P34
P25=P13+P14+P34
P35=P12+P14+P24
P45=P12+P13+P23
Q2=P12+P13+P14+P23+P24+P34
P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1)
P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1)
P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1)
P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1)
KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1)
LEIE = 0
LEIE = LEIE + P1K*P2K * ( - 2*Q2*P13*P15*P25 + Q2*
+ P23*P25**2 + 2*P13*P14*P15*P25 - 4*P14*P15*P25*P45 - P14*P23*
+ P25**2 + 2*P14**2*P15*P45 - 2*P15*P23*P25*P34 + 2*P15*P25**2*
+ P45 - P23*P25**2*P45 )/(P13*P25*P45)
LEIE = LEIE + P1K*P3K * ( - 2*Q2*P13*P15 + Q2*P23*
+ P25 - 2*P12*P14*P15 - P14*P23*P25 - 2*P15*P23*P34 + 2*P15*P25
+ *P45 - P23*P25*P45 )/(P13*P45)
LEIE = LEIE + P1K*P4K * ( 2*Q2*P25 + 2*P13*P15 - 2*P14*
+ P25 - 2*P15*P25 - 2*P15*P34 - 2*P25*P45 )*P23/(P13*P45)
LEIE = LEIE + P1K**2 * ( - 2*P15*P23*P34 )/(P13*P45)
LEIE = LEIE + P2K*P3K * ( - 2*P12*P15*P25 + 2*P13*P15*
+ P25 + 2*P14*P15*P45 - 2*P23*P25**2 )*P14/(P13*P25*P45)
LEIE = LEIE + P2K*P4K * ( Q2*P23*P25**2 + 2*P12*P15*P25
+ *P34 - 2*P13*P15*P24*P25 + 2*P14*P15*P23*P25 - 2*P14*P15*P25*
+ P45 - P14*P23*P25**2 + 2*P14**2*P15*P45 - P15*P23*P25**2 )
+ /(P13*P25*P45)
LEIE = LEIE + P2K**2 * ( 2*P13*P15 - P23*P25 )*P14/(P13*P45)
LEIE = LEIE + P3K*P4K * ( - 2*Q2*P13*P15 + Q2*P23*
+ P25 - 2*P12*P14*P15 + 2*P13*P15*P23 - P14*P23*P25 - 3*P15*P23
+ *P25 + 2*P15*P25*P45 )/(P13*P45)
LEIE = LEIE + P3K**2 * ( - 2*P12*P15 - P23*P25 )*P14/(P13*P45)
LEIE = LEIE + P4K**2 * ( 2*P15*P23 )/P45
c$$$ LEIE = LEIE + KK * ( Q2*P13*P15*P24*P25 + Q2*P14*P23*
c$$$ + P25**2 - Q2*P14**2*P15*P45 + Q2*P15*P23*P25**2 + Q2*P23
c$$$ + *P25**2*P45 - Q2**2*P23*P25**2 - P12*P13*P15*P25*P34 + P12
c$$$ + *P14*P15*P23*P25 - P12*P15*P24*P25*P34 - P12*P15*P25*P34**2 -
c$$$ + P12**2*P15*P25*P34 - P13*P14*P15*P23*P25 - P13*P14*P15*P24*
c$$$ + P25 - P13*P15*P23*P24*P25 - P14*P15*P23*P24*P25 + P14*P15*P23
c$$$ + *P25*P34 + P14*P15*P25*P35*P45 - P15*P23*P25**2*P45 )
c$$$ + /(2*P13*P25*P45)
LEIE = LEIE/(P15*P25*P34)
END
C-----------------------------------------------------------------------
SUBROUTINE WARNER(SUBR,TYPE,*)
IMPLICIT NONE
C---GIVE WARNING MESSAGES ABOUT NUMERICAL ERRORS
INTEGER TYPE
LOGICAL FIRST
CHARACTER*(*) SUBR
INTEGER METYPE,NF
DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF
COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF,
$ NF,METYPE
DATA FIRST/.TRUE./
WRITE (*,'(A,I2,3A)')
$' Numerical error of type',TYPE,' in routine ',SUBR,'!'
IF (FIRST) THEN
WRITE (*,'(/A)')
$' This should be preventable by increasing the internal variable'
WRITE (*,'(A,1PE8.1,A)')
$' CUTOFF, whose present value is',CUTOFF,'. If this fails to'
WRITE (*,'(A)')
$' reduce the frequency of errors, something is probably more'
WRITE (*,'(A)')
$' seriously wrong, so please contact seymour@surya11.cern.ch'
WRITE (*,'(A)')
$' giving the date of this version and the latest value of ISEED.'
WRITE (*,'(A/)') ' Event generation has not been affected.'
FIRST=.FALSE.
ENDIF
RETURN 1
END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
subroutine convert(p)
implicit none
integer i,j
double precision p(4,7),plab
common /mom/ plab(4,10)
do j=1,4
do i=1,4
plab(j,i+2)=p(j,i)
enddo
plab(j,1)=p(j,7)
plab(j,2)=p(j,6)
enddo
end
C-----------------------------------------------------------------------
subroutine setup
implicit double precision(a-h,o-z)
parameter(pi=3.141592653589793238D0)
common /inppar/ njets,nloop,itmax1,itmax2,nshot1,nshot2,nshot3
common /zcup/czur,czul,czdr,czdl,czlr,czll,cznr,cznl
common /gcup/cgur,cgul,cgdr,cgdl,cglr,cgll,cgnr,cgnl
common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon
common /strong/ as,gs,renscl,qcdl,b0,nf
common /sighad/ sig0
common /order/iorder
common /tcuts/ ymin,dly
common /ecuts/ ycut,evis,jetalg
common /polar/ Peleft,Ppright
character*4 alg
logical done
data done/.false./
if (done) return
done=.true.
c---set the general parameters
rmz=91.17D0
rgz=2.516D0
w=0
c---set the polarisations
Peleft=0.5
Ppright=0.5
c---and a few other parameters
ymin=1D-8
as=2*pi
c
c iorder determines whether we keep all terms, or leading N, subleading N
c or nf pieces
c iorder=0 all, iorder=1, leading N, iorder=2 1/N**2, iorder=3 nf
c
iorder=0
c
c icon determines which interference terms to use
c icon=0 all, icon=1 gives gg, icon=2 gives gz, icon=3 gives zz
c
icon=1
c
c check that matrix elements vanish when pz is contracted in : momcon=1
c
momcon=0
c
c physics parameters:
c
c nutf : number of up type flavours
c ndtf : number of down type flavours
nutf=2
ndtf=3
nf=nutf+ndtf
if(iorder.eq.0)b0=(33D0-2D0*nf)/12D0/pi
if(iorder.eq.1)b0=33D0/12D0/pi
if(iorder.eq.2)b0=0D0
if(iorder.eq.3)b0=-2D0*nf/12D0/pi
c
c renscl : scale at which alpha_s is evaluated = scl* rmz
c can adjust alpha_s by changing scl or qcdl below
c
scl=1D0
renscl=scl*rmz
dly=log(ymin)
c
c qcdl : two loop lambda qcd for 4 flavours
c
qcdl=0.200D0
alfa=alfas(renscl,qcdl,2)
gs=sqrt(4.D0*pi*as)
am=1.D0/128D0
em=sqrt(4.D0*pi*am)
sw2=0.23D0
c determine the width of the z
scw=sqrt(sw2)*sqrt(1D0-sw2)
cgll=-1D0
cglr=-1D0
cgnl=0D0
cgnr=0D0
cgur=+2D0/3D0
cgdr=-1D0/3D0
cgul=+2D0/3D0
cgdl=-1D0/3D0
czll=(-0.5D0+sw2)/scw
czlr=sw2/scw
cznl=0.5D0/scw
cznr=0D0
czul=(+0.5D0-sw2*2D0/3D0)/scw
czur=(-sw2*2D0/3D0)/scw
czdl=(-0.5D0+sw2*1D0/3D0)/scw
czdr=(+sw2*1D0/3D0)/scw
c calculate the width of the z ( no top quark! )
rgee= em**2*(czll**2+czlr**2)*rmz/24.D0/pi
rgnn= em**2*(cznl**2+cznr**2)*rmz/24.D0/pi
rguu= 3D0*em**2*(czul**2+czur**2)*rmz/24.D0/pi
rgdd= 3D0*em**2*(czdl**2+czdr**2)*rmz/24.D0/pi
rglep= 3D0*rgee
rginv= 3D0*rgnn
rghad= (1D0+(alfa/pi)+1.41*(alfa/pi)**2)*(2D0*rguu+3D0*rgdd)
rghad= (2D0*rguu+3D0*rgdd)
rgz0 = rglep + rginv + rghad
rgz1 = rglep + rginv + rghad*(1D0+(alfa/pi))
rgz2 = rglep + rginv + rghad*(1D0+(alfa/pi)+1.41*(alfa/pi)**2)
* test=abs((rgztest-rgz)/rgz*100D0)
* write(*,40)rgz0
40 format(/,'lowest order z boson width = ',f6.3,' GeV ')
* write(*,41)rgz1
41 format(/,' first order z boson width = ',f6.3,' GeV ')
* write(*,42)rgz2
42 format(/,'second order z boson width = ',f6.3,' GeV ')
wtnano=389385.73D0
c
c hadronic cross section
c
sig0=wtnano*rghad*rgee*12D0*pi/rmz**2/rgz**2
if(jetalg.eq.1)alg=' E '
if(jetalg.eq.2)alg=' E0 '
if(jetalg.eq.3)alg=' P '
if(jetalg.eq.4)alg=' P0 '
if(jetalg.eq.5)alg=' D '
if(jetalg.eq.6)alg=' DE '
if(jetalg.eq.7)alg=' DP '
if(jetalg.eq.8)alg=' D1 '
write(*,*) '*************************************************',
. '*****************'
write(*,*)' ',
. ' prelim 13/6/94 '
write(*,11)njets,(njets-2+nloop),w
11 format(/,' e+e- to ',i2,' jets at O(alphas**',i1,') at ',
. ' root s = ',f7.2, ' GeV ')
if(iorder.eq.0)write(*,*)' including all colour pieces'
if(iorder.eq.1)write(*,*)' including only leading colour pieces'
if(iorder.eq.2)write(*,*)' including only subleading ',
. 'colour pieces'
if(iorder.eq.3)write(*,*)' including only nf pieces'
if(icon.eq.0)write(*,*)' including both photon and z exchange'
if(icon.eq.1)write(*,*)' including only photon exchange'
if(icon.eq.2)write(*,*)' including photon - z interference'
if(icon.eq.3)write(*,*)' including only z exchange'
if(momcon.eq.1)write(*,*)' checking momentum conservation => 0'
write(*,12)as
12 format(/,' alphas(Mz) = ',f8.4)
write(*,13)4D0*pi/em**2,sw2
13 format(/,' 1/alpha = ',f7.2,'; sin**2 theta = ',f6.3)
write(*,14)sig0
14 format(/,' sigma(hadronic) = ',f8.3,' nb')
write(*,15)alg,ycut
15 format(/,' with jets defined in',a3,'-scheme with ycut = ',f8.4)
write(*,16)rmz,rgz
16 format(/,' Z mass = ',f6.3,' GeV, Z width = ',f6.3,' GeV')
write(*,17)Peleft*100D0
17 format(/,' electron beam = ',f4.1,'% left handed')
write(*,18)Ppright*100D0
18 format(/,' positron beam = ',f4.1,'% right handed')
write(*,19)ymin
19 format(/,' parton resolution parameter = ',g12.5)
write(*,20)itmax1,nshot1/10,nshot2/10,nshot3/10
20 format(/,' warm up with ',i6,' sweeps of ',3i9,' shots ')
write(*,21)itmax2,nshot1,nshot2,nshot3
21 format(/,' final run of ',i6,' sweeps of ',3i9,' shots ')
write(*,*) '*************************************************',
. '*****************'
return
end
c***********************************************************************
c
c compute matrix element weight using library
c
subroutine sfmee(npar,nloop,eetot)
c***********************************************************************
implicit double precision (a-h,o-y)
implicit complex*16(z)
common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon
common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr
common /polar/ Peleft,Ppright
common /mom/ plab(4,10)
common /kinema/ p(4,7)
common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
common /linpr/ rin(7,7)
common /curent/ zrcur(2,2),zlcur(2,2)
common /check/ momcon
c
c transform vectors - all momenta outgoing
do 10 mu=1,4
do 5 n=1,npar
p(mu,n)=plab(mu,n+2)
5 continue
p(mu,6)=-plab(mu,1)
p(mu,7)=-plab(mu,2)
10 continue
c
c initialise the spinor and dot products
call setspn(npar)
c all quarks have positive energy so
negaen=1
c
c initialise the lepton current
do 20 ia=1,2
do 20 ib=1,2
if(momcon.eq.0)then
zrcur(ia,ib)=zkod(6,ia)*zko(7,ib)
elseif(momcon.eq.1)then
zrcur(ia,ib)=(0D0,0D0)
do 15 ipar=1,npar+1
zrcur(ia,ib)=zrcur(ia,ib)+zkod(ipar,ia)*zko(ipar,ib)
15 continue
endif
zlcur(ib,ia)=conjg(zrcur(ia,ib))
20 continue
c
if(npar.eq.2) call sfee2j(nloop,eetot)
if(npar.eq.3) call sfee3j(nloop,eetot)
if(npar.eq.4) call sfee4j(nloop,eetot)
c
if(icon.eq.1)eetot=Peleft*Ppright*eeggl
. +(1D0-Peleft)*(1D0-Ppright)*eeggr
if(icon.eq.2)eetot=Peleft*Ppright*eegzl
. +(1D0-Peleft)*(1D0-Ppright)*eegzr
if(icon.eq.3)eetot=Peleft*Ppright*eezzl
. +(1D0-Peleft)*(1D0-Ppright)*eezzr
c
return
end
***********************************************************************
c sfee2j calculates the process e+ e- --> g,z --> q + qb
c includes: no spin averages, final state statistics, couplings
c**********************************************************************
subroutine sfee2j(nloop,eetot)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon
common /strong/ as,gs,renscl,qcdl,b0,nf
common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr
common /polar/ Peleft,Ppright
dimension rl(3,2),rr(3,2)
dimension zhl(3,2),zhr(3,2)
c
call eeq2gn(0,nloop,rl(1,1),rr(1,1),zhl,zhr,rmz,rgz)
c
factor=em**4
c
eeggl=(nutf*rl(1,1)+ndtf*rl(1,2))*factor
eegzl=(nutf*rl(2,1)+ndtf*rl(2,2))*factor
eezzl=(nutf*rl(3,1)+ndtf*rl(3,2))*factor
eeggr=(nutf*rr(1,1)+ndtf*rr(1,2))*factor
eegzr=(nutf*rr(2,1)+ndtf*rr(2,2))*factor
eezzr=(nutf*rr(3,1)+ndtf*rr(3,2))*factor
eetotl=eeggl+eegzl+eezzl
eetotr=eeggr+eegzr+eezzr
eetot=Peleft*Ppright*eetotl
. +(1D0-Peleft)*(1D0-Ppright)*eetotr
end
c
c**********************************************************************
c sfee3j calculates the process e+ e- --> g,z --> q + qb + g
c includes: no spin averages, final state statistics, couplings
c**********************************************************************
subroutine sfee3j(nloop,eetot)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon
common /strong/ as,gs,renscl,qcdl,b0,nf
common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr
common /reshel/zeegghl,zeegzhl,zeezzhl,zeegghr,zeegzhr,
$ zeezzhr,zeetohl,zeetohr,zeetoh
common /polar/ Peleft,Ppright
dimension rl(3,2),rr(3,2)
dimension zhl(3,2),zhr(3,2)
c
call eeq2gn(1,nloop,rl(1,1),rr(1,1),zhl,zhr,rmz,rgz)
c
factor=em**4*gs**2
c
eeggl=(nutf*rl(1,1)+ndtf*rl(1,2))*factor
eegzl=(nutf*rl(2,1)+ndtf*rl(2,2))*factor
eezzl=(nutf*rl(3,1)+ndtf*rl(3,2))*factor
eeggr=(nutf*rr(1,1)+ndtf*rr(1,2))*factor
eegzr=(nutf*rr(2,1)+ndtf*rr(2,2))*factor
eezzr=(nutf*rr(3,1)+ndtf*rr(3,2))*factor
eetotl=eeggl+eegzl+eezzl
eetotr=eeggr+eegzr+eezzr
eetot=Peleft*Ppright*eetotl
. +(1D0-Peleft)*(1D0-Ppright)*eetotr
c---also apply the appropriate factors to the gluon helicity amplitudes
if (nloop.eq.0) then
zeegghl=(nutf*zhl(1,1)+ndtf*zhl(1,2))*factor
zeegzhl=(nutf*zhl(2,1)+ndtf*zhl(2,2))*factor
zeezzhl=(nutf*zhl(3,1)+ndtf*zhl(3,2))*factor
zeegghr=(nutf*zhr(1,1)+ndtf*zhr(1,2))*factor
zeegzhr=(nutf*zhr(2,1)+ndtf*zhr(2,2))*factor
zeezzhr=(nutf*zhr(3,1)+ndtf*zhr(3,2))*factor
zeetohl=zeegghl+zeegzhl+zeezzhl
zeetohr=zeegghr+zeegzhr+zeezzhr
if(icon.eq.1) then
zeetoh=Peleft*Ppright*zeegghl
. +(1D0-Peleft)*(1D0-Ppright)*zeegghr
elseif(icon.eq.2) then
zeetoh=Peleft*Ppright*zeegzhl
. +(1D0-Peleft)*(1D0-Ppright)*zeegzhr
elseif(icon.eq.3) then
zeetoh=Peleft*Ppright*zeezzhl
. +(1D0-Peleft)*(1D0-Ppright)*zeezzhr
else
zeetoh=Peleft*Ppright*zeetohl
. +(1D0-Peleft)*(1D0-Ppright)*zeetohr
endif
endif
end
c
c**********************************************************************
c sfee4j calculates the process e+ e- --> g,z --> q + qb + ( gg or qqb)
c includes: no spin averages, final state statistics, couplings
c**********************************************************************
subroutine sfee4j(nloop,eetot)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon
common /strong/ as,gs,renscl,qcdl,b0,nf
common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr
common /polar/ Peleft,Ppright
dimension rggl(3,2),rggr(3,2),rqql(3,5),rqqr(3,5)
dimension zhl(3,2),zhr(3,2)
common /order/iorder
common /perm/ips
c
do i=1,3
do j=1,5
rqql(i,j)=0D0
rqqr(i,j)=0D0
enddo
do k=1,2
rggl(i,k)=0D0
rggr(i,k)=0D0
enddo
enddo
sum=1
if(ips.eq.1)then
if(iorder.eq.0)then
call eeq2gn(2,nloop,rggl(1,1),rggr(1,1),zhl,zhr,rmz,rgz)
call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,1,3,2,4)
call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,2,4,1,3)
sum=2
elseif(iorder.eq.1)then
call eeq2gn(2,nloop,rggl(1,1),rggr(1,1),zhl,zhr,rmz,rgz)
elseif(iorder.eq.3)then
call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,1,3,2,4)
call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,2,4,1,3)
sum=2
endif
elseif(ips.eq.2)then
if(iorder.eq.0.or.iorder.eq.2)then
call eeq2gn(2,nloop,rggl(1,1),rggr(1,1),zhl,zhr,rmz,rgz)
endif
endif
c
factor=em**4*gs**4
c
eeggl=(+nutf*rggl(1,1)/2D0+ndtf*rggl(1,2)/2D0
. +nutf*rqql(1,1)/4D0/sum+ndtf*rqql(1,2)/4D0/sum
. +0.5D0*nutf*(nutf-1)*rqql(1,3)/sum
. +0.5D0*ndtf*(ndtf-1)*rqql(1,4)/sum
. +nutf*ndtf*rqql(1,5)/sum)*factor
eegzl=(+nutf*rggl(2,1)/2D0+ndtf*rggl(2,2)/2D0
. +nutf*rqql(2,1)/4D0/sum+ndtf*rqql(2,2)/4D0/sum
. +0.5D0*nutf*(nutf-1)*rqql(2,3)/sum
. +0.5D0*ndtf*(ndtf-1)*rqql(2,4)/sum
. +nutf*ndtf*rqql(2,5)/sum)*factor
eezzl=(+nutf*rggl(3,1)/2D0+ndtf*rggl(3,2)/2D0
. +nutf*rqql(3,1)/4D0/sum+ndtf*rqql(3,2)/4D0/sum
. +0.5D0*nutf*(nutf-1)*rqql(3,3)/sum
. +0.5D0*ndtf*(ndtf-1)*rqql(3,4)/sum
. +nutf*ndtf*rqql(3,5)/sum)*factor
eeggr=(+nutf*rggr(1,1)/2D0+ndtf*rggr(1,2)/2D0
. +nutf*rqqr(1,1)/4D0/sum+ndtf*rqqr(1,2)/4D0/sum
. +0.5D0*nutf*(nutf-1)*rqqr(1,3)/sum
. +0.5D0*ndtf*(ndtf-1)*rqqr(1,4)/sum
. +nutf*ndtf*rqqr(1,5)/sum)*factor
eegzr=(+nutf*rggr(2,1)/2D0+ndtf*rggr(2,2)/2D0
. +nutf*rqqr(2,1)/4D0/sum+ndtf*rqqr(2,2)/4D0/sum
. +0.5D0*nutf*(nutf-1)*rqqr(2,3)/sum
. +0.5D0*ndtf*(ndtf-1)*rqqr(2,4)/sum
. +nutf*ndtf*rqqr(2,5)/sum)*factor
eezzr=(+nutf*rggr(3,1)/2D0+ndtf*rggr(3,2)/2D0
. +nutf*rqqr(3,1)/4D0/sum+ndtf*rqqr(3,2)/4D0/sum
. +0.5D0*nutf*(nutf-1)*rqqr(3,3)/sum
. +0.5D0*ndtf*(ndtf-1)*rqqr(3,4)/sum
. +nutf*ndtf*rqqr(3,5)/sum)*factor
eetotl=eeggl+eegzl+eezzl
eetotr=eeggr+eegzr+eezzr
eetot=Peleft*Ppright*eetotl
. +(1D0-Peleft)*(1D0-Ppright)*eetotr
end
c
c**********************************************************************
c matrix element for: e+ e- --> g,z --> q qb + n gluons
c**********************************************************************
c information: the particles in plab are ordered as follows:
c 1,2 incoming particles 1: e+ and 2: e-
c 3,4 quark and anti-quark
c 5,.., 4+n gluons in final state.
c**********************************************************************
c output: the matrix elements for the subprocesses returned in
c left-handed contribution
c rl(3,2): first index are mixings ==> 1: gg 2: gz 3: zz
c second index are flavours ==> 1: u-type, 2: d-type.
c right-handed contribution
c rr(3,2): first index are mixings ==> 1: gg 2: gz 3: zz
c second index are flavours ==> 1: u-type, 2: d-type.
c not included: strong coupling constants!
c statistics, spin averages.
c**********************************************************************
c literature:
c lowest order : berends, giele & kuijf, np 321, 595, 1989
c next to lowest order : giele & glover, prd 46, 1980, 1992
c**********************************************************************
c function eeq2gn(n,nloop,rl,rr,rmz,rgz) calculates the matrix element
c squared of: e+ e- --> g,z --> q + qb + n gluons
c the input vectors are in /kinema/ p(4,7), where the sixth
c and seventh vectors are the e- and the e+ outgoing!.
c the first and second are the quark and the anti-quark.
c**********************************************************************
subroutine eeq2gn(n,nloop,rl,rr,zhl,zhr,rmz,rgz)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /esadb/ zsadb(2,2,2,4)
common /curent/ zrcur(2,2),zlcur(2,2)
common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
common /linpr/ rin(7,7)
common /tcuts/ ymin,dly
common /strong/ as,gs,renscl,qcdl,b0,nf
common /zcup/czur,czul,czdr,czdl,czlr,czll,cznr,cznl
common /gcup/cgur,cgul,cgdr,cgdl,cglr,cgll,cgnr,cgnl
common /order/iorder
common /perm/ips
parameter(pi=3.141592653589793238D0)
dimension col(2,2,0:2)
dimension rl(3,2),rr(3,2)
dimension zhl(3,2),zhr(3,2)
dimension zrr(2,4),zrl(2,4),zlr(2,4),zll(2,4)
dimension z1rr(2,4),z1rl(2,4),z1lr(2,4),z1ll(2,4)
dimension zrsadb(2,2,2,4),zlsadb(2,2,2,4)
dimension zr1adb(2,2,2,4),zl1adb(2,2,2,4)
col(1,1,0)=3D0
col(1,2,0)=0D0
col(2,1,0)=0D0
col(2,2,0)=0D0
col(1,1,1)=4D0
col(1,2,1)=0D0
col(2,1,1)=0D0
col(2,2,1)=0D0
*
if(ips.eq.1)then
if(iorder.eq.0.or.iorder.eq.1)then
c$$$ c2a=12D0
C---THIS MODIFICATION SENT BY NIGEL GLOVER TO AVOID SUMMING OVER PERMS
c2a=6D0
c2b=0D0
endif
elseif(ips.eq.2)then
if(iorder.eq.0.or.iorder.eq.2)then
c2a=-2D0/3D0
c2b=c2a
endif
endif
c$$$ col(1,1,2)=c2b
C---THIS MODIFICATION SENT BY NIGEL GLOVER TO AVOID SUMMING OVER PERMS
col(1,1,2)=c2a
col(1,2,2)=c2b
col(2,1,2)=c2b
col(2,2,2)=c2a
c
s=2D0*rin(6,7)
zprg=(0D0,-1D0)/s
zprz=(0D0,-1D0)/(s-rmz**2+(0D0,1D0)*rgz*rmz)
c
if (n.eq.0) then
factor=1D0
ihel=1
iperm=1
do 6 ia=1,2
do 6 ib=1,2
zrsadb(ia,ib,1,1)=zkod(1,ia)*zko(2,ib)
zlsadb(ib,ia,1,1)=conjg(zrsadb(ia,ib,1,1))
6 continue
endif
c
if (n.eq.1) then
factor=2D0
ihel=2
iperm=1
do 7 ia=1,2
do 7 ib=1,2
zrsadb(ia,ib,1,1)=+(zkod(1,ia)*zud(1,2)
. +zkod(3,ia)*zud(3,2))*zko(2,ib)/(zud(1,3)*zud(3,2))
zrsadb(ia,ib,1,2)=-zkod(1,ia)*(zko(3,ib)*zd(3,1)
. +zko(2,ib)*zd(2,1))/(zd(1,3)*zd(3,2))
zlsadb(ib,ia,1,1)=conjg(zrsadb(ia,ib,1,1))
zlsadb(ib,ia,1,2)=conjg(zrsadb(ia,ib,1,2))
7 continue
if(nloop.eq.1)then
call ert(1,2,3,alp,bep,gap,dep)
call ert(2,1,3,alm,bem,gam,dem)
do 17 ia=1,2
do 17 ib=1,2
zr1adb(ia,ib,1,1)=
. +alp*zkod(1,ia)*zud(1,2)*zko(2,ib)/(zud(1,3)*zud(3,2))
. +bep*zkod(3,ia)*zud(3,2)*zko(2,ib)/(zud(1,3)*zud(3,2))
. +gap*zkod(3,ia)*zd(1,3)*zko(1,ib)/(zud(1,3)*zd(3,2))
. +dep*zd(1,3)/(zud(1,3)*zd(2,1))
. *(zkod(1,ia)*zko(1,ib)+zkod(2,ia)*zko(2,ib)+zkod(3,ia)*zko(3,ib))
zr1adb(ia,ib,1,2)=
. -alm*zkod(1,ia)*zd(2,1)*zko(2,ib)/(zd(1,3)*zd(3,2))
. -bem*zkod(1,ia)*zd(3,1)*zko(3,ib)/(zd(1,3)*zd(3,2))
. -gam*zkod(2,ia)*zud(2,3)*zko(3,ib)/(zud(1,3)*zd(3,2))
. -dem*zud(2,3)/(zd(2,3)*zud(1,2))
. *(zkod(1,ia)*zko(1,ib)+zkod(2,ia)*zko(2,ib)+zkod(3,ia)*zko(3,ib))
zl1adb(ib,ia,1,1)=conjg(zr1adb(ia,ib,1,1))
zl1adb(ib,ia,1,2)=conjg(zr1adb(ia,ib,1,2))
17 continue
endif
endif
c
if (n.eq.2) then
factor=4D0
ihel=4
iperm=2
call f2sadb(3,4,1,1,2,3,4)
call f2sadb(4,3,2,1,3,2,4)
do 8 ih=1,ihel
do 8 ip=1,iperm
do 8 ia=1,2
do 8 ib=1,2
zrsadb(ia,ib,ip,ih)=zsadb(ia,ib,ip,ih)
zlsadb(ib,ia,ip,ih)=conjg(zsadb(ia,ib,ip,ih))
8 continue
endif
c
do 10 ih=1,ihel
do 10 ip=1,iperm
zrr(ip,ih)=zspv(zrsadb(1,1,ip,ih),zrcur(1,1))
zrl(ip,ih)=zspv(zlsadb(1,1,ip,ih),zrcur(1,1))
zlr(ip,ih)=zspv(zrsadb(1,1,ip,ih),zlcur(1,1))
zll(ip,ih)=zspv(zlsadb(1,1,ip,ih),zlcur(1,1))
c
c non-tree level loop amplitudes contracted into lepton current
c
if(nloop.eq.1.and.n.eq.1)then
z1rr(ip,ih)=zspv(zr1adb(1,1,ip,ih),zrcur(1,1))
z1rl(ip,ih)=zspv(zl1adb(1,1,ip,ih),zrcur(1,1))
z1lr(ip,ih)=zspv(zr1adb(1,1,ip,ih),zlcur(1,1))
z1ll(ip,ih)=zspv(zl1adb(1,1,ip,ih),zlcur(1,1))
endif
10 continue
c
do 15 i=1,3
do 15 j=1,2
rl(i,j)=0D0
rr(i,j)=0D0
15 continue
c
if(nloop.eq.0)then
tloop=1D0
elseif(nloop.eq.1)then
c
c one-loop contribution proportional to tree level
c
if(n.eq.0)then
y12=rin(1,2)/rin(6,7)
c$$$ dl12=log(y12/ymin)
c$$$ e0=-dl12**2+3D0*dl12/2D0+pi**2/6D0-1D0/2D0
c$$$ tloop=1D0+4D0*as/3D0/pi*e0
C---THIS MODIFICATION MADE BY MIKE BY ANALOGY WITH THE ONE BELOW
e0=pi**2/2-4
tloop=0D0+4D0*as/3D0/pi*e0
elseif(n.eq.1)then
y12=rin(1,2)/rin(6,7)
y13=rin(1,3)/rin(6,7)
y23=rin(2,3)/rin(6,7)
c$$$ dl12=log(y12/ymin)
c$$$ dl13=log(y13/ymin)
c$$$ dl23=log(y23/ymin)
c$$$ dlm2=log(renscl**2/rmz**2/ymin)
c$$$ e0=29D0/9D0+pi**2/3D0+3D0/4D0*dl13+3D0/4D0*dl23
c$$$ . -dl13**2-dl23**2
c$$$ e2=-dl12**2+3D0*dl12/2D0+pi**2/6D0-1D0/2D0
c$$$ e1=5D0/9D0
C---THIS MODIFICATION SENT BY NIGEL GLOVER TO GIVE ERT FINITE TERMS
dly12=log(y12)
dly13=log(y13)
dly23=log(y23)
dlm2=0D0
C---PI**2 COEFFICIENT CHANGED BY MIKE
e0=pi**2-4D0-dly13**2/2D0-dly23**2/2D0
e2=pi**2/2-4D0-dly12**2/2D0
e1=0D0
tloop=0D0+as*3D0/2D0/pi*(e0-e2/9D0-nf*e1/3D0)+as*b0*dlm2
if(iorder.eq.1)tloop=0D0+as*3D0/2D0/pi*e0+as*b0*dlm2
if(iorder.eq.2)tloop=0D0-as*3D0/2D0/pi*e2/9D0+as*b0*dlm2
if(iorder.eq.3)tloop=0D0-as*3D0/2D0/pi*nf*e1/3D0+as*b0*dlm2
elseif(n.eq.2)then
tloop=1D0
endif
endif
fac=factor*tloop
zprgg=zprg*conjg(zprg)
zprzg=zprz*conjg(zprg)
zprzz=zprz*conjg(zprz)
do 20 i=1,ihel
do 20 j1=1,iperm
do 20 j2=1,iperm
zrrrr=zrr(j1,i)*conjg(zrr(j2,i))
zrlrl=zrl(j1,i)*conjg(zrl(j2,i))
zlrlr=zlr(j1,i)*conjg(zlr(j2,i))
zllll=zll(j1,i)*conjg(zll(j2,i))
*
* left-handed contributions
*
rl(1,1)=rl(1,1)+fac*col(j1,j2,n)*zprgg*(
. +(zlrlr*cgur*cgur+zllll*cgul*cgul)*cgll*cgll )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zlrlr*czur*cgur+zllll*czul*cgul)*czll*cgll )
rl(2,1)=rl(2,1)+2D0*real(zz)
rl(3,1)=rl(3,1)+fac*col(j1,j2,n)*zprzz*(
. +(zlrlr*czur*czur+zllll*czul*czul)*czll*czll )
rl(1,2)=rl(1,2)+fac*col(j1,j2,n)*zprgg*(
. +(zlrlr*cgdr*cgdr+zllll*cgdl*cgdl)*cgll*cgll )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zlrlr*czdr*cgdr+zllll*czdl*cgdl)*czll*cgll )
rl(2,2)=rl(2,2)+2D0*real(zz)
rl(3,2)=rl(3,2)+fac*col(j1,j2,n)*zprzz*(
. +(zlrlr*czdr*czdr+zllll*czdl*czdl)*czll*czll )
*
* right-handed contributions
*
rr(1,1)=rr(1,1)+fac*col(j1,j2,n)*zprgg*(
. +(zrrrr*cgur*cgur+zrlrl*cgul*cgul)*cglr*cglr )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zrrrr*czur*cgur+zrlrl*czul*cgul)*czlr*cglr )
rr(2,1)=rr(2,1)+2D0*real(zz)
rr(3,1)=rr(3,1)+fac*col(j1,j2,n)*zprzz*(
. +(zrrrr*czur*czur+zrlrl*czul*czul)*czlr*czlr )
rr(1,2)=rr(1,2)+fac*col(j1,j2,n)*zprgg*(
. +(zrrrr*cgdr*cgdr+zrlrl*cgdl*cgdl)*cglr*cglr )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zrrrr*czdr*cgdr+zrlrl*czdl*cgdl)*czlr*cglr )
rr(2,2)=rr(2,2)+2D0*real(zz)
rr(3,2)=rr(3,2)+fac*col(j1,j2,n)*zprzz*(
. +(zrrrr*czdr*czdr+zrlrl*czdl*czdl)*czlr*czlr )
20 continue
*
* Equivalent expressions for the gluon interference term: (+-)+(-+)
*
if (n.eq.1.and.nloop.eq.0) then
j1=1
j2=1
zrrrr=zrr(j1,1)*conjg(zrr(j2,2))
zrlrl=zrl(j1,1)*conjg(zrl(j2,2))
zlrlr=zlr(j1,1)*conjg(zlr(j2,2))
zllll=zll(j1,1)*conjg(zll(j2,2))
*
* left-handed contributions
*
zhl(1,1)=fac*col(j1,j2,n)*zprgg*(
. +(zlrlr*cgur*cgur+zllll*cgul*cgul)*cgll*cgll )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zlrlr*czur*cgur+zllll*czul*cgul)*czll*cgll )
zhl(2,1)=2D0*real(zz)
zhl(3,1)=fac*col(j1,j2,n)*zprzz*(
. +(zlrlr*czur*czur+zllll*czul*czul)*czll*czll )
zhl(1,2)=fac*col(j1,j2,n)*zprgg*(
. +(zlrlr*cgdr*cgdr+zllll*cgdl*cgdl)*cgll*cgll )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zlrlr*czdr*cgdr+zllll*czdl*cgdl)*czll*cgll )
zhl(2,2)=2D0*real(zz)
zhl(3,2)=fac*col(j1,j2,n)*zprzz*(
. +(zlrlr*czdr*czdr+zllll*czdl*czdl)*czll*czll )
*
* right-handed contributions
*
zhr(1,1)=fac*col(j1,j2,n)*zprgg*(
. +(zrrrr*cgur*cgur+zrlrl*cgul*cgul)*cglr*cglr )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zrrrr*czur*cgur+zrlrl*czul*cgul)*czlr*cglr )
zhr(2,1)=2D0*real(zz)
zhr(3,1)=fac*col(j1,j2,n)*zprzz*(
. +(zrrrr*czur*czur+zrlrl*czul*czul)*czlr*czlr )
zhr(1,2)=fac*col(j1,j2,n)*zprgg*(
. +(zrrrr*cgdr*cgdr+zrlrl*cgdl*cgdl)*cglr*cglr )
zz=fac*col(j1,j2,n)*zprzg*(
. +(zrrrr*czdr*cgdr+zrlrl*czdl*cgdl)*czlr*cglr )
zhr(2,2)=2D0*real(zz)
zhr(3,2)=fac*col(j1,j2,n)*zprzz*(
. +(zrrrr*czdr*czdr+zrlrl*czdl*czdl)*czlr*czlr )
endif
c
c higher order contributions
c
if(nloop.eq.1)then
c
c interference between non-tree level structures and tree level
c
if(n.eq.1)then
fac=factor*3D0*as/4D0/pi
do 30 i=1,ihel
do 30 j1=1,iperm
do 30 j2=1,iperm
z1rrrr=z1rr(j1,i)*conjg(zrr(j2,i))
z1rlrl=z1rl(j1,i)*conjg(zrl(j2,i))
z1lrlr=z1lr(j1,i)*conjg(zlr(j2,i))
z1llll=z1ll(j1,i)*conjg(zll(j2,i))
*
* left-handed contributions
*
zz=fac*col(j1,j2,n)*zprgg*(
. +(z1lrlr*cgur*cgur+z1llll*cgul*cgul)*cgll*cgll )
rl(1,1)=rl(1,1)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzg*(
. +(z1lrlr*czur*cgur+z1llll*czul*cgul)*czll*cgll )
rl(2,1)=rl(2,1)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzz*(
. +(z1lrlr*czur*czur+z1llll*czul*czul)*czll*czll )
rl(3,1)=rl(3,1)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprgg*(
. +(z1lrlr*cgdr*cgdr+z1llll*cgdl*cgdl)*cgll*cgll )
rl(1,2)=rl(1,2)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzg*(
. +(z1lrlr*czdr*cgdr+z1llll*czdl*cgdl)*czll*cgll )
rl(2,2)=rl(2,2)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzz*(
. +(z1lrlr*czdr*czdr+z1llll*czdl*czdl)*czll*czll )
rl(3,2)=rl(3,2)+2D0*real(zz)
*
* right-handed contributions
*
zz=fac*col(j1,j2,n)*zprgg*(
. +(z1rrrr*cgur*cgur+z1rlrl*cgul*cgul)*cglr*cglr )
rr(1,1)=rr(1,1)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzg*(
. +(z1rrrr*czur*cgur+z1rlrl*czul*cgul)*czlr*cglr )
rr(2,1)=rr(2,1)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzz*(
. +(z1rrrr*czur*czur+z1rlrl*czul*czul)*czlr*czlr )
rr(3,1)=rr(3,1)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprgg*(
. +(z1rrrr*cgdr*cgdr+z1rlrl*cgdl*cgdl)*cglr*cglr )
rr(1,2)=rr(1,2)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzg*(
. +(z1rrrr*czdr*cgdr+z1rlrl*czdl*cgdl)*czlr*cglr )
rr(2,2)=rr(2,2)+2D0*real(zz)
zz=fac*col(j1,j2,n)*zprzz*(
. +(z1rrrr*czdr*czdr+z1rlrl*czdl*czdl)*czlr*czlr )
rr(3,2)=rr(3,2)+2D0*real(zz)
30 continue
endif
endif
c
end
c
c***********************************************************************
c subroutine ert(iq,ip,ik,al,be,ga,de) computes the functions multiplyin
c the different one-loop helicity structures
c***********************************************************************
c call ert(iq,ip,ik,al,be,ga,de)
c iq,ip and ik point to the momenta of the quark, antiquark and gluon
c al, be and ga are the functions multiplying the helicity structures
c***********************************************************************
subroutine ert(iq,ip,ik,al,be,ga,de)
implicit double precision(a-h,o-y)
parameter(pi=3.141592653589793238D0)
common /linpr/ rin(7,7)
common /order/iorder
c
yqp=rin(iq,ip)/rin(6,7)
yqk=rin(iq,ik)/rin(6,7)
ykp=rin(ik,ip)/rin(6,7)
omyqp=1D0-yqp
omyqk=1D0-yqk
omykp=1D0-ykp
dlqp=log(yqp)
dlqk=log(yqk)
dlkp=log(ykp)
dlomqp=log(omyqp)
dlomqk=log(omyqk)
dlomkp=log(omykp)
spqp=rsp(yqp)
spqk=rsp(yqk)
spkp=rsp(ykp)
rqpkp=dlqp*dlkp-dlqp*dlomqp-dlkp*dlomkp+pi**2/6D0-spqp-spkp
rqpqk=dlqp*dlqk-dlqp*dlomqp-dlqk*dlomqk+pi**2/6D0-spqp-spqk
rqkkp=dlqk*dlkp-dlqk*dlomqk-dlkp*dlomkp+pi**2/6D0-spqk-spkp
c
al0=-rqkkp
. -ykp*(4D0-3D0*ykp)/2D0/omykp**2*dlkp
. -ykp/2D0/omykp
be0=-rqkkp
. +(4D0-3D0*ykp)/2D0/omykp*dlkp
. +1D0
ga0=+ykp/2D0/omykp**2*dlkp
. +ykp/2D0/omykp
de0=+(yqp*ykp/2D0/omykp**2+yqp/omykp)*dlkp
. +yqp/2D0/omykp
c
al2=-rqpqk
. -omyqp**2/yqk**2*rqpkp
. -ykp/yqk*dlqp
. -(ykp*(4D0-3D0*ykp)/2D0/omykp**2+ykp**2/yqk/omykp)*dlkp
. -ykp/2D0/omykp
be2=-rqpqk
. +(yqp*omyqp/yqk**2+1D0/yqk)*rqpkp
. +(yqk/omyqp**2+omykp/yqk)*dlqp
. +((4D0-3D0*ykp)/2D0/omykp+ykp/yqk)*dlkp
. +yqk/omyqp
ga2=+ykp/yqk**2*rqpkp
. -(yqk/omyqp**2-1D0/yqk)*dlqp
. +(ykp/2D0/omykp**2+ykp/yqk/omykp)*dlkp
. +ykp/omyqp+ykp/2D0/omykp
de2=+yqp*omyqp/yqk**2*rqpkp
. +yqp/yqk*dlqp
. +(yqp*ykp/2D0/omykp**2+yqp*omyqp/yqk/omykp)*dlkp
. +yqp/2D0/omykp
c
C---THESE LINES REMOVED BY MIKE FOR ERT DEFINITION OF COLLINEAR PIECES
c$$$ al0=al0-3D0/4D0*dlqk-3D0/4D0*dlkp
c$$$ be0=be0-3D0/4D0*dlqk-3D0/4D0*dlkp
c$$$ al2=al2-3D0/2D0*dlqp
c$$$ be2=be2-3D0/2D0*dlqp
c
if(iorder.eq.0)then
al=al0-al2/9D0
be=be0-be2/9D0
ga=ga0-ga2/9D0
de=de0-de2/9D0
elseif(iorder.eq.1)then
al=al0
be=be0
ga=ga0
de=de0
elseif(iorder.eq.2)then
al=-al2/9D0
be=-be2/9D0
ga=-ga2/9D0
de=-de2/9D0
endif
c check : t0 and t2 should be zero
c t0=-al0+be0-ga0-2D0/yqp*de0
c t2=-al2+be2-ga2-2D0/yqp*de2
c write(*,*)t0,t2
end
c**********************************************************************
c subroutine f2sadb(i1,i2,index,j1,j2,j3,j4) calculates the four
c helicity amplitudes for the permutation of the gluons (i1,i2).
c the quark(q) having + helicity and the anti-quark(p) - helictity.
c the calculation is organized as follows:
c after calculating the various propagators, all the spinor
c currents from the form: ( propagator contracted with a spinor)
c are composed. the currents containing the q are dotted. those with
c the qb are undotted.
c example: zD12m2 is the spinor (q+g(i1)+g(i2))adb * g(i2)b
c after initializing various constants belonging to a certain helicity
c combination, e.g. zppx belongs to the combination where both
c gluons have helicity +, the currents zsadb are calculated.
c**********************************************************************
c calling conventions:
c i1,i2 = permutation of the gluons. because they are used
c as indices we added two.
c this way we have (i1,i2)= (3,4) or (4,3)
c index = permutation counter
c j1,j2,j3,j4 = indices to get the helicity currents on the
c right spot in the data structures sadb.
c also needed are commons /spinor/ and /dotpr/
c**********************************************************************
c output conventions:
c put all the helicity combinations of a certain
c permutation (i1,i2) in the data-structure zsadb(2,2,2,4)
c**********************************************************************
subroutine f2sadb(i1,i2,index,j1,j2,j3,j4)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /esadb/ zsadb(2,2,2,4)
common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
common /linpr/ rin(7,7)
dimension zD12m1(2),zD12m2(2),zD12mp(2),
. z12pm1(2),z12pm2(2),z12pmq(2)
c
pro12=2D0*rin(i1,i2)
proD12=pro12+2D0*(rin(1,i1)+rin(1,i2))
pro12p=pro12+2D0*(rin(2,i1)+rin(2,i2))
c
do 10 i=1,2
zD12m1(i)=zkod(1,i)*zud(1,i1)+zkod(i2,i)*zud(i2,i1)
zD12m2(i)=zkod(1,i)*zud(1,i2)+zkod(i1,i)*zud(i1,i2)
zD12mp(i)=zkod(1,i)*zud(1,2)+zkod(i1,i)*zud(i1,2)
. +zkod(i2,i)*zud(i2,2)
z12pm1(i)=zko(2,i)*zd(2,i1)+zko(i2,i)*zd(i2,i1)
z12pm2(i)=zko(2,i)*zd(2,i2)+zko(i1,i)*zd(i1,i2)
z12pmq(i)=zko(2,i)*zd(2,1)+zko(i1,i)*zd(i1,1)+zko(i2,i)*zd(i2,1)
10 continue
c
zppa=-1D0/(zud(1,i1)*zud(i1,i2)*zud(i2,2))
zpma=zd(1,i1)*zud(1,i2)/(pro12*proD12*zud(1,i1))
zpmb=zd(2,i1)*zud(2,i2)/(pro12*pro12p*zd(2,i2))
zpmc=1D0/(pro12*zud(1,i1)*zd(i2,2))
zmpa=zud(2,i1)**2/(pro12*pro12p*zud(2,i2))
zmpb=zd(1,i2)**2/(pro12*proD12*zd(1,i1))
zmpc=-zd(1,i2)*zud(2,i1)/(pro12*zd(1,i1)*zud(2,i2))
zmma=1D0/(zd(1,i1)*zd(i1,i2)*zd(i2,2))
c
do 20 ia=1,2
do 20 ib=1,2
zsadb(ia,ib,index,j1)=-zppa*zD12mp(ia)*zko(2,ib)
zsadb(ia,ib,index,j2)=-zpma*zD12m2(ia)*zko(2,ib)
. -zpmb*zkod(1,ia)*z12pm1(ib)
. -zpmc*zD12m2(ia)*z12pm1(ib)
zsadb(ia,ib,index,j3)=-zmpa*zkod(1,ia)*z12pm2(ib)
. -zmpb*zD12m1(ia)*zko(2,ib)
. -zmpc*zkod(1,ia)*zko(2,ib)
zsadb(ia,ib,index,j4)=-zmma*zkod(1,ia)*z12pmq(ib)
20 continue
c
end
c
c**********************************************************************
c matrix element for: e+ e- --> g,z --> q qb q qb
c**********************************************************************
c information: the particles in plab are ordered as follows:
c 1,2 : incoming e+ and e-
c 3,5 : outgoing quarks
c 4,6 : outgoing anti-quarks
c sw2: square of sin(weak angle)
c rmz: mass of z
c rgz: width of z
c
c output: the matrix elements for the subprocesses returned in
c rl(3,5) for left handed contribution
c rr(3,5) for right handed contribution
c first index 1 = gg
c 2 = gz
c 3 = zz
c second index 1 = same quarks up-type
c 2 = same quarks down-type
c 3 = different quarks up-type
c 4 = different quarks down-type
c 5 = different quarks up - down pairs
c**********************************************************************
c not included: strong coupling constants!
c statistics, spin averages.
c**********************************************************************
c calling procedure for subroutines:
c 1 : anti-particle
c 2 : particle
c 3,5: labels of quarks in plab array
c 4,6: labels of anti-quarks in plab array
c particles 1 and 2 are always flipped: we consider all particle
c to be outgoing.
c**********************************************************************
c most of the formulae in this program can be found in:
c exact expressions for processes involving a vector boson
c and up to five partons. (berends, giele, kuijf)
c**********************************************************************
C---MODIFIED BY MIKE TO ENABLE SIMPLE SUMMING OVER PERMS
C QUARK 1 HAS BEEN REPLACED BY VARIABLE I1 EVERYWHERE, ETC
subroutine eeq4(rl,rr,rmz,rgz,I1,I2,I3,I4)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /mom/ plab(4,10)
common /kinema/ p(4,7)
common /zqqqq/ zr(4,6),zl(4,6)
common /curent/ zrcur(2,2),zlcur(2,2)
dimension za(2,2),zb(2,2)
dimension rl(3,5),rr(3,5)
c
c fill h functions
c the problem of an odd number of quark-particles with negative
c energy is handled within the h-functions.
c handle ppmm and mmpp
call hpmmp(I1,I4,I3,I2,-1D0,za(1,1),zb(1,1))
zr(1,1)=zspv(zrcur(1,1),za(1,1))
zl(1,1)=zspv(zlcur(1,1),za(1,1))
zr(1,4)=zspv(zrcur(1,1),zb(1,1))
zl(1,4)=zspv(zlcur(1,1),zb(1,1))
call hpmmp(I3,I2,I1,I4,-1D0,za(1,1),zb(1,1))
zr(2,4)=zspv(zrcur(1,1),za(1,1))
zl(2,4)=zspv(zlcur(1,1),za(1,1))
zr(2,1)=zspv(zrcur(1,1),zb(1,1))
zl(2,1)=zspv(zlcur(1,1),zb(1,1))
c
c handle pmmp and mppm
call hpmmp(I1,I2,I3,I4,1D0,za(1,1),zb(1,1))
zr(1,2)=zspv(zrcur(1,1),za(1,1))
zl(1,2)=zspv(zlcur(1,1),za(1,1))
zr(1,3)=zspv(zrcur(1,1),zb(1,1))
zl(1,3)=zspv(zlcur(1,1),zb(1,1))
call hpmmp(I3,I4,I1,I2,1D0,za(1,1),zb(1,1))
zr(2,3)=zspv(zrcur(1,1),za(1,1))
zl(2,3)=zspv(zlcur(1,1),za(1,1))
zr(2,2)=zspv(zrcur(1,1),zb(1,1))
zl(2,2)=zspv(zlcur(1,1),zb(1,1))
c
c handle pmpm and mpmp
call hpmpm(I1,I2,I3,I4,1D0,za(1,1),zb(1,1))
zr(1,5)=zspv(zrcur(1,1),za(1,1))
zl(1,5)=zspv(zlcur(1,1),za(1,1))
zr(1,6)=zspv(zrcur(1,1),zb(1,1))
zl(1,6)=zspv(zlcur(1,1),zb(1,1))
call hpmpm(I1,I4,I3,I2,-1D0,za(1,1),zb(1,1))
zr(2,5)=zspv(zrcur(1,1),za(1,1))
zl(2,5)=zspv(zlcur(1,1),za(1,1))
zr(2,6)=zspv(zrcur(1,1),zb(1,1))
zl(2,6)=zspv(zlcur(1,1),zb(1,1))
call hpmpm(I3,I4,I1,I2,1D0,za(1,1),zb(1,1))
zr(3,5)=zspv(zrcur(1,1),za(1,1))
zl(3,5)=zspv(zlcur(1,1),za(1,1))
zr(3,6)=zspv(zrcur(1,1),zb(1,1))
zl(3,6)=zspv(zlcur(1,1),zb(1,1))
call hpmpm(I3,I2,I1,I4,-1D0,za(1,1),zb(1,1))
zr(4,5)=zspv(zrcur(1,1),za(1,1))
zl(4,5)=zspv(zlcur(1,1),za(1,1))
zr(4,6)=zspv(zrcur(1,1),zb(1,1))
zl(4,6)=zspv(zlcur(1,1),zb(1,1))
c
c insert different couplings to obtain subprocesses
c
call squar(1,1,rmz,rgz,rl(1,1),rr(1,1))
call squar(2,2,rmz,rgz,rl(1,2),rr(1,2))
call squar(1,3,rmz,rgz,rl(1,3),rr(1,3))
call squar(2,4,rmz,rgz,rl(1,4),rr(1,4))
call squar(1,2,rmz,rgz,rl(1,5),rr(1,5))
c
end
c
c**********************************************************************
c function squar calculates the matrix element squared of the
c e+ e- --> g,z --> q + qb + q + qb
c**********************************************************************
c desciption of the common zqqqq
c the first is the index for the squared matrix
c the second index are the different subamplitudes (6)
c the first and the second are indices in the squared matrix
c index 3=1 --> ++--, =2 --> +--+, =3 --> -++-, =4 --> --++
c index 3=5 --> +-+-, =6 --> -+-+
c**********************************************************************
c output : rl(3) 1: gg 2: gz 3: zz
c rr(3) 1: gg 2: gz 3: zz
c**********************************************************************
subroutine squar(ityp1,ityp2,rmz,rgz,rl,rr)
implicit double precision (a-h,o-y)
implicit complex*16(z)
parameter(c1=8D0,c2=-8D0/3D0)
common /zqqqq/ zr(4,6),zl(4,6)
dimension s2g(4,2),s2z(4,2),s4g(2,4),s4z(2,4)
dimension col4(4,4),col2(2,2),rl(3),rr(3)
common /zcup/czur,czul,czdr,czdl,czlr,czll,cznr,cznl
common /gcup/cgur,cgul,cgdr,cgdl,cglr,cgll,cgnr,cgnl
data col4/c1,c2,c1,c2,c2,c1,c2,c1,c1,c2,c1,c2,c2,c1,c2,c1/
data col2/c1,c1,c1,c1/
C---THESE LINES REMOVED BY MIKE TO FACILITATE SUMMATION OVER PERMS
C (THEY WERE SUPERFLUOUS ANYWAY, SINCE THE ARRAYS WERE EMPTY ALREADY)
c$$$c
c$$$c init answer
c$$$ rl(1)=0D0
c$$$ rl(2)=0D0
c$$$ rl(3)=0D0
c$$$ rr(1)=0D0
c$$$ rr(2)=0D0
c$$$ rr(3)=0D0
c
c couplings and all that (like energy dependent width for z propagator).
s=2D0*rinpro(6,7)
zprg=(0D0,-1D0)/s
zprz=(0D0,-1D0)/(s-rmz**2+(0D0,1D0)*s*rgz/rmz)
if(ityp1.eq.1)then
cg1r=cgur
cg1l=cgul
cz1r=czur
cz1l=czul
elseif(ityp1.eq.2)then
cg1r=cgdr
cg1l=cgdl
cz1r=czdr
cz1l=czdl
endif
if(ityp2.eq.1.or.ityp2.eq.3)then
cg2r=cgur
cg2l=cgul
cz2r=czur
cz2l=czul
elseif(ityp2.eq.2.or.ityp2.eq.4)then
cg2r=cgdr
cg2l=cgdl
cz2r=czdr
cz2l=czdl
endif
c
if(ityp1.eq.ityp2)then
delta=1D0
else
delta=0D0
endif
c
c s2(1,.) is ppmm
s2g(1,1)=delta*cg1r
s2z(1,1)=delta*cz1r
s2g(1,2)=delta*cg1l
s2z(1,2)=delta*cz1l
c s2(2,.) is pmmp
s2g(2,1)=cg1r
s2z(2,1)=cz1r
s2g(2,2)=cg2l
s2z(2,2)=cz2l
c s2(3,.) is mppm
s2g(3,1)=cg1l
s2z(3,1)=cz1l
s2g(3,2)=cg2r
s2z(3,2)=cz2r
c s2(4,.) is mmpp
s2g(4,1)=delta*cg1l
s2z(4,1)=delta*cz1l
s2g(4,2)=delta*cg1r
s2z(4,2)=delta*cz1r
c
zprgg=zprg*conjg(zprg)
zprzg=zprz*conjg(zprg)
zprzz=zprz*conjg(zprz)
do 10 i=1,4
ij=i
do 11 ix=1,2
do 12 iy=1,2
zrr=zr(ix,i)*conjg(zr(iy,i))
zll=zl(ix,i)*conjg(zl(iy,i))
c
c left-handed contributions
c
rl(1)=rl(1)+col2(ix,iy)*s2g(ij,ix)*s2g(ij,iy)*
. zprgg*cgll*cgll*zll
zz=col2(ix,iy)*s2g(ij,ix)*s2z(ij,iy)*
. zprzg*czll*cgll*zll
rl(2)=rl(2)+2D0*real(zz)
rl(3)=rl(3)+col2(ix,iy)*s2z(ij,ix)*s2z(ij,iy)*
. zprzz*czll*czll*zll
c
c right-handed contributions
c
rr(1)=rr(1)+col2(ix,iy)*s2g(ij,ix)*s2g(ij,iy)*
. zprgg*cglr*cglr*zrr
zz=col2(ix,iy)*s2g(ij,ix)*s2z(ij,iy)*
. zprzg*czlr*cglr*zrr
rr(2)=rr(2)+2D0*real(zz)
rr(3)=rr(3)+col2(ix,iy)*s2z(ij,ix)*s2z(ij,iy)*
. zprzz*czlr*czlr*zrr
12 continue
11 continue
10 continue
c
c s4(1,.) is pmpm
s4g(1,1)=cg1r
s4z(1,1)=cz1r
s4g(1,2)=delta*cg1r
s4z(1,2)=delta*cz1r
s4g(1,3)=cg2r
s4z(1,3)=cz2r
s4g(1,4)=delta*cg1r
s4z(1,4)=delta*cz1r
c s4(2,.) is mpmp
s4g(2,1)=cg1l
s4z(2,1)=cz1l
s4g(2,2)=delta*cg1l
s4z(2,2)=delta*cz1l
s4g(2,3)=cg2l
s4z(2,3)=cz2l
s4g(2,4)=delta*cg1l
s4z(2,4)=delta*cz1l
c
do 20 i=5,6
ij=(i-4)
do 21 ix=1,4
do 22 iy=1,4
zrr=zr(ix,i)*conjg(zr(iy,i))
zll=zl(ix,i)*conjg(zl(iy,i))
c
c left-handed contributions
c
rl(1)=rl(1)+col4(ix,iy)*s4g(ij,ix)*s4g(ij,iy)*
. zprgg*cgll*cgll*zll
zz=col4(ix,iy)*s4g(ij,ix)*s4z(ij,iy)*
. zprzg*czll*cgll*zll
rl(2)=rl(2)+2D0*real(zz)
rl(3)=rl(3)+col4(ix,iy)*s4z(ij,ix)*s4z(ij,iy)*
. zprzz*czll*czll*zll
c
c right-handed contributions
c
rr(1)=rr(1)+col4(ix,iy)*s4g(ij,ix)*s4g(ij,iy)*
. zprgg*cglr*cglr*zrr
zz=col4(ix,iy)*s4g(ij,ix)*s4z(ij,iy)*
. zprzg*czlr*cglr*zrr
rr(2)=rr(2)+2D0*real(zz)
rr(3)=rr(3)+col4(ix,iy)*s4z(ij,ix)*s4z(ij,iy)*
. zprzz*czlr*czlr*zrr
22 continue
21 continue
20 continue
c
end
c the problem of an odd number of quark-particles with negative
c energy is handled within the h-functions.
c**********************************************************************
c subroutine hpmpm calculates the h-function when the helicities of
c the quarks are +-+-.
c call the h-functions with:
c argument 1-4 : quark permutation
c argument 5 is the sign of the h-function (exchange of two
c identical fermions)
c argument 6 is where to put it in zessen
c argument 7 is where to put its complex conjugate.
c z1 and z2 are the 2*2 complex matrices which are the h-vectors
c in spinor language. z2 is always the c.c. of z1. when there
c are an odd number of quarks with negative energy we have to add
c an extra minus sign when complex conjugating. this sign is
c in the common /dotpr/ negaen and is calculated before through
c the call setspn(4,4) in squar
c**********************************************************************
subroutine hpmpm(i1,i2,i3,i4,sign,z1,z2)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
common /linpr/ rin(7,7)
dimension z1(2,2),z2(2,2)
c
pro34=2D0*rin(i3,i4)
pro134=pro34+2D0*(rin(i1,i3)+rin(i1,i4))
pro234=pro34+2D0*(rin(i2,i3)+rin(i2,i4))
c
zf1=-sign*zd(i1,i3)/(pro34*pro134)
zf2=-sign*zud(i2,i4)/(pro34*pro234)
do 10 ix=1,2
do 10 iy=1,2
z1(ix,iy)=+zf1*zko(i2,iy)*(zud(i1,i4)*zkod(i1,ix)+
. zud(i3,i4)*zkod(i3,ix))
. -zf2*zkod(i1,ix)*(zd(i2,i3)*zko(i2,iy)+
. zd(i4,i3)*zko(i4,iy))
10 continue
call congat(z1,z2)
end
c
c**********************************************************************
c subroutine hpmmp calculates the h-function when the helicities of
c the quarks are +--+. see also subroutine hpmpm
c**********************************************************************
subroutine hpmmp(i1,i2,i3,i4,sign,z1,z2)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
common /linpr/ rin(7,7)
dimension z1(2,2),z2(2,2)
c
pro34=2D0*rin(i3,i4)
pro134=pro34+2D0*(rin(i1,i3)+rin(i1,i4))
pro234=pro34+2D0*(rin(i2,i3)+rin(i2,i4))
c
zf1=-sign*zd(i1,i4)/(pro34*pro134)
zf2=-sign*zud(i2,i3)/(pro34*pro234)
do 10 ix=1,2
do 10 iy=1,2
z1(ix,iy)=+zf1*zko(i2,iy)*(zud(i1,i3)*zkod(i1,ix)+
. zud(i4,i3)*zkod(i4,ix))
. -zf2*zkod(i1,ix)*(zd(i2,i4)*zko(i2,iy)+
. zd(i3,i4)*zko(i3,iy))
10 continue
call congat(z1,z2)
end
c
c**********************************************************************
c subroutine congat(z1,z2) takes the complex conjugate
c spinor-tensor of z1 into z2. use negaen to correct for
c an odd number of quarks with negative energy.
c**********************************************************************
subroutine congat(z1,z2)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
dimension z1(2,2),z2(2,2)
c
z2(1,1)=negaen*conjg(z1(1,1))
z2(2,1)=negaen*conjg(z1(1,2))
z2(1,2)=negaen*conjg(z1(2,1))
z2(2,2)=negaen*conjg(z1(2,2))
end
c
c**********************************************************************
c subroutine setspn(n) determines the co-variant weyl-vd waerden
c spinors and all the possible spinor-inproducts.
c special treatment for particles 6 and 7 (leptons)
c particles with negative energy have an artificial minus
c sign in their dotted spinors.
c**********************************************************************
c calling conventions:
c n = number of of vectors (n must be <= 7 )
c common /kinema/ p = array of covariant vectors (-,-,-,+) metric
c note: none of the vectors should be parallel to the y-axis.
c**********************************************************************
c output conventions:
c common /dotpr/zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
c input parameters remain unchanged.
c with:
c zud(i,j) = undotted spinor inproduct
c zd(i,j) = dotted spinor inproduct
c zko(i,1:2) = undotted co-variant spinors
c zkod(i,1:2) = dotted co-variant spinors
c negaen = number of quarks with negative
c energy. needed when complex-conjugating.
c negaen=1,-1 when that number is even,odd
c**********************************************************************
subroutine setspn(n)
implicit double precision(a-h,o-y)
implicit complex*16(z)
common /kinema/ p(4,7)
common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen
common /linpr/ rin(7,7)
c
do 10 i=1,7
if((i.le.n).or.(i.ge.6))then
zko(i,1)=(p(1,i)-(0D0,1D0)*p(3,i))
. /sqrt(dcmplx(p(4,i)-p(2,i)))
zko(i,2)=sqrt(dcmplx(p(4,i)-p(2,i)))
zkod(i,1)=conjg(zko(i,1))
zkod(i,2)=conjg(zko(i,2))
if (p(4,i).lt.0D0) then
zkod(i,1)=-zkod(i,1)
zkod(i,2)=-zkod(i,2)
endif
endif
10 continue
c
do 20 i=1,n-1
do 20 j=i,n
zud(i,j)=-zko(i,1)*zko(j,2)+zko(i,2)*zko(j,1)
zud(j,i)=-zud(i,j)
zd(i,j)=-zkod(i,1)*zkod(j,2)+zkod(i,2)*zkod(j,1)
zd(j,i)=-zd(i,j)
20 continue
c
do 30 i=1,n-1
do 30 j=i,n
rin(i,j)=rinpro(i,j)
rin(j,i)=rin(i,j)
30 continue
rin(1,6)=rinpro(1,6)
rin(1,7)=rinpro(1,7)
rin(2,6)=rinpro(2,6)
rin(2,7)=rinpro(2,7)
rin(6,7)=rinpro(6,7)
rin(7,6)=rin(6,7)
end
c
c**********************************************************************
c function rinpro(i,j) determines the minkowski inproduct.
c the parameters i and j refer to the momentum array /kinema/ p.
c both vectors are contra-variant.
c**********************************************************************
function rinpro(i,j)
implicit double precision(a-h,o-y)
common /kinema/ p(4,7)
c
rinpro=p(4,i)*p(4,j)-p(1,i)*p(1,j)-p(2,i)*p(2,j)-p(3,i)*p(3,j)
end
c
c**********************************************************************
c
c two loop strong coupling constant at scale rq
c
function alfas(mu,lam,nloops)
c alfas in terms of lambda of four flavors for one or two loops.
c matching acheived using renorm group eqn. approximately
c above mu=mb,mu=mt
implicit none
integer nloops
double precision mu,lam,mc,mb,mt
double precision b4,b4p,b5,b5p,b6,b6p,one,two
double precision alfas,atinv,abinv,atinv1,abinv1,asinv,xqc,xqb,
$ xqt,xb,xt
parameter(one=1.D0,two=2.D0)
parameter(b4=1.326291192D0,b5=1.220187897D0,b6=1.114084601D0)
parameter(b4p=0.490197225D0,b5p=0.401347248D0,b6p=0.295573466D0)
c b=(33.D0-2.D0*nf)/6.D0/pi
c bp=(153.D0-19.D0*nf)/(2.D0*pi*(33.D0-2.D0*nf))
parameter(mc=1.5D0,mb=5.0D0,mt=140.D0)
if (mu.lt.mc) then
write(6,*) 'unimplemented, mu too low',mu
alfas=0.D0
return
endif
xb=log(mb/lam)
abinv=b4*xb
if (mu.lt.mb) then
xqc=log(mu/lam)
asinv=b4*xqc
alfas=one/asinv
elseif (mu.lt.mt) then
xqb=log(mu/mb)
asinv=abinv+b5*xqb
alfas=one/asinv
else
xqt=log(mu/mt)
xt=log(mt/mb)
atinv=abinv+b5*xt
asinv=atinv+b6*xqt
alfas=one/asinv
endif
if (nloops.eq.1) return
abinv1=abinv/(one-b4p*log(two*xb)/abinv)
if (mu.lt.mb) then
asinv=asinv/(one-b4p*log(two*xqc)/asinv)
alfas=one/asinv
elseif (mu.lt.mt) then
asinv=abinv1+b5*xqb+b5p*log((asinv+b5p)/(abinv+b5p))
alfas=one/asinv
else
atinv1=abinv1+b5*xt+b5p*log((b5p+atinv)/(b5p+abinv))
asinv=atinv1+b6*xqt+b6p*log((b6p+asinv)/(atinv+b6p))
alfas=one/asinv
endif
return
end
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C THE ROUTINES FROM HERE ON ARE MORE GENERAL UTILITIES
C - AN APPROXIMATION TO THE DILOGARITHM
C - A RANDOM NUMBER GENERATOR
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
FUNCTION DILOG(X)
IMPLICIT NONE
c$$$C---RETURN THE DILOGARITHM (MAXIMUM ERROR = 0.0058, MAX FRAC ERROR = 1%)
c$$$ DOUBLE PRECI\ION DILOG,X,PISQO6,OX,XX,L2
c$$$ DATA PISQO6,L2/2*0/
c$$$ IF (PISQO6.EQ.0) PISQO6=(ATAN(1D0)*4)**2/6
c$$$ IF (L2.EQ.0) L2=LOG(2D0)
c$$$ IF (X.LT.-1) THEN
c$$$ XX=1/X
c$$$ ELSE
c$$$ XX=X
c$$$ ENDIF
c$$$ IF (XX.LT.-0.5) THEN
c$$$ OX=1+XX
c$$$ DILOG=-PISQO6/2-L2*LOG(-XX)
c$$$ $ -OX**2/4-5*OX**3/24-OX**4/6-131*OX**5/960
c$$$ ELSEIF (XX.LT.0.5) THEN
c$$$ DILOG=XX+XX**2/4+XX**3/9
c$$$ ELSEIF (XX.LT.1) THEN
c$$$ OX=1-XX
c$$$ DILOG=PISQO6-LOG(OX)*LOG(XX)-OX-OX**2/4-OX**3/9
c$$$ ELSEIF (XX.EQ.1) THEN
c$$$ DILOG=PISQO6
c$$$ ELSE
c$$$ WRITE (*,*) 'DILOG CALLED FOR X=',X
c$$$ DILOG=0
c$$$ ENDIF
c$$$ IF (X.LT.-1) DILOG=-DILOG-PISQO6-LOG(-X)**2/2
double precision dilog,rsp,x
dilog=rsp(x)
END
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
SUBROUTINE RANGEN(N,R)
IMPLICIT NONE
C---RANDOM NUMBER GENERATOR
C USES METHOD OF l'Ecuyer, (VIA F.JAMES, COMP PHYS COMM 60(1990)329)
C RETURNS A VECTOR OF N RANDOM VALUES
C IF (N.EQ.0) THE FIRST TWO VALUES IN R SET THE SEEDS
C IF (N.LT.0) PRINT THE CURRENT VALUES OF THE SEEDS
DOUBLE PRECISION R(*)
INTEGER N,I,ISEED(2),K,IZ
DATA ISEED/123450,678900/
IF (N.LT.0) WRITE (*,'(I10,A,I10,I11)') -N-1,', ISEED=',ISEED
IF (N.GT.0) THEN
DO I=1,N
K=ISEED(1)/53668
ISEED(1)=40014*(ISEED(1)-K*53668)-K*12211
IF (ISEED(1).LT.0) ISEED(1)=ISEED(1)+2147483563
K=ISEED(2)/52774
ISEED(2)=40692*(ISEED(2)-K*52774)-K*3791
IF (ISEED(2).LT.0) ISEED(2)=ISEED(2)+2147483399
IZ=ISEED(1)-ISEED(2)
IF (IZ.LT.1) IZ=IZ+2147483562
R(I)=DBLE(IZ)*4.656613D-10
ENDDO
ELSEIF (N.EQ.0) THEN
ISEED(1)=NINT(R(1))
ISEED(2)=NINT(R(2))
ENDIF
END
c function zspv(smu,snu) calculates a smu*snu inproduct in
c spinor-tensor language. both smu and snu are co-variant tensors.
c the factor 2 is due to spinor transformation algebra.
c**********************************************************************
function zspv(smu,snu)
complex*16 zspv,smu(1:2,1:2),snu(1:2,1:2)
zspv=2D0*(+smu(1,1)*snu(2,2)
. +smu(2,2)*snu(1,1)
. -smu(1,2)*snu(2,1)
. -smu(2,1)*snu(1,2) )
end
c
c
c spence function taking only real arguments
c
function rsp(x)
implicit double precision(a-h,o-z)
common/spint/a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,zeta2
x2=x*x
if(x.gt.1.D0)then
write(*,*)' argument greater than 1 passed to spence function'
rsp=0.D0
return
endif
if(x2.gt.1.D0.and.x.gt.0.5D0)then
y=(x-1.D0)/x
z=-log(1.D0-y)
z2=z*z
rsp=z*(1.D0+a1*z*(1.D0+a2*z*(1.D0+a3*z2*(1.D0+a4*z2*
1 (1.D0+a5*z2*(1.D0+a6*z2*(1.D0+a7*z2*(1.D0+a8*z2*(1.D0+a9*z2*
2 (1.D0+a10*z2))))))))))
3 +zeta2-log(x)*log(1.D0-x)+0.5D0*log(x)**2
return
elseif(x2.gt.1.D0.and.x.le.0.5D0)then
y=1.D0/x
z=-log(1.D0-y)
z2=z*z
rsp=-z*(1.D0+a1*z*(1.D0+a2*z*(1.D0+a3*z2*(1.D0+a4*z2*
1 (1.D0+a5*z2*(1.D0+a6*z2*(1.D0+a7*z2*(1.D0+a8*z2*(1.D0+a9*z2*
2 (1.D0+a10*z2))))))))))
3 -zeta2-0.5D0*log(-x)**2
return
elseif(x2.eq.1.D0)then
rsp=zeta2
return
elseif(x2.le.1.D0.and.x.gt.0.5D0)then
y=1.D0-x
z=-log(1.D0-y)
z2=z*z
rsp=-z*(1.D0+a1*z*(1.D0+a2*z*(1.D0+a3*z2*(1.D0+a4*z2*
1 (1.D0+a5*z2*(1.D0+a6*z2*(1.D0+a7*z2*(1.D0+a8*z2*(1.D0+a9*z2*
2 (1.D0+a10*z2))))))))))
3 +zeta2-log(x)*log(1.D0-x)
return
elseif(x2.le.1.D0.and.x.le.0.5D0)then
y=x
z=-log(1.D0-y)
z2=z*z
rsp=z*(1.D0+a1*z*(1.D0+a2*z*(1.D0+a3*z2*(1.D0+a4*z2*
1 (1.D0+a5*z2*(1.D0+a6*z2*(1.D0+a7*z2*(1.D0+a8*z2*(1.D0+a9*z2*
2 (1.D0+a10*z2))))))))))
return
else
write(*,*)' illegal x value in spence function'
rsp=0.D0
endif
return
end
c spence function
c
block data splint
implicit double precision (a-h,o-z)
common/spint/a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,zeta2
data a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,zeta2/
1 -0.250000000000000D0,
2 -0.111111111111111D0,
3 -0.010000000000000D0,
4 -0.017006802721088D0,
5 -0.019444444444444D0,
6 -0.020661157024793D0,
7 -0.021417300648069D0,
8 -0.021948866377231D0,
9 -0.022349233811171D0,
1 -0.022663689135191D0,
2 1.644934066848226D0/
end block data splint
c
c***********************************************************************
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
Reference in New Issue View Git Blame Copy Permalink