PROGRAM MAIN IMPLICIT NONE INTEGER NFL,NEV DOUBLE PRECISION EM COMMON /PARAMS/ NFL,NEV,EM EXTERNAL DEMO CALL READCARD CALL DEMOIN CALL EVENT2(NEV,EM,NFL,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 VERSION 0.4 - 5th May 2026 - separated hist ouput C - added runcard input C - allowed different bin sizes C - seeds can now be externally set C - added log-binning capability 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 CALL PRINTHEADER 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 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 PRINTHEADER IMPLICIT NONE 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)') ' Edited by Giorgio Chiurato, April 2026' WRITE (*,'(A/)') ' Version 0.4, May 2026' END C----------------------------------------------------------------------- SUBROUTINE READCARD implicit none integer :: stat,n,nl,i,iline,ichar character(3) :: flag(1:5) character(40) :: val(1:5) character(20) :: fname character(5) :: RID character(72) :: line character(12) :: keys(20),settings(20) logical fexists INTEGER NFL,NEV DOUBLE PRECISION EM INTEGER METYPE,NF DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP, $ CQ,FPAL2,ONF INTEGER CB1,CB2,DB,TB,YB,EB DOUBLE PRECISION CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH DOUBLE PRECISION CBS1,CBS2,DBS,TBS,YBS,EBS INTEGER LHI DOUBLE PRECISION ISEEDIN(2) CHARACTER*40 PATH DATA fname/'card.input'/ DATA RID /'00000'/ DATA ISEEDIN/123450,678900/ COMMON /PARAMS/ NFL,NEV,EM COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP, $ CQ,FPAL2,ONF,NF,METYPE COMMON /GBINS/CB1,CB2,DB,TB,YB,EB COMMON /GBINS/CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH COMMON /GBINS/CBS1,CBS2,DBS,TBS,YBS,EBS COMMON /GBINS/LHI COMMON /ID/ RID common /runcard/keys,settings PATH = './results/' n = iargc() nl = -1 if (n.ge.2) then call getarg(1,flag(1)) call getarg(2,val(1)) nl = 1 endif if (n.ge.4)then call getarg(3,flag(2)) call getarg(4,val(2)) nl = 2 endif if (n.ge.6)then call getarg(5,flag(3)) call getarg(6,val(3)) nl = 3 endif if (n.ge.8)then call getarg(7,flag(4)) call getarg(8,val(4)) nl = 4 endif if (n.ge.10)then call getarg(9,flag(5)) call getarg(10,val(5)) nl = 5 endif C TODO: need to add a proper info page and maybe change the input parameter handling. do i=1,nl if (flag(i).eq.'-i')then fname = val(i) elseif (flag(i).eq.'-o')then PATH = val(i) elseif (flag(i).eq.'-sa') then read(val(i), *) ISEEDIN(1) elseif (flag(i).eq.'-sb') then read(val(i), *) ISEEDIN(2) elseif (flag(i).eq.'-c') then RID = val(i) else CALL PRINTHEADER print *, 'Option ', flag(i), $ ' not recognized.' STOP "Exiting program..." endif enddo call RANGEN(0, ISEEDIN) call CREATEOUT(PATH) inquire(file=fname, exist=fexists) if (.not.fexists)then CALL PRINTHEADER stop 'run card ' // trim(fname) // ' not found' endif open(9,file=fname) iline = 1 do call getline(9, line, stat) if (line.ne.'')then ichar = index(line, "=") if (ichar.gt.0)then keys(iline) = line(:ichar-1) settings(iline) = line(ichar+1:) endif iline = iline+1 endif if (stat.lt.0) exit end do call readparm('EM ', EM, real(91.187, 8)) call readparm('CUTOFF ', CUTOFF, 1D-8) call readparm('CUTUP ', CUTUP , 1D0) call readmode('NF ', NFL, 5) call readmode('NEV ', NEV, 5 000 000) call readmode('CB1 ', CB1, 50) call readmode('CB2 ', CB2, 37) call readmode('DB ', DB, 50) call readmode('TB ', TB, 50) call readmode('YB ', YB, 40) call readmode('EB ', EB, 50) call readparm('CL1 ', CL1 , 0.00D0) call readparm('CH1 ', CH1 , 1.00D0) call readparm('CL2 ', CL2 , 0.00D0) call readparm('CH2 ', CH2 , 0.74D0) call readparm('DL ', DL , 0.00D0) call readparm('DH ', DH , 1.00D0) call readparm('TL ', TL , 0.00D0) call readparm('TH ', TH , 0.50D0) call readparm('YL ', YL , 0.00D0) call readparm('YH ', YH , 0.40D0) call readparm('EL ', EL , -1.0D0) call readparm('EH ', EH , 1.00D0) call readmode('LH ', LHI , 0) C Close file. close(9) END C----------------------------------------------------------------------- SUBROUTINE DEMOIN IMPLICIT NONE CHARACTER(5) RID INTEGER I INTEGER CB1,CB2,DB,TB,YB,EB DOUBLE PRECISION CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH DOUBLE PRECISION CBS1,CBS2,DBS,TBS,YBS,EBS LOGICAL LH DOUBLE PRECISION CSUM,CSQR,BSUM,BSQR COMMON /DEMCOM/CSUM(8),CSQR(8),BSUM(2,6),BSQR(2,6) COMMON /GBINS/CB1,CB2,DB,TB,YB,EB COMMON /GBINS/CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH COMMON /GBINS/CBS1,CBS2,DBS,TBS,YBS,EBS COMMON /GBINS/LH COMMON /ID/ RID CBS1=(CH1-CL1)/CB1 CBS2=(CH2-CL2)/CB2 DBS=(DH-DL)/DB TBS=(TH-TL)/TB YBS=(YH-YL)/YB EBS=(EH-EL)/EB 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_'//RID,CB1,CL1,CH1,LH) CALL GBOOK1(101,'C-PARAMETER_'//RID,CB1,CL1,CH1,LH) CALL GBOOK1(201,'C-PARAMETER_'//RID,CB1,CL1,CH1,LH) CALL GBOOK1( 11,'C-PARAMETER_'//RID,CB2,CL2,CH2,LH) CALL GBOOK1(111,'C-PARAMETER_'//RID,CB2,CL2,CH2,LH) CALL GBOOK1(211,'C-PARAMETER_'//RID,CB2,CL2,CH2,LH) CALL GBOOK1( 2,'D-PARAMETER_'//RID,DB,DL,DH,LH) CALL GBOOK1(102,'D-PARAMETER_'//RID,DB,DL,DH,LH) CALL GBOOK1(202,'D-PARAMETER_'//RID,DB,DL,DH,LH) CALL GBOOK1( 12,'D-PARAMETER_'//RID,DB,DL,DH,LH) CALL GBOOK1(112,'D-PARAMETER_'//RID,DB,DL,DH,LH) CALL GBOOK1(212,'D-PARAMETER_'//RID,DB,DL,DH,LH) CALL GBOOK1( 3,'THRUST_'//RID,TB,TL,TH,LH) CALL GBOOK1(103,'THRUST_'//RID,TB,TL,TH,LH) CALL GBOOK1(203,'THRUST_'//RID,TB,TL,TH,LH) CALL GBOOK1( 13,'THRUST_'//RID,TB,TL,TH,LH) CALL GBOOK1(113,'THRUST_'//RID,TB,TL,TH,LH) CALL GBOOK1(213,'THRUST_'//RID,TB,TL,TH,LH) CALL GBOOK1( 4,'y3(JADE,P)_'//RID,YB,YL,YH,LH) CALL GBOOK1(104,'y3(JADE,P)_'//RID,YB,YL,YH,LH) CALL GBOOK1(204,'y3(JADE,P)_'//RID,YB,YL,YH,LH) CALL GBOOK1( 14,'y3(JADE,P)_'//RID,YB,YL,YH,LH) CALL GBOOK1(114,'y3(JADE,P)_'//RID,YB,YL,YH,LH) CALL GBOOK1(214,'y3(JADE,P)_'//RID,YB,YL,YH,LH) CALL GBOOK1( 5,'EEC_'//RID,EB,EL,EH,LH) CALL GBOOK1(105,'EEC_'//RID,EB,EL,EH,LH) CALL GBOOK1(205,'EEC_'//RID,EB,EL,EH,LH) CALL GBOOK1( 15,'EEC_'//RID,EB,EL,EH,LH) CALL GBOOK1(115,'EEC_'//RID,EB,EL,EH,LH) CALL GBOOK1(215,'EEC_'//RID,EB,EL,EH,LH) 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,TTL,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 INTEGER CB1,CB2,DB,TB,YB,EB DOUBLE PRECISION CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH DOUBLE PRECISION CBS1,CBS2,DBS,TBS,YBS,EBS LOGICAL LH COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB COMMON /DEMCOM/CSUM,CSQR,BSUM,BSQR COMMON /GBINS/CB1,CB2,DB,TB,YB,EB COMMON /GBINS/CL1,CL2,DL,TL,YL,EL,CH1,CH2,DH,TH,YH,EH COMMON /GBINS/CBS1,CBS2,DBS,TBS,YBS,EBS COMMON /GBINS/LH 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 GFILLSC1(201+ORD,C,WEIGHT) IF (ORD.EQ.0) CS(3)=CS(3)+WEIGHT 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 GFILLSC1(202+ORD,D,WEIGHT) IF (ORD.EQ.0) CS(4)=CS(4)+WEIGHT 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 TTL=PC(1,I)**2+PC(2,I)**2+PC(3,I)**2 IF (TTL.GT.T) T=TTL ENDDO ENDDO ENDDO T=SQRT(T)*ORS ELSE T=MAX(E(1),E(2),E(3))*ORS*2 ENDIF TAU = 1 - T CALL GFILLSC1(203+ORD,TAU,WEIGHT) IF (ORD.EQ.0) CS(5)=CS(5)+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 GFILLSC1(204+ORD,Y3,WEIGHT) IF (ORD.EQ.0) CS(6)=CS(6)+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 GFILLSC1(205+ORD,COSANG,BEEC) 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-----------------------------------------------------------------------