commit 82f7abd4ecb107a74225494d1b860f591e208bb7 Author: Giorgio Chiurato Date: Tue Mar 31 18:04:12 2026 +0200 Initial commit. diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c9a4638 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +bin/ + +gtopdraw.top diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..cc878a2 --- /dev/null +++ b/Makefile @@ -0,0 +1,23 @@ +# Compiler and flags +FC = gfortran +FFLAGS = -O2 + +# Output folder +OUT = bin + +# Target executable +TARGET = event2 + +# Source files +SRCS = event2_03.f gbook.f + +# Default rule +all: $(TARGET) + +# Build rule +$(TARGET): $(SRCS) + $(FC) $(FFLAGS) -o $(OUT)/$(TARGET) $(SRCS) + +# Clean rule +clean: + rm -f $(TARGET) *.o diff --git a/event2_03.f b/event2_03.f new file mode 100644 index 0000000..20f1eaa --- /dev/null +++ b/event2_03.f @@ -0,0 +1,3618 @@ + IMPLICIT NONE + INTEGER NF,NEV + DOUBLE PRECISION EM + EXTERNAL DEMO + NF=5 + EM=91.1876 + NEV=5 000 000 + CALL DEMOIN + CALL EVENT2(NEV,EM,NF,DEMO) + CALL DEMOUT(NEV) + END +C----------------------------------------------------------------------- + SUBROUTINE EVENT2(NEV,EM,NFL,USER) + IMPLICIT NONE +C---CALCULATE E+E- EVENT FEATURES TO NEXT-TO-LEADING ORDER +C ACCORDING TO THE METHOD OF CATANI AND SEYMOUR PLB378 (1996) 287 +C +C VERSION 0.0 - 18th January 1996 +C VERSION 0.1 - 26th February 1996 - bug fix in event orientation +C VERSION 0.2 - 15th March 1996 - removed need for Lorentz boosts +C VERSION 0.3 - 10th November 1997 - improved numerical convergence +C and safety against arithmetic errors +C - added CUTUP variable +C +C - NEV IS THE NUMBER OF EVENTS TO GENERATE +C - EM IS THE TOTAL CENTRE-OF-MASS ENERGY +C - NFL IS THE NUMBER OF FLAVOURS +C - USER IS THE NAME OF ROUTINE TO ANALYSE THE EVENTS +C IT TAKES FIVE ARGUMENTS: +C SUBROUTINE USER(N,NA,ITYPE,P,WEIGHT) +C N = THE NUMBER OF PARTONS +C NA = THE ORDER IN ALPHA +C ITYPE = THE TYPE: 0=TREE LEVEL, 1=SUBTRACTION, 2=FINITE VIRTUAL +C P(4,7) = THEIR MOMENTA (P(1-4,I) IS 4-MOMENTUM OF PARTICLE I, +C IN ORDER Q,QBAR, Q,QBAR,G, Q,QBAR,G,G OR Q',Q,Q'BAR,QBAR +C P(1-4,5) IS TOTAL 4-MOM OF THE N PARTICLES, +C P(1-4,6/7) IS THE 4-MOM OF THE ELECTRON/POSITRON) +C WEIGHT = THE WEIGHT IN THE INTEGRAL +C IT IS ALSO CALLED AT THE END OF EACH FULL EVENT WITH N=0 +C A SIMPLE EXAMPLE (CALLED DEMO) IS GIVEN BELOW +C +C IT IS IMPORTANT TO THINK ABOUT THE EFFICIENCY OF THE USER ROUTINE, +C AS IT IS CALLED 17 TIMES PER EVENT +C (THIS IS 1 4-PARTON CALL, 12 3-PARTON CALLS AND 4 2-PARTON CALLS) +C + INTEGER NEV,NFL,NPERM3,NPERM4,NU,ND,I,J,K + LOGICAL KEEP + PARAMETER (NPERM3=2,NPERM4=10) + DOUBLE PRECISION EM,WTWO,WTHR,WFOR,JTWO,JTHR,JFOR,JTMP, + $ MTWO,MTHR,MFOR,VTWO,VTHR,STHR(NPERM3),SFOR(NPERM4),WEIGHT, + $ ZERO,ONE,P(4,7),Q(4,7,NPERM4),NRM,QD,QU,AEM,SUBS(6,NPERM4) + PARAMETER (ZERO=0,ONE=1) + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + INTEGER NPOW1,NPOW2 + DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ + COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2 + DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + EXTERNAL USER +C---PRINT OPENING MESSAGE + WRITE (*,'(/2A)') ' This is EVENT2, a program for calculating', + $ ' jet quantities' + WRITE (*,'(2A)') ' in e+e- annihilation to next-to-leading', + $ ' order in alpha_s' + WRITE (*,'(A)') ' If you use this program, please reference:' + WRITE (*,'(2A)') ' S.Catani & M.H.Seymour,', + $ ' Phys. Lett. B378 (1996) 287.' + WRITE (*,'(/A)') ' Written by Mike Seymour, January 1996' + WRITE (*,'(A/)') ' Version 0.3, November 1997' + WRITE (*,'(A,I10)')' NEV=',NEV +C---INITIALIZE CONSTANTS + PI=ATAN(ONE)*4 + PISQ=PI**2 + HF=0.5 + CA=3 + CF=(CA**2-1)/(2*CA) + TR=HF + NF=NFL + ONF=1 + IF (NF.NE.0) ONF=ONE/NF + ND=INT((NF+1)/2) + NU=INT(NF/2) + QD=-ONE/3 + QU=QD+1 + CQ=ND*QD**2+NU*QU**2 + AEM=ONE/128 + FPAL2=(4*PI*AEM)**2 + CUTOFF=1D-8 + CUTUP=1 +C---METYPE IS 0 FOR ERT, 1 FOR LEIDEN AND 2 FOR GIELE-GLOVER +C (CROSS-SECTION NORMALIZATION IS DIFFERENT IN THE DIFFERENT CASES. +C FOR ERT, 2-PARTON CROSS-SECTION IS 1, FOR THE OTHERS IT IS THE +C BORN CROSS-SECTION FOR E+E- --> QQBAR IN NANOBARNS.) + METYPE=1 + IF (METYPE.EQ.0) THEN + NRM=16*PI*EM**2 + ELSE + NRM=389385 + ENDIF +C---PARAMETERS RELATED TO IMPORTANCE SAMPLING + NPOW1=2 + NPOW2=2 + XPOW1=1-ONE/NPOW1 + XPOW2=1-ONE/NPOW2 + COSP=1 + COSM=1 + COSZ=0 +C---START MAIN LOOP + DO I=1,NEV + IF (MOD(I,100 000).EQ.1) CALL RANGEN(-I,SFOR) +C---GENERATE A TWO-PARTON STATE + CALL GENTWO(EM,P,WTWO,*1000) +C---CALCULATE THE JACOBIAN FACTOR + CALL JACTWO(P,JTWO) +C---EVALUATE THE TWO-PARTON TREE-LEVEL MATRIX ELEMENT + CALL MATTWO(P,MTWO) +C---GIVE IT TO THE USER + WEIGHT=NRM*MTWO*WTWO/JTWO/NEV + CALL USER(2,0,0,P,WEIGHT) + IF (CUTUP.GE.1) THEN +C---EVALUATE THE TWO-PARTON ONE-LOOP MATRIX ELEMENT + CALL VIRTWO(P,VTWO) +C---GIVE IT TO THE USER + WEIGHT=NRM*VTWO*WTWO/JTWO/NEV + CALL USER(2,1,2,P,WEIGHT) + ENDIF +C---GENERATE A THREE-PARTON STATE + KEEP=.FALSE. + CALL GENTHR(P,WTHR,KEEP,*1000) +C---CALCULATE THE JACOBIAN FACTOR AND SUBTRACTION CONFIGURATIONS + JTHR=0 + DO J=1,NPERM3 + CALL SUBTHR(J,P,Q(1,1,J),STHR(J),JTMP,*1000) + JTHR=JTHR+JTMP + ENDDO + IF (KEEP) THEN +C---EVALUATE THE THREE-PARTON TREE-LEVEL MATRIX ELEMENT + CALL MATTHR(P,MTHR) +C---GIVE IT TO THE USER + WEIGHT=NRM*MTHR*WTWO*WTHR/JTHR/NEV + CALL USER(3,1,0,P,WEIGHT) +C---GIVE THE SUBTRACTION CONFIGURATIONS TO THE USER + DO J=1,NPERM3 + WEIGHT=-NRM*STHR(J)*WTWO*WTHR/JTHR/NEV + CALL USER(2,1,1,Q(1,1,J),WEIGHT) + ENDDO +C---EVALUATE THE THREE-PARTON ONE-LOOP MATRIX ELEMENT + CALL VIRTHR(P,VTHR) +C---GIVE IT TO THE USER + WEIGHT=NRM*VTHR*WTWO*WTHR/JTHR/NEV + CALL USER(3,2,2,P,WEIGHT) + ENDIF +C---GENERATE A FOUR-PARTON STATE + CALL GENFOR(P,WFOR,KEEP,*1000) +C---CALCULATE THE JACOBIAN FACTOR AND SUBTRACTION CONFIGURATIONS + JFOR=0 + DO J=1,NPERM4 + CALL SUBFOR(J,P,Q(1,1,J),SFOR(J),JTMP,KEEP,*1000) + JFOR=JFOR+JTMP + SUBS(1,J)=CFSUB + SUBS(2,J)=CASUB + SUBS(3,J)=TFSUB + SUBS(4,J)=GGSUB + SUBS(5,J)=QQSUB + SUBS(6,J)=QPSUB + ENDDO + IF (KEEP) THEN +C---EVALUATE THE FOUR-PARTON TREE-LEVEL MATRIX ELEMENT + CALL MATFOR(P,MFOR) +C---GIVE IT TO THE USER + WEIGHT=NRM*MFOR*WTWO*WTHR*WFOR/JFOR/NEV + CALL USER(4,2,0,P,WEIGHT) +C---GIVE THE SUBTRACTION CONFIGURATIONS TO THE USER + DO J=1,NPERM4 + CFSUB=SUBS(1,J) + CASUB=SUBS(2,J) + TFSUB=SUBS(3,J) + GGSUB=SUBS(4,J) + QQSUB=SUBS(5,J) + QPSUB=SUBS(6,J) + WEIGHT=-NRM*SFOR(J)*WTWO*WTHR*WFOR/JFOR/NEV + CALL USER(3,2,1,Q(1,1,J),WEIGHT) + ENDDO + ENDIF +C---TELL THE USER THAT THE EVENT IS COMPLETE + 1000 CALL USER(0,0,0,P,ZERO) + ENDDO + END +C----------------------------------------------------------------------- + SUBROUTINE DEMOIN + IMPLICIT NONE + INTEGER I + DOUBLE PRECISION CSUM,CSQR,BSUM,BSQR + COMMON /DEMCOM/CSUM(8),CSQR(8),BSUM(2,6),BSQR(2,6) + DO I=1,8 + CSUM(I)=0 + CSQR(I)=0 + IF (I.LE.6) THEN + BSUM(1,I)=0 + BSUM(2,I)=0 + BSQR(1,I)=0 + BSQR(2,I)=0 + ENDIF + ENDDO + CALL GBOOK1( 1,'C-PARAMETER',50,0D0,1D0) + CALL GBOOK1(101,'C-PARAMETER',50,0D0,1D0) + CALL GBOOK1(201,'C-PARAMETER',50,0D0,1D0) + CALL GBOOK1( 11,'C-PARAMETER',37,0D0,0.74D0) + CALL GBOOK1(111,'C-PARAMETER',37,0D0,0.74D0) + CALL GBOOK1(211,'C-PARAMETER',37,0D0,0.74D0) + CALL GBOOK1( 2,'D-PARAMETER',50,0D0,1D0) + CALL GBOOK1(102,'D-PARAMETER',50,0D0,1D0) + CALL GBOOK1(202,'D-PARAMETER',50,0D0,1D0) + CALL GBOOK1( 12,'D-PARAMETER',50,0D0,1D0) + CALL GBOOK1(112,'D-PARAMETER',50,0D0,1D0) + CALL GBOOK1(212,'D-PARAMETER',50,0D0,1D0) + CALL GBOOK1( 3,'THRUST',50,0.5D0,1D0) + CALL GBOOK1(103,'THRUST',50,0.5D0,1D0) + CALL GBOOK1(203,'THRUST',50,0.5D0,1D0) + CALL GBOOK1( 13,'THRUST',50,0.5D0,1D0) + CALL GBOOK1(113,'THRUST',50,0.5D0,1D0) + CALL GBOOK1(213,'THRUST',50,0.5D0,1D0) + CALL GBOOK1( 4,'y3(JADE,P)',40,0D0,0.4D0) + CALL GBOOK1(104,'y3(JADE,P)',40,0D0,0.4D0) + CALL GBOOK1(204,'y3(JADE,P)',40,0D0,0.4D0) + CALL GBOOK1( 14,'y3(JADE,P)',40,0D0,0.4D0) + CALL GBOOK1(114,'y3(JADE,P)',40,0D0,0.4D0) + CALL GBOOK1(214,'y3(JADE,P)',40,0D0,0.4D0) + CALL GBOOK1( 5,'EEC',50,-1D0,1D0) + CALL GBOOK1(105,'EEC',50,-1D0,1D0) + CALL GBOOK1(205,'EEC',50,-1D0,1D0) + CALL GBOOK1( 15,'EEC',50,-1D0,1D0) + CALL GBOOK1(115,'EEC',50,-1D0,1D0) + CALL GBOOK1(215,'EEC',50,-1D0,1D0) + END +C----------------------------------------------------------------------- + SUBROUTINE DEMOUT(NEV) + IMPLICIT NONE + INTEGER I,J,NEV + DOUBLE PRECISION CSUM,CSQR,BSUM,BSQR,SF,SQR + COMMON /DEMCOM/CSUM(8),CSQR(8),BSUM(2,6),BSQR(2,6) + CHARACTER*10 NAME(8) + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + SQR(SF)=SIGN(SQRT(ABS(SF)),SF) + DATA NAME + $/' sigma_0',' sigma_1',' C-PARAM',' D-PARAM', + $ ' 1-THRUST','y3(JADE,p)','EEC(90deg)','EEC(total)'/ + SF=CSUM(1) + IF (SF.NE.0) THEN + DO I=1,5 + CALL GSCALE( I,1/SF) + CALL GSCALE(100+I,1/SF**2) + CALL GSCALE( 10+I,1/SF) + CALL GSCALE(110+I,1/SF**2) + ENDDO + DO I=1,7 + CSUM(I+1)=CSUM(I+1)/SF + CSQR(I+1)=CSQR(I+1)/SF**2 + ENDDO + DO I=1,6 + BSUM(1,I)=BSUM(1,I)/SF + BSQR(1,I)=BSQR(1,I)/SF**2 + BSUM(2,I)=BSUM(2,I)/SF + BSQR(2,I)=BSQR(2,I)/SF**2 + ENDDO + ENDIF + DO I=1,5 + CALL GOPERA( I,'*', I,200+I,1D0,1D0) + CALL GOPERA(100+I,'-',200+I,100+I,1D0,1D0/NEV) + CALL GOPERA(100+I,'S', 0,100+I,1D0,0D0) + CALL GOPERA( 10+I,'*', 10+I,210+I,1D0,1D0) + CALL GOPERA(110+I,'-',210+I,110+I,1D0,1D0/NEV) + CALL GOPERA(110+I,'S', 0,110+I,1D0,0D0) + CALL GTOPER( I,1,1,0,100+I) + CALL GTOPER(10+I,0,0,0,110+I) + ENDDO + WRITE (*,'(8(1X,2(A,F15.5)/))') + $ (NAME(I),CSUM(I),' +-',SQR(CSQR(I)-CSUM(I)**2/NEV),I=1,8) + WRITE (*,'(2(/A))') ' Moments of EEC function:',' m n' + WRITE (*,'(2I4,F15.5,A,F15.5)') + $ ((I-1,J-1,BSUM(I,J),' +-',SQR(BSQR(I,J)-BSUM(I,J)**2/NEV) + $ ,J=1,6),I=1,2) + END +C----------------------------------------------------------------------- + SUBROUTINE DEMO(N,NA,ITYPE,P,WEIGHT) + IMPLICIT NONE + INTEGER N,NA,ITYPE,I,J,K,L,ORD + DOUBLE PRECISION P(4,7),WEIGHT,DOT,C,D,OS,ORS,E(7),OE(7),COSANG, + $ COSLIM,BEEC,T,TL,TA(3),PT(3),PC(3,4),Y3,Y,Q(4,7), + $ CSUM(8),CSQR(8),CS(8),BSUM(2,6),BSQR(2,6),BS(2,6) + DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /DEMCOM/CSUM,CSQR,BSUM,BSQR + DATA CS,BS/20*0/ +C---FILL THE HISTOGRAMS AT THE END OF EACH EVENT +C (THIS IS SO THAT PAIRS OF LARGE POSITIVE AND NEGATIVE WEIGHTS +C GET CANCELLED BEFORE CALCULATING THE VARIANCE) + IF (N.EQ.0) THEN + DO I=0,10,10 + DO J=1,5 + CALL GOPERA(200+I+J,'+', I+J, I+J,1D0,1D0) + CALL GOPERA(200+I+J,'*',200+I+J,200+I+J,1D0,1D0) + CALL GOPERA(200+I+J,'+',100+I+J,100+I+J,1D0,1D0) + CALL GRESET(200+I+J) + ENDDO + ENDDO + DO I=1,8 + CSUM(I)=CSUM(I)+CS(I) + CSQR(I)=CSQR(I)+CS(I)**2 + CS(I)=0 + IF (I.LE.6) THEN + BSUM(1,I)=BSUM(1,I)+BS(1,I) + BSUM(2,I)=BSUM(2,I)+BS(2,I) + BSQR(1,I)=BSQR(1,I)+BS(1,I)**2 + BSQR(2,I)=BSQR(2,I)+BS(2,I)**2 + BS(1,I)=0 + BS(2,I)=0 + ENDIF + ENDDO + RETURN + ENDIF +C---CALCULATE THE ENERGY OF EACH PARTICLE + OS=1/DOT(P,5,5) + ORS=SQRT(OS) + DO I=1,7 + IF (I.LE.N.OR.I.GE.5) THEN + E(I)=DOT(P,I,5)*ORS + OE(I)=1/E(I) + ENDIF + ENDDO +C---CALCULATE THE TOTAL CROSS-SECTION + IF (NA.EQ.0) CS(1)=CS(1)+WEIGHT + IF (NA.EQ.1) CS(2)=CS(2)+WEIGHT + IF (N.LT.3) RETURN + ORD=20-10*NA +c$$$C---ONLY INCLUDE ONE OF THE COLOUR PIECES IN THE ALPHA_S**2 TERMS +c$$$ IF (NA.EQ.2) WEIGHT=WEIGHT*CASUB +C---CALCULATE THE C-PARAMETER + C=3 + DO I=2,N + DO J=1,I-1 + C=C-3*DOT(P,I,J)**2*OS/(E(I)*E(J)) + ENDDO + ENDDO + CALL GFILL1(201+ORD,C,C*WEIGHT*50) + IF (ORD.EQ.0) CS(3)=CS(3)+WEIGHT*C +C---CALCULATE THE D-PARAMETER + IF (N.EQ.4) THEN + D=27*(2*(DOT(P,1,2)*DOT(P,1,3)*DOT(P,2,4)*DOT(P,3,4)+ + $ DOT(P,1,2)*DOT(P,1,4)*DOT(P,2,3)*DOT(P,3,4)+ + $ DOT(P,1,3)*DOT(P,1,4)*DOT(P,2,3)*DOT(P,2,4))-( + $ DOT(P,1,2)**2*DOT(P,3,4)**2+DOT(P,1,3)**2*DOT(P,2,4)**2+ + $ DOT(P,1,4)**2*DOT(P,2,3)**2))/( + $ (DOT(P,1,2)+DOT(P,1,3)+DOT(P,1,4))* + $ (DOT(P,1,2)+DOT(P,2,3)+DOT(P,2,4))* + $ (DOT(P,1,3)+DOT(P,2,3)+DOT(P,3,4))* + $ (DOT(P,1,4)+DOT(P,2,4)+DOT(P,3,4))) + ELSE + D=0 + ENDIF + CALL GFILL1(202+ORD,D,D*WEIGHT*50) + IF (ORD.EQ.0) CS(4)=CS(4)+WEIGHT*D +C---CALCULATE THE THRUST + IF (N.EQ.4) THEN + T=0 + DO K=2,N + DO J=1,K-1 + TA(1)=P(2,J)*P(3,K)-P(3,J)*P(2,K) + TA(2)=P(3,J)*P(1,K)-P(1,J)*P(3,K) + TA(3)=P(1,J)*P(2,K)-P(2,J)*P(1,K) + PT(1)=0 + PT(2)=0 + PT(3)=0 + DO I=1,N + IF (I.NE.J.AND.I.NE.K) THEN + IF (P(1,I)*TA(1)+P(2,I)*TA(2)+P(3,I)*TA(3).GE.0.) THEN + PT(1)=PT(1)+P(1,I) + PT(2)=PT(2)+P(2,I) + PT(3)=PT(3)+P(3,I) + ELSE + PT(1)=PT(1)-P(1,I) + PT(2)=PT(2)-P(2,I) + PT(3)=PT(3)-P(3,I) + ENDIF + ENDIF + ENDDO + DO I=1,3 + PC(I,1)=PT(I)+P(I,J)+P(I,K) + PC(I,2)=PT(I)+P(I,J)-P(I,K) + PC(I,3)=PT(I)-P(I,J)+P(I,K) + PC(I,4)=PT(I)-P(I,J)-P(I,K) + ENDDO + DO I=1,4 + TL=PC(1,I)**2+PC(2,I)**2+PC(3,I)**2 + IF (TL.GT.T) T=TL + ENDDO + ENDDO + ENDDO + T=SQRT(T)*ORS + ELSE + T=MAX(E(1),E(2),E(3))*ORS*2 + ENDIF + CALL GFILL1(203+ORD,T,(1-T)*WEIGHT*100) + IF (ORD.EQ.0) CS(5)=CS(5)+(1-T)*WEIGHT +C---CALCULATE THE Y3 VALUE (USES P SCHEME FOR NO PARTICULAR REASON) + IF (N.EQ.4) THEN + Y3=1D6 + DO I=2,4 + DO J=1,I-1 + Y=2*DOT(P,I,J)*OS + IF (Y.LT.Y3) THEN + Y3=Y + K=I + L=J + ENDIF + ENDDO + ENDDO + DO I=1,3 + Q(I,L)=P(I,L)+P(I,K) + ENDDO + Q(4,L)=SQRT(Q(1,L)**2+Q(2,L)**2+Q(3,L)**2) + IF (K.LT.4) THEN + DO I=1,4 + Q(I,K)=P(I,4) + ENDDO + ENDIF + DO I=1,3 + IF (I.NE.K.AND.I.NE.L) THEN + DO J=1,4 + Q(J,I)=P(J,I) + ENDDO + ENDIF + ENDDO + Y3=1D6 + DO I=2,3 + DO J=1,I-1 + Y=2*DOT(Q,I,J)*OS + IF (Y.LT.Y3) Y3=Y + ENDDO + ENDDO + ELSE + Y3=1-T + ENDIF + CALL GFILL1(204+ORD,Y3,Y3*WEIGHT*100) + IF (ORD.EQ.0) CS(6)=CS(6)+Y3*WEIGHT +C---CALCULATE THE ENERGY-ENERGY CORRELATION + DO I=2,N + DO J=1,I-1 + COSANG=1-DOT(P,I,J)*OE(I)*OE(J) + IF (ABS(COSANG).GE.1) COSANG=SIGN(1-1D-12,COSANG) + BEEC=(1-COSANG**2)*2*WEIGHT*E(I)*E(J)*OS + CALL GFILL1(205+ORD,COSANG,BEEC*50/2) + COSLIM=0.1D0 + IF (ORD.EQ.0) THEN + IF (ABS(COSANG).LT.COSLIM) + & CS(7)=CS(7)+BEEC/(2*COSLIM) + CS(8)=CS(8)+BEEC + DO K=1,2 + DO L=1,6 + BS(K,L)=BS(K,L)+ + $ BEEC*SQRT(1-COSANG**2)**(L-1)*COSANG**(K-1) + ENDDO + ENDDO + ENDIF + ENDDO + ENDDO + END +C----------------------------------------------------------------------- + FUNCTION DOT(P,I,J) + IMPLICIT NONE +C---RETURN THE DOT PRODUCT OF P(*,I) AND P(*,J) + INTEGER I,J + DOUBLE PRECISION DOT,P(4,7) + DOT=P(4,I)*P(4,J)-P(3,I)*P(3,J)-P(2,I)*P(2,J)-P(1,I)*P(1,J) + END +C----------------------------------------------------------------------- + SUBROUTINE GENTWO(EM,P,W,*) + IMPLICIT NONE +C---GENERATE A TWO-PARTON CONFIGURATION +C THE WEIGHT GIVES THE TOTAL VOLUME OF PHASE-SPACE + DOUBLE PRECISION EM,P(4,7),W,R(5),E,C,S,A + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + INTEGER NPOW1,NPOW2 + DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ + COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2 +C---GENERATE C=COS(THETA) ACCORDING TO COSZ+COSP*(1+C)**2+COSM*(1-C)**2 + CALL RANGEN(5,R) + E=EM/2 + IF (R(1)*(4*COSP+4*COSM+3*COSZ).LT.4*COSP) THEN + C=2*MAX(R(2),R(3),R(4))-1 + ELSEIF (R(1)*(4*COSP+4*COSM+3*COSZ).LT.4*COSP+4*COSM) THEN + C=2*MIN(R(2),R(3),R(4))-1 + ELSE + C=2*R(2)-1 + ENDIF + S=1-C**2 + IF (S.LT.0) CALL WARNER('GENTWO',1,*999) + S=SQRT(S) + A=2*PI*R(5) + P(1,1)= E*COS(A)*S + P(2,1)= E*SIN(A)*S + P(3,1)= E*C + P(4,1)= E + P(1,2)= -E*COS(A)*S + P(2,2)= -E*SIN(A)*S + P(3,2)= -E*C + P(4,2)= E + P(1,5)= 0 + P(2,5)= 0 + P(3,5)= 0 + P(4,5)=2*E + P(1,6)= 0 + P(2,6)= 0 + P(3,6)= E + P(4,6)= E + P(1,7)= 0 + P(2,7)= 0 + P(3,7)= -E + P(4,7)= E + W=1/(16*PI*EM**2) + RETURN + 999 RETURN 1 + END +C----------------------------------------------------------------------- + SUBROUTINE GENTHR(P,W,KEEP,*) + IMPLICIT NONE +C---GENERATE A THREE-PARTON CONFIGURATION FROM A GIVEN TWO-PARTON ONE +C IF ANY DOT PRODUCTS ARE BELOW CUTOFF, EVENT IS REJECTED, +C BY JUMPING TO THE SPECIFIED LABEL +C IF ANY DOT PRODUCTS ARE BELOW CUTUP, EVENT IS MARKED BY KEEP=.TRUE. + INTEGER EMIT,OMIT + LOGICAL KEEP + DOUBLE PRECISION P(4,7),W,X(2),R(4),EMSQ,DOT,T + INTEGER I + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + INTEGER NPOW1,NPOW2 + DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ + COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2 +C---GENERATE X1,X2 VALUES + CALL RANGEN(4,R) + X(1)=1-(R(1)*R(2))**NPOW1 + IF (X(1).GT.1) CALL WARNER('GENTHR',1,*999) + X(2)=2-X(1)-(1-X(1))**R(3) +C---ENFORCE INVARIANT MASS CUTOFFS + IF (1-X(1).LT.CUTOFF.OR.1-X(2).LT.CUTOFF.OR.X(1)+X(2)-1.LT.CUTOFF) + $ RETURN 1 + IF (1-X(1).LT.CUTUP.OR.1-X(2).LT.CUTUP.OR.X(1)+X(2)-1.LT.CUTUP) + $ KEEP=.TRUE. +C---CHOOSE WHICH ONE EMITTED + EMIT=INT(R(4)*2)+1 + OMIT=3-EMIT +C---GENERATE THEIR MOMENTA + CALL GENDIP(P,EMIT,OMIT,3,X,*999) +C---NOW SWAP EMITTER AND 3 HALF THE TIME + IF (MOD(INT(4*R(4)),2).EQ.1) THEN + DO I=1,4 + T=P(I,EMIT) + P(I,EMIT)=P(I,3) + P(I,3)=T + ENDDO + ENDIF +C---CALCULATE WEIGHT + EMSQ=DOT(P,5,5) + W=HF*EMSQ/(16*PISQ) + RETURN + 999 RETURN 1 + END +C----------------------------------------------------------------------- + SUBROUTINE GENFOR(P,W,KEEP,*) + IMPLICIT NONE +C---GENERATE A FOUR-BODY CONFIGURATION FROM A GIVEN THREE-BODY ONE +C IF ANY DOT PRODUCTS ARE BELOW CUTOFF, EVENT IS REJECTED, +C BY JUMPING TO THE SPECIFIED LABEL +C IF ANY DOT PRODUCTS ARE BELOW CUTUP, EVENT IS MARKED BY KEEP=.TRUE. + INTEGER I,J,EMIT,OMIT + LOGICAL KEEP + DOUBLE PRECISION P(4,7),W,X(2),R(6),T,EMSQ,DOT + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + INTEGER NPOW1,NPOW2 + DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ + COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2 +C---GENERATE X1,X2 VALUES + CALL RANGEN(6,R) + X(1)=1-(R(1)*R(2))**NPOW2 + IF (X(1).GT.1) CALL WARNER('GENFOR',1,*999) + X(2)=2-X(1)-(1-X(1))**R(3) +C---ENFORCE INVARIANT MASS CUTOFFS + W=0 + IF (1-X(1).LT.CUTOFF.OR.1-X(2).LT.CUTOFF.OR.X(1)+X(2)-1.LT.CUTOFF) + $ RETURN 1 +C---CHOOSE WHICH ONE TO SPLIT UNIFORMLY + EMIT=INT(R(4)*3)+1 +C---AND WHICH OF THE OTHERS TO BALANCE ITS MOMENTUM WITH + OMIT=INT(R(5)*2)+1 + IF (OMIT.EQ.EMIT) OMIT=3 +C---GENERATE THEIR MOMENTA + CALL GENDIP(P,EMIT,OMIT,4,X,*999) +C---HALF THE TIME, SWAP PARTICLES 3 AND 4 SO THEY ARE SYMMETRIC + IF (R(6).GT.0.5) THEN + DO I=1,4 + T=P(I,3) + P(I,3)=P(I,4) + P(I,4)=T + ENDDO + ENDIF +C---REJECT THE EVENT IF ANY PAIR MASSES ARE BELOW CUTOFF*Q**2 + EMSQ=DOT(P,5,5) + DO I=1,3 + DO J=I+1,4 + IF (2*DOT(P,I,J).LT.CUTOFF*EMSQ) RETURN 1 + IF (2*DOT(P,I,J).LT.CUTUP*EMSQ) KEEP=.TRUE. + ENDDO + ENDDO +C---CALCULATE WEIGHT (NOTE: THIS IS NORMALIZED TO THE TOTAL CMF MASS, +C RATHER THAN THE SUBSYSTEM MASS. THE EXTRA FACTOR GETS INSERTED IN +C THE JACOBIAN FACTOR CALCULATED IN SUBFOR) + W=HF*EMSQ/(16*PISQ) + RETURN + 999 RETURN 1 + END +C----------------------------------------------------------------------- + SUBROUTINE GENDIP(P,EMIT,OMIT,THRD,X,*) + IMPLICIT NONE +C---GENERATE A 2->3 DIPOLE EMISSION + DOUBLE PRECISION P(4,7),X(2),R(1),Q(4,2),PTDQ,C(4),D(4) + INTEGER EMIT,OMIT,THRD,I + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + CALL RANGEN(1,R) + PTDQ=(X(1)+X(2)-1)*(1-X(1))*(1-X(2))/X(1)**2 + IF (PTDQ.LT.0) CALL WARNER('GENDIP',1,*999) + PTDQ=SQRT(PTDQ) + CALL GTPERP(PTDQ,P,EMIT,OMIT,6,C,D,*999) + DO I=1,4 + Q(I,1)=X(1)*P(I,OMIT) + Q(I,2)=(1-(1-X(2))/X(1))*P(I,EMIT) + $ +(1-X(1))*(1-X(2))/X(1)*P(I,OMIT) + $ +COS(R(1)*2*PI)*C(I)+SIN(R(1)*2*PI)*D(I) + ENDDO + DO I=1,4 + P(I,THRD)=P(I,EMIT)+P(I,OMIT)-Q(I,1)-Q(I,2) + P(I,OMIT)=Q(I,1) + P(I,EMIT)=Q(I,2) + ENDDO + RETURN + 999 RETURN 1 + END +C----------------------------------------------------------------------- + SUBROUTINE GTPERP(PTDQ,P,I,J,K,C,D,*) + IMPLICIT NONE +C---FIND THE VECTORS PERPENDICULAR TO P(I) AND P(J) +C C AND D ARE PURELY SPACE-LIKE VECTORS IN THE P(I)+P(J) CMF, +C WITH C IN THE SAME PLANE AS P(K) AND D PERPENDICULAR TO IT, +C BOTH HAVING LENGTH PTDQ*SQRT(2*DOT(P,I,J)) + DOUBLE PRECISION PTDQ,P(4,7),C(4),D(4),PTF,DIJ,DIK,DJK,DOT,EPS4 + INTEGER I,J,K,L + DIJ=DOT(P,I,J) + DIK=DOT(P,I,K) + DJK=DOT(P,J,K) + IF (DIK.LE.0.OR.DJK.LE.0) CALL WARNER('GTPERP',1,*999) + PTF=PTDQ/SQRT(DIK*DJK) + DO L=1,4 + C(L)=PTF*(DIJ*P(L,K)-DJK*P(L,I)-DIK*P(L,J)) + ENDDO + DO L=1,4 + D(L)=EPS4(L,P(1,I),P(1,J),C)/DIJ + ENDDO + RETURN + 999 RETURN 1 + END +C----------------------------------------------------------------------- + FUNCTION EPS4(I,A,B,C) + IMPLICIT NONE + DOUBLE PRECISION EPS4,EPS3,A(4),B(4),C(4),AA(3),BB(3),CC(3) + INTEGER I,J,K,S(4) + DATA S/+1,-1,+1,+1/ + J=1 + DO K=1,3 + IF (I.EQ.J) J=J+1 + AA(K)=A(J) + BB(K)=B(J) + CC(K)=C(J) + J=J+1 + ENDDO + EPS4=0 + DO J=1,3 + EPS4=EPS4+CC(J)*EPS3(J,AA,BB) + ENDDO + EPS4=S(I)*EPS4 + END +C----------------------------------------------------------------------- + FUNCTION EPS3(I,A,B) + IMPLICIT NONE + DOUBLE PRECISION EPS3,A(3),B(3),AA(2),BB(2) + INTEGER I,J,K,S(3) + DATA S/+1,-1,+1/ + J=1 + DO K=1,2 + IF (I.EQ.J) J=J+1 + AA(K)=A(J) + BB(K)=B(J) + J=J+1 + ENDDO + EPS3=S(I)*(AA(1)*BB(2)-AA(2)*BB(1)) + END +C----------------------------------------------------------------------- + SUBROUTINE JACTWO(P,J) + IMPLICIT NONE +C---CALCULATE THE TWO-PARTON JACOBIAN FACTOR + DOUBLE PRECISION P(4,7),J,CP,CM,DOT + INTEGER NPOW1,NPOW2 + DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ + COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2 + CP=DOT(P,1,7)*DOT(P,5,5)/(DOT(P,1,5)*DOT(P,7,5)) + CM=DOT(P,1,6)*DOT(P,5,5)/(DOT(P,1,5)*DOT(P,6,5)) + J=(COSZ+COSP*CP**2+COSM*CM**2)*3/(3*COSZ+4*COSP+4*COSM) + END +C----------------------------------------------------------------------- + SUBROUTINE MATTWO(P,M) + IMPLICIT NONE +C---EVALUATE THE TWO-PARTON MATRIX ELEMENT SQUARED FOR THE GIVEN +C CONFIGURATION. USING ERT NORMALIZATION, THIS IS JUST UNITY. + DOUBLE PRECISION P(4,7),M,EMSQ,DOT + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + IF (METYPE.EQ.0) THEN + M=1 + ELSEIF (METYPE.EQ.1) THEN + EMSQ=DOT(P,5,5) + M=2*(-8*DOT(P,1,6)*DOT(P,2,6))+EMSQ/2*(4*EMSQ) +C---INCLUDE EXTERNAL FACTORS + M=M/4/EMSQ + M=M*FPAL2*CA*CQ*4/EMSQ + ELSEIF (METYPE.EQ.2) THEN + CALL CONVERT(P) + CALL SETUP + CALL SFMEE(2,0,M) + ELSE + STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN MATTWO!' + ENDIF + END +C----------------------------------------------------------------------- + SUBROUTINE MATTHR(P,M) + IMPLICIT NONE +C---EVALUATE THE THREE-PARTON MATRIX ELEMENT SQUARED FOR THE GIVEN +C CONFIGURATION. INCLUDES COLOUR FACTORS, BUT NOT ALPHA_S/2PI. + DOUBLE PRECISION P(4,7),M,EMSQ,E1,E2,S13,S23,DOT,K1,K2,K3 + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + IF (METYPE.EQ.0) THEN + EMSQ=DOT(P,5,5) + S13=2*DOT(P,1,3) + S23=2*DOT(P,2,3) + E1=EMSQ-S23 + E2=EMSQ-S13 + M=CF*16*PISQ*(E1**2+E2**2)/(EMSQ*S13*S23) + ELSEIF (METYPE.EQ.1) THEN + EMSQ=DOT(P,5,5) + S13=2*DOT(P,1,3) + S23=2*DOT(P,2,3) + E1=EMSQ-S23 + E2=EMSQ-S13 + K1=2*DOT(P,1,6) + K2=2*DOT(P,2,6) + K3=2*DOT(P,3,6) + M=CF*16*PISQ/(S13*S23)*( + $ 2*(2*(K1**2*S23+K2**2*S13 + $ -K1*K2*(E1+E2)-K2*K3*E2-K3*K1*E1)) + $ +EMSQ/2*(4*(E1**2+E2**2))) +C---INCLUDE EXTERNAL FACTORS + M=M/4/EMSQ + M=M*FPAL2*CA*CQ*4/EMSQ + ELSEIF (METYPE.EQ.2) THEN + CALL CONVERT(P) + CALL SETUP + CALL SFMEE(3,0,M) + ELSE + STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN MATTHR!' + ENDIF + END +C----------------------------------------------------------------------- + SUBROUTINE CONTHR(P,V,VV,M) + IMPLICIT NONE +C---EVALUATE THE CONTRACTION OF THE THREE-PARTON MATRIX ELEMENT SQUARED +C FOR THE GIVEN CONFIGURATION WITH THE VECTOR V. IT IS ASSUMED THAT +C V IS PERPENDICULAR TO THE GLUON, V.P3.EQ.0, AND THAT V.V.NE.0 . +C THE NORMALIZATION IS SUCH THAT THE AVERAGE OVER V'S AZIMUTH AROUND +C P3 IS EQUAL TO HALF OF MATTHR + INTEGER I,J + DOUBLE PRECISION P(4,7),V(4),M,DOT,VDOT,EMSQ,OMSQ,VV,Y1,Y2,L1,L2, + $ DOTS,T1,T2,T3,T4,T +C---EERAD COMMON BLOCKS + DOUBLE PRECISION CZUR,CZUL,CZDR,CZDL,CZLR,CZLL,CZNR,CZNL, + $ CGUR,CGUL,CGDR,CGDL,CGLR,CGLL,CGNR,CGNL, + $ W,EM,SW2,RMZ,RGZ, + $ CGGVV,CGZVV,CZZVV,CVV,CGGAA,CGZAA,CZZAA,CAA + INTEGER NUTF,NDTF,ICON + COMMON /ZCUP/CZUR,CZUL,CZDR,CZDL,CZLR,CZLL,CZNR,CZNL + COMMON /GCUP/CGUR,CGUL,CGDR,CGDL,CGLR,CGLL,CGNR,CGNL + COMMON /GENERA/ W,EM,SW2,RMZ,RGZ,NUTF,NDTF,ICON +C---END EERAD COMMON BLOCKS + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + VDOT(I)=P(4,I)*V(4)-P(3,I)*V(3)-P(2,I)*V(2)-P(1,I)*V(1) + IF (METYPE.EQ.0) THEN + EMSQ=DOT(P,5,5) + OMSQ=1/EMSQ + Y1=2*DOT(P,2,3) + Y2=2*DOT(P,1,3) + T=CF*16*PISQ/(Y1*Y2)**2 + M=-T*(2*(VDOT(1)*Y1-VDOT(2)*Y2)**2/VV + $ -HF*(Y1**2+Y2**2)*Y1*Y2*OMSQ) + ELSEIF (METYPE.EQ.1.OR.METYPE.EQ.2) THEN + EMSQ=DOT(P,5,5) + OMSQ=1/EMSQ + Y1=2*DOT(P,2,3) + Y2=2*DOT(P,1,3) + L1=2*DOT(P,1,7) + L2=2*DOT(P,2,7) + DOTS=1/(VV*(Y1+Y2)**2) +C---INCLUDE COUPLING STRENGTHS + IF (METYPE.EQ.1) THEN + CVV=CQ*FPAL2 + CAA=0 + ELSE + cggvv=nutf*0.25*(cgur**2+cgul**2)*(cglr**2+cgll**2) + $ +ndtf*0.25*(cgdr**2+cgdl**2)*(cglr**2+cgll**2) + cggaa=nutf*0.25*(cgur**2-cgul**2)*(cglr**2-cgll**2) + $ +ndtf*0.25*(cgdr**2-cgdl**2)*(cglr**2-cgll**2) + cgzvv=nutf*0.5*(cgur*czur+cgul*czul)*(cglr*czlr+cgll*czll) + $ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2) + $ +ndtf*0.5*(cgdr*czdr+cgdl*czdl)*(cglr*czlr+cgll*czll) + $ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2) + cgzaa=nutf*0.5*(cgur*czur-cgul*czul)*(cglr*czlr-cgll*czll) + $ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2) + $ +ndtf*0.5*(cgdr*czdr-cgdl*czdl)*(cglr*czlr-cgll*czll) + $ *emsq*(emsq-rmz**2)/((emsq-rmz**2)**2+(rgz*rmz)**2) + czzvv=nutf*0.25*(czur**2+czul**2)*(czlr**2+czll**2) + $ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2) + $ +ndtf*0.25*(czdr**2+czdl**2)*(czlr**2+czll**2) + $ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2) + czzaa=nutf*0.25*(czur**2-czul**2)*(czlr**2-czll**2) + $ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2) + $ +ndtf*0.25*(czdr**2-czdl**2)*(czlr**2-czll**2) + $ *emsq**2/((emsq-rmz**2)**2+(rgz*rmz)**2) + if (icon.eq.1) then + cvv=cggvv + caa=cggaa + elseif (icon.eq.2) then + cvv=cgzvv + caa=cgzaa + elseif (icon.eq.3) then + cvv=czzvv + caa=czzaa + else + cvv=cggvv+cgzvv+czzvv + caa=cggaa+cgzaa+czzaa + endif + cvv=cvv*em**4 + caa=caa*em**4 + ENDIF +C---THE VECTOR-VECTOR PART + T=CF*16*PISQ*CA*CVV/(Y1*Y2)**2 + T1=T*4*(Y1**2*(L2**2+(EMSQ-L2)**2)+Y2**2*(L1**2+(EMSQ-L1)**2) + $ +2*Y1*Y2*(L1*EMSQ+L2*EMSQ-2*L1*L2))*OMSQ**2*DOTS + T2=T*4*(Y2*(EMSQ-2*L1)-Y1*(EMSQ-2*L2))*Y1*Y2*OMSQ**2*DOTS + T3=T*8*(Y1*Y2*OMSQ)**2*DOTS + T4=T*((EMSQ-Y1-L1-L2)**2+(EMSQ-Y2-L1-L2)**2) + $ *Y1*Y2*OMSQ +C---THE AXIAL-AXIAL PART + T=CF*16*PISQ*CA*CAA/(Y1*Y2)**2 + T1=T1+T*4*(Y1*(EMSQ-2*L2)-Y2*(EMSQ-2*L1))*(Y1+Y2)*OMSQ*DOTS + T2=T2-T*4*(Y1+Y2)*Y1*Y2*OMSQ*DOTS + T3=T3 + T4=T4-T*(Y1+Y2+2*(L1+L2-EMSQ))*(Y1-Y2)*Y1*Y2*OMSQ + M=-T1*(VDOT(1)*Y1-VDOT(2)*Y2)**2 + $ -T2*2*(VDOT(1)*Y1-VDOT(2)*Y2)* + $ (VDOT(7)*(Y2+Y1-EMSQ+L1+L2)-(EMSQ-L1-L2)*VDOT(6)) + $ -T3*(VDOT(7)*(Y2+Y1-EMSQ+L1+L2)-(EMSQ-L1-L2)*VDOT(6))**2 + $ +T4 + ELSE + STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN CONTHR!' + ENDIF + END +C----------------------------------------------------------------------- + SUBROUTINE MATFOR(P,M) + IMPLICIT NONE +C---EVALUATE THE FOUR-PARTON MATRIX ELEMENT SQUARED FOR THE GIVEN +C CONFIGURATION. INCLUDES COLOUR FACTORS, BUT NOT (ALPHA_S/2PI)**2. + DOUBLE PRECISION P(4,7),M,ERTA,ERTB,ERTC,ERTD,ERTE, + $ LEIA,LEIB,LEIC,LEID,LEIE, + $ A,B,C,DS,DD,E,EMSQ,DOT,T + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB +C---EERAD COMMON BLOCKS + INTEGER IPS + COMMON /PERM/IPS +C---END EERAD COMMON BLOCKS + IF (METYPE.EQ.0) THEN + A=HF*( + $ ERTA(P,1,2,3,4)+ERTA(P,2,1,3,4)+ERTA(P,1,2,4,3)+ERTA(P,2,1,4,3)) + B=HF*( + $ ERTB(P,1,2,3,4)+ERTB(P,2,1,3,4)+ERTB(P,1,2,4,3)+ERTB(P,2,1,4,3)) + C=HF*( + $ ERTC(P,1,2,3,4)+ERTC(P,2,1,3,4)+ERTC(P,1,2,4,3)+ERTC(P,2,1,4,3)) + DD=HF*( + $ ERTD(P,1,2,3,4)+ERTD(P,2,1,4,3)+ERTD(P,3,4,1,2)+ERTD(P,4,3,2,1)) + DS=HF*HF*( + $ ERTD(P,2,1,3,4)+ERTD(P,1,2,4,3)+ERTD(P,3,4,2,1)+ERTD(P,4,3,1,2)) + E=HF*HF*( + $ ERTE(P,1,2,3,4)+ERTE(P,2,1,4,3)+ERTE(P,3,4,1,2)+ERTE(P,4,3,2,1)+ + $ ERTE(P,2,1,3,4)+ERTE(P,1,2,4,3)+ERTE(P,3,4,2,1)+ERTE(P,4,3,1,2)) + CFSUB=A+B+E + CASUB=C-(B+E)/2 + TFSUB=DD+(DS-HF*DD)*ONF + M=CF*CFSUB+CA*CASUB+NF*TR*TFSUB + CFSUB=CFSUB/(M*CF) + CASUB=CASUB/(M*CF) + TFSUB=TFSUB/(M*CF) + GGSUB=(CF*A+(CF-CA/2)*B+CA*C)/M + QQSUB=(TR*(DD*HF+DS)+(CF-CA/2)*E)/M + QPSUB=TR*(NF-1)*DD/M +C---INCLUDE EXTERNAL FACTORS + EMSQ=DOT(P,5,5) + M=M*256*PI**4*CF/EMSQ + ELSEIF (METYPE.EQ.1) THEN + EMSQ=DOT(P,5,5) + A=2*(LEIA(P,P(1,6),1,2,3,4)+LEIA(P,P(1,6),2,1,3,4) + $ +LEIA(P,P(1,6),1,2,4,3)+LEIA(P,P(1,6),2,1,4,3))+EMSQ/2*( + $ ERTA(P,1,2,3,4)+ERTA(P,2,1,3,4)+ERTA(P,1,2,4,3)+ERTA(P,2,1,4,3)) + B=2*(LEIB(P,P(1,6),1,2,3,4)+LEIB(P,P(1,6),2,1,3,4) + $ +LEIB(P,P(1,6),1,2,4,3)+LEIB(P,P(1,6),2,1,4,3))+EMSQ/2*( + $ ERTB(P,1,2,3,4)+ERTB(P,2,1,3,4)+ERTB(P,1,2,4,3)+ERTB(P,2,1,4,3)) + C=2*(LEIC(P,P(1,6),1,2,3,4)+LEIC(P,P(1,6),2,1,3,4) + $ +LEIC(P,P(1,6),1,2,4,3)+LEIC(P,P(1,6),2,1,4,3))+EMSQ/2*( + $ ERTC(P,1,2,3,4)+ERTC(P,2,1,3,4)+ERTC(P,1,2,4,3)+ERTC(P,2,1,4,3)) + DD=2*(LEID(P,P(1,6),1,2,3,4)+LEID(P,P(1,6),2,1,4,3) + $ +LEID(P,P(1,6),3,4,1,2)+LEID(P,P(1,6),4,3,2,1))+EMSQ/2*( + $ ERTD(P,1,2,3,4)+ERTD(P,2,1,4,3)+ERTD(P,3,4,1,2)+ERTD(P,4,3,2,1)) + DS=2*(LEID(P,P(1,6),2,1,3,4)+LEID(P,P(1,6),1,2,4,3) + $ +LEID(P,P(1,6),3,4,2,1)+LEID(P,P(1,6),4,3,1,2))+EMSQ/2*( + $ ERTD(P,2,1,3,4)+ERTD(P,1,2,4,3)+ERTD(P,3,4,2,1)+ERTD(P,4,3,1,2)) + E=2*(LEIE(P,P(1,6),1,2,3,4)+LEIE(P,P(1,6),2,1,4,3) + $ +LEIE(P,P(1,6),3,4,1,2)+LEIE(P,P(1,6),4,3,2,1) + $ +LEIE(P,P(1,6),2,1,3,4)+LEIE(P,P(1,6),1,2,4,3) + $ +LEIE(P,P(1,6),3,4,2,1)+LEIE(P,P(1,6),4,3,1,2))+EMSQ/2*( + $ ERTE(P,1,2,3,4)+ERTE(P,2,1,4,3)+ERTE(P,3,4,1,2)+ERTE(P,4,3,2,1)+ + $ ERTE(P,2,1,3,4)+ERTE(P,1,2,4,3)+ERTE(P,3,4,2,1)+ERTE(P,4,3,1,2)) + A=HF*A + B=HF*B + C=HF*C + DD=HF*DD + DS=HF*HF*DS + E=HF*HF*E + CFSUB=A+B+E + CASUB=C-(B+E)/2 + TFSUB=DD+(DS-HF*DD)*ONF + M=CF*CFSUB+CA*CASUB+NF*TR*TFSUB + CFSUB=CFSUB/(M*CF) + CASUB=CASUB/(M*CF) + TFSUB=TFSUB/(M*CF) + GGSUB=(CF*A+(CF-CA/2)*B+CA*C)/M + QQSUB=(TR*(DD*HF+DS)+(CF-CA/2)*E)/M + QPSUB=TR*(NF-1)*DD/M +C---INCLUDE EXTERNAL FACTORS + M=M*256*PI**4*CF/EMSQ + M=M*FPAL2*CA*CQ*4/EMSQ + ELSEIF (METYPE.EQ.2) THEN + CALL CONVERT(P) + CALL SETUP + IPS=1 + CALL SFMEE(4,0,M) + IPS=2 + CALL SFMEE(4,0,T) + M=M+T + ELSE + STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN MATFOR!' + ENDIF + END +C----------------------------------------------------------------------- + SUBROUTINE VIRTWO(P,V) + IMPLICIT NONE +C---CALCULATE THE TWO-PARTON MATRIX-ELEMENT AT NEXT-TO-LEADING ORDER + DOUBLE PRECISION P(4,7),V,M,SUB + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE +C---CALCULATE THE LOWEST-ORDER MATRIX-ELEMENT + CALL MATTWO(P,M) + IF (METYPE.EQ.0.OR.METYPE.EQ.1) THEN +C---VIRTUAL CROSS-SECTION, WITH THE POLES REMOVED AND THE LIMIT TAKEN + V=M*CF*(PISQ-8) + ELSEIF (METYPE.EQ.2) THEN + CALL CONVERT(P) + CALL SETUP + CALL SFMEE(2,1,V) + ELSE + STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN VIRTWO!' + ENDIF +C---THE SUBTRACTION + SUB=CF*(10-PISQ) +C---THE TOTAL + V=V + M*SUB + END +C----------------------------------------------------------------------- + SUBROUTINE VIRTHR(P,V) + IMPLICIT NONE +C---CALCULATE THE THREE-PARTON MATRIX-ELEMENT AT NEXT-TO-LEADING ORDER + DOUBLE PRECISION P(4,7),V,DOT,EMSQ,Y12,Y13,Y23,L12,L13,L23, + $ M,ERTV,LEIV,CFTOT,CATOT,TFTOT,TOT + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + EMSQ=DOT(P,5,5) + Y12=2*DOT(P,1,2)/EMSQ + Y13=2*DOT(P,1,3)/EMSQ + Y23=2*DOT(P,2,3)/EMSQ + L12=LOG(Y12) + L13=LOG(Y13) + L23=LOG(Y23) +C---CALCULATE THE LOWEST-ORDER MATRIX-ELEMENT + CALL MATTHR(P,M) + IF (METYPE.EQ.0) THEN +C---THE FIRST ERT EXPRESSION, WITH THE POLES REMOVED AND THE LIMIT TAKEN + CFTOT=M*(HF*PISQ*2 - 8 - L12**2) + CATOT=M*(HF*PISQ + HF*(L12**2-L13**2-L23**2)) + TFTOT=0 +C---THE FINITE TERMS NOT PROPORTIONAL TO TREE-LEVEL RESULT + TOT=ERTV(P) + CFTOT=CFTOT+CF*CFSUB*TOT + CATOT=CATOT+CF*CASUB*TOT + TFTOT=TFTOT+CF*TFSUB*TOT + ELSEIF (METYPE.EQ.1) THEN +C---THE FIRST ERT EXPRESSION, WITH THE POLES REMOVED AND THE LIMIT TAKEN + CFTOT=M*(HF*PISQ*2 - 8 - L12**2) + CATOT=M*(HF*PISQ + HF*(L12**2-L13**2-L23**2)) + TFTOT=0 +C---THE FINITE TERMS NOT PROPORTIONAL TO TREE-LEVEL RESULT + TOT=(FPAL2*CA*CQ*4/EMSQ)*(2*LEIV(P,P(1,6))) + CFTOT=CFTOT+CF*CFSUB*TOT + CATOT=CATOT+CF*CASUB*TOT + TFTOT=TFTOT+CF*TFSUB*TOT + TOT=(FPAL2*CA*CQ*4/EMSQ)*(EMSQ/2*ERTV(P)) + CFTOT=CFTOT+CF*CFSUB*TOT + CATOT=CATOT+CF*CASUB*TOT + TFTOT=TFTOT+CF*TFSUB*TOT + ELSEIF (METYPE.EQ.2) THEN + CALL CONVERT(P) + CALL SETUP + CALL SFMEE(3,1,V) + CFTOT=V/CF + CATOT=0 + TFTOT=0 + ELSE + STOP 'UNKNOWN MATRIX-ELEMENT TYPE IN VIRTHR!' + ENDIF +C---THE SUBTRACTION + CFTOT=CFTOT+M*(-HF*PISQ*2 + 10 + L12**2 - 3D0/2*(2*L12)) + CATOT=CATOT+M*(-HF*PISQ + 50D0/9 - HF*(L12**2-L13**2-L23**2) + $ - HF*(11D0/6)*(L13+L23) - 3D0/2*(HF*(L13+L23-2*L12))) + TFTOT=TFTOT+M*( - 16D0/9 - HF*(-2D0/3)*(L13+L23)) +C---THE TOTAL + V=CF*CFTOT+CA*CATOT+TR*NF*TFTOT + CFSUB=CFTOT/(CF*V) + CASUB=CATOT/(CF*V) + TFSUB=TFTOT/(CF*V) + GGSUB=CF*(CF*CFSUB+CA*CASUB) + QQSUB=CF*TR*TFSUB + QPSUB=CF*TR*(NF-1)*TFSUB + END +C----------------------------------------------------------------------- + DOUBLE PRECISION FUNCTION ERTV(P) + IMPLICIT NONE +C---RETURN THE ERT F FUNCTION FOR VIRTUAL TERMS THAT ARE NOT +C TRIVIALLY PROPORTIONAL TO TREE-LEVEL + DOUBLE PRECISION P(4,7),DOT,R,X,Y,DILOG,EMSQ, + $ Y12,Y13,Y23,L12,L13,L23 + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + R(X,Y)=LOG(X)*LOG(Y)-LOG(X)*LOG(1-X)-LOG(Y)*LOG(1-Y)+PISQ/6 + $ -DILOG(X)-DILOG(Y) + EMSQ=DOT(P,5,5) + Y12=2*DOT(P,1,2)/EMSQ + Y13=2*DOT(P,1,3)/EMSQ + Y23=2*DOT(P,2,3)/EMSQ + L12=LOG(Y12) + L13=LOG(Y13) + L23=LOG(Y23) + CFSUB=Y12/(Y12+Y13)+Y12/(Y12+Y23)+(Y12+Y23)/Y13+(Y12+Y13)/Y23 + $ +L13*(4*Y12**2+2*Y12*Y13+4*Y12*Y23+Y13*Y23)/(Y12+Y23)**2 + $ +L23*(4*Y12**2+2*Y12*Y23+4*Y12*Y13+Y13*Y23)/(Y12+Y13)**2 + $ -2*((Y12**2+(Y12+Y13)**2)/(Y13*Y23)*R(Y12,Y23) + $ +(Y12**2+(Y12+Y23)**2)/(Y13*Y23)*R(Y12,Y13) + $ +(Y13**2+Y23**2)/(Y13*Y23*(Y13+Y23)) + $ -2*L12*(Y12**2/(Y13+Y23)**2+2*Y12/(Y13+Y23))) + CASUB=L13*Y13/(Y12+Y23)+L23*Y23/(Y12+Y13) + $ +((Y12**2+(Y12+Y13)**2)/(Y13*Y23)*R(Y12,Y23) + $ +(Y12**2+(Y12+Y23)**2)/(Y13*Y23)*R(Y12,Y13) + $ +(Y13**2+Y23**2)/(Y13*Y23*(Y13+Y23)) + $ -2*L12*(Y12**2/(Y13+Y23)**2+2*Y12/(Y13+Y23))) + $ -R(Y13,Y23)*(Y13/Y23+Y23/Y13+2*Y12/(Y13*Y23)) + TFSUB=0 + ERTV=CF*(CF*CFSUB+CA*CASUB) + CFSUB=CFSUB/ERTV + CASUB=CASUB/ERTV + ERTV=16*PISQ/EMSQ*ERTV + END +C----------------------------------------------------------------------- + DOUBLE PRECISION FUNCTION LEIV(P,K) + IMPLICIT NONE +C---RETURN THE PART OF THE ONE-LOOP HADRONIC TENSOR NOT TRIVIALLY +C PROPORTIONAL TO TREE-LEVEL ONE CONTRACTED WITH A TENSOR -K(MU)K(NU) + DOUBLE PRECISION P(4,7),K(4),DOT,R,X,Y,Z,DILOG,EMSQ,Y12,Y13,Y23, + $ L12,L13,L23,A,B,C,D,P1K,P2K,P3K,KK, + $ CF1,CF2,CF3,CF4,CA1,CA2,CA3,CA4 + PARAMETER (Z=0) + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + R(X,Y)=LOG(X)*LOG(Y)-LOG(X)*LOG(1-X)-LOG(Y)*LOG(1-Y)+PISQ/6 + $ -DILOG(X)-DILOG(Y) + EMSQ=DOT(P,5,5) + Y12=2*DOT(P,1,2)/EMSQ + Y13=2*DOT(P,1,3)/EMSQ + Y23=2*DOT(P,2,3)/EMSQ + L12=LOG(Y12) + L13=LOG(Y13) + L23=LOG(Y23) +C---FIRST DECOMPOSE K AS A*P1+B*P2+C*P3+D*E_(MU,P1,P2,P3)/Q^2 + P1K=P(4,1)*K(4)-P(3,1)*K(3)-P(2,1)*K(2)-P(1,1)*K(1) + P2K=P(4,2)*K(4)-P(3,2)*K(3)-P(2,2)*K(2)-P(1,2)*K(1) + P3K=P(4,3)*K(4)-P(3,3)*K(3)-P(2,3)*K(2)-P(1,3)*K(1) + KK=K(4)*K(4)-K(3)*K(3)-K(2)*K(2)-K(1)*K(1) + A=(P3K*Y12+P2K*Y13-P1K*Y23)/(Y12*Y13*EMSQ) + B=(P1K*Y23+P3K*Y12-P2K*Y13)/(Y12*Y23*EMSQ) + C=(P2K*Y13+P1K*Y23-P3K*Y12)/(Y23*Y13*EMSQ) + D=SQRT(MAX(Z,4*(KK/EMSQ-A*B*Y12-A*C*Y13-B*C*Y23)/(-Y12*Y13*Y23))) +C---THEN CALCULATE THE CONTRACTIONS WITH -P1(MU)P1(NU) ETC + CF1=R(Y13,Y12)*Y12-R(Y23,Y12)*Y12**2/Y13 + $ -L12*Y12*Y13/(Y13+Y23)**2-L23*Y12+Y12/2*(Y23-Y13)/(Y23+Y13) + CA1=-R(Y13,Y12)*Y12/2+R(Y23,Y12)*Y12**2/Y13/2+R(Y23,Y13)*Y12/2 + $ +L12*Y12*Y13/(Y13+Y23)**2/2-Y12*Y23/(Y23+Y13)/2 + CF2=R(Y23,Y12)*Y12-R(Y13,Y12)*Y12**2/Y23 + $ -L12*Y12*Y23/(Y13+Y23)**2-L13*Y12+Y12/2*(Y13-Y23)/(Y13+Y23) + CA2=-R(Y23,Y12)*Y12/2+R(Y13,Y12)*Y12**2/Y23/2+R(Y23,Y13)*Y12/2 + $ +L12*Y12*Y23/(Y13+Y23)**2/2-Y12*Y13/(Y13+Y23)/2 + CF3=R(Y13,Y12)*Y12+R(Y23,Y12)*Y12-2*L12*Y12 + $ -L13*Y12*Y23/2/(Y12+Y23)**2*(1+Y12+Y23) + $ -L23*Y12*Y13/2/(Y12+Y13)**2*(1+Y12+Y13) + $ -Y12/2*(Y13/(Y12+Y13)+Y23/(Y12+Y23)) + CA3=-R(Y13,Y12)*Y12/2-R(Y23,Y12)*Y12/2+R(Y23,Y13)*Y12 + $ +L12*Y12+L13*Y12*Y13/(Y12+Y23)/2+L23*Y12*Y23/(Y12+Y13)/2 + CF4=R(Y13,Y12)*Y12/4/Y23*(Y12-Y12**2-Y12*Y23+2*Y12**2*Y23 + $ +2*Y12*Y23**2+Y23**3) + $ +R(Y23,Y12)*Y12/4/Y13*(Y12-Y12**2-Y12*Y13+2*Y12**2*Y13 + $ +2*Y12*Y13**2+Y13**3) + $ +L12*Y12**2/4/(Y13+Y23)*(Y13+Y23-2*Y13*Y23) + $ -L13*Y12/8/(Y12+Y23)*Y13*(Y12+Y13)*(Y23-2*Y12) + $ -L23*Y12/8/(Y12+Y13)*Y23*(Y12+Y23)*(Y13-2*Y12) + $ +Y12/8*(Y13+Y23-2*Y13*Y23) + CA4=-R(Y13,Y12)*Y12/8/Y23*(Y12-Y12**2-Y12*Y23+2*Y12**2*Y23 + $ +2*Y12*Y23**2+Y23**3) + $ -R(Y23,Y12)*Y12/8/Y13*(Y12-Y12**2-Y12*Y13+2*Y12**2*Y13 + $ +2*Y12*Y13**2+Y13**3) + $ +R(Y13,Y23)*Y12/8*((1-Y13)**2+(1-Y23)**2) + $ -L12*Y12**2/8/(Y13+Y23)*(Y13+Y23-2*Y13*Y23) + $ -L13*Y12/8*Y13*(Y12+Y13)-L23*Y12/8*Y23*(Y12+Y23) + $ -Y12/8*(Y13+Y23-2*Y13*Y23) +C---AND COMBINE THEM WITH THE APPROPRIATE WEIGHTS + CFSUB=(A-B)*(A-C)*CF1+(B-C)*(B-A)*CF2+(C-A)*(C-B)*CF3+D*D*CF4 + CASUB=(A-B)*(A-C)*CA1+(B-C)*(B-A)*CA2+(C-A)*(C-B)*CA3+D*D*CA4 + TFSUB=0 + LEIV=CF*(CF*CFSUB+CA*CASUB) + CFSUB=CFSUB/LEIV + CASUB=CASUB/LEIV + LEIV=16*PISQ*LEIV + END +C----------------------------------------------------------------------- + SUBROUTINE SUBTHR(PERM,P,Q,S,JAC,*) + IMPLICIT NONE +C---GENERATE A TWO-PARTON STATE FROM A THREE-PARTON STATE, +C CALCULATE THE JACOBIAN FACTOR FOR THE CORRESPONDING CHANNEL, +C AND (IF PERM.GT.0) THE APPROXIMATE MATRIX-ELEMENT. + INTEGER NPERM,PERM,IPERM,IJF,KF,I,J,K,M + DOUBLE PRECISION P(4,7),Q(4,7),S,JAC,DOT,EMSQ,DEN,Y,ZI,ZJ,OY,ZTI + PARAMETER (NPERM=2) + DIMENSION IPERM(3,NPERM),IJF(NPERM),KF(NPERM) + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + INTEGER NPOW1,NPOW2 + DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ + COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2 + DATA IPERM/1,3,2, 2,3,1/ + DATA IJF/1,2/ + DATA KF/2,1/ + IF (ABS(PERM).GT.NPERM.OR.PERM.EQ.0) + $ STOP 'PERM TOO BIG IN SUBTHR!' +C---FIND WHICH PARTONS TO TREAT + I=IPERM(1,ABS(PERM)) + J=IPERM(2,ABS(PERM)) + K=IPERM(3,ABS(PERM)) +C---FIND THE KINEMATIC VARIABLES + EMSQ=2*(DOT(P,I,J)+DOT(P,J,K)+DOT(P,K,I)) + IF (EMSQ.LE.0) CALL WARNER('SUBTHR',1,*999) + DEN=2/EMSQ + Y=DOT(P,I,J)*DEN + IF (Y.LE.0.OR.Y.GE.1) CALL WARNER('SUBTHR',2,*999) + ZI=DOT(P,I,K)*DEN + IF (ZI.LE.0.OR.ZI.GE.1) CALL WARNER('SUBTHR',3,*999) + ZJ=DOT(P,J,K)*DEN + IF (ZJ.LE.0.OR.ZJ.GE.1) CALL WARNER('SUBTHR',4,*999) + OY=1/(1-Y) + ZTI=ZI*OY +C---COPY INTO Q, REPLACING K BY KTILDE AND I BY IJTILDE + DO M=1,4 + Q(M,IJF(ABS(PERM)))=P(M,I)+P(M,J)-Y*OY*P(M,K) + Q(M,KF(ABS(PERM)))=OY*P(M,K) + Q(M,5)=P(M,5) + Q(M,6)=P(M,6) + Q(M,7)=P(M,7) + ENDDO +C---CALCULATE THE CORRESPONDING JACOBIAN FACTOR + CALL JACTWO(Q,JAC) + JAC=JAC/(4*NPOW1**2*(1-ZI)*(1-ZJ)*Y**XPOW1)*(2-ZI-ZJ) +C---INCLUDE A PRIORI CHANNEL WEIGHTS + JAC=JAC/2 + IF (PERM.LT.0) RETURN +C---CALCULATE WEIGHT FOR QUARK-GLUON SPLITTING FUNCTION + CALL MATTWO(Q,S) + S=CF*S*(2/(1-ZTI*(1-Y))-(1+ZTI)) +C---INCLUDE OVERALL FACTOR + DEN=1/(Y*EMSQ) + S=S*16*PISQ*DEN + RETURN + 999 RETURN 1 + END +C----------------------------------------------------------------------- + SUBROUTINE SUBFOR(PERM,P,Q,S,JAC,KEEP,*) + IMPLICIT NONE +C---GENERATE A THREE-PARTON STATE FROM A FOUR-PARTON STATE, +C CALCULATE THE JACOBIAN FACTOR FOR THE CORRESPONDING CHANNEL, +C AND (IF PERM.GT.0) THE APPROXIMATE MATRIX-ELEMENT. +C IF ANY DOT PRODUCTS ARE BELOW CUTOFF, EVENT IS REJECTED, +C BY JUMPING TO THE SPECIFIED LABEL +C IF ANY DOT PRODUCTS ARE BELOW CUTUP, EVENT IS MARKED BY KEEP=.TRUE. + INTEGER NPERM,PERM,IPERM,IJF,KF,LF,NPERM3,I,J,K,L,M,N + LOGICAL KEEP + DOUBLE PRECISION P(4,7),Q(4,7),S,S1,S2,S3,JAC,DOT,EMSQ,QUSQ,T, + $ DEN,Y,ZI,ZJ,OY,ZTI,ZTJ,OX1,OX2,JTMP,STMP,QTMP(4,7),V(4),VV + PARAMETER (NPERM=10,NPERM3=2) + DIMENSION IPERM(4,NPERM),IJF(NPERM),KF(NPERM),LF(NPERM) + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + INTEGER NPOW1,NPOW2 + DOUBLE PRECISION XPOW1,XPOW2,COSP,COSM,COSZ + COMMON /SAMCOM/ XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2 + DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB + DATA IPERM/1,3,2,4, 1,3,4,2, 1,4,2,3, 1,4,3,2, + $ 2,4,1,3, 2,4,3,1, 2,3,1,4, 2,3,4,1, + $ 3,4,1,2, 3,4,2,1/ + DATA IJF/1,1,1,1,2,2,2,2,3,3/ + DATA KF/2,3,2,3,1,3,1,3,1,2/ + DATA LF/3,2,3,2,3,1,3,1,2,1/ + IF (ABS(PERM).GT.NPERM.OR.PERM.EQ.0) + $ STOP 'PERM TOO BIG IN SUBFOR!' +C---FIND WHICH PARTONS TO TREAT + I=IPERM(1,ABS(PERM)) + J=IPERM(2,ABS(PERM)) + K=IPERM(3,ABS(PERM)) + L=IPERM(4,ABS(PERM)) +C---FIND THE KINEMATIC VARIABLES + EMSQ=2*(DOT(P,I,J)+DOT(P,J,K)+DOT(P,K,I)) + IF (EMSQ.LE.0) CALL WARNER('SUBFOR',1,*999) + DEN=2/EMSQ + Y=DOT(P,I,J)*DEN + IF (Y.LE.0.OR.Y.GE.1) CALL WARNER('SUBFOR',2,*999) + ZI=DOT(P,I,K)*DEN + IF (ZI.LE.0.OR.ZI.GE.1) CALL WARNER('SUBFOR',3,*999) + ZJ=DOT(P,J,K)*DEN + IF (ZJ.LE.0.OR.ZJ.GE.1) CALL WARNER('SUBFOR',4,*999) + OY=1/(1-Y) + ZTI=ZI*OY + ZTJ=ZJ*OY +C---COPY INTO Q, REPLACING K BY KTILDE AND I BY IJTILDE + DO M=1,4 + Q(M,IJF(ABS(PERM)))=P(M,I)+P(M,J)-Y*OY*P(M,K) + Q(M,KF(ABS(PERM)))=OY*P(M,K) + Q(M,LF(ABS(PERM)))=P(M,L) + Q(M,5)=P(M,5) + Q(M,6)=P(M,6) + Q(M,7)=P(M,7) + ENDDO +C---REENFORCE THE MINIMUM CUTOFF ON ALL PAIR MASSES + QUSQ=DOT(Q,5,5) + OX1=2*DOT(Q,2,3)/QUSQ + OX2=2*DOT(Q,1,3)/QUSQ + IF (OX1.LT.CUTOFF.OR.OX2.LT.CUTOFF.OR.1-OX1-OX2.LT.CUTOFF) + $ RETURN 1 + IF (OX1.LT.CUTUP.OR.OX2.LT.CUTUP.OR.1-OX1-OX2.LT.CUTUP) + $ KEEP=.TRUE. +C---CALCULATE THE CORRESPONDING JACOBIAN FACTOR + JAC=0 + DO M=1,NPERM3 + CALL SUBTHR(-M,Q,QTMP,STMP,JTMP,*999) + JAC=JAC+JTMP + ENDDO + IF (ABS(PERM).LE.8) THEN + JAC=JAC/(2*NPOW2**2*(1-ZI)*Y**XPOW2) + ELSE + JAC=JAC/(4*NPOW2**2*(1-ZI)*(1-ZJ)*Y**XPOW2)*(2-ZI-ZJ) + ENDIF + JAC=JAC*QUSQ/EMSQ +C---INCLUDE A PRIORI CHANNEL WEIGHTS + JAC=JAC/12 + IF (ABS(PERM).GT.8) JAC=JAC*2 + IF (PERM.LT.0) RETURN +C---CALCULATE WEIGHT FOR QUARK-GLUON SPLITTING FUNCTION + IF (PERM.LE.8) THEN + CALL MATTHR(Q,S1) + S1=S1*(2/(1-ZTI*(1-Y))-(1+ZTI)) + IF (KF(PERM).EQ.3) THEN + S1=S1*HF*CA + ELSE + S1=S1*(CF-HF*CA) + ENDIF +C---THEN PERMUTE Q + IF (PERM.LE.4) THEN + DO M=1,4 + T=Q(M,1) + Q(M,1)=Q(M,2) + Q(M,2)=Q(M,3) + Q(M,3)=T + ENDDO + ELSE + DO M=1,4 + T=Q(M,3) + Q(M,3)=Q(M,2) + Q(M,2)=T + ENDDO + ENDIF + ELSE + S1=0 + ENDIF +C---CALCULATE WEIGHT FOR GLUON DECAY + DEN=1/(Y*EMSQ) + CALL MATTHR(Q,S2) + DO M=1,4 + V(M)=ZTI*P(M,I)-ZTJ*P(M,J) + ENDDO + VV=-ZTI*ZTJ*Y*EMSQ + CALL CONTHR(Q,V,VV,S3) + S3=4*ZTI*ZTJ*S3 + IF (PERM.LE.8) THEN + S2=S2-S3 + IF (MOD(PERM-1,4).LE.1) THEN + S2=S2*HF*TR*(NF-HF)*2 + ELSE + S2=S2*HF*TR + ENDIF + ELSE + S2=S2*(2/(1-ZTI*(1-Y))+2/(1-ZTJ*(1-Y))-4)+S3 + S2=S2*HF*CA + ENDIF +C---INCLUDE SYMMETRY FACTORS + S1=S1/2 + IF (PERM.LE.8) THEN + S2=S2/4 + ELSE + S2=S2/2 + ENDIF +C---INCLUDE OVERALL FACTOR + S=(S1+S2)*16*PISQ*DEN + IF (PERM.LE.8) THEN + CASUB=S1/(S1+S2) + TFSUB=(1-CASUB)/(TR*CF)*ONF + IF (KF(PERM).EQ.3) THEN + CFSUB=0 + CASUB=CASUB/(CA*CF) + ELSE + CFSUB=CASUB/((CF-CA/2)*CF) + CASUB=CASUB/((CA-2*CF)*CF) + ENDIF + ELSE + CASUB=1/(CA*CF) + CFSUB=0 + TFSUB=0 + ENDIF + RETURN + 999 RETURN 1 + END +C----------------------------------------------------------------------- + FUNCTION ERTA(P,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE ERT A FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION ERTA,P(4,7),DOT, + $ S12,S13,S14,S23,S24,S34,S134,S234,S + S12=2*DOT(P,I,J) + S13=2*DOT(P,I,K) + S14=2*DOT(P,I,L) + S23=2*DOT(P,J,K) + S24=2*DOT(P,J,L) + S34=2*DOT(P,K,L) + S134=S13+S14+S34 + S234=S23+S24+S34 + S=S12+S13+S14+S23+S24+S34 + ERTA=(S12*S34**2-S13*S24*S34+S14*S23*S34+3*S12*S23*S34+ + $ 3*S12*S14*S34+4*S12**2*S34-S13*S23*S24+2*S12*S23*S24- + $ S13*S14*S24-2*S12*S13*S24+2*S12**2*S24+S14*S23**2+ + $ 2*S12*S23**2+S14**2*S23+4*S12*S14*S23+4*S12**2*S23+ + $ 2*S12*S14**2+2*S12*S13*S14+4*S12**2*S14+2*S12**2*S13+ + $ 2*S12**3)/(2*S13*S134*S234*S24)+ + $ (S24*S34+S12*S34+S13*S24-S14*S23+S12*S13)/(S13*S134**2)+ + $ 2*S23*(S-S13)/(S13*S134*S24)+ + $ S34/(2*S13*S24) + END +C----------------------------------------------------------------------- + FUNCTION ERTB(P,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE ERT B FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION ERTB,P(4,7),DOT, + $ S12,S13,S14,S23,S24,S34,S123,S124,S134,S234,S + S12=2*DOT(P,I,J) + S13=2*DOT(P,I,K) + S14=2*DOT(P,I,L) + S23=2*DOT(P,J,K) + S24=2*DOT(P,J,L) + S34=2*DOT(P,K,L) + S123=S12+S13+S23 + S124=S12+S14+S24 + S134=S13+S14+S34 + S234=S23+S24+S34 + S=S12+S13+S14+S23+S24+S34 + ERTB=(S12*S24*S34+S12*S14*S34-S13*S24**2+S13*S14*S24+ + $ 2*S12*S14*S24)/(S13*S134*S23*S14)+ + $ S12*(S+S34)*S124/(S134*S234*S14*S24)- + $ (2*S13*S24+S14**2+S13*S23+2*S12*S13)/(S13*S134*S14)+ + $ S12*S123*S124/(2*S13*S14*S23*S24) + END +C----------------------------------------------------------------------- + FUNCTION ERTC(P,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE ERT C FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION ERTC,P(4,7),DOT, + $ S12,S13,S14,S23,S24,S34,S134,S234 + S12=2*DOT(P,I,J) + S13=2*DOT(P,I,K) + S14=2*DOT(P,I,L) + S23=2*DOT(P,J,K) + S24=2*DOT(P,J,L) + S34=2*DOT(P,K,L) + S134=S13+S14+S34 + S234=S23+S24+S34 + ERTC=-(5*S12*S34**2+2*S12*S24*S34+2*S12*S23*S34+2*S12*S14*S34+ + $ 2*S12*S13*S34+4*S12**2*S34-S13*S24**2+S14*S23*S24+ + $ S13*S23*S24+S13*S14*S24-S12*S14*S24-S13**2*S24- + $ 3*S12*S13*S24-S14*S23**2-S14**2*S23+S13*S14*S23- + $ 3*S12*S14*S23-S12*S13*S23)/(4*S134*S234*S34**2)+ + $ (3*S12*S34**2-3*S13*S24*S34+3*S12*S24*S34+3*S14*S23*S34- + $ S13*S24**2-S12*S23*S34+6*S12*S14*S34+2*S12*S13*S34- + $ 2*S12**2*S34+S14*S23*S24-3*S13*S23*S24-2*S13*S14*S24+ + $ 4*S12*S14*S24+2*S12*S13*S24+3*S14*S23**2+2*S14**2*S23+ + $ 2*S14**2*S12+2*S12**2*S14+6*S12*S14*S23-2*S12*S13**2- + $ 2*S12**2*S13)/(4*S13*S134*S234*S34)+ + $ (2*S12*S34**2-2*S13*S24*S34+S12*S24*S34+4*S13*S23*S34+ + $ 4*S12*S14*S34+2*S12*S13*S34+2*S12**2*S34-S13*S24**2+ + $ 3*S14*S23*S24+4*S13*S23*S24-2*S13*S14*S24+4*S12*S14*S24+ + $ 2*S12*S13*S24+2*S14*S23**2+4*S13*S23**2+2*S13*S14*S23+ + $ 2*S12*S14*S23+4*S12*S13*S23+2*S12*S14**2+4*S12**2*S13+ + $ 4*S12*S13*S14+2*S12**2*S14)/(4*S13*S134*S24*S34) + ERTC=ERTC- + $ (S12*S34**2-2*S14*S24*S34-2*S13*S24*S34-S14*S23*S34+ + $ S13*S23*S34+S12*S14*S34+2*S12*S13*S34-2*S14**2*S24- + $ 4*S13*S14*S24-4*S13**2*S24-S14**2*S23-S13**2*S23+ + $ S12*S13*S14-S12*S13**2)/(2*S13*S34*S134**2)+ + $ (S12*S34**2-4*S14*S24*S34-2*S13*S24*S34-2*S14*S23*S34- + $ 4*S13*S23*S34-4*S12*S14*S34-4*S12*S13*S34-2*S13*S14*S24+ + $ 2*S13**2*S24+2*S14**2*S23-2*S13*S14*S23-S12*S14**2- + $ 6*S12*S13*S14-S12*S13**2)/(4*S34**2*S134**2) + END +C----------------------------------------------------------------------- + FUNCTION ERTD(P,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE ERT D FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION ERTD,P(4,7),DOT, + $ S12,S13,S14,S23,S24,S34,S123,S134 + S12=2*DOT(P,I,J) + S13=2*DOT(P,I,K) + S14=2*DOT(P,I,L) + S23=2*DOT(P,J,K) + S24=2*DOT(P,J,L) + S34=2*DOT(P,K,L) + S123=S12+S13+S23 + S134=S13+S14+S34 + ERTD=(S13*S23*S34+S12*S23*S34-S12**2*S34+S13*S23*S24+ + $ 2*S12*S23*S24-S14*S23**2+S12*S13*S24+S12*S14*S23+ + $ S12*S13*S14)/(S13**2*S123**2)- + $ (S12*S34**2-S13*S24*S34+S12*S24*S34-S14*S23*S34- + $ S12*S23*S34-S13*S24**2+S14*S23*S24-S13*S23*S24- + $ S13**2*S24+S14*S23**2)/(S13**2*S123*S134) + END +C----------------------------------------------------------------------- + FUNCTION ERTE(P,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE ERT E FUNCTION WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION ERTE,P(4,7),DOT, + $ S12,S13,S14,S23,S24,S34,S123,S124,S134,S234 + S12=2*DOT(P,I,J) + S13=2*DOT(P,I,K) + S14=2*DOT(P,I,L) + S23=2*DOT(P,J,K) + S24=2*DOT(P,J,L) + S34=2*DOT(P,K,L) + S123=S12+S13+S23 + S124=S12+S14+S24 + S134=S13+S14+S34 + S234=S23+S24+S34 + ERTE=(S12*S23*S34-S12*S24*S34+S12*S14*S34+S12*S13*S34+S13*S24**2- + $ S14*S23*S24+S13*S23*S24+S13*S14*S24+S13**2*S24-S14*S23**2- + $ S14**2*S23-S13*S14*S23)/(S13*S23*S123*S134)- + $ S12*(S12*S34-S23*S24-S13*S24-S14*S23- + $ S14*S13)/(S13*S23*S123**2)- + $ (S14+S13)*(S24+S23)*S34/(S13*S23*S134*S234) + END +C----------------------------------------------------------------------- + FUNCTION LEIA(P,Q,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU) +C WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION LEIA,P(4,7),Q(4),DOT, + $ P12,P13,P14,P23,P24,P34,P15,P25,P35,Q2, + $ P1K,P2K,P3K,P4K,KK + P12=2*DOT(P,I,J) + P13=2*DOT(P,I,K) + P14=2*DOT(P,I,L) + P23=2*DOT(P,J,K) + P24=2*DOT(P,J,L) + P34=2*DOT(P,K,L) + P15=P23+P24+P34 + P25=P13+P14+P34 + P35=P12+P14+P24 + Q2=P12+P13+P14+P23+P24+P34 + P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1) + P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1) + P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1) + P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1) + KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1) + LEIA = 0 + LEIA = LEIA + P1K*P2K * ( - 2*Q2*P12*P25 - Q2*P14* + + P25 - Q2*P23*P25 + P12*P13*P25 + P12*P24*P25 + 2*P12*P25* + + P34 - P13*P23*P25 - 2*P13*P24*P25 - P13*P25*P34 + 2*P14*P15* + + P24 - 2*P14*P23*P25 - P14*P24*P25 - 2*P15*P23*P25 - 2*P15*P24 + + *P25 - 2*P15*P25*P34 - P24*P25*P34 ) + LEIA = LEIA + P1K*P3K * ( - Q2*P12 + Q2*P24 + P12* + + P13 + P14*P24 + 2*P15*P24 - P24*P34 )*P25 + LEIA = LEIA + P1K*P4K * ( - Q2*P12 - Q2*P23 + P12* + + P13 + P12*P34 - P14*P23 )*P25 + LEIA = LEIA + P1K**2 * ( Q2*P23 + Q2*P24 - P13*P23 + + - P13*P24 - P24*P34 )*P25 + LEIA = LEIA + P2K*P3K * ( - Q2*P12*P25 - Q2*P14* + + P25 - 2*P12*P15*P25 + P12*P24*P25 + P12*P25*P34 + 2*P14*P15* + + P24 - P14*P23*P25 - 4*P15*P23*P25 - 2*P15*P24*P25 - 2*P15*P25 + + *P34 ) + LEIA = LEIA + P2K*P4K * ( - Q2*P12*P25 + Q2*P13* + + P25 + P12*P24*P25 + 2*P13*P15*P25 + P13*P23*P25 - P13*P25*P34 + + + 2*P14*P15*P24 - 2*P15*P23*P25 - 2*P15*P24*P25 ) + LEIA = LEIA + P2K**2 * ( Q2*P13 + Q2*P14 + 2*P13*P15 + + - P13*P24 - P13*P34 - P14*P24 + 2*P15*P34 )*P25 + LEIA = LEIA + P3K*P4K * ( - 2*P12*P15 + P12*P34 - P13* + + P24 - P14*P23 - 2*P15*P23 - P15*P25 )*P25 + LEIA = LEIA + P3K**2 * ( P14 + 2*P15 )*P24*P25 + LEIA = LEIA + P4K**2 * ( P13*P23*P25 ) +c$$$ LEIA = LEIA + KK * ( 2*Q2*P12*P14*P25 + 2*Q2*P12*P23 +c$$$ + *P25 + 2*Q2*P12**2*P25 - 2*Q2*P14*P15*P24 + 2*Q2*P14* +c$$$ + P23*P25 + P12*P13*P25*P34 + 4*P12*P15*P23*P25 + 2*P12*P15*P25 +c$$$ + *P34 + P12*P24*P25*P34 - P13*P14*P23*P25 - 2*P13*P14*P24*P25 +c$$$ + - 2*P13*P15*P24*P25 - 2*P13*P23*P24*P25 - P13*P24**2*P25 - +c$$$ + P13**2*P24*P25 + 2*P14*P15*P23*P25 - P14*P23*P24*P25 + 4*P15* +c$$$ + P23*P24*P25 + 4*P15*P23*P25*P34 + 4*P15*P23**2*P25 + 2*P15* +c$$$ + P24*P25*P35 + P15*P25**2*P34 )/4 + LEIA = LEIA/(P15*P25**2*P13*P24) + END +C----------------------------------------------------------------------- + FUNCTION LEIB(P,Q,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU) +C WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION LEIB,P(4,7),Q(4),DOT, + $ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2, + $ P1K,P2K,P3K,P4K,KK + P12=2*DOT(P,I,J) + P13=2*DOT(P,I,K) + P14=2*DOT(P,I,L) + P23=2*DOT(P,J,K) + P24=2*DOT(P,J,L) + P34=2*DOT(P,K,L) + P15=P23+P24+P34 + P25=P13+P14+P34 + P35=P12+P14+P24 + P45=P12+P13+P23 + Q2=P12+P13+P14+P23+P24+P34 + P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1) + P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1) + P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1) + P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1) + KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1) + LEIB = 0 + LEIB = LEIB + P1K*P2K * ( - 4*Q2*P14*P24*P45 - 2*P12 + + *P15*P25*P34 - P12*P15*P25*P35 - P12*P15*P25*P45 + 8*P13*P15* + + P23*P24 + 4*P13*P15*P24*P34 + 4*P13*P15*P24**2 + 4*P14*P15* + + P23*P24 + 2*P14*P15*P24*P45 + 2*P14*P24*P25*P45 - 4*P15*P23* + + P24*P25 + 4*P15*P23*P24*P34 ) + LEIB = LEIB + P1K*P3K * ( - P12*P14*P15*P25 - 2*P12* + + P14*P24*P45 + P12*P15*P24*P25 + 4*P12*P15*P24*P34 - P12**2* + + P15*P25 - 4*P13*P15*P24**2 + 4*P14*P15*P23*P24 ) + LEIB = LEIB + P1K*P4K * ( - P12*P13*P15*P25 - 2*P12* + + P14*P24*P45 + P12*P15*P23*P25 - P12**2*P15*P25 + 4*P13*P15* + + P23*P24 + 4*P14*P15*P23*P24 ) + LEIB = LEIB + P1K**2 * ( 2*Q2*P14*P24*P45 - 2*P12*P14* + + P24*P45 + P12*P15**2*P25 - 2*P14*P24*P25*P45 - 4*P15*P23*P24* + + P34 ) + LEIB = LEIB + P2K*P3K * ( 4*P12*P13*P15*P24 + P12*P14* + + P15*P25 - 2*P12*P14*P24*P45 - P12*P15*P24*P25 - P12**2*P15* + + P25 + 4*P13*P15*P23*P24 + 4*P13*P15*P24**2 ) + LEIB = LEIB + P2K*P4K * ( - 4*P12*P13*P15*P24 + P12* + + P13*P15*P25 - 8*P12*P14*P15*P24 - 2*P12*P14*P24*P45 - P12*P15 + + *P23*P25 - 4*P12*P15*P24*P34 - P12**2*P15*P25 - 4*P13*P14*P15 + + *P24 + 4*P13*P15*P23*P24 + 4*P13*P15*P24**2 - 4*P13**2*P15* + + P24 ) + LEIB = LEIB + P2K**2 * ( 2*Q2*P14*P24*P45 - 2*P12*P14* + + P24*P45 + P12*P15*P25**2 - 4*P13*P14*P15*P24 - 4*P13*P15*P24* + + P34 - 4*P13**2*P15*P24 - 2*P14*P15*P24*P45 ) + LEIB = LEIB + P3K*P4K * ( - 2*P12*P25 - 4*P14*P24 ) + + *P12*P15 + LEIB = LEIB + P3K**2 * ( - 4*P12*P14*P15*P24 ) +c$$$ LEIB = LEIB + KK * ( 2*Q2*P12*P14*P24*P45 + Q2* +c$$$ + P12**2*P15*P25 - 2*Q2*P13*P15*P23*P24 - 2*Q2*P14*P15*P23* +c$$$ + P24 + 2*Q2*P15*P23*P24*P25 + P12*P13*P14*P15*P25 - 4*P12* +c$$$ + P13*P15*P23*P24 - 2*P12*P13*P15*P24*P34 + 2*P12*P14*P15*P24* +c$$$ + P34 + 4*P12*P14*P15*P24**2 + P12*P15*P23*P24*P25 - 2*P12*P15* +c$$$ + P23*P24*P34 + 2*P12*P15*P24**2*P34 - 2*P13*P14*P15*P23*P24 + +c$$$ + 2*P13*P14*P15*P24**2 - 2*P13*P15*P23*P24**2 - 2*P13*P15* +c$$$ + P24**3 + 2*P13**2*P15*P24**2 + 2*P14*P15*P23*P24**2 + 2*P14* +c$$$ + P15*P23**2*P24 - 2*P14**2*P15*P23*P24 - 2*P15**2*P23*P24*P25 +c$$$ + )/2 + LEIB = LEIB/(2*P13*P14*P15*P23*P24*P25) + END +C----------------------------------------------------------------------- + FUNCTION LEIC(P,Q,I,J,K,L) + IMPLICIT NONE +C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU) +C WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L + DOUBLE PRECISION LEIC,P(4,7),Q(4),DOT, + $ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2, + $ P1K,P2K,P3K,P4K,KK + P12=2*DOT(P,I,J) + P13=2*DOT(P,I,K) + P14=2*DOT(P,I,L) + P23=2*DOT(P,J,K) + P24=2*DOT(P,J,L) + P34=2*DOT(P,K,L) + P15=P23+P24+P34 + P25=P13+P14+P34 + P35=P12+P14+P24 + P45=P12+P13+P23 + Q2=P12+P13+P14+P23+P24+P34 + P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1) + P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1) + P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1) + P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1) + KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1) + LEIC = 0 + LEIC = LEIC + P1K*P2K * ( - 4*P12*P13*P15*P25*P34 + 6* + + P12*P13*P24*P25*P34 - 2*P12*P14*P15*P25*P34 - 2*P12*P14*P24* + + P25*P34 - 2*P12*P15*P25*P34**2 + 2*P12*P24*P25*P34**2 + 4*P13 + + *P14*P15*P24*P34 - 2*P13*P14*P15*P25*P34 - 3*P13*P14*P23*P24* + + P25 - P13*P14*P24**2*P25 + P13*P14**2*P15*P24 - 2*P13*P15*P23 + + *P25*P34 + P13*P15*P24*P25*P34 - P13*P15*P24*P34**2 - P13*P15 + + *P25*P34**2 + P13*P23*P24*P25*P34 - 2*P13*P24**2*P25*P34 + 6* + + P13**2*P14*P15*P24 - P13**2*P23*P24*P25 + P13**2*P24*P25*P34 + + - 3*P13**2*P24**2*P25 + P13**3*P15*P24 - P14*P15*P23*P25*P34 + + - 2*P14*P15*P24*P25*P34 - 2*P14*P15*P25*P34**2 - 3*P14*P23* + + P24*P25*P34 + P14*P24*P25*P34**2 - 2*P14*P24**2*P25*P34 - + + P14**2*P15*P25*P34 - P14**2*P24*P25*P34 + P15*P23*P25*P34**2 + + + P15*P24*P25*P34**2 + P15*P25*P34**3 + P24*P25*P34**3 - + + P24**2*P25*P34**2 ) + LEIC = LEIC + P1K*P3K * ( 3*P12*P13*P34 - P12*P14*P34 + + - P12*P34**2 + P13*P14*P24 + 2*P13*P15*P34 - P13**2*P24 - + + P14*P15*P34 + 3*P14*P24*P34 + P15*P34**2 - 3*P24*P34**2 ) + + *P24*P25 + LEIC = LEIC + P1K*P4K * ( - 2*P12*P13*P15*P34 + 3*P12* + + P13*P24*P34 - P12*P14*P15*P34 - P12*P14*P24*P34 - 3*P12*P15* + + P34**2 + 3*P12*P24*P34**2 - P13*P14*P23*P24 - P13*P15*P23*P34 + + - 2*P13*P15*P24*P34 - 2*P13*P23*P24*P34 + P13**2*P23*P24 + 2 + + *P14*P15*P23*P34 - 3*P14*P23*P24*P34 + P23*P24*P34**2 )*P25 + LEIC = LEIC + P1K**2 * ( 2*P13*P15 - 3*P13*P23 - 3*P13* + + P24 + P14*P15 + P14*P23 + P14*P24 + 3*P15*P34 + P23*P34 - 3* + + P24*P34 )*P24*P25*P34 + LEIC = LEIC + P2K*P3K * ( - 2*P12*P13*P15*P25*P34 + + + P12*P13*P24*P25*P34 - P12*P14*P15*P25*P34 - 3*P12*P14*P24*P25 + + *P34 - P12*P15*P25*P34**2 + P12*P24*P25*P34**2 + 2*P13*P14* + + P15*P24*P34 - P13*P14*P15*P25*P34 - P13*P14*P23*P24*P25 - 2* + + P13*P14*P24*P25*P34 + P13*P14*P24**2*P25 - 2*P13*P14**2*P15* + + P24 - 4*P13*P15*P23*P25*P34 - 2*P13*P15*P25*P34**2 + 2*P13**2 + + *P14*P15*P24 - 2*P14*P15*P23*P25*P34 - 3*P14*P15*P24*P25*P34 + + - 2*P14*P15*P25*P34**2 - 3*P14*P23*P24*P25*P34 - P14*P24**2* + + P25*P34 - P14**2*P24*P25*P34 ) + LEIC = LEIC + P2K*P4K * ( - P12*P13*P15*P25*P34 + 2* + + P12*P13*P24*P25*P34 - 2*P12*P14*P15*P25*P34 - 2*P12*P14*P24* + + P25*P34 - 2*P12*P24*P25*P34**2 + 4*P13*P14*P15*P24*P34 - P13* + + P14*P15*P25*P34 + P13*P14*P24*P25*P34 - 2*P13*P15*P23*P25*P34 + + - 3*P13*P15*P24*P25*P34 - P13*P15*P25*P34**2 + 3*P13*P23*P24 + + *P25*P34 - 2*P13*P24*P25*P34**2 + P13*P24**2*P25*P34 + 2* + + P13**2*P14*P15*P24 - 2*P13**2*P15*P24*P34 + 2*P13**2*P15*P25* + + P34 + P13**2*P23*P24*P25 - P13**2*P24**2*P25 - 2*P13**3*P15* + + P24 - 2*P14*P15*P23*P25*P34 - 4*P14*P15*P24*P25*P34 ) + LEIC = LEIC + P2K**2 * ( 2*P13*P14*P15 + P13*P14*P24 + + + P13*P15*P34 - P13*P24*P34 + 2*P13**2*P15 - P13**2*P24 + 2*P14 + + *P15*P34 + 2*P14*P24*P34 + 2*P14**2*P15 + 2*P14**2*P24 ) + + *P25*P34 + LEIC = LEIC + P3K*P4K * ( - P12*P13*P15*P34 + 2*P12* + + P13*P24*P34 - P12*P14*P15*P34 - P12*P15*P34**2 + 2*P12*P24* + + P34**2 - 2*P13*P14*P23*P24 - 2*P13*P15*P23*P34 - 2*P13*P24**2 + + *P34 - 2*P13**2*P24**2 - 2*P14*P23*P24*P34 )*P25 + LEIC = LEIC + P3K**2 * ( 2*P13*P14*P24 + 2*P13*P15*P34 + + + 2*P14*P24*P34 )*P24*P25 + LEIC = LEIC + P4K**2 * ( 2*P12*P15*P34 + 2*P13*P23*P24 + + + 2*P23*P24*P34 )*P13*P25 +c$$$ LEIC = LEIC + KK * ( - 4*P12*P13*P14*P15*P24*P34 + 4* +c$$$ + P12*P13*P14*P15*P25*P34 + 3*P12*P13*P14*P23*P24*P25 - 2*P12* +c$$$ + P13*P14*P24*P25*P34 + P12*P13*P14*P24**2*P25 - P12*P13*P14**2 +c$$$ + *P15*P24 + 4*P12*P13*P15*P23*P25*P34 + P12*P13*P15*P24*P34**2 +c$$$ + + 2*P12*P13*P15*P25*P34**2 - 2*P12*P13*P23*P24*P25*P34 - 2* +c$$$ + P12*P13*P24*P25*P34**2 - 6*P12*P13**2*P14*P15*P24 + P12* +c$$$ + P13**2*P23*P24*P25 - 4*P12*P13**2*P24*P25*P34 + 3*P12*P13**2* +c$$$ + P24**2*P25 - P12*P13**3*P15*P24 + 2*P12*P14*P15*P23*P25*P34 +c$$$ + + 4*P12*P14*P15*P24*P25*P34 + 5*P12*P14*P15*P25*P34**2 + 6* +c$$$ + P12*P14*P23*P24*P25*P34 + P12*P14*P24*P25*P34**2 + 4*P12*P14* +c$$$ + P24**2*P25*P34 + 2*P12*P14**2*P15*P25*P34 + 2*P12*P14**2*P24* +c$$$ + P25*P34 - P12*P15*P24*P25*P34**2 + P12*P15*P25*P34**3 - P12* +c$$$ + P23*P24*P25*P34**2 + 3*P12*P24**2*P25*P34**2 + 4*P12**2*P13* +c$$$ + P15*P25*P34 - 6*P12**2*P13*P24*P25*P34 + 2*P12**2*P14*P15*P25 +c$$$ + *P34 + 2*P12**2*P14*P24*P25*P34 + 2*P12**2*P15*P25*P34**2 - 2 +c$$$ + *P12**2*P24*P25*P34**2 )/4 +c$$$ LEIC = LEIC + KK * ( - 2*P13*P14*P15*P23*P24*P34 + 2* +c$$$ + P13*P14*P15*P23*P25*P34 + P13*P14*P15*P24*P25*P34 - 4*P13*P14 +c$$$ + *P15*P24**2*P34 + 3*P13*P14*P23*P24*P25*P34 - P13*P14*P23* +c$$$ + P24**2*P25 + P13*P14*P23**2*P24*P25 - 3*P13*P14*P24**2*P25* +c$$$ + P34 + 2*P13*P14**2*P15*P23*P24 + P13*P14**2*P23*P24*P25 + 2* +c$$$ + P13*P15*P23*P24*P25*P34 + 4*P13*P15*P23*P25*P34**2 + 4*P13* +c$$$ + P15*P23**2*P25*P34 - P13*P15*P24*P25*P34**2 + 3*P13*P15* +c$$$ + P24**2*P25*P34 - 3*P13*P23*P24**2*P25*P34 + 2*P13*P24**2*P25* +c$$$ + P34**2 - P13*P24**3*P25*P34 - 2*P13**2*P14*P15*P23*P24 - 2* +c$$$ + P13**2*P14*P15*P24**2 - P13**2*P14*P23*P24*P25 - P13**2*P14* +c$$$ + P24**2*P25 - 2*P13**2*P15*P24*P25*P34 + 2*P13**2*P15*P24**2* +c$$$ + P34 - P13**2*P23*P24**2*P25 + P13**2*P24**2*P25*P34 + P13**2* +c$$$ + P24**3*P25 + 2*P13**3*P15*P24**2 + P13**3*P24**2*P25 + 5*P14* +c$$$ + P15*P23*P24*P25*P34 + P14*P15*P23*P25*P34**2 + 2*P14*P15* +c$$$ + P23**2*P25*P34 + 4*P14*P15*P24**2*P25*P34 + 2*P14*P23*P24*P25 +c$$$ + *P34**2 + )/4 +c$$$ LEIC = LEIC + KK * ( P14*P23*P24**2*P25*P34 + 3*P14* +c$$$ + P23**2*P24*P25*P34 - P14**2*P15*P23*P25*P34 + 3*P14**2*P23* +c$$$ + P24*P25*P34 )/4 + LEIC = LEIC/(2*P15*P25**2*P13*P34**2*P24) + END +C----------------------------------------------------------------------- + FUNCTION LEID(P,Q,II,JJ,KKK,LL) + IMPLICIT NONE +C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU) +C WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L,II,JJ,KKK,LL + DOUBLE PRECISION LEID,P(4,7),Q(4),DOT, + $ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2, + $ P1K,P2K,P3K,P4K,KK +C---CONVERT FROM ERT NOTATION (Q'QQ'BARQBAR) TO LEIDEN (QQBARQ'BARQ') + I=JJ + J=LL + K=KKK + L=II + P12=2*DOT(P,I,J) + P13=2*DOT(P,I,K) + P14=2*DOT(P,I,L) + P23=2*DOT(P,J,K) + P24=2*DOT(P,J,L) + P34=2*DOT(P,K,L) + P15=P23+P24+P34 + P25=P13+P14+P34 + P35=P12+P14+P24 + P45=P12+P13+P23 + Q2=P12+P13+P14+P23+P24+P34 + P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1) + P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1) + P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1) + P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1) + KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1) + LEID = 0 + LEID = LEID + P1K*P2K * ( - 2*P12*P25*P34 - 4*P13*P14* + + P15 - 2*P13*P15*P34 + 2*P13*P24*P25 - 2*P14*P15*P34 + 2*P14* + + P23*P25 )/P25 + LEID = LEID + P1K*P3K * ( - P12*P34 - 2*P14*P24 + P24*P25 ) + LEID = LEID + P1K*P4K * ( - P12*P34 - 2*P13*P23 + P23*P25 ) + LEID = LEID + P1K**2 * ( P23 + P24 )*P34 + LEID = LEID + P2K*P3K * ( - P12*P25*P34 - 4*P13*P14*P15 + + - 2*P13*P15*P34 + 3*P14*P15*P25 - 2*P14*P15*P34 - 2*P14*P24* + + P25 )/P25 + LEID = LEID + P2K*P4K * ( - P12*P25*P34 - 4*P13*P14*P15 + + + 3*P13*P15*P25 - 2*P13*P15*P34 - 2*P13*P23*P25 - 2*P14*P15* + + P34 )/P25 + LEID = LEID + P2K**2 * ( P13 + P14 )*P34 + LEID = LEID + P3K*P4K * ( - 2*P12*P34 + 2*P13*P24 + 2* + + P14*P23 ) + LEID = LEID + P3K**2 * ( - 2*P14*P24 ) + LEID = LEID + P4K**2 * ( - 2*P13*P23 ) +c$$$ LEID = LEID + KK * ( Q2*P12*P25*P34 - Q2*P13*P24*P25 +c$$$ + - Q2*P14*P23*P25 + 2*P12*P13*P15*P25 - P12*P13*P24*P25 - 2 +c$$$ + *P12*P13**2*P15 + 2*P12*P14*P15*P25 - P12*P14*P23*P25 - 2*P12 +c$$$ + *P14**2*P15 + P12*P25*P34**2 + P12**2*P25*P34 + 2*P13*P14*P23 +c$$$ + *P25 + 2*P13*P14*P24*P25 + 2*P13*P15*P23*P25 + 2*P13*P23*P24* +c$$$ + P25 - P13*P24*P25*P34 - 2*P13**2*P15*P23 - 2*P13**2*P15*P24 +c$$$ + + 2*P14*P15*P24*P25 + 2*P14*P23*P24*P25 - P14*P23*P25*P34 - +c$$$ + 2*P14**2*P15*P23 - 2*P14**2*P15*P24 )/(4*P25) + LEID = LEID/(P15*P25*P34**2) + END +C----------------------------------------------------------------------- + FUNCTION LEIE(P,Q,II,JJ,KKK,LL) + IMPLICIT NONE +C---EVALUATE THE MATRIX ELEMENT CONTRACTED WITH A TENSOR -Q(MU)Q(NU) +C WITH P(*,I)=P1, P(*,J)=P2, ETC + INTEGER I,J,K,L,II,JJ,KKK,LL + DOUBLE PRECISION LEIE,P(4,7),Q(4),DOT, + $ P12,P13,P14,P23,P24,P34,P15,P25,P35,P45,Q2, + $ P1K,P2K,P3K,P4K,KK +C---CONVERT FROM ERT NOTATION (Q'QQ'BARQBAR) TO LEIDEN (QQBARQ'BARQ') + I=JJ + J=LL + K=KKK + L=II + P12=2*DOT(P,I,J) + P13=2*DOT(P,I,K) + P14=2*DOT(P,I,L) + P23=2*DOT(P,J,K) + P24=2*DOT(P,J,L) + P34=2*DOT(P,K,L) + P15=P23+P24+P34 + P25=P13+P14+P34 + P35=P12+P14+P24 + P45=P12+P13+P23 + Q2=P12+P13+P14+P23+P24+P34 + P1K=P(4,I)*Q(4)-P(3,I)*Q(3)-P(2,I)*Q(2)-P(1,I)*Q(1) + P2K=P(4,J)*Q(4)-P(3,J)*Q(3)-P(2,J)*Q(2)-P(1,J)*Q(1) + P3K=P(4,K)*Q(4)-P(3,K)*Q(3)-P(2,K)*Q(2)-P(1,K)*Q(1) + P4K=P(4,L)*Q(4)-P(3,L)*Q(3)-P(2,L)*Q(2)-P(1,L)*Q(1) + KK=Q(4)*Q(4)-Q(3)*Q(3)-Q(2)*Q(2)-Q(1)*Q(1) + LEIE = 0 + LEIE = LEIE + P1K*P2K * ( - 2*Q2*P13*P15*P25 + Q2* + + P23*P25**2 + 2*P13*P14*P15*P25 - 4*P14*P15*P25*P45 - P14*P23* + + P25**2 + 2*P14**2*P15*P45 - 2*P15*P23*P25*P34 + 2*P15*P25**2* + + P45 - P23*P25**2*P45 )/(P13*P25*P45) + LEIE = LEIE + P1K*P3K * ( - 2*Q2*P13*P15 + Q2*P23* + + P25 - 2*P12*P14*P15 - P14*P23*P25 - 2*P15*P23*P34 + 2*P15*P25 + + *P45 - P23*P25*P45 )/(P13*P45) + LEIE = LEIE + P1K*P4K * ( 2*Q2*P25 + 2*P13*P15 - 2*P14* + + P25 - 2*P15*P25 - 2*P15*P34 - 2*P25*P45 )*P23/(P13*P45) + LEIE = LEIE + P1K**2 * ( - 2*P15*P23*P34 )/(P13*P45) + LEIE = LEIE + P2K*P3K * ( - 2*P12*P15*P25 + 2*P13*P15* + + P25 + 2*P14*P15*P45 - 2*P23*P25**2 )*P14/(P13*P25*P45) + LEIE = LEIE + P2K*P4K * ( Q2*P23*P25**2 + 2*P12*P15*P25 + + *P34 - 2*P13*P15*P24*P25 + 2*P14*P15*P23*P25 - 2*P14*P15*P25* + + P45 - P14*P23*P25**2 + 2*P14**2*P15*P45 - P15*P23*P25**2 ) + + /(P13*P25*P45) + LEIE = LEIE + P2K**2 * ( 2*P13*P15 - P23*P25 )*P14/(P13*P45) + LEIE = LEIE + P3K*P4K * ( - 2*Q2*P13*P15 + Q2*P23* + + P25 - 2*P12*P14*P15 + 2*P13*P15*P23 - P14*P23*P25 - 3*P15*P23 + + *P25 + 2*P15*P25*P45 )/(P13*P45) + LEIE = LEIE + P3K**2 * ( - 2*P12*P15 - P23*P25 )*P14/(P13*P45) + LEIE = LEIE + P4K**2 * ( 2*P15*P23 )/P45 +c$$$ LEIE = LEIE + KK * ( Q2*P13*P15*P24*P25 + Q2*P14*P23* +c$$$ + P25**2 - Q2*P14**2*P15*P45 + Q2*P15*P23*P25**2 + Q2*P23 +c$$$ + *P25**2*P45 - Q2**2*P23*P25**2 - P12*P13*P15*P25*P34 + P12 +c$$$ + *P14*P15*P23*P25 - P12*P15*P24*P25*P34 - P12*P15*P25*P34**2 - +c$$$ + P12**2*P15*P25*P34 - P13*P14*P15*P23*P25 - P13*P14*P15*P24* +c$$$ + P25 - P13*P15*P23*P24*P25 - P14*P15*P23*P24*P25 + P14*P15*P23 +c$$$ + *P25*P34 + P14*P15*P25*P35*P45 - P15*P23*P25**2*P45 ) +c$$$ + /(2*P13*P25*P45) + LEIE = LEIE/(P15*P25*P34) + END +C----------------------------------------------------------------------- + SUBROUTINE WARNER(SUBR,TYPE,*) + IMPLICIT NONE +C---GIVE WARNING MESSAGES ABOUT NUMERICAL ERRORS + INTEGER TYPE + LOGICAL FIRST + CHARACTER*(*) SUBR + INTEGER METYPE,NF + DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF + COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CUTUP,CQ,FPAL2,ONF, + $ NF,METYPE + DATA FIRST/.TRUE./ + WRITE (*,'(A,I2,3A)') + $' Numerical error of type',TYPE,' in routine ',SUBR,'!' + IF (FIRST) THEN + WRITE (*,'(/A)') + $' This should be preventable by increasing the internal variable' + WRITE (*,'(A,1PE8.1,A)') + $' CUTOFF, whose present value is',CUTOFF,'. If this fails to' + WRITE (*,'(A)') + $' reduce the frequency of errors, something is probably more' + WRITE (*,'(A)') + $' seriously wrong, so please contact seymour@surya11.cern.ch' + WRITE (*,'(A)') + $' giving the date of this version and the latest value of ISEED.' + WRITE (*,'(A/)') ' Event generation has not been affected.' + FIRST=.FALSE. + ENDIF + RETURN 1 + END +C----------------------------------------------------------------------- +C----------------------------------------------------------------------- +C----------------------------------------------------------------------- + subroutine convert(p) + implicit none + integer i,j + double precision p(4,7),plab + common /mom/ plab(4,10) + do j=1,4 + do i=1,4 + plab(j,i+2)=p(j,i) + enddo + plab(j,1)=p(j,7) + plab(j,2)=p(j,6) + enddo + end +C----------------------------------------------------------------------- + subroutine setup + implicit double precision(a-h,o-z) + parameter(pi=3.141592653589793238D0) + common /inppar/ njets,nloop,itmax1,itmax2,nshot1,nshot2,nshot3 + common /zcup/czur,czul,czdr,czdl,czlr,czll,cznr,cznl + common /gcup/cgur,cgul,cgdr,cgdl,cglr,cgll,cgnr,cgnl + common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon + common /strong/ as,gs,renscl,qcdl,b0,nf + common /sighad/ sig0 + + common /order/iorder + common /tcuts/ ymin,dly + common /ecuts/ ycut,evis,jetalg + common /polar/ Peleft,Ppright + character*4 alg + logical done + data done/.false./ + if (done) return + done=.true. + +c---set the general parameters + rmz=91.17D0 + rgz=2.516D0 + w=0 + +c---set the polarisations + Peleft=0.5 + Ppright=0.5 + +c---and a few other parameters + ymin=1D-8 + as=2*pi + +c +c iorder determines whether we keep all terms, or leading N, subleading N +c or nf pieces +c iorder=0 all, iorder=1, leading N, iorder=2 1/N**2, iorder=3 nf +c + iorder=0 + +c +c icon determines which interference terms to use +c icon=0 all, icon=1 gives gg, icon=2 gives gz, icon=3 gives zz +c + icon=1 +c +c check that matrix elements vanish when pz is contracted in : momcon=1 +c + momcon=0 +c +c physics parameters: +c +c nutf : number of up type flavours +c ndtf : number of down type flavours + nutf=2 + ndtf=3 + nf=nutf+ndtf + if(iorder.eq.0)b0=(33D0-2D0*nf)/12D0/pi + if(iorder.eq.1)b0=33D0/12D0/pi + if(iorder.eq.2)b0=0D0 + if(iorder.eq.3)b0=-2D0*nf/12D0/pi +c +c renscl : scale at which alpha_s is evaluated = scl* rmz +c can adjust alpha_s by changing scl or qcdl below +c + scl=1D0 + renscl=scl*rmz + dly=log(ymin) +c +c qcdl : two loop lambda qcd for 4 flavours +c + qcdl=0.200D0 + alfa=alfas(renscl,qcdl,2) + gs=sqrt(4.D0*pi*as) + am=1.D0/128D0 + em=sqrt(4.D0*pi*am) + sw2=0.23D0 +c determine the width of the z + scw=sqrt(sw2)*sqrt(1D0-sw2) + cgll=-1D0 + cglr=-1D0 + cgnl=0D0 + cgnr=0D0 + cgur=+2D0/3D0 + cgdr=-1D0/3D0 + cgul=+2D0/3D0 + cgdl=-1D0/3D0 + czll=(-0.5D0+sw2)/scw + czlr=sw2/scw + cznl=0.5D0/scw + cznr=0D0 + czul=(+0.5D0-sw2*2D0/3D0)/scw + czur=(-sw2*2D0/3D0)/scw + czdl=(-0.5D0+sw2*1D0/3D0)/scw + czdr=(+sw2*1D0/3D0)/scw +c calculate the width of the z ( no top quark! ) + rgee= em**2*(czll**2+czlr**2)*rmz/24.D0/pi + rgnn= em**2*(cznl**2+cznr**2)*rmz/24.D0/pi + rguu= 3D0*em**2*(czul**2+czur**2)*rmz/24.D0/pi + rgdd= 3D0*em**2*(czdl**2+czdr**2)*rmz/24.D0/pi + rglep= 3D0*rgee + rginv= 3D0*rgnn + rghad= (1D0+(alfa/pi)+1.41*(alfa/pi)**2)*(2D0*rguu+3D0*rgdd) + rghad= (2D0*rguu+3D0*rgdd) + rgz0 = rglep + rginv + rghad + rgz1 = rglep + rginv + rghad*(1D0+(alfa/pi)) + rgz2 = rglep + rginv + rghad*(1D0+(alfa/pi)+1.41*(alfa/pi)**2) +* test=abs((rgztest-rgz)/rgz*100D0) +* write(*,40)rgz0 + 40 format(/,'lowest order z boson width = ',f6.3,' GeV ') +* write(*,41)rgz1 + 41 format(/,' first order z boson width = ',f6.3,' GeV ') +* write(*,42)rgz2 + 42 format(/,'second order z boson width = ',f6.3,' GeV ') + wtnano=389385.73D0 +c +c hadronic cross section +c + sig0=wtnano*rghad*rgee*12D0*pi/rmz**2/rgz**2 + if(jetalg.eq.1)alg=' E ' + if(jetalg.eq.2)alg=' E0 ' + if(jetalg.eq.3)alg=' P ' + if(jetalg.eq.4)alg=' P0 ' + if(jetalg.eq.5)alg=' D ' + + if(jetalg.eq.6)alg=' DE ' + + if(jetalg.eq.7)alg=' DP ' + if(jetalg.eq.8)alg=' D1 ' + write(*,*) '*************************************************', + . '*****************' + + write(*,*)' ', + + . ' prelim 13/6/94 ' + write(*,11)njets,(njets-2+nloop),w + 11 format(/,' e+e- to ',i2,' jets at O(alphas**',i1,') at ', + . ' root s = ',f7.2, ' GeV ') + if(iorder.eq.0)write(*,*)' including all colour pieces' + if(iorder.eq.1)write(*,*)' including only leading colour pieces' + if(iorder.eq.2)write(*,*)' including only subleading ', + . 'colour pieces' + if(iorder.eq.3)write(*,*)' including only nf pieces' + if(icon.eq.0)write(*,*)' including both photon and z exchange' + if(icon.eq.1)write(*,*)' including only photon exchange' + if(icon.eq.2)write(*,*)' including photon - z interference' + if(icon.eq.3)write(*,*)' including only z exchange' + if(momcon.eq.1)write(*,*)' checking momentum conservation => 0' + write(*,12)as + 12 format(/,' alphas(Mz) = ',f8.4) + write(*,13)4D0*pi/em**2,sw2 + 13 format(/,' 1/alpha = ',f7.2,'; sin**2 theta = ',f6.3) + write(*,14)sig0 + + 14 format(/,' sigma(hadronic) = ',f8.3,' nb') + write(*,15)alg,ycut + 15 format(/,' with jets defined in',a3,'-scheme with ycut = ',f8.4) + write(*,16)rmz,rgz + 16 format(/,' Z mass = ',f6.3,' GeV, Z width = ',f6.3,' GeV') + write(*,17)Peleft*100D0 + 17 format(/,' electron beam = ',f4.1,'% left handed') + write(*,18)Ppright*100D0 + 18 format(/,' positron beam = ',f4.1,'% right handed') + write(*,19)ymin + 19 format(/,' parton resolution parameter = ',g12.5) + write(*,20)itmax1,nshot1/10,nshot2/10,nshot3/10 + 20 format(/,' warm up with ',i6,' sweeps of ',3i9,' shots ') + write(*,21)itmax2,nshot1,nshot2,nshot3 + 21 format(/,' final run of ',i6,' sweeps of ',3i9,' shots ') + write(*,*) '*************************************************', + . '*****************' + + return + end +c*********************************************************************** +c +c compute matrix element weight using library +c + subroutine sfmee(npar,nloop,eetot) +c*********************************************************************** + implicit double precision (a-h,o-y) + implicit complex*16(z) + common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon + common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr + common /polar/ Peleft,Ppright + common /mom/ plab(4,10) + common /kinema/ p(4,7) + common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen + common /linpr/ rin(7,7) + common /curent/ zrcur(2,2),zlcur(2,2) + common /check/ momcon +c +c transform vectors - all momenta outgoing + do 10 mu=1,4 + do 5 n=1,npar + + p(mu,n)=plab(mu,n+2) + 5 continue + p(mu,6)=-plab(mu,1) + p(mu,7)=-plab(mu,2) + 10 continue +c +c initialise the spinor and dot products + call setspn(npar) +c all quarks have positive energy so + negaen=1 +c +c initialise the lepton current + do 20 ia=1,2 + do 20 ib=1,2 + if(momcon.eq.0)then + zrcur(ia,ib)=zkod(6,ia)*zko(7,ib) + elseif(momcon.eq.1)then + zrcur(ia,ib)=(0D0,0D0) + do 15 ipar=1,npar+1 + zrcur(ia,ib)=zrcur(ia,ib)+zkod(ipar,ia)*zko(ipar,ib) + 15 continue + endif + zlcur(ib,ia)=conjg(zrcur(ia,ib)) + 20 continue +c + if(npar.eq.2) call sfee2j(nloop,eetot) + if(npar.eq.3) call sfee3j(nloop,eetot) + if(npar.eq.4) call sfee4j(nloop,eetot) +c + if(icon.eq.1)eetot=Peleft*Ppright*eeggl + . +(1D0-Peleft)*(1D0-Ppright)*eeggr + if(icon.eq.2)eetot=Peleft*Ppright*eegzl + . +(1D0-Peleft)*(1D0-Ppright)*eegzr + if(icon.eq.3)eetot=Peleft*Ppright*eezzl + . +(1D0-Peleft)*(1D0-Ppright)*eezzr +c + return + end +*********************************************************************** +c sfee2j calculates the process e+ e- --> g,z --> q + qb +c includes: no spin averages, final state statistics, couplings +c********************************************************************** + subroutine sfee2j(nloop,eetot) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon + common /strong/ as,gs,renscl,qcdl,b0,nf + common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr + common /polar/ Peleft,Ppright + dimension rl(3,2),rr(3,2) + dimension zhl(3,2),zhr(3,2) +c + call eeq2gn(0,nloop,rl(1,1),rr(1,1),zhl,zhr,rmz,rgz) +c + factor=em**4 + +c + eeggl=(nutf*rl(1,1)+ndtf*rl(1,2))*factor + eegzl=(nutf*rl(2,1)+ndtf*rl(2,2))*factor + eezzl=(nutf*rl(3,1)+ndtf*rl(3,2))*factor + eeggr=(nutf*rr(1,1)+ndtf*rr(1,2))*factor + eegzr=(nutf*rr(2,1)+ndtf*rr(2,2))*factor + eezzr=(nutf*rr(3,1)+ndtf*rr(3,2))*factor + eetotl=eeggl+eegzl+eezzl + eetotr=eeggr+eegzr+eezzr + eetot=Peleft*Ppright*eetotl + . +(1D0-Peleft)*(1D0-Ppright)*eetotr + end +c +c********************************************************************** +c sfee3j calculates the process e+ e- --> g,z --> q + qb + g +c includes: no spin averages, final state statistics, couplings +c********************************************************************** + subroutine sfee3j(nloop,eetot) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon + common /strong/ as,gs,renscl,qcdl,b0,nf + common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr + common /reshel/zeegghl,zeegzhl,zeezzhl,zeegghr,zeegzhr, + $ zeezzhr,zeetohl,zeetohr,zeetoh + common /polar/ Peleft,Ppright + dimension rl(3,2),rr(3,2) + dimension zhl(3,2),zhr(3,2) +c + call eeq2gn(1,nloop,rl(1,1),rr(1,1),zhl,zhr,rmz,rgz) +c + factor=em**4*gs**2 + +c + eeggl=(nutf*rl(1,1)+ndtf*rl(1,2))*factor + eegzl=(nutf*rl(2,1)+ndtf*rl(2,2))*factor + eezzl=(nutf*rl(3,1)+ndtf*rl(3,2))*factor + eeggr=(nutf*rr(1,1)+ndtf*rr(1,2))*factor + eegzr=(nutf*rr(2,1)+ndtf*rr(2,2))*factor + eezzr=(nutf*rr(3,1)+ndtf*rr(3,2))*factor + eetotl=eeggl+eegzl+eezzl + eetotr=eeggr+eegzr+eezzr + eetot=Peleft*Ppright*eetotl + . +(1D0-Peleft)*(1D0-Ppright)*eetotr +c---also apply the appropriate factors to the gluon helicity amplitudes + if (nloop.eq.0) then + zeegghl=(nutf*zhl(1,1)+ndtf*zhl(1,2))*factor + zeegzhl=(nutf*zhl(2,1)+ndtf*zhl(2,2))*factor + zeezzhl=(nutf*zhl(3,1)+ndtf*zhl(3,2))*factor + zeegghr=(nutf*zhr(1,1)+ndtf*zhr(1,2))*factor + zeegzhr=(nutf*zhr(2,1)+ndtf*zhr(2,2))*factor + zeezzhr=(nutf*zhr(3,1)+ndtf*zhr(3,2))*factor + zeetohl=zeegghl+zeegzhl+zeezzhl + zeetohr=zeegghr+zeegzhr+zeezzhr + if(icon.eq.1) then + zeetoh=Peleft*Ppright*zeegghl + . +(1D0-Peleft)*(1D0-Ppright)*zeegghr + elseif(icon.eq.2) then + zeetoh=Peleft*Ppright*zeegzhl + . +(1D0-Peleft)*(1D0-Ppright)*zeegzhr + elseif(icon.eq.3) then + zeetoh=Peleft*Ppright*zeezzhl + . +(1D0-Peleft)*(1D0-Ppright)*zeezzhr + else + zeetoh=Peleft*Ppright*zeetohl + . +(1D0-Peleft)*(1D0-Ppright)*zeetohr + endif + endif + end +c +c********************************************************************** +c sfee4j calculates the process e+ e- --> g,z --> q + qb + ( gg or qqb) +c includes: no spin averages, final state statistics, couplings +c********************************************************************** + subroutine sfee4j(nloop,eetot) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /genera/ w,em,sw2,rmz,rgz,nutf,ndtf,icon + common /strong/ as,gs,renscl,qcdl,b0,nf + common /result/eeggl,eegzl,eezzl,eeggr,eegzr,eezzr,eetotl,eetotr + common /polar/ Peleft,Ppright + dimension rggl(3,2),rggr(3,2),rqql(3,5),rqqr(3,5) + dimension zhl(3,2),zhr(3,2) + common /order/iorder + common /perm/ips +c + do i=1,3 + do j=1,5 + rqql(i,j)=0D0 + rqqr(i,j)=0D0 + enddo + do k=1,2 + rggl(i,k)=0D0 + rggr(i,k)=0D0 + enddo + enddo + sum=1 + if(ips.eq.1)then + if(iorder.eq.0)then + call eeq2gn(2,nloop,rggl(1,1),rggr(1,1),zhl,zhr,rmz,rgz) + call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,1,3,2,4) + call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,2,4,1,3) + sum=2 + elseif(iorder.eq.1)then + call eeq2gn(2,nloop,rggl(1,1),rggr(1,1),zhl,zhr,rmz,rgz) + elseif(iorder.eq.3)then + call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,1,3,2,4) + call eeq4(rqql(1,1),rqqr(1,1),rmz,rgz,2,4,1,3) + sum=2 + endif + elseif(ips.eq.2)then + if(iorder.eq.0.or.iorder.eq.2)then + call eeq2gn(2,nloop,rggl(1,1),rggr(1,1),zhl,zhr,rmz,rgz) + endif + endif +c + factor=em**4*gs**4 + +c + eeggl=(+nutf*rggl(1,1)/2D0+ndtf*rggl(1,2)/2D0 + . +nutf*rqql(1,1)/4D0/sum+ndtf*rqql(1,2)/4D0/sum + . +0.5D0*nutf*(nutf-1)*rqql(1,3)/sum + . +0.5D0*ndtf*(ndtf-1)*rqql(1,4)/sum + . +nutf*ndtf*rqql(1,5)/sum)*factor + eegzl=(+nutf*rggl(2,1)/2D0+ndtf*rggl(2,2)/2D0 + . +nutf*rqql(2,1)/4D0/sum+ndtf*rqql(2,2)/4D0/sum + . +0.5D0*nutf*(nutf-1)*rqql(2,3)/sum + . +0.5D0*ndtf*(ndtf-1)*rqql(2,4)/sum + . +nutf*ndtf*rqql(2,5)/sum)*factor + eezzl=(+nutf*rggl(3,1)/2D0+ndtf*rggl(3,2)/2D0 + . +nutf*rqql(3,1)/4D0/sum+ndtf*rqql(3,2)/4D0/sum + . +0.5D0*nutf*(nutf-1)*rqql(3,3)/sum + . +0.5D0*ndtf*(ndtf-1)*rqql(3,4)/sum + . +nutf*ndtf*rqql(3,5)/sum)*factor + eeggr=(+nutf*rggr(1,1)/2D0+ndtf*rggr(1,2)/2D0 + . +nutf*rqqr(1,1)/4D0/sum+ndtf*rqqr(1,2)/4D0/sum + . +0.5D0*nutf*(nutf-1)*rqqr(1,3)/sum + . +0.5D0*ndtf*(ndtf-1)*rqqr(1,4)/sum + . +nutf*ndtf*rqqr(1,5)/sum)*factor + eegzr=(+nutf*rggr(2,1)/2D0+ndtf*rggr(2,2)/2D0 + . +nutf*rqqr(2,1)/4D0/sum+ndtf*rqqr(2,2)/4D0/sum + . +0.5D0*nutf*(nutf-1)*rqqr(2,3)/sum + . +0.5D0*ndtf*(ndtf-1)*rqqr(2,4)/sum + . +nutf*ndtf*rqqr(2,5)/sum)*factor + eezzr=(+nutf*rggr(3,1)/2D0+ndtf*rggr(3,2)/2D0 + . +nutf*rqqr(3,1)/4D0/sum+ndtf*rqqr(3,2)/4D0/sum + . +0.5D0*nutf*(nutf-1)*rqqr(3,3)/sum + . +0.5D0*ndtf*(ndtf-1)*rqqr(3,4)/sum + . +nutf*ndtf*rqqr(3,5)/sum)*factor + eetotl=eeggl+eegzl+eezzl + eetotr=eeggr+eegzr+eezzr + eetot=Peleft*Ppright*eetotl + . +(1D0-Peleft)*(1D0-Ppright)*eetotr + end +c +c********************************************************************** +c matrix element for: e+ e- --> g,z --> q qb + n gluons +c********************************************************************** +c information: the particles in plab are ordered as follows: +c 1,2 incoming particles 1: e+ and 2: e- +c 3,4 quark and anti-quark +c 5,.., 4+n gluons in final state. +c********************************************************************** +c output: the matrix elements for the subprocesses returned in +c left-handed contribution +c rl(3,2): first index are mixings ==> 1: gg 2: gz 3: zz +c second index are flavours ==> 1: u-type, 2: d-type. +c right-handed contribution +c rr(3,2): first index are mixings ==> 1: gg 2: gz 3: zz +c second index are flavours ==> 1: u-type, 2: d-type. +c not included: strong coupling constants! +c statistics, spin averages. +c********************************************************************** +c literature: +c lowest order : berends, giele & kuijf, np 321, 595, 1989 +c next to lowest order : giele & glover, prd 46, 1980, 1992 +c********************************************************************** +c function eeq2gn(n,nloop,rl,rr,rmz,rgz) calculates the matrix element +c squared of: e+ e- --> g,z --> q + qb + n gluons +c the input vectors are in /kinema/ p(4,7), where the sixth +c and seventh vectors are the e- and the e+ outgoing!. +c the first and second are the quark and the anti-quark. +c********************************************************************** + subroutine eeq2gn(n,nloop,rl,rr,zhl,zhr,rmz,rgz) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /esadb/ zsadb(2,2,2,4) + common /curent/ zrcur(2,2),zlcur(2,2) + common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen + common /linpr/ rin(7,7) + common /tcuts/ ymin,dly + common /strong/ as,gs,renscl,qcdl,b0,nf + common /zcup/czur,czul,czdr,czdl,czlr,czll,cznr,cznl + common /gcup/cgur,cgul,cgdr,cgdl,cglr,cgll,cgnr,cgnl + common /order/iorder + common /perm/ips + parameter(pi=3.141592653589793238D0) + dimension col(2,2,0:2) + dimension rl(3,2),rr(3,2) + dimension zhl(3,2),zhr(3,2) + dimension zrr(2,4),zrl(2,4),zlr(2,4),zll(2,4) + dimension z1rr(2,4),z1rl(2,4),z1lr(2,4),z1ll(2,4) + dimension zrsadb(2,2,2,4),zlsadb(2,2,2,4) + dimension zr1adb(2,2,2,4),zl1adb(2,2,2,4) + col(1,1,0)=3D0 + col(1,2,0)=0D0 + col(2,1,0)=0D0 + col(2,2,0)=0D0 + col(1,1,1)=4D0 + col(1,2,1)=0D0 + col(2,1,1)=0D0 + col(2,2,1)=0D0 +* + if(ips.eq.1)then + if(iorder.eq.0.or.iorder.eq.1)then +c$$$ c2a=12D0 +C---THIS MODIFICATION SENT BY NIGEL GLOVER TO AVOID SUMMING OVER PERMS + c2a=6D0 + c2b=0D0 + endif + elseif(ips.eq.2)then + if(iorder.eq.0.or.iorder.eq.2)then + c2a=-2D0/3D0 + c2b=c2a + endif + endif +c$$$ col(1,1,2)=c2b +C---THIS MODIFICATION SENT BY NIGEL GLOVER TO AVOID SUMMING OVER PERMS + col(1,1,2)=c2a + col(1,2,2)=c2b + col(2,1,2)=c2b + col(2,2,2)=c2a +c + s=2D0*rin(6,7) + zprg=(0D0,-1D0)/s + zprz=(0D0,-1D0)/(s-rmz**2+(0D0,1D0)*rgz*rmz) +c + if (n.eq.0) then + factor=1D0 + ihel=1 + iperm=1 + do 6 ia=1,2 + do 6 ib=1,2 + zrsadb(ia,ib,1,1)=zkod(1,ia)*zko(2,ib) + zlsadb(ib,ia,1,1)=conjg(zrsadb(ia,ib,1,1)) + 6 continue + endif +c + if (n.eq.1) then + factor=2D0 + ihel=2 + iperm=1 + do 7 ia=1,2 + do 7 ib=1,2 + zrsadb(ia,ib,1,1)=+(zkod(1,ia)*zud(1,2) + . +zkod(3,ia)*zud(3,2))*zko(2,ib)/(zud(1,3)*zud(3,2)) + zrsadb(ia,ib,1,2)=-zkod(1,ia)*(zko(3,ib)*zd(3,1) + . +zko(2,ib)*zd(2,1))/(zd(1,3)*zd(3,2)) + zlsadb(ib,ia,1,1)=conjg(zrsadb(ia,ib,1,1)) + zlsadb(ib,ia,1,2)=conjg(zrsadb(ia,ib,1,2)) + 7 continue + if(nloop.eq.1)then + call ert(1,2,3,alp,bep,gap,dep) + call ert(2,1,3,alm,bem,gam,dem) + do 17 ia=1,2 + do 17 ib=1,2 + zr1adb(ia,ib,1,1)= + . +alp*zkod(1,ia)*zud(1,2)*zko(2,ib)/(zud(1,3)*zud(3,2)) + . +bep*zkod(3,ia)*zud(3,2)*zko(2,ib)/(zud(1,3)*zud(3,2)) + . +gap*zkod(3,ia)*zd(1,3)*zko(1,ib)/(zud(1,3)*zd(3,2)) + . +dep*zd(1,3)/(zud(1,3)*zd(2,1)) + . *(zkod(1,ia)*zko(1,ib)+zkod(2,ia)*zko(2,ib)+zkod(3,ia)*zko(3,ib)) + zr1adb(ia,ib,1,2)= + . -alm*zkod(1,ia)*zd(2,1)*zko(2,ib)/(zd(1,3)*zd(3,2)) + . -bem*zkod(1,ia)*zd(3,1)*zko(3,ib)/(zd(1,3)*zd(3,2)) + . -gam*zkod(2,ia)*zud(2,3)*zko(3,ib)/(zud(1,3)*zd(3,2)) + . -dem*zud(2,3)/(zd(2,3)*zud(1,2)) + . *(zkod(1,ia)*zko(1,ib)+zkod(2,ia)*zko(2,ib)+zkod(3,ia)*zko(3,ib)) + zl1adb(ib,ia,1,1)=conjg(zr1adb(ia,ib,1,1)) + zl1adb(ib,ia,1,2)=conjg(zr1adb(ia,ib,1,2)) + 17 continue + endif + endif +c + if (n.eq.2) then + factor=4D0 + ihel=4 + iperm=2 + call f2sadb(3,4,1,1,2,3,4) + call f2sadb(4,3,2,1,3,2,4) + do 8 ih=1,ihel + do 8 ip=1,iperm + do 8 ia=1,2 + do 8 ib=1,2 + zrsadb(ia,ib,ip,ih)=zsadb(ia,ib,ip,ih) + zlsadb(ib,ia,ip,ih)=conjg(zsadb(ia,ib,ip,ih)) + 8 continue + endif +c + do 10 ih=1,ihel + do 10 ip=1,iperm + zrr(ip,ih)=zspv(zrsadb(1,1,ip,ih),zrcur(1,1)) + zrl(ip,ih)=zspv(zlsadb(1,1,ip,ih),zrcur(1,1)) + zlr(ip,ih)=zspv(zrsadb(1,1,ip,ih),zlcur(1,1)) + zll(ip,ih)=zspv(zlsadb(1,1,ip,ih),zlcur(1,1)) +c +c non-tree level loop amplitudes contracted into lepton current +c + if(nloop.eq.1.and.n.eq.1)then + z1rr(ip,ih)=zspv(zr1adb(1,1,ip,ih),zrcur(1,1)) + z1rl(ip,ih)=zspv(zl1adb(1,1,ip,ih),zrcur(1,1)) + z1lr(ip,ih)=zspv(zr1adb(1,1,ip,ih),zlcur(1,1)) + z1ll(ip,ih)=zspv(zl1adb(1,1,ip,ih),zlcur(1,1)) + endif + 10 continue +c + do 15 i=1,3 + do 15 j=1,2 + rl(i,j)=0D0 + rr(i,j)=0D0 + 15 continue +c + if(nloop.eq.0)then + tloop=1D0 + elseif(nloop.eq.1)then +c +c one-loop contribution proportional to tree level +c + if(n.eq.0)then + y12=rin(1,2)/rin(6,7) +c$$$ dl12=log(y12/ymin) +c$$$ e0=-dl12**2+3D0*dl12/2D0+pi**2/6D0-1D0/2D0 +c$$$ tloop=1D0+4D0*as/3D0/pi*e0 +C---THIS MODIFICATION MADE BY MIKE BY ANALOGY WITH THE ONE BELOW + e0=pi**2/2-4 + tloop=0D0+4D0*as/3D0/pi*e0 + elseif(n.eq.1)then + y12=rin(1,2)/rin(6,7) + y13=rin(1,3)/rin(6,7) + y23=rin(2,3)/rin(6,7) +c$$$ dl12=log(y12/ymin) +c$$$ dl13=log(y13/ymin) +c$$$ dl23=log(y23/ymin) +c$$$ dlm2=log(renscl**2/rmz**2/ymin) +c$$$ e0=29D0/9D0+pi**2/3D0+3D0/4D0*dl13+3D0/4D0*dl23 +c$$$ . -dl13**2-dl23**2 +c$$$ e2=-dl12**2+3D0*dl12/2D0+pi**2/6D0-1D0/2D0 +c$$$ e1=5D0/9D0 +C---THIS MODIFICATION SENT BY NIGEL GLOVER TO GIVE ERT FINITE TERMS + dly12=log(y12) + dly13=log(y13) + dly23=log(y23) + dlm2=0D0 +C---PI**2 COEFFICIENT CHANGED BY MIKE + e0=pi**2-4D0-dly13**2/2D0-dly23**2/2D0 + e2=pi**2/2-4D0-dly12**2/2D0 + e1=0D0 + tloop=0D0+as*3D0/2D0/pi*(e0-e2/9D0-nf*e1/3D0)+as*b0*dlm2 + if(iorder.eq.1)tloop=0D0+as*3D0/2D0/pi*e0+as*b0*dlm2 + if(iorder.eq.2)tloop=0D0-as*3D0/2D0/pi*e2/9D0+as*b0*dlm2 + if(iorder.eq.3)tloop=0D0-as*3D0/2D0/pi*nf*e1/3D0+as*b0*dlm2 + elseif(n.eq.2)then + tloop=1D0 + endif + endif + fac=factor*tloop + zprgg=zprg*conjg(zprg) + zprzg=zprz*conjg(zprg) + zprzz=zprz*conjg(zprz) + do 20 i=1,ihel + do 20 j1=1,iperm + do 20 j2=1,iperm + zrrrr=zrr(j1,i)*conjg(zrr(j2,i)) + zrlrl=zrl(j1,i)*conjg(zrl(j2,i)) + zlrlr=zlr(j1,i)*conjg(zlr(j2,i)) + zllll=zll(j1,i)*conjg(zll(j2,i)) +* +* left-handed contributions +* + rl(1,1)=rl(1,1)+fac*col(j1,j2,n)*zprgg*( + . +(zlrlr*cgur*cgur+zllll*cgul*cgul)*cgll*cgll ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zlrlr*czur*cgur+zllll*czul*cgul)*czll*cgll ) + rl(2,1)=rl(2,1)+2D0*real(zz) + rl(3,1)=rl(3,1)+fac*col(j1,j2,n)*zprzz*( + . +(zlrlr*czur*czur+zllll*czul*czul)*czll*czll ) + rl(1,2)=rl(1,2)+fac*col(j1,j2,n)*zprgg*( + . +(zlrlr*cgdr*cgdr+zllll*cgdl*cgdl)*cgll*cgll ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zlrlr*czdr*cgdr+zllll*czdl*cgdl)*czll*cgll ) + rl(2,2)=rl(2,2)+2D0*real(zz) + rl(3,2)=rl(3,2)+fac*col(j1,j2,n)*zprzz*( + . +(zlrlr*czdr*czdr+zllll*czdl*czdl)*czll*czll ) +* +* right-handed contributions +* + rr(1,1)=rr(1,1)+fac*col(j1,j2,n)*zprgg*( + . +(zrrrr*cgur*cgur+zrlrl*cgul*cgul)*cglr*cglr ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zrrrr*czur*cgur+zrlrl*czul*cgul)*czlr*cglr ) + rr(2,1)=rr(2,1)+2D0*real(zz) + rr(3,1)=rr(3,1)+fac*col(j1,j2,n)*zprzz*( + . +(zrrrr*czur*czur+zrlrl*czul*czul)*czlr*czlr ) + rr(1,2)=rr(1,2)+fac*col(j1,j2,n)*zprgg*( + . +(zrrrr*cgdr*cgdr+zrlrl*cgdl*cgdl)*cglr*cglr ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zrrrr*czdr*cgdr+zrlrl*czdl*cgdl)*czlr*cglr ) + rr(2,2)=rr(2,2)+2D0*real(zz) + rr(3,2)=rr(3,2)+fac*col(j1,j2,n)*zprzz*( + . +(zrrrr*czdr*czdr+zrlrl*czdl*czdl)*czlr*czlr ) + 20 continue +* +* Equivalent expressions for the gluon interference term: (+-)+(-+) +* + if (n.eq.1.and.nloop.eq.0) then + j1=1 + j2=1 + zrrrr=zrr(j1,1)*conjg(zrr(j2,2)) + zrlrl=zrl(j1,1)*conjg(zrl(j2,2)) + zlrlr=zlr(j1,1)*conjg(zlr(j2,2)) + zllll=zll(j1,1)*conjg(zll(j2,2)) +* +* left-handed contributions +* + zhl(1,1)=fac*col(j1,j2,n)*zprgg*( + . +(zlrlr*cgur*cgur+zllll*cgul*cgul)*cgll*cgll ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zlrlr*czur*cgur+zllll*czul*cgul)*czll*cgll ) + zhl(2,1)=2D0*real(zz) + zhl(3,1)=fac*col(j1,j2,n)*zprzz*( + . +(zlrlr*czur*czur+zllll*czul*czul)*czll*czll ) + zhl(1,2)=fac*col(j1,j2,n)*zprgg*( + . +(zlrlr*cgdr*cgdr+zllll*cgdl*cgdl)*cgll*cgll ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zlrlr*czdr*cgdr+zllll*czdl*cgdl)*czll*cgll ) + zhl(2,2)=2D0*real(zz) + zhl(3,2)=fac*col(j1,j2,n)*zprzz*( + . +(zlrlr*czdr*czdr+zllll*czdl*czdl)*czll*czll ) +* +* right-handed contributions +* + zhr(1,1)=fac*col(j1,j2,n)*zprgg*( + . +(zrrrr*cgur*cgur+zrlrl*cgul*cgul)*cglr*cglr ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zrrrr*czur*cgur+zrlrl*czul*cgul)*czlr*cglr ) + zhr(2,1)=2D0*real(zz) + zhr(3,1)=fac*col(j1,j2,n)*zprzz*( + . +(zrrrr*czur*czur+zrlrl*czul*czul)*czlr*czlr ) + zhr(1,2)=fac*col(j1,j2,n)*zprgg*( + . +(zrrrr*cgdr*cgdr+zrlrl*cgdl*cgdl)*cglr*cglr ) + zz=fac*col(j1,j2,n)*zprzg*( + . +(zrrrr*czdr*cgdr+zrlrl*czdl*cgdl)*czlr*cglr ) + zhr(2,2)=2D0*real(zz) + zhr(3,2)=fac*col(j1,j2,n)*zprzz*( + . +(zrrrr*czdr*czdr+zrlrl*czdl*czdl)*czlr*czlr ) + endif +c +c higher order contributions +c + if(nloop.eq.1)then +c +c interference between non-tree level structures and tree level +c + if(n.eq.1)then + fac=factor*3D0*as/4D0/pi + do 30 i=1,ihel + do 30 j1=1,iperm + do 30 j2=1,iperm + z1rrrr=z1rr(j1,i)*conjg(zrr(j2,i)) + z1rlrl=z1rl(j1,i)*conjg(zrl(j2,i)) + z1lrlr=z1lr(j1,i)*conjg(zlr(j2,i)) + z1llll=z1ll(j1,i)*conjg(zll(j2,i)) +* +* left-handed contributions +* + zz=fac*col(j1,j2,n)*zprgg*( + . +(z1lrlr*cgur*cgur+z1llll*cgul*cgul)*cgll*cgll ) + rl(1,1)=rl(1,1)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzg*( + . +(z1lrlr*czur*cgur+z1llll*czul*cgul)*czll*cgll ) + rl(2,1)=rl(2,1)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzz*( + . +(z1lrlr*czur*czur+z1llll*czul*czul)*czll*czll ) + rl(3,1)=rl(3,1)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprgg*( + . +(z1lrlr*cgdr*cgdr+z1llll*cgdl*cgdl)*cgll*cgll ) + rl(1,2)=rl(1,2)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzg*( + . +(z1lrlr*czdr*cgdr+z1llll*czdl*cgdl)*czll*cgll ) + rl(2,2)=rl(2,2)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzz*( + . +(z1lrlr*czdr*czdr+z1llll*czdl*czdl)*czll*czll ) + rl(3,2)=rl(3,2)+2D0*real(zz) +* +* right-handed contributions +* + zz=fac*col(j1,j2,n)*zprgg*( + . +(z1rrrr*cgur*cgur+z1rlrl*cgul*cgul)*cglr*cglr ) + rr(1,1)=rr(1,1)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzg*( + . +(z1rrrr*czur*cgur+z1rlrl*czul*cgul)*czlr*cglr ) + rr(2,1)=rr(2,1)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzz*( + . +(z1rrrr*czur*czur+z1rlrl*czul*czul)*czlr*czlr ) + rr(3,1)=rr(3,1)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprgg*( + . +(z1rrrr*cgdr*cgdr+z1rlrl*cgdl*cgdl)*cglr*cglr ) + rr(1,2)=rr(1,2)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzg*( + . +(z1rrrr*czdr*cgdr+z1rlrl*czdl*cgdl)*czlr*cglr ) + rr(2,2)=rr(2,2)+2D0*real(zz) + zz=fac*col(j1,j2,n)*zprzz*( + . +(z1rrrr*czdr*czdr+z1rlrl*czdl*czdl)*czlr*czlr ) + rr(3,2)=rr(3,2)+2D0*real(zz) + 30 continue + endif + endif +c + end +c +c*********************************************************************** +c subroutine ert(iq,ip,ik,al,be,ga,de) computes the functions multiplyin +c the different one-loop helicity structures +c*********************************************************************** +c call ert(iq,ip,ik,al,be,ga,de) +c iq,ip and ik point to the momenta of the quark, antiquark and gluon +c al, be and ga are the functions multiplying the helicity structures +c*********************************************************************** + subroutine ert(iq,ip,ik,al,be,ga,de) + implicit double precision(a-h,o-y) + parameter(pi=3.141592653589793238D0) + common /linpr/ rin(7,7) + common /order/iorder +c + yqp=rin(iq,ip)/rin(6,7) + yqk=rin(iq,ik)/rin(6,7) + ykp=rin(ik,ip)/rin(6,7) + omyqp=1D0-yqp + omyqk=1D0-yqk + omykp=1D0-ykp + dlqp=log(yqp) + dlqk=log(yqk) + dlkp=log(ykp) + dlomqp=log(omyqp) + dlomqk=log(omyqk) + dlomkp=log(omykp) + spqp=rsp(yqp) + spqk=rsp(yqk) + spkp=rsp(ykp) + rqpkp=dlqp*dlkp-dlqp*dlomqp-dlkp*dlomkp+pi**2/6D0-spqp-spkp + rqpqk=dlqp*dlqk-dlqp*dlomqp-dlqk*dlomqk+pi**2/6D0-spqp-spqk + rqkkp=dlqk*dlkp-dlqk*dlomqk-dlkp*dlomkp+pi**2/6D0-spqk-spkp +c + al0=-rqkkp + . -ykp*(4D0-3D0*ykp)/2D0/omykp**2*dlkp + . -ykp/2D0/omykp + be0=-rqkkp + . +(4D0-3D0*ykp)/2D0/omykp*dlkp + . +1D0 + ga0=+ykp/2D0/omykp**2*dlkp + . +ykp/2D0/omykp + de0=+(yqp*ykp/2D0/omykp**2+yqp/omykp)*dlkp + . +yqp/2D0/omykp +c + al2=-rqpqk + . -omyqp**2/yqk**2*rqpkp + . -ykp/yqk*dlqp + . -(ykp*(4D0-3D0*ykp)/2D0/omykp**2+ykp**2/yqk/omykp)*dlkp + . -ykp/2D0/omykp + be2=-rqpqk + . +(yqp*omyqp/yqk**2+1D0/yqk)*rqpkp + . +(yqk/omyqp**2+omykp/yqk)*dlqp + . +((4D0-3D0*ykp)/2D0/omykp+ykp/yqk)*dlkp + . +yqk/omyqp + ga2=+ykp/yqk**2*rqpkp + . -(yqk/omyqp**2-1D0/yqk)*dlqp + . +(ykp/2D0/omykp**2+ykp/yqk/omykp)*dlkp + . +ykp/omyqp+ykp/2D0/omykp + de2=+yqp*omyqp/yqk**2*rqpkp + . +yqp/yqk*dlqp + . +(yqp*ykp/2D0/omykp**2+yqp*omyqp/yqk/omykp)*dlkp + . +yqp/2D0/omykp +c +C---THESE LINES REMOVED BY MIKE FOR ERT DEFINITION OF COLLINEAR PIECES +c$$$ al0=al0-3D0/4D0*dlqk-3D0/4D0*dlkp +c$$$ be0=be0-3D0/4D0*dlqk-3D0/4D0*dlkp +c$$$ al2=al2-3D0/2D0*dlqp +c$$$ be2=be2-3D0/2D0*dlqp +c + if(iorder.eq.0)then + al=al0-al2/9D0 + be=be0-be2/9D0 + ga=ga0-ga2/9D0 + de=de0-de2/9D0 + elseif(iorder.eq.1)then + al=al0 + + be=be0 + + ga=ga0 + + de=de0 + + elseif(iorder.eq.2)then + al=-al2/9D0 + be=-be2/9D0 + ga=-ga2/9D0 + de=-de2/9D0 + endif + +c check : t0 and t2 should be zero +c t0=-al0+be0-ga0-2D0/yqp*de0 +c t2=-al2+be2-ga2-2D0/yqp*de2 +c write(*,*)t0,t2 + end +c********************************************************************** +c subroutine f2sadb(i1,i2,index,j1,j2,j3,j4) calculates the four +c helicity amplitudes for the permutation of the gluons (i1,i2). +c the quark(q) having + helicity and the anti-quark(p) - helictity. +c the calculation is organized as follows: +c after calculating the various propagators, all the spinor +c currents from the form: ( propagator contracted with a spinor) +c are composed. the currents containing the q are dotted. those with +c the qb are undotted. +c example: zD12m2 is the spinor (q+g(i1)+g(i2))adb * g(i2)b +c after initializing various constants belonging to a certain helicity +c combination, e.g. zppx belongs to the combination where both +c gluons have helicity +, the currents zsadb are calculated. +c********************************************************************** +c calling conventions: +c i1,i2 = permutation of the gluons. because they are used +c as indices we added two. +c this way we have (i1,i2)= (3,4) or (4,3) +c index = permutation counter +c j1,j2,j3,j4 = indices to get the helicity currents on the +c right spot in the data structures sadb. +c also needed are commons /spinor/ and /dotpr/ +c********************************************************************** +c output conventions: +c put all the helicity combinations of a certain +c permutation (i1,i2) in the data-structure zsadb(2,2,2,4) +c********************************************************************** + subroutine f2sadb(i1,i2,index,j1,j2,j3,j4) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /esadb/ zsadb(2,2,2,4) + common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen + common /linpr/ rin(7,7) + dimension zD12m1(2),zD12m2(2),zD12mp(2), + . z12pm1(2),z12pm2(2),z12pmq(2) +c + pro12=2D0*rin(i1,i2) + proD12=pro12+2D0*(rin(1,i1)+rin(1,i2)) + pro12p=pro12+2D0*(rin(2,i1)+rin(2,i2)) +c + do 10 i=1,2 + zD12m1(i)=zkod(1,i)*zud(1,i1)+zkod(i2,i)*zud(i2,i1) + zD12m2(i)=zkod(1,i)*zud(1,i2)+zkod(i1,i)*zud(i1,i2) + zD12mp(i)=zkod(1,i)*zud(1,2)+zkod(i1,i)*zud(i1,2) + . +zkod(i2,i)*zud(i2,2) + z12pm1(i)=zko(2,i)*zd(2,i1)+zko(i2,i)*zd(i2,i1) + z12pm2(i)=zko(2,i)*zd(2,i2)+zko(i1,i)*zd(i1,i2) + z12pmq(i)=zko(2,i)*zd(2,1)+zko(i1,i)*zd(i1,1)+zko(i2,i)*zd(i2,1) + 10 continue +c + zppa=-1D0/(zud(1,i1)*zud(i1,i2)*zud(i2,2)) + zpma=zd(1,i1)*zud(1,i2)/(pro12*proD12*zud(1,i1)) + zpmb=zd(2,i1)*zud(2,i2)/(pro12*pro12p*zd(2,i2)) + zpmc=1D0/(pro12*zud(1,i1)*zd(i2,2)) + zmpa=zud(2,i1)**2/(pro12*pro12p*zud(2,i2)) + zmpb=zd(1,i2)**2/(pro12*proD12*zd(1,i1)) + zmpc=-zd(1,i2)*zud(2,i1)/(pro12*zd(1,i1)*zud(2,i2)) + zmma=1D0/(zd(1,i1)*zd(i1,i2)*zd(i2,2)) +c + do 20 ia=1,2 + do 20 ib=1,2 + zsadb(ia,ib,index,j1)=-zppa*zD12mp(ia)*zko(2,ib) + zsadb(ia,ib,index,j2)=-zpma*zD12m2(ia)*zko(2,ib) + . -zpmb*zkod(1,ia)*z12pm1(ib) + . -zpmc*zD12m2(ia)*z12pm1(ib) + zsadb(ia,ib,index,j3)=-zmpa*zkod(1,ia)*z12pm2(ib) + . -zmpb*zD12m1(ia)*zko(2,ib) + . -zmpc*zkod(1,ia)*zko(2,ib) + zsadb(ia,ib,index,j4)=-zmma*zkod(1,ia)*z12pmq(ib) + 20 continue +c + end +c +c********************************************************************** +c matrix element for: e+ e- --> g,z --> q qb q qb +c********************************************************************** +c information: the particles in plab are ordered as follows: +c 1,2 : incoming e+ and e- +c 3,5 : outgoing quarks +c 4,6 : outgoing anti-quarks +c sw2: square of sin(weak angle) +c rmz: mass of z +c rgz: width of z +c +c output: the matrix elements for the subprocesses returned in +c rl(3,5) for left handed contribution +c rr(3,5) for right handed contribution +c first index 1 = gg +c 2 = gz +c 3 = zz +c second index 1 = same quarks up-type +c 2 = same quarks down-type +c 3 = different quarks up-type +c 4 = different quarks down-type +c 5 = different quarks up - down pairs +c********************************************************************** +c not included: strong coupling constants! +c statistics, spin averages. +c********************************************************************** +c calling procedure for subroutines: +c 1 : anti-particle +c 2 : particle +c 3,5: labels of quarks in plab array +c 4,6: labels of anti-quarks in plab array +c particles 1 and 2 are always flipped: we consider all particle +c to be outgoing. +c********************************************************************** +c most of the formulae in this program can be found in: +c exact expressions for processes involving a vector boson +c and up to five partons. (berends, giele, kuijf) +c********************************************************************** +C---MODIFIED BY MIKE TO ENABLE SIMPLE SUMMING OVER PERMS +C QUARK 1 HAS BEEN REPLACED BY VARIABLE I1 EVERYWHERE, ETC + subroutine eeq4(rl,rr,rmz,rgz,I1,I2,I3,I4) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /mom/ plab(4,10) + common /kinema/ p(4,7) + common /zqqqq/ zr(4,6),zl(4,6) + common /curent/ zrcur(2,2),zlcur(2,2) + dimension za(2,2),zb(2,2) + dimension rl(3,5),rr(3,5) +c +c fill h functions +c the problem of an odd number of quark-particles with negative +c energy is handled within the h-functions. +c handle ppmm and mmpp + call hpmmp(I1,I4,I3,I2,-1D0,za(1,1),zb(1,1)) + zr(1,1)=zspv(zrcur(1,1),za(1,1)) + zl(1,1)=zspv(zlcur(1,1),za(1,1)) + zr(1,4)=zspv(zrcur(1,1),zb(1,1)) + zl(1,4)=zspv(zlcur(1,1),zb(1,1)) + call hpmmp(I3,I2,I1,I4,-1D0,za(1,1),zb(1,1)) + zr(2,4)=zspv(zrcur(1,1),za(1,1)) + zl(2,4)=zspv(zlcur(1,1),za(1,1)) + zr(2,1)=zspv(zrcur(1,1),zb(1,1)) + zl(2,1)=zspv(zlcur(1,1),zb(1,1)) +c +c handle pmmp and mppm + call hpmmp(I1,I2,I3,I4,1D0,za(1,1),zb(1,1)) + zr(1,2)=zspv(zrcur(1,1),za(1,1)) + zl(1,2)=zspv(zlcur(1,1),za(1,1)) + zr(1,3)=zspv(zrcur(1,1),zb(1,1)) + zl(1,3)=zspv(zlcur(1,1),zb(1,1)) + call hpmmp(I3,I4,I1,I2,1D0,za(1,1),zb(1,1)) + zr(2,3)=zspv(zrcur(1,1),za(1,1)) + zl(2,3)=zspv(zlcur(1,1),za(1,1)) + zr(2,2)=zspv(zrcur(1,1),zb(1,1)) + zl(2,2)=zspv(zlcur(1,1),zb(1,1)) +c +c handle pmpm and mpmp + call hpmpm(I1,I2,I3,I4,1D0,za(1,1),zb(1,1)) + zr(1,5)=zspv(zrcur(1,1),za(1,1)) + zl(1,5)=zspv(zlcur(1,1),za(1,1)) + zr(1,6)=zspv(zrcur(1,1),zb(1,1)) + zl(1,6)=zspv(zlcur(1,1),zb(1,1)) + call hpmpm(I1,I4,I3,I2,-1D0,za(1,1),zb(1,1)) + zr(2,5)=zspv(zrcur(1,1),za(1,1)) + zl(2,5)=zspv(zlcur(1,1),za(1,1)) + zr(2,6)=zspv(zrcur(1,1),zb(1,1)) + zl(2,6)=zspv(zlcur(1,1),zb(1,1)) + call hpmpm(I3,I4,I1,I2,1D0,za(1,1),zb(1,1)) + zr(3,5)=zspv(zrcur(1,1),za(1,1)) + zl(3,5)=zspv(zlcur(1,1),za(1,1)) + zr(3,6)=zspv(zrcur(1,1),zb(1,1)) + zl(3,6)=zspv(zlcur(1,1),zb(1,1)) + call hpmpm(I3,I2,I1,I4,-1D0,za(1,1),zb(1,1)) + zr(4,5)=zspv(zrcur(1,1),za(1,1)) + zl(4,5)=zspv(zlcur(1,1),za(1,1)) + zr(4,6)=zspv(zrcur(1,1),zb(1,1)) + zl(4,6)=zspv(zlcur(1,1),zb(1,1)) +c +c insert different couplings to obtain subprocesses +c + call squar(1,1,rmz,rgz,rl(1,1),rr(1,1)) + call squar(2,2,rmz,rgz,rl(1,2),rr(1,2)) + call squar(1,3,rmz,rgz,rl(1,3),rr(1,3)) + call squar(2,4,rmz,rgz,rl(1,4),rr(1,4)) + call squar(1,2,rmz,rgz,rl(1,5),rr(1,5)) +c + end +c +c********************************************************************** +c function squar calculates the matrix element squared of the +c e+ e- --> g,z --> q + qb + q + qb +c********************************************************************** +c desciption of the common zqqqq +c the first is the index for the squared matrix +c the second index are the different subamplitudes (6) +c the first and the second are indices in the squared matrix +c index 3=1 --> ++--, =2 --> +--+, =3 --> -++-, =4 --> --++ +c index 3=5 --> +-+-, =6 --> -+-+ +c********************************************************************** +c output : rl(3) 1: gg 2: gz 3: zz +c rr(3) 1: gg 2: gz 3: zz +c********************************************************************** + subroutine squar(ityp1,ityp2,rmz,rgz,rl,rr) + implicit double precision (a-h,o-y) + implicit complex*16(z) + parameter(c1=8D0,c2=-8D0/3D0) + common /zqqqq/ zr(4,6),zl(4,6) + dimension s2g(4,2),s2z(4,2),s4g(2,4),s4z(2,4) + dimension col4(4,4),col2(2,2),rl(3),rr(3) + common /zcup/czur,czul,czdr,czdl,czlr,czll,cznr,cznl + common /gcup/cgur,cgul,cgdr,cgdl,cglr,cgll,cgnr,cgnl + data col4/c1,c2,c1,c2,c2,c1,c2,c1,c1,c2,c1,c2,c2,c1,c2,c1/ + data col2/c1,c1,c1,c1/ +C---THESE LINES REMOVED BY MIKE TO FACILITATE SUMMATION OVER PERMS +C (THEY WERE SUPERFLUOUS ANYWAY, SINCE THE ARRAYS WERE EMPTY ALREADY) +c$$$c +c$$$c init answer +c$$$ rl(1)=0D0 +c$$$ rl(2)=0D0 +c$$$ rl(3)=0D0 +c$$$ rr(1)=0D0 +c$$$ rr(2)=0D0 +c$$$ rr(3)=0D0 +c +c couplings and all that (like energy dependent width for z propagator). + s=2D0*rinpro(6,7) + zprg=(0D0,-1D0)/s + zprz=(0D0,-1D0)/(s-rmz**2+(0D0,1D0)*s*rgz/rmz) + if(ityp1.eq.1)then + cg1r=cgur + cg1l=cgul + cz1r=czur + cz1l=czul + elseif(ityp1.eq.2)then + cg1r=cgdr + cg1l=cgdl + cz1r=czdr + cz1l=czdl + endif + if(ityp2.eq.1.or.ityp2.eq.3)then + cg2r=cgur + cg2l=cgul + cz2r=czur + cz2l=czul + elseif(ityp2.eq.2.or.ityp2.eq.4)then + cg2r=cgdr + cg2l=cgdl + cz2r=czdr + cz2l=czdl + endif +c + if(ityp1.eq.ityp2)then + delta=1D0 + else + delta=0D0 + endif +c +c s2(1,.) is ppmm + s2g(1,1)=delta*cg1r + s2z(1,1)=delta*cz1r + s2g(1,2)=delta*cg1l + s2z(1,2)=delta*cz1l +c s2(2,.) is pmmp + s2g(2,1)=cg1r + s2z(2,1)=cz1r + s2g(2,2)=cg2l + s2z(2,2)=cz2l +c s2(3,.) is mppm + s2g(3,1)=cg1l + s2z(3,1)=cz1l + s2g(3,2)=cg2r + s2z(3,2)=cz2r +c s2(4,.) is mmpp + s2g(4,1)=delta*cg1l + s2z(4,1)=delta*cz1l + s2g(4,2)=delta*cg1r + s2z(4,2)=delta*cz1r +c + zprgg=zprg*conjg(zprg) + zprzg=zprz*conjg(zprg) + zprzz=zprz*conjg(zprz) + do 10 i=1,4 + ij=i + do 11 ix=1,2 + do 12 iy=1,2 + zrr=zr(ix,i)*conjg(zr(iy,i)) + zll=zl(ix,i)*conjg(zl(iy,i)) +c +c left-handed contributions +c + rl(1)=rl(1)+col2(ix,iy)*s2g(ij,ix)*s2g(ij,iy)* + . zprgg*cgll*cgll*zll + zz=col2(ix,iy)*s2g(ij,ix)*s2z(ij,iy)* + . zprzg*czll*cgll*zll + rl(2)=rl(2)+2D0*real(zz) + rl(3)=rl(3)+col2(ix,iy)*s2z(ij,ix)*s2z(ij,iy)* + . zprzz*czll*czll*zll +c +c right-handed contributions +c + rr(1)=rr(1)+col2(ix,iy)*s2g(ij,ix)*s2g(ij,iy)* + . zprgg*cglr*cglr*zrr + zz=col2(ix,iy)*s2g(ij,ix)*s2z(ij,iy)* + . zprzg*czlr*cglr*zrr + rr(2)=rr(2)+2D0*real(zz) + rr(3)=rr(3)+col2(ix,iy)*s2z(ij,ix)*s2z(ij,iy)* + . zprzz*czlr*czlr*zrr + + 12 continue + 11 continue + 10 continue +c +c s4(1,.) is pmpm + s4g(1,1)=cg1r + s4z(1,1)=cz1r + s4g(1,2)=delta*cg1r + s4z(1,2)=delta*cz1r + s4g(1,3)=cg2r + s4z(1,3)=cz2r + s4g(1,4)=delta*cg1r + s4z(1,4)=delta*cz1r +c s4(2,.) is mpmp + s4g(2,1)=cg1l + s4z(2,1)=cz1l + s4g(2,2)=delta*cg1l + s4z(2,2)=delta*cz1l + s4g(2,3)=cg2l + s4z(2,3)=cz2l + s4g(2,4)=delta*cg1l + s4z(2,4)=delta*cz1l +c + do 20 i=5,6 + ij=(i-4) + do 21 ix=1,4 + do 22 iy=1,4 + zrr=zr(ix,i)*conjg(zr(iy,i)) + zll=zl(ix,i)*conjg(zl(iy,i)) +c +c left-handed contributions +c + rl(1)=rl(1)+col4(ix,iy)*s4g(ij,ix)*s4g(ij,iy)* + . zprgg*cgll*cgll*zll + zz=col4(ix,iy)*s4g(ij,ix)*s4z(ij,iy)* + . zprzg*czll*cgll*zll + rl(2)=rl(2)+2D0*real(zz) + rl(3)=rl(3)+col4(ix,iy)*s4z(ij,ix)*s4z(ij,iy)* + . zprzz*czll*czll*zll +c +c right-handed contributions +c + rr(1)=rr(1)+col4(ix,iy)*s4g(ij,ix)*s4g(ij,iy)* + . zprgg*cglr*cglr*zrr + + zz=col4(ix,iy)*s4g(ij,ix)*s4z(ij,iy)* + . zprzg*czlr*cglr*zrr + + rr(2)=rr(2)+2D0*real(zz) + rr(3)=rr(3)+col4(ix,iy)*s4z(ij,ix)*s4z(ij,iy)* + . zprzz*czlr*czlr*zrr + + 22 continue + 21 continue + 20 continue +c + end +c the problem of an odd number of quark-particles with negative +c energy is handled within the h-functions. +c********************************************************************** +c subroutine hpmpm calculates the h-function when the helicities of +c the quarks are +-+-. +c call the h-functions with: +c argument 1-4 : quark permutation +c argument 5 is the sign of the h-function (exchange of two +c identical fermions) +c argument 6 is where to put it in zessen +c argument 7 is where to put its complex conjugate. +c z1 and z2 are the 2*2 complex matrices which are the h-vectors +c in spinor language. z2 is always the c.c. of z1. when there +c are an odd number of quarks with negative energy we have to add +c an extra minus sign when complex conjugating. this sign is +c in the common /dotpr/ negaen and is calculated before through +c the call setspn(4,4) in squar +c********************************************************************** + subroutine hpmpm(i1,i2,i3,i4,sign,z1,z2) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen + common /linpr/ rin(7,7) + dimension z1(2,2),z2(2,2) +c + pro34=2D0*rin(i3,i4) + pro134=pro34+2D0*(rin(i1,i3)+rin(i1,i4)) + pro234=pro34+2D0*(rin(i2,i3)+rin(i2,i4)) +c + zf1=-sign*zd(i1,i3)/(pro34*pro134) + zf2=-sign*zud(i2,i4)/(pro34*pro234) + do 10 ix=1,2 + do 10 iy=1,2 + z1(ix,iy)=+zf1*zko(i2,iy)*(zud(i1,i4)*zkod(i1,ix)+ + . zud(i3,i4)*zkod(i3,ix)) + . -zf2*zkod(i1,ix)*(zd(i2,i3)*zko(i2,iy)+ + . zd(i4,i3)*zko(i4,iy)) + 10 continue + call congat(z1,z2) + end +c +c********************************************************************** +c subroutine hpmmp calculates the h-function when the helicities of +c the quarks are +--+. see also subroutine hpmpm +c********************************************************************** + subroutine hpmmp(i1,i2,i3,i4,sign,z1,z2) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen + common /linpr/ rin(7,7) + dimension z1(2,2),z2(2,2) +c + pro34=2D0*rin(i3,i4) + pro134=pro34+2D0*(rin(i1,i3)+rin(i1,i4)) + pro234=pro34+2D0*(rin(i2,i3)+rin(i2,i4)) +c + zf1=-sign*zd(i1,i4)/(pro34*pro134) + zf2=-sign*zud(i2,i3)/(pro34*pro234) + do 10 ix=1,2 + do 10 iy=1,2 + z1(ix,iy)=+zf1*zko(i2,iy)*(zud(i1,i3)*zkod(i1,ix)+ + . zud(i4,i3)*zkod(i4,ix)) + . -zf2*zkod(i1,ix)*(zd(i2,i4)*zko(i2,iy)+ + . zd(i3,i4)*zko(i3,iy)) + 10 continue + call congat(z1,z2) + end +c +c********************************************************************** +c subroutine congat(z1,z2) takes the complex conjugate +c spinor-tensor of z1 into z2. use negaen to correct for +c an odd number of quarks with negative energy. +c********************************************************************** + subroutine congat(z1,z2) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen + dimension z1(2,2),z2(2,2) +c + z2(1,1)=negaen*conjg(z1(1,1)) + z2(2,1)=negaen*conjg(z1(1,2)) + z2(1,2)=negaen*conjg(z1(2,1)) + z2(2,2)=negaen*conjg(z1(2,2)) + end +c +c********************************************************************** +c subroutine setspn(n) determines the co-variant weyl-vd waerden +c spinors and all the possible spinor-inproducts. +c special treatment for particles 6 and 7 (leptons) +c particles with negative energy have an artificial minus +c sign in their dotted spinors. +c********************************************************************** +c calling conventions: +c n = number of of vectors (n must be <= 7 ) +c common /kinema/ p = array of covariant vectors (-,-,-,+) metric +c note: none of the vectors should be parallel to the y-axis. +c********************************************************************** +c output conventions: +c common /dotpr/zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen +c input parameters remain unchanged. +c with: +c zud(i,j) = undotted spinor inproduct +c zd(i,j) = dotted spinor inproduct +c zko(i,1:2) = undotted co-variant spinors +c zkod(i,1:2) = dotted co-variant spinors +c negaen = number of quarks with negative +c energy. needed when complex-conjugating. +c negaen=1,-1 when that number is even,odd +c********************************************************************** + subroutine setspn(n) + implicit double precision(a-h,o-y) + implicit complex*16(z) + common /kinema/ p(4,7) + common /dotpr/ zud(7,7),zd(7,7),zko(7,2),zkod(7,2),negaen + common /linpr/ rin(7,7) +c + do 10 i=1,7 + if((i.le.n).or.(i.ge.6))then + zko(i,1)=(p(1,i)-(0D0,1D0)*p(3,i)) + . /sqrt(dcmplx(p(4,i)-p(2,i))) + zko(i,2)=sqrt(dcmplx(p(4,i)-p(2,i))) + zkod(i,1)=conjg(zko(i,1)) + zkod(i,2)=conjg(zko(i,2)) + if (p(4,i).lt.0D0) then + zkod(i,1)=-zkod(i,1) + zkod(i,2)=-zkod(i,2) + endif + endif + 10 continue +c + do 20 i=1,n-1 + do 20 j=i,n + zud(i,j)=-zko(i,1)*zko(j,2)+zko(i,2)*zko(j,1) + zud(j,i)=-zud(i,j) + zd(i,j)=-zkod(i,1)*zkod(j,2)+zkod(i,2)*zkod(j,1) + zd(j,i)=-zd(i,j) + 20 continue +c + do 30 i=1,n-1 + do 30 j=i,n + rin(i,j)=rinpro(i,j) + rin(j,i)=rin(i,j) + 30 continue + rin(1,6)=rinpro(1,6) + rin(1,7)=rinpro(1,7) + rin(2,6)=rinpro(2,6) + rin(2,7)=rinpro(2,7) + rin(6,7)=rinpro(6,7) + rin(7,6)=rin(6,7) + end +c +c********************************************************************** +c function rinpro(i,j) determines the minkowski inproduct. +c the parameters i and j refer to the momentum array /kinema/ p. +c both vectors are contra-variant. +c********************************************************************** + function rinpro(i,j) + implicit double precision(a-h,o-y) + common /kinema/ p(4,7) +c + rinpro=p(4,i)*p(4,j)-p(1,i)*p(1,j)-p(2,i)*p(2,j)-p(3,i)*p(3,j) + end +c +c********************************************************************** +c 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 +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 +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*********************************************************************** +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----------------------------------------------------------------------- +C----------------------------------------------------------------------- +C----------------------------------------------------------------------- diff --git a/gbook.f b/gbook.f new file mode 100644 index 0000000..5bb9f79 --- /dev/null +++ b/gbook.f @@ -0,0 +1,955 @@ +C********************************************************************* +C*** HEADER ******************************************************** +C +C GBOOK - A SMALL HISTOGRAM PACKAGE +C WRITTEN DECEMBER 1979, LAST CHANGED MARCH 1983 +C AUTHOR: TORBJORN SJOSTRAND, DEPT. OF THEORETICAL PHYSICS, +C UNIVERSITY OF LUND, SOLVEGATAN 14A, S-223 62 LUND, SWEDEN +C PRESENT ADDRESS: FNAL, THEORY, TEL. 312-840-3753 +C PLEASE REPORT ANY ERRORS TO THE AUTHOR +C +C*** INTRODUCTION ************************************************** +C +C GBOOK IS A SMALL SUBROUTINE PACKAGE FOR GETTING ONE- OR TWO- +C DIMENSIONAL HISTOGRAMS ON AN ORDINARY LINE PRINTER OR TERMINAL. +C IT CAN BE USED IN A MANNER VERY SIMILAR TO HBOOK, BUT HAS BEEN +C WRITTEN COMPLETELY INDEPENDENTLY. THE CODE IS FULLY IN FORTRAN77, +C WITH A MINIMUM OF MACHINE DEPENDENCE. THE USAGE CAN BE DIVIDED INTO +C FOUR STEPS: BOOKING, FILLING, EDITING AND PRINTING. ALL SUBROUTINES +C HAVE NAMES SIX CHARACTERS LONG AND BEGINNING WITH G. +C +C*** BOOKING ******************************************************* +C +C THERE ARE NMAX HISTOGRAMS AT THE DISPOSAL OF THE USER, EACH GIVEN BY +C A NUMBER BETWEEN 1 AND NMAX. BEFORE A HISTOGRAM CAN BE USED, SPACE +C MUST BE RESERVED FOR IT. NMAX IS A COMPILE-TIME PARAMETER ADDED BY +C MIKE SEYMOUR. +C +C CALL GBOOK1(ID,TITLE,NX,XL,XU) +C PURPOSE: BOOK A ONE-DIMENSIONAL HISTOGRAM. +C ID: HISTOGRAM NUMBER, INTEGER BETWEEN 1 AND NMAX. +C TITLE: HISTOGRAM TITLE, CAN BE GIVEN EITHER AS A CHARACTER STRING +C OF AT MOST 60 CHARACTERS OR AS A CHARACTER*60 VARIABLE. +C NX: NUMBER OF BINS (IN X DIRECTION) IN HISTOGRAM, INTEGER BETWEEN +C 1 AND 100. +C XL, XU: LOWER AND UPPER BOUND, RESPECTIVELY, ON THE (X) RANGE +C COVERED BY THE HISTOGRAM. +C +C CALL GBOOK2(ID,TITLE,NX,XL,XU,NY,YL,YU) +C PURPOSE: BOOK A TWO-DIMENSIONAL HISTOGRAM. +C ID: HISTOGRAM NUMBER, INTEGER BETWEEN 1 AND NMAX. +C TITLE: HISTOGRAM TITLE, SEE GBOOK1. +C NX: NUMBER OF BINS IN X DIRECTION, INTEGER BETWEEN 1 AND 50. +C XL, XU: LOWER AND UPPER BOUND ON X RANGE OF HISTOGRAM. +C NY: NUMBER OF BINS IN Y DIRECTION, ARBITRARY POSITIVE INTEGER. +C YL, YU: LOWER AND UPPER BOUND OF Y RANGE OF HISTOGRAM. +C +C*** FILLING ******************************************************* +C +C FOR BOOKED HISTOGRAMS WEIGHTS CAN BE FILLED AT GIVEN COORDINATES. +C +C CALL GFILL1(ID,X,W) +C PURPOSE: FILL IN A ONE-DIMENSIONAL HISTOGRAM. +C ID: HISTOGRAM NUMBER. +C X: X COORDINATE OF POINT. +C W: WEIGHT TO BE ADDED IN THIS POINT. +C +C CALL GFILL2(ID,X,Y,W) +C PURPOSE: FILL IN A TWO-DIMENSIONAL HISTOGRAM. +C ID: HISTOGRAM NUMBER. +C X: X COORDINATE OF POINT. +C Y: Y COORDINATE OF POINT. +C W: WEIGHT TO BE ADDED IN THIS POINT. +C +C*** EDITING ******************************************************* +C +C FOR EDITING OF HISTOGRAMS BEFORE PRINTOUT TWO ROUTINES ARE AVAILABLE. +C +C CALL GSCALE(ID,F) +C PURPOSE: RESCALE THE CONTENTS OF A HISTOGRAM. +C ID: HISTOGRAM NUMBER. +C F: RESCALING FACTOR, I.E. FACTOR THAT ALL BIN CONTENTS (INCLUDING +C OVERFLOW ETC.) ARE MULTIPLIED BY. +C REMARK: A TYPICAL RESCALING FACTOR FOR A ONE-DIMENSIONAL HISTOGRAM +C COULD BE F = 1/(BIN SIZE * NUMBER OF EVENTS) = +C = NX/(XU-XL) * 1/(NUMBER OF EVENTS). +C +C CALL GOPERA(ID1,OPER,ID2,ID3,F1,F2) +C PURPOSE: THIS IS A GENERAL PURPOSE ROUTINE FOR EDITING ONE OR SEVERAL +C HISTOGRAMS, WHICH ALL ARE ASSUMED TO HAVE THE SAME NUMBER OF +C BINS. OPERATIONS ARE CARRIED OUT BIN BY BIN, INCLUDING OVERFLOW +C BINS ETC. +C OPER: GIVES THE TYPE OF OPERATION TO BE CARRIED OUT, A ONE-CHARACTER +C STRING OR A CHARACTER*1 VARIABLE. +C OPER= '+', '-', '*', '/': ADD, SUBTRACT, MULTIPLY OR DIVIDE THE +C CONTENTS IN ID1 AND ID2 AND PUT THE RESULT IN ID3. F1 AND F2, IF +C NOT 1., GIVE FACTORS BY WHICH THE ID1 AND ID2 BIN CONTENTS ARE +C MULTIPLIED BEFORE THE INDICATED OPERATION. (DIVISION WITH A +C VANISHING BIN CONTENT WILL GIVE 0.) +C OPER= 'A', 'S', 'L': FOR 'S' THE SQUARE ROOT OF THE CONTENT IN ID1 +C IS TAKEN (RESULT 0 FOR NEGATIVE BIN CONTENTS) AND FOR 'L' THE +C 10-LOGARITHM IS TAKEN (A NONPOSITIVE BIN CONTENT IS BEFORE THAT +C REPLACED BY 0.8 TIMES THE SMALLEST POSITIVE BIN CONTENT). +C THEREAFTER, IN ALL THREE CASES, THE CONTENT IS MULTIPLIED BY F1 +C AND ADDED WITH F2, AND THE RESULT IS PLACED IN ID3. THUS ID2 +C IS DUMMY IN THESE CASES. +C OPER= 'M': INTENDED FOR STATISTICAL ANALYSIS, BIN-BY-BIN MEAN AND +C STANDARD DEVIATION OF A VARIABLE, ASSUMING THAT ID1 CONTAINS +C ACCUMULATED WEIGHTS, ID2 ACCUMULATED WEIGHT*VARIABLE AND +C ID3 ACCUMULATED WEIGHT*VARIABLE-SQUARED. AFTERWARDS ID2 WILL +C CONTAIN THE MEAN VALUES (=ID2/ID1) AND ID3 THE STANDARD +C DEVIATIONS (=SQRT(ID3/ID1-(ID2/ID1)**2)). IN THE END, F1 +C MULTIPLIES ID1 (FOR NORMALIZATION PURPOSES), WHILE F2 IS DUMMY. +C +C*** PRINTING ****************************************************** +C +C AT PRINTING AXES ARE CHOSEN SUCH THAT HISTOGRAMS ARE SUPPOSED +C TO FIT INTO ONE PAGE. FOR ONE-DIMENSIONAL HISTOGRAMS NUMBERS +C SMALLER IN MAGNITUDE THAN 10**(-10) ARE CONSIDERED TO BE 0, +C OTHERWISE SCALES ARE CHOSEN AUTOMATICALLY FOR MAXIMUM RESOLUTION. +C TWO-DIMENSIONAL HISTOGRAMS ARE STRONGLY ORIENTED TOWARDS HAVING +C WEIGHTS IN THE ORDER OF OR BIGGER THAN UNITY : SIGNS 1 - 9 ARE +C CHOSEN IN A LINEAR SCALE UP TO WEIGHTS 9.5, WHEREAS A - Z WILL +C EITHER BE CHOSEN IN A LINEAR SCALE WITH STEP 1 OR, IF MAXIMUM +C WEIGHT IS LARGER THAN 36.5, IN A LOGARITHMICALLY EVEN SCALE. +C NEGATIVE NUMBERS ARE ALLOWED BOTH IN ONE- AND TWO-DIMENSIONAL +C HISTOGRAMS, AND ARE PRECEDED BY A - SIGN. +C +C CALL GCLEAR +C PURPOSE: PRINT OUT ALL HISTOGRAMS THAT HAVE BEEN FILLED, AND +C RESET THEM THEREAFTER TO 0. +C +C CALL GPRINT(ID) +C PURPOSE: PRINT OUT A SINGLE HISTOGRAM. +C ID: HISTOGRAM TO BE PRINTED. +C +C CALL GRESET(ID) +C PURPOSE: RESET ALL BIN CONTENTS, INCLUDING OVERFLOW ETC., TO 0. +C ID: HISTOGRAM TO BE RESET. +C +C*** COMMON BLOCK AND SPACE REQUIREMENTS *************************** +C +C A COMMONBLOCK +C PARAMETER (NSIZE=200000,NMAX=2000) +C COMMON /GBOOK/ A(NSIZE) +C IS USED TO STORE HISTOGRAM INFORMATION. THE HISTOGRAM INDEX TAKES +C NMAX+2 POSITIONS. EACH BOOKED ONE- (TWO-) DIMENSIONAL HISTOGRAM TAKES +C AN ADDITIONAL 38+NX (38+NX*NY) POSITIONS. THE PROGRAM HAS TO BE +C RECOMPILED WITH CHANGED COMMONBLOCK IF MORE IS REQUIRED. +C +C*** END DESCRIPTION ********************************************** + SUBROUTINE IDATE_MOD(iyr, imo, iday) +C Drop-in replacement for old 'idate' subroutine +C Returns year, month, day as integers + + INTEGER iyr, imo, iday + INTEGER values(8) + + CALL DATE_AND_TIME(values=values) + + iyr = values(1) + imo = values(2) + iday = values(3) + + RETURN + END +C********************************************************************* + SUBROUTINE GBOOK1(ID,TITLE,NX,XL,XU) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + CHARACTER TITLE*(*),TITFX*60 + EQUIVALENCE (REQ,IEQ) + IF (ID.GT.NMAX .OR. A(1)+A(2)+38+NX.GT.NSIZE) THEN + WRITE (6,200) ID,NMAX,INT(A(1)+A(2)+38+NX+0.5),NSIZE + RETURN + ENDIF + A(ID+2)=A(1)+A(2) + A(2)=A(2)+38+NX + IS=A(ID+2)+0.5 + A(IS+1)=NX + A(IS+2)=XL + A(IS+3)=XU + A(IS+4)=(XU-XL)/NX + A(IS+5)=1 + CALL GRESET(ID) + TITFX=TITLE//' ' + DO 100 IT=1,20 + IEQ=256**2*ICHAR(TITFX(3*IT-2:3*IT-2))+256*ICHAR(TITFX(3*IT-1: + &3*IT-1))+ICHAR(TITFX(3*IT:3*IT)) + 100 A(IS+18+NX+IT)=REQ + RETURN + 200 FORMAT (' ERROR: Too much space requested in GBOOK1!'/ + & ' Requested ID=',I4,', Maximum ID=',I4/ + & ' Requested space=',I6,', Maximum space=',I6/ + & ' Recompile with larger NMAX and/or NSIZE') + END +C********************************************************************* + SUBROUTINE GBOOK2(ID,TITLE,NX,XL,XU,NY,YL,YU) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + CHARACTER TITLE*(*),TITFX*60 + EQUIVALENCE (REQ,IEQ) + IF (ID.GT.NMAX .OR. A(1)+A(2)+38+NX*NY.GT.NSIZE) THEN + WRITE (6,200) ID,NMAX,INT(A(1)+A(2)+38+NX*NY+0.5),NSIZE + RETURN + ENDIF + A(ID+2)=A(1)+A(2) + A(2)=A(2)+38+NX*NY + IS=A(ID+2)+0.5 + A(IS+1)=NX + A(IS+2)=XL + A(IS+3)=XU + A(IS+4)=(XU-XL)/NX + A(IS+5)=NY + A(IS+6)=YL + A(IS+7)=YU + A(IS+8)=(YU-YL)/NY + CALL GRESET(ID) + TITFX=TITLE//' ' + DO 100 IT=1,20 + IEQ=256**2*ICHAR(TITFX(3*IT-2:3*IT-2))+256*ICHAR(TITFX(3*IT-1: + &3*IT-1))+ICHAR(TITFX(3*IT:3*IT)) + 100 A(IS+18+NX*NY+IT)=REQ + RETURN + 200 FORMAT (' ERROR: Too much space requested in GBOOK2!'/ + & ' Requested ID=',I4,', Maximum ID=',I4/ + & ' Requested space=',I6,', Maximum space=',I6/ + & ' Recompile with larger NMAX and/or NSIZE') + END +C********************************************************************* + SUBROUTINE GFILL1(ID,X,W) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + IS=A(ID+2)+0.5 + A(IS+9)=A(IS+9)+1. + IOX=2 + IF(X.LT.A(IS+2)) IOX=1 + IF(X.GE.A(IS+3)) IOX=3 + A(IS+12+IOX)=A(IS+12+IOX)+W + IF(IOX.NE.2) RETURN + IX=(X-A(IS+2))/A(IS+4) + A(IS+19+IX)=A(IS+19+IX)+W + RETURN + END +C********************************************************************* + SUBROUTINE GFILL2(ID,X,Y,W) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + IS=A(ID+2)+0.5 + A(IS+9)=A(IS+9)+1. + IOX=2 + IF(X.LT.A(IS+2)) IOX=1 + IF(X.GE.A(IS+3)) IOX=3 + IOY=2 + IF(Y.LT.A(IS+6)) IOY=1 + IF(Y.GE.A(IS+7)) IOY=3 + A(IS+6+3*IOY+IOX)=A(IS+6+3*IOY+IOX)+W + IF(IOX.NE.2.OR.IOY.NE.2) RETURN + IX=(X-A(IS+2))/A(IS+4) + IY=(Y-A(IS+6))/A(IS+8) + IC=INT(A(IS+1)+0.5)*IY+IX + A(IS+19+IC)=A(IS+19+IC)+W + RETURN + END +C********************************************************************* + SUBROUTINE GSCALE(ID,F) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + IS=A(ID+2)+0.5 + DO 100 IC=IS+10,IS+18+INT(A(IS+1)+0.5)*INT(A(IS+5)+0.5) + 100 A(IC)=F*A(IC) + RETURN + END +C********************************************************************* + SUBROUTINE GOPERA(ID1,OPER,ID2,ID3,F1,F2) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + CHARACTER OPER*1 + IF (ID1.GT.NMAX .OR. ID2.GT.NMAX .OR. ID3.GT.NMAX) RETURN + IS1=A(ID1+2)+0.5 + IS2=A(ID2+2)+0.5 + IS3=A(ID3+2)+0.5 + NC=INT(A(IS3+1)+0.5)*INT(A(IS3+5)+0.5) + IF(OPER.EQ.'+'.OR.OPER.EQ.'-'.OR.OPER.EQ.'*'.OR.OPER.EQ.'/') + &A(IS3+9)=A(IS1+9)+A(IS2+9) + IF(OPER.EQ.'A'.OR.OPER.EQ.'S'.OR.OPER.EQ.'L') A(IS3+9)=A(IS1+9) + IF(OPER.EQ.'+') THEN + DO 100 IC=10,18+NC + 100 A(IS3+IC)=F1*A(IS1+IC)+F2*A(IS2+IC) + ELSEIF(OPER.EQ.'-') THEN + DO 110 IC=10,18+NC + 110 A(IS3+IC)=F1*A(IS1+IC)-F2*A(IS2+IC) + ELSEIF(OPER.EQ.'*') THEN + DO 120 IC=10,18+NC + 120 A(IS3+IC)=F1*A(IS1+IC)*F2*A(IS2+IC) + ELSEIF(OPER.EQ.'/') THEN + DO 130 IC=10,18+NC + FA2=F2*A(IS2+IC) + IF(ABS(FA2).LE.1E-10) A(IS3+IC)=0. + 130 IF(ABS(FA2).GT.1E-10) A(IS3+IC)=F1*A(IS1+IC)/FA2 + ELSEIF(OPER.EQ.'A') THEN + DO 140 IC=10,18+NC + 140 A(IS3+IC)=F1*A(IS1+IC)+F2 + ELSEIF(OPER.EQ.'S') THEN + ZERO=0 + DO 150 IC=10,18+NC + 150 A(IS3+IC)=F1*SQRT(MAX(ZERO,A(IS1+IC)))+F2 + ELSEIF(OPER.EQ.'L') THEN + ZMIN=1E30 + DO 160 IC=19,18+NC + 160 IF(A(IS1+IC).LT.ZMIN.AND.A(IS1+IC).GT.1E-20) ZMIN=0.8*A(IS1+IC) + DO 170 IC=10,18+NC + 170 A(IS3+IC)=F1*LOG10(MAX(A(IS1+IC),ZMIN))+F2 + ELSEIF(OPER.EQ.'M') THEN + DO 180 IC=10,18+NC + IF(ABS(A(IS1+IC)).LE.1E-10) A(IS2+IC)=0. + IF(ABS(A(IS1+IC)).GT.1E-10) A(IS2+IC)=A(IS2+IC)/A(IS1+IC) + IF(ID3.NE.0.AND.ABS(A(IS1+IC)).LE.1E-10) A(IS3+IC)=0. + ZERO=0 + IF(ID3.NE.0.AND.ABS(A(IS1+IC)).GT.1E-10) A(IS3+IC)= + &SQRT(MAX(A(IS3+IC)/A(IS1+IC)-A(IS2+IC)**2,ZERO)) + 180 A(IS1+IC)=F1*A(IS1+IC) + ENDIF + RETURN + END +C********************************************************************* + SUBROUTINE GCLEAR + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + DO 100 ID=1,INT(A(1)+0.5) + IS=A(ID+2)+0.5 + IF(IS.EQ.0.OR.A(IS+9).LT.0.5) GOTO 100 + CALL GPRINT(ID) + CALL GRESET(ID) + 100 CONTINUE + RETURN + END +C********************************************************************* + SUBROUTINE GPRINT(ID) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + INTEGER GLAST + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + CHARACTER TITLE*60, CTIME*8, OUT*108, CHA(40)*1 + DIMENSION IROW(100), IFRA(100), DYAC(10), EV(20) + EQUIVALENCE (REQ,IEQ) + DATA DYAC/.04,.05,.06,.08,.10,.12,.15,.20,.25,.30/, LIN/41/ + DATA CHA/' ','0','1','2','3','4','5','6','7','8','9','A','B', + &'C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q', + &'R','S','T','U','V','W','X','Y','Z','*','-','!'/ + IF (ID.GT.NMAX) RETURN + IS=A(ID+2)+0.5 + IF(A(IS+9).LT.0.5) WRITE(6,400) ID + IF(A(IS+9).LT.0.5) RETURN + NX=INT(A(IS+1)+0.5) + NY=INT(A(IS+5)+0.5) + DO 100 IT=1,20 + EV(IT)=0. + REQ=A(IS+18+NX*NY+IT) + 100 TITLE(3*IT-2:3*IT)=CHAR(IEQ/256**2)//CHAR(MOD(IEQ,256**2)/256) + &//CHAR(MOD(IEQ,256)) + CALL IDATE_MOD(IMON,IDAY,IYEAR) + CALL DATE_AND_TIME(CTIME) + WRITE(6,410) ID, TITLE, IYEAR, IMON, IDAY, CTIME(1:5) + IF(NY.EQ.1) THEN + YMIN=A(IS+19) + YMAX=A(IS+19) + DO 120 IX=IS+20,IS+18+NX + IF(A(IX).LT.YMIN) YMIN=A(IX) + 120 IF(A(IX).GT.YMAX) YMAX=A(IX) + IF(YMAX-YMIN.GT.LIN*DYAC(1)*1E-9) THEN + IF(YMIN.GT.0..AND.YMIN.LT.YMAX/10.) YMIN=0. + IF(YMAX.LT.0..AND.YMAX.GT.YMIN/10.) YMAX=0. + IPOT=INT(LOG10(YMAX-YMIN)+10.)-10 + IF(YMAX-YMIN.LT.LIN*DYAC(1)*10.**IPOT) IPOT=IPOT-1 + IF(YMAX-YMIN.GT.LIN*DYAC(10)*10.**IPOT) IPOT=IPOT+1 + DELY=DYAC(1) + DO 130 IDEL=1,9 + 130 IF(YMAX-YMIN.GE.LIN*DYAC(IDEL)*10.**IPOT) DELY=DYAC(IDEL+1) + DY=DELY*10.**IPOT + DO 140 IX=1,NX + CTA=ABS(A(IS+18+IX))/DY + IROW(IX)=SIGN(CTA+0.95,A(IS+18+IX)) + 140 IFRA(IX)=10.*(CTA+1.25-AINT(CTA+0.95)) + IRMI=SIGN(ABS(YMIN)/DY+0.95,YMIN) + IRMA=SIGN(ABS(YMAX)/DY+0.95,YMAX) + DO 160 IR=IRMA,IRMI,-1 + IF(IR.EQ.0) GOTO 160 + OUT=' ' + DO 150 IX=1,NX + IF(IR.EQ.IROW(IX)) OUT(IX:IX)=CHA(IFRA(IX)+23*(IFRA(IX)/12)) + 150 IF(IR*(IROW(IX)-IR).GT.0) OUT(IX:IX)=CHA(35) + WRITE(6,420) IR*DELY, IPOT, OUT(1:GLAST(OUT)) + 160 CONTINUE + IPOT=INT(LOG10(MAX(YMAX,-YMIN))+10.0001)-10 + OUT=' ' + DO 170 IX=1,NX + IF(A(IS+18+IX).LT.-10.**(IPOT-4)) OUT(IX:IX)=CHA(39) + 170 IROW(IX)=10.**(3-IPOT)*ABS(A(IS+18+IX))+0.5 + WRITE(6,430) OUT(1:GLAST(OUT)) + DO 190 IR=4,1,-1 + DO 180 IX=1,NX + 180 OUT(IX:IX)=CHA(2+MOD(IROW(IX),10**IR)/10**(IR-1)) + 190 WRITE(6,440) IPOT+IR-4, OUT(1:GLAST(OUT)) + IPOT=INT(LOG10(MAX(-A(IS+2),A(IS+3)-A(IS+4)))+10.0001)-10 + OUT=' ' + DO 200 IX=1,NX + IF(A(IS+2)+(IX-1)*A(IS+4).LT.-10**(IPOT-3)) OUT(IX:IX)=CHA(39) + 200 IROW(IX)=10.**(2-IPOT)*ABS(A(IS+2)+(IX-1)*A(IS+4))+0.5 + WRITE(6,450) OUT(1:GLAST(OUT)) + DO 220 IR=3,1,-1 + DO 210 IX=1,NX + 210 OUT(IX:IX)=CHA(2+MOD(IROW(IX),10**IR)/10**(IR-1)) + 220 WRITE(6,440) IPOT+IR-3, OUT(1:GLAST(OUT)) + ENDIF + DO 230 IX=1,NX + CTA=ABS(A(IS+18+IX)) + X=A(IS+2)+(IX-0.5)*A(IS+4) + EV(1)=EV(1)+CTA + EV(2)=EV(2)+CTA*X + 230 EV(3)=EV(3)+CTA*X**2 + SMALL=1E-20 + XMEAN=EV(2)/MAX(EV(1),SMALL) + ZERO=0 + XRMS=SQRT(MAX(ZERO,EV(3)/MAX(EV(1),SMALL)-XMEAN**2)) + WRITE(6,460) INT(A(IS+9)+0.5),XMEAN,A(IS+13),A(IS+2),A(IS+14), + &XRMS,A(IS+15),A(IS+3) + ELSE + ZMAX=35.4 + DO 240 IC=IS+19,IS+NX*NY+18 + 240 IF(ABS(A(IC)).GT.ZMAX) ZMAX=1.0001*ABS(A(IC)) + ZFAC=26./LOG(ZMAX/9.5) + DO 290 IY=1,MAX(NY+4,38) + OUT=' ' + IF(IY.EQ.1.OR.IY.EQ.NY+2) THEN + OUT(4:4)=CHA(38) + DO 250 IX=5,2*NX+5 + 250 OUT(IX:IX)=CHA(39) + OUT(2*NX+6:2*NX+6)=CHA(38) + ELSEIF(IY.EQ.NY+3) THEN + DO 260 IX=6,2*NX+4,2 + 260 OUT(IX:IX)=CHA(2+(IX-3)/20) + ELSEIF(IY.EQ.NY+4) THEN + DO 270 IX=6,2*NX+4,2 + 270 OUT(IX:IX)=CHA(IX/2-10*((IX-3)/20)) + ELSEIF(IY.LE.NY+1) THEN + OUT(1:1)=CHA(2+(NY+2-IY)/10) + OUT(2:2)=CHA(4+NY-IY-10*((NY+2-IY)/10)) + OUT(4:4)=CHA(40) + DO 280 IX=6,2*NX+4,2 + CT=A(IS+16+NX*(NY+1-IY)+IX/2) + IF(CT.LE.-0.05) OUT(IX-1:IX-1)=CHA(39) + CTA=ABS(CT) + IF(CTA.GE.0.1) OUT(IX:IX)=CHA(2) + IF(CTA.GE.0.5.AND.CTA.LT.36.) OUT(IX:IX)=CHA(2+INT(CTA+0.5)) + 280 IF(ZMAX.GT.35.5.AND.CTA.GE.9.5) OUT(IX:IX)= + &CHA(12+INT(ZFAC*LOG(CTA/9.5))) + OUT(2*NX+6:2*NX+6)=CHA(40) + ENDIF + IF(IY.GE.2.AND.IY.LE.38) THEN + WMIN=0.1 + IF(IY.GE.3) WMIN=IY-2.5 + IF(ZMAX.GT.35.5.AND.IY.GE.12) WMIN=9.5*EXP((IY-12)/ZFAC) + ENDIF + IF(IY.EQ.1) WRITE(6,470) OUT(1:GLAST(OUT)) + IF(IY.GT.1.AND.IY.LE.38)WRITE(6,480)CHA(IY),WMIN,OUT(1:GLAST(OUT)) + 290 IF(IY.GE.39) WRITE(6,490) OUT(1:GLAST(OUT)) + DO 300 IY=1,NY + Y=A(IS+6)+(IY-0.5)*A(IS+8) + DO 300 IX=1,NX + X=A(IS+2)+(IX-0.5)*A(IS+4) + CTA=ABS(A(IS+18+NX*(IY-1)+IX)) + EV(1)=EV(1)+CTA + EV(2)=EV(2)+CTA*X + EV(3)=EV(3)+CTA*X**2 + EV(4)=EV(4)+CTA*Y + EV(5)=EV(5)+CTA*Y**2 + 300 EV(6)=EV(6)+CTA*X*Y + SMALL=1E-20 + XMEAN=EV(2)/MAX(EV(1),SMALL) + ZERO=0 + XRMS=SQRT(MAX(ZERO,EV(3)/MAX(EV(1),SMALL)-XMEAN**2)) + YMEAN=EV(4)/MAX(EV(1),SMALL) + ZERO=0 + YRMS=SQRT(MAX(ZERO,EV(5)/MAX(EV(1),SMALL)-YMEAN**2)) + XYCOR=(EV(6)/MAX(EV(1),SMALL)-XMEAN*YMEAN)/MAX(SMALL,XRMS*YRMS) + WRITE(6,500) INT(A(IS+9)+0.5),(A(IS+J), J=16,18), XMEAN, A(IS+2), + &XRMS, A(IS+3), (A(IS+J), J=13,15), YMEAN, A(IS+6), YRMS, A(IS+7), + &(A(IS+J), J=10,12), XYCOR + ENDIF + RETURN + 400 FORMAT(/5X,'HISTOGRAM NO',I4,' : NO ENTRIES') + 410 FORMAT('1'/5X,'HISTOGRAM NO',I4,4X,A60,'19',I2,'-',I2,'-',I2, + &1X,A5/) + 420 FORMAT(3X,F7.2,'*10**',I2,4X,A) + 430 FORMAT(/9X,'CONTENTS',4X,A) + 440 FORMAT(10X,'*10**',I2,4X,A) + 450 FORMAT(/9X,'LOW EDGE',4X,A) + 460 FORMAT(/5X,'ENTRIES =',I10,1P,6X,'MEAN =',E10.3,6X,'UNDERFLOW =' + &,E10.3,6X,'LOW EDGE =',E10.3/5X,'ALL CHAN =',E10.3,6X, + &'RMS =',E10.3,6X,'OVERFLOW =',E10.3,6X,'HIGH EDGE =',E10.3) + 470 FORMAT(8X,'SCALE',6X,A) + 480 FORMAT(6X,A1,1X,F7.1,'-',3X,A) + 490 FORMAT(19X,A) + 500 FORMAT(/5X,'ENTRIES =',I10,1P,17X,E10.3,' ? ',E10.3,' ? ',E10.3, + &7X,'XMEAN =',E10.3/5X,'XMIN =',E10.3,17X,36('-'),7X,'XRMS =', + &E10.3/5X,'XMAX =',E10.3,7X,'OVERFLOW',2X,E10.3,' ? ',E10.3, + &' ? ',E10.3,7X,'YMEAN =',E10.3/5X,'YMIN =',E10.3,17X,36('-'), + &7X,'YRMS =',E10.3/5X,'YMAX =',E10.3,17X,E10.3,' ? ',E10.3, + &' ? ',E10.3,7X,'XYCOR =',E10.3) + END +C********************************************************************* + INTEGER FUNCTION GLAST(STRING) +C LOOK FOR LAST NON-SPACE IN STRING + CHARACTER STRING*108 + I=108 + 10 IF (STRING(I:I).EQ.' ') THEN + I=I-1 + IF (I.GT.1) GOTO 10 + ENDIF + GLAST=I + END +C********************************************************************* + SUBROUTINE GRESET(ID) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + IS=A(ID+2)+0.5 + DO 110 IC=IS+9,IS+18+INT(A(IS+1)+0.5)*INT(A(IS+5)+0.5) + 110 A(IC)=0. + RETURN + END +C********************************************************************* + BLOCK DATA GDATA + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + DATA (A(J), J=1,NMAX+2)/NMAX,2,NMAX*0/ + END +C*** END GBOOK ***************************************************** +C----------------------------------------------------------------------- +C SOME NEW ROUTINES WRITTEN BY MIKE SEYMOUR: +C GSCAT2(ID,X,Y,N) +C EXACTLY THE SAME AS GFILL2(ID,X,Y,1.0) EXCEPT THAT ONLY THE +C FIRST N CALLS ARE USED. VERY USEFUL FOR MAKING SCATTER PLOTS +C THAT ALL HAVE THE SAME NUMBER OF ENTRIES. +C +C GAREA(ID,AREA) +C RESCALE HISTOGRAM SO THAT AREA UNDER IT BECOMES AREA +C +C GACCUM(ID) +C REPLACE HISTOGRAM BY A CUMULATIVE SUM OF ITS ENTRIES +C----------------------------------------------------------------------- + SUBROUTINE GSCAT2(ID,X,Y,N) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + ONE=1 + IS=A(ID+2)+0.5 + IF (INT(A(IS+9)+0.5).LT.N) CALL GFILL2(ID,X,Y,ONE) + END +C----------------------------------------------------------------------- + SUBROUTINE GAREA(ID,AREA) +C +C SCALES HISTOGRAM ID SO THAT THE AREA UNDER IT BECOMES AREA +C + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + IS=NINT(A(ID+2)) + BNWDTH=A(IS+4) + AREAOL=A(IS+14) + IF (AREAOL*BNWDTH.EQ.0.0) RETURN + AREASC=AREA/(AREAOL*BNWDTH) + CALL GSCALE (ID,AREASC) + END +C----------------------------------------------------------------------- + SUBROUTINE GACCUM(ID) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=400) + COMMON /GBOOK/ A(NSIZE) + IF (ID.GT.NMAX) RETURN + TOTAL=0.0 + IS=NINT(A(ID+2)) + IF (IS.EQ.0) RETURN + DO 50 IX=1,NINT(A(IS+1)) + TOTAL=TOTAL+A(IS+18+IX) + A(IS+18+IX)=TOTAL + 50 CONTINUE + END +C----------------------------------------------------------------------- + SUBROUTINE GTOPDR(ID,NEW,HIST,LOG) +C +C M SEYMOUR'S GBOOK->TOPDRAWER INTERFACE +C WRITES HISTOGRAM WITH TOPDRAWER PACKAGE +C +C NEW IS .TRUE. TO PUT HISTOGRAM ON A NEW PAGE +C .FALSE. TO PUT IT ON TOP OF THE LAST ONE +C +C HIST IS .TRUE. FOR HISTOGRAM +C .FALSE. FOR GRAPH +C +C LOG IS .TRUE. FOR LOGARITHMIC Y-SCALE +C .FALSE. FOR LINEAR Y-SCALE +C +C .TRUE. MAY BE ABBREVIATED TO 1 WHEN CALLING, AND +C .FALSE. TO 0 +C +C IF LOG IS SELECTED, THEN ALL ZERO OR NEGATIVE BINS ARE SET TO +C ONE THOUSANTH OF THE LOWEST OTHER BIN +C + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + LOGICAL NEW,HIST,LOG + CALL GTOPER(ID,NEW,HIST,LOG,0) + END +C----------------------------------------------------------------------- + SUBROUTINE GTOPER(ID,NEW,HIST,LOG,IDERR) +C +C IDENTICAL TO GTOPDR EXCEPT THAT IF IDERR IS NON-ZERO IT IS USED FOR +C ERROR BARS FOR 1-DIM HISTOGRAMS +C + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + LOGICAL NEW,HIST,LOG + IF (ID.GT.NMAX) RETURN + IS=NINT(A(ID+2)) + IF (IS .EQ. 0) RETURN + IF (A(IS+14) .EQ. 0.0) THEN + WRITE (6,1000) ID + ELSE + WRITE (6,1010) ID + IF (A(IS+5) .LT. 1.5) THEN + CALL GTOP1(ID,NEW,HIST,LOG,IDERR) + ELSE + CALL GTOP2(ID,NEW) + ENDIF + ENDIF + 1000 FORMAT (/5X,'HISTOGRAM NO',I4,' : NO ENTRIES FOUND') + 1010 FORMAT (/5X,'HISTOGRAM NO',I4,' : OUTPUTTING TO TOPDRAWER FILE') + END +C----------------------------------------------------------------------- + SUBROUTINE GTOP1(ID,NEW,HIST,LOG,IDERR) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + LOGICAL NEW,HIST,LOG, DIV + CHARACTER*8 TYPE(4),PLOT + COMMON /GTCOM/ TYPE,NTYPE + PARAMETER (NXMAX=1000) + DIMENSION X(NXMAX),Y(NXMAX),E(NXMAX) + CHARACTER TITLE*60 + IF (ID.GT.NMAX) RETURN + CALL GOPENF + WEIGHT=0.0 + AVGE=0.0 + VARIAN=0.0 + IS=NINT(A(ID+2)) + NX=NINT(A(IS+1)) + IF (NX.GT.NXMAX) THEN + WRITE (6,2000) NXMAX + NX=NXMAX + ENDIF + ISERR=0 + IF (IDERR.NE.0) THEN + ISERR=NINT(A(IDERR+2)) + IF (ISERR .NE. 0) THEN + IF (A(ISERR+5) .GT. 1.5) THEN + WRITE (6,1000) IDERR + ISERR=0 + ELSEIF (NX.NE.NINT(A(ISERR+1))) THEN + WRITE (6,1010) IDERR + ISERR=0 + ELSE + WRITE (6,1020) IDERR + ENDIF + ENDIF + ENDIF + YMIN=1D30 + DO 50 IX=1, NX + 50 IF (A(IS+18+IX).GT.0.0) YMIN=MIN(YMIN,A(IS+18+IX)) + DIV=.FALSE. + DO 100 IX=1, NX + X(IX)=(IX-0.5)*A(IS+4)+A(IS+2) + Y(IX)=A(IS+18+IX) + E(IX)=0 + IF (ISERR.NE.0) E(IX)=A(ISERR+18+IX) + WEIGHT=WEIGHT+Y(IX) + AVGE=AVGE+X(IX)*Y(IX) + VARIAN=VARIAN+X(IX)**2*Y(IX) + IF (LOG .AND. Y(IX).LE.0.0) THEN + WRITE (6,2010) Y(IX),YMIN/1000 + Y(IX)=YMIN/1000 + DIV=.TRUE. + ENDIF + 100 CONTINUE + IF (WEIGHT.EQ.0) WEIGHT=1 + AVGE=AVGE/WEIGHT + VARIAN=VARIAN/WEIGHT-AVGE**2 + CALL GTITLE(ID,TITLE,LENTIT) +C SET UP PAGE... + IF (NEW) THEN + WRITE (21,*) 'NEW FRAME' + WRITE (21,*) 'SET WINDOW X 4.9 9.5 Y 2 9' + WRITE (21,*) 'SET FONT DUPLEX' + IF (LOG) WRITE (21,*) 'SET SCALE Y LOG' + IF (DIV) WRITE (21,*) 'SET LIMIT YMIN',YMIN/10 + WRITE (21,*) 'TITLE TOP ''',TITLE(1:LENTIT),'''' + IF (NINT(A(IS+2)*10).EQ.0.AND.NINT(A(IS+3)*10).EQ.10) THEN + WRITE (21,*) 'SET LIMITS X 0 1' + ELSE + WRITE (21,*) 'SET LIMITS X',SNGL(A(IS+2)),SNGL(A(IS+3)) + ENDIF + WRITE (21,*) 'SET PATT 0.02 0.08' + WRITE (21,*) 'SET ORDER X Y DY' + NTYPE=0 + ENDIF + NTYPE=MOD(NTYPE,4)+1 +C PLOT... + DO 200 I=1, NX + IF ((ABS(X(I)).GE.1E-3.AND.ABS(X(I)).LT.1E5.OR.X(I).EQ.0).AND. + & (ABS(Y(I)).GE.1E-5.AND.ABS(Y(I)).LT.1E3.OR.Y(I).EQ.0)) THEN + WRITE (21,'(F14.5,2F14.7)') X(I),Y(I),E(I) + ELSE + WRITE (21,'(3(1PE14.4))') X(I),Y(I),E(I) + ENDIF + 200 CONTINUE + PLOT=' ' + IF (ISERR.NE.0) PLOT=' ; PLOT' + IF (HIST) THEN + WRITE (21,*) 'HIST ',TYPE(NTYPE),PLOT + ELSE + WRITE (21,*) 'JOIN ',TYPE(NTYPE),PLOT + ENDIF + WRITE (6,1100) A(IS+13) + WRITE (6,1110) A(IS+14) + WRITE (6,1120) A(IS+15) + WRITE (6,1130) AVGE + IF (VARIAN.GE.0.0) THEN + WRITE (6,1140) SQRT(VARIAN) + ELSE + WRITE (6,1150) VARIAN + ENDIF + 1000 FORMAT (5X,'HISTOGRAM NO',I4,' : SHOULD BE 1-DIM BUT ISN''T') + 1010 FORMAT (5X,'HISTOGRAM NO',I4,' : SHOULD BE SAME SIZE BUT ISN''T') + 1020 FORMAT (5X,'HISTOGRAM NO',I4,' : BEING USED FOR ERROR BARS') + 1100 FORMAT (5X,'Total underflow=',G10.3) + 1110 FORMAT (5X,'Total entry =',G10.3) + 1120 FORMAT (5X,'Total overflow =',G10.3) + 1130 FORMAT (5X,'Mean value =',G10.3) + 1140 FORMAT (5X,'Standard deviation=',G10.3) + 1150 FORMAT (5X,'Variance =',G10.3) + 2000 FORMAT (5X,'WARNING: Histogram too big, using first',I4,' bins') + 2010 FORMAT (5X,'WARNING:',G10.3,' cannot be logged. Now using',G10.3) + END +C----------------------------------------------------------------------- + SUBROUTINE GTOP2(ID,NEW) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + LOGICAL NEW + PARAMETER (NZ=250) + DIMENSION Z(NZ,NZ) + CHARACTER TITLE*60 + REAL RAN + DATA ISEED/12345/ + IF (ID.GT.NMAX) RETURN + CALL GOPENF + WEIGHT=0.0 + XAVGE=0.0 + XVARIN=0.0 + YAVGE=0.0 + YVARIN=0.0 + ZMAX=0 + IS=NINT(A(ID+2)) + NX=NINT(A(IS+1)) + NY=NINT(A(IS+5)) + IF (NX.GT.NZ) THEN + WRITE (6,2000) NZ + NX=NZ + ENDIF + IF (NY.GT.NZ) THEN + WRITE (6,2000) NZ + NY=NZ + ENDIF + DO 100 IX=1, NX + DO 100 IY=1, NY + X=(IX-1)*A(IS+4)+A(IS+2) + Y=(IY-1)*A(IS+8)+A(IS+6) + Z(IX,IY)=A(IS+18+NX*(IY-1)+IX) + SMALL=1E-20 + IF (ABS(Z(IX,IY)).LT.SMALL) Z(IX,IY)=0.0 + IF (Z(IX,IY).GT.ZMAX) ZMAX=Z(IX,IY) + WEIGHT=WEIGHT+Z(IX,IY) + XAVGE=XAVGE+X*Z(IX,IY) + XVARIN=XVARIN+X**2*Z(IX,IY) + YAVGE=YAVGE+Y*Z(IX,IY) + YVARIN=YVARIN+Y**2*Z(IX,IY) + 100 CONTINUE + IF (ZMAX.GE.1 .AND. ZMAX.LE.20) THEN + ZMAX=20 + ELSE + WRITE (6,2010) ZMAX/20 + ENDIF + XAVGE=XAVGE/WEIGHT + XVARIN=XVARIN/WEIGHT-XAVGE**2 + YAVGE=YAVGE/WEIGHT + YVARIN=YVARIN/WEIGHT-YAVGE**2 + CALL GTITLE(ID,TITLE,LENTIT) +C SET UP PAGE... + IF (NEW) THEN + WRITE (21,*) 'NEW FRAME' + IF (A(IS+2).EQ.A(IS+6) .AND. A(IS+3).EQ.A(IS+7)) + & WRITE (21,*) 'SET WINDOW X 3.7 10.7 Y 2 9' + WRITE (21,*) 'SET FONT DUPLEX' + WRITE (21,*) 'TITLE TOP ''',TITLE(1:LENTIT),'''' + WRITE (21,*) 'SET LIMITS X',SNGL(A(IS+2)),SNGL(A(IS+3)), + & ' Y',SNGL(A(IS+6)),SNGL(A(IS+7)) + ENDIF +C PLOT... + NN=0 + DO 210 IX=1, NX + DO 210 IY=1, NY + DO 200 IZ=1, NINT(20*Z(IX,IY)/ZMAX) + X=(IX-RAN(ISEED))*A(IS+4)+A(IS+2) + Y=(IY-RAN(ISEED))*A(IS+8)+A(IS+6) + IF ((ABS(X).GE.1E-5.AND.ABS(X).LT.1E3.OR.X.EQ.0).AND. + & (ABS(Y).GE.1E-5.AND.ABS(Y).LT.1E3.OR.Y.EQ.0)) THEN + WRITE (21,'(2F14.7)') X,Y + ELSE + WRITE (21,'(2(1PE14.4))') X,Y + ENDIF + NN=NN+1 + 200 CONTINUE + IF (NN.GT.500) THEN + WRITE (21,*) 'PLOT' + NN=0 + ENDIF + 210 CONTINUE + IF (NN.GT.0) WRITE (21,*) 'PLOT' + WRITE (6,1000) + WRITE (6,1030) (A(IS+J), J=16,18) + WRITE (6,1020) (A(IS+J), J=13,15) + WRITE (6,1010) (A(IS+J), J=10,12) + WRITE (6,1040) XAVGE + IF (XVARIN.GT.0) WRITE (6,1050) SQRT(XVARIN) + WRITE (6,1060) YAVGE + IF (YVARIN.GT.0) WRITE (6,1070) SQRT(YVARIN) + 1000 FORMAT (20X,' X under',5X,' X on grid',5X,' X over') + 1010 FORMAT (5X,' Y under',3(5X,G10.3)) + 1020 FORMAT (5X,' Y on grid',3(5X,G10.3)) + 1030 FORMAT (5X,' Y over',3(5X,G10.3)) + 1040 FORMAT (5X,'Mean x value =',G10.3) + 1050 FORMAT (5X,'X standard deviation=',G10.3) + 1060 FORMAT (5X,'Mean y value =',G10.3) + 1070 FORMAT (5X,'Y standard deviation=',G10.3) + 2000 FORMAT (5X,'WARNING: Histogram too big, using first',I4,' bins') + 2010 FORMAT (5X,'WARNING: Dividing by ',G10.3) + END +C----------------------------------------------------------------------- + SUBROUTINE GTOP +C---OUTPUT ALL HISTOGRAMS TO TOPDRAWER + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + DO 100 ID=1,INT(A(1)+0.5) + IS=A(ID+2)+0.5 + IF(IS.EQ.0.OR.A(IS+9).LT.0.5) GOTO 100 + CALL GTOPDR(ID,.TRUE.,.TRUE.,.FALSE.) + 100 CONTINUE + END +C----------------------------------------------------------------------- + SUBROUTINE GTITLE(ID,TITLE,LENTIT) + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + PARAMETER (NSIZE=200000,NMAX=2000) + COMMON /GBOOK/ A(NSIZE) + CHARACTER TITLE*60 + EQUIVALENCE (REQ,IEQ) + IF (ID.GT.NMAX) RETURN + IS=NINT(A(ID+2)) + NX=NINT(A(IS+1)) + NY=NINT(A(IS+5)) + DO 100 IT=1, 20 + REQ=A(IS+18+NX*NY+IT) + TITLE(3*IT-2:3*IT)=CHAR(IEQ/256**2)//CHAR(MOD(IEQ,256**2) + & /256)//CHAR(MOD(IEQ,256)) + 100 CONTINUE + LENTIT=60 + 200 IF (LENTIT.GT.1 .AND. TITLE(LENTIT:LENTIT).EQ.' ') THEN + LENTIT=LENTIT-1 + GOTO 200 + ENDIF + END +C----------------------------------------------------------------------- + SUBROUTINE GOPENF +C---OPEN TopDrawer FILE + CHARACTER*14 TITLE + LOGICAL OPEN + DATA OPEN/.FALSE./ + IF (OPEN) RETURN + OPEN=.TRUE. + N=0 + TITLE='gtopdraw.top' + OPEN (21,FILE='gtopdraw.top',STATUS='NEW',ERR=10) + WRITE (6,*) 'Using file "gtopdraw.top"' + RETURN + 10 N=N+1 + IF (N.GE.100) STOP + WRITE (6,*) 'Could not open file "',TITLE,'"' + WRITE (TITLE,'(A8,I2.2,A4)') 'gtopdraw',N,'.top' + OPEN (21,FILE=TITLE,STATUS='NEW',ERR=10) + WRITE (6,*) 'Using file "',TITLE,'" instead' + END +C----------------------------------------------------------------------- + BLOCK DATA GTDAT + IMPLICIT INTEGER (I-N) + IMPLICIT DOUBLE PRECISION (A-H,O-Z) + CHARACTER*8 TYPE(4) + COMMON /GTCOM/ TYPE,NTYPE + DATA TYPE,NTYPE/' ','DASH ','DOT-DASH','PATTERN ',0/ + END +C----------------------------------------------------------------------- + double precision function ran(iseed) + implicit none + integer iseed + double precision r(1) + call rangen(1,r) + iseed=iseed+1 + ran=r(1) + end +C-----------------------------------------------------------------------