C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** C*********************************************************************** C*********************************************************************** C C JETSUBS C Version 3.0, 23 February 1991 C Version 3.01, 14 March 1991: C - misc small changes to enhance portability C Version 3.1, 26 February 1992: C - changes in soft integrals to match paper C - a few other changes to match paper and simplify C - common block /ARRAYS/ enlarged C - switch for proton-proton collisions C Version 3.4, 13 June 1996 C - only the name is changed, to indicate by C the name that this is the version to use C with jet3_4.f C - name change is incorportated in PRINTVERSION; C version date left unchanged. C Version 3.4.01, 2 September 1997 C - check for divide by zero incorporated in C variable SINGULARITY. C C*********************************************************************** C C Subroutines for JET. Those that one is likely to want to modify, C which are consequently kept with the main program, are designated C with a (*). C C SUBROUTINE RENO(NRENO) C C SUBROUTINE INTEGRATE2TO2(Y1,Y2,PHI2,P2,JACOBIAN) C Virtual: C FUNCTION PSITILDE4(PROCESS,K) C FUNCTION PSITILDE6NS(PROCESS,K) C Collinear: C FUNCTION EFFPARTON(A,X,SCALE) C FUNCTION ALTARELLI(AOUT,AIN,Z) C FUNCTION ALTRLIPRIME(AOUT,AIN,Z) C FUNCTION GAMMA(NPARTON) C FUNCTION GAMMAPRIME(NPARTON) C Soft: C FUNCTION PSICTILDE(PROCESS,I1,I2,K) C FUNCTION LI2(X) C C SUBROUTINE INTEGRATEA(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF) C FUNCTION FA(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO) C SUBROUTINE INTEGRATEB(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF) C FUNCTION FB(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO) C SUBROUTINE INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF) C FUNCTION F1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUCO) C SUBROUTINE INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF) C FUNCTION F2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUCO) C C Used in INTEGRATE2TO2,INTEGRATEA,INTEGRATEB,INTEGRATE1,INTEGRATE2. C SUBROUTINE JETDEFBORN(Y1,P2,Y2,PHI2, C > OK,PJ,YJ,UNUSED,SVALUE) (*) C SUBROUTINE JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3, C > OK,PJ,YJ,UNUSED,SVALUE) (*) C SUBROUTINE SETMU(PJ,YJ,UNUSED,MUUV,MUCO) (*) C SUBROUTINE BINIT(PJ,YJ,UNUSED,INTEGRAND) (*) C FUNCTION TRIAL(PJ,YJ) (*) C The JETDEF subroutines calculate the 'measured' variables, like C PJ,YJ, given the parton momenta. Then SETMU uses the measured C variables to set the mu values. At this point the physics functions C can be calculated. Finally, the results are smeared a bit over C the physics variables and stored by BINIT. C C Used in FA,FB,F1,F2: C FUNCTION CROSSINGSIGN(PROCESS,PERMINV) C FUNCTION CONVERT(PHI) C FUNCTION RESIDUE(PROCESS,J,N,M,KIN) C SUBROUTINE RESIDUEINIT C LOGICAL FUNCTION MATCH(J,N,M,JJ,NN,MM) C SUBROUTINE PERMUTE C SUBROUTINE PARTONSIN(LMNSTY,XA,XB,MUCO) C FUNCTION LUMINOSITY(LMNSTY,PROCESS,PERMINV) C C Auxiliary functions: C FUNCTION PARTON(NPARTON,X,SCALE) Parton distributions (*) C FUNCTION ALPHAS(Q) Alpha-s C FUNCTION WNUM(NPARTON) Spin and color weight C FUNCTION COLOR(NPARTON) Color charge C FUNCTION THETA(BOOL) Theta function C FUNCTION FEXP(Y) Exponential that doesn't crash for large y C SUBROUTINE PRINTVERSION C C Smearing: C FUNCTION UNSMEAR(M,J,DELTA) C FUNCTION DELTA(M,J) C FUNCTION FACTORIAL(N) C Fitting: C SUBROUTINE FITJET(VALARRAY,ERRARRAY,PJGRID,YJGRID,NPJ,NYJ,NRENO, C > LNPJSMEAR,YJSMEAR,JMAX,KMAX,CTILDE,C) (*) C SUBROUTINE FITVAL(FUN,VAL,DEV,NPT,NPAR,PARAMS,CHISQ) C SUBROUTINE MMIN(FUNGRD,FUNHES,DELVEC,N) C Random numbers: C FUNCTION RANDOM(DUMMY) C SUBROUTINE RANDOMINIT(IRAN) C SUBROUTINE NEWRAN C C*********************************************************************** C*********************************************************************** C*********************************************************************** C*********************************************************************** C SUBROUTINE PRINTVERSION C C INTEGER NIN,NOUT,NWRT COMMON /IOUNIT/ NIN,NOUT,NWRT C WRITE(NOUT,1) 1 FORMAT(' JETSUBS version 3.4, 26 Feb 1992 '/, > ' Version 3.4.01 2 September 1997 ',/50('-')) RETURN END C C*********************************************************************** C SUBROUTINE RENO(NRENO) C INTEGER NRENO C C (Azimuthal angles are defined to lie between - pi and pi) C DOUBLE PRECISION RANDOM DOUBLE PRECISION X1,X2,X3,X4,X5,X6,X7 INTEGER ICALL1,ICALL2 C DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C DOUBLE PRECISION YSCALE,PSCALE COMMON /RENOSCALES/ YSCALE,PSCALE C LOGICAL S2TO2 LOGICAL SAFINITE,SBFINITE,S1FINITE,S2FINITE COMMON /SW2TO3/S2TO2,SAFINITE,SBFINITE,S1FINITE,S2FINITE C DOUBLE PRECISION TINY DOUBLE PRECISION Y1,Y2,Y3,PHI1,PHI2,PHI3,P1,P2,P3,W,XI DOUBLE PRECISION COS0,PHI0,SIN0,X0,Y0,Z0 DOUBLE PRECISION DY1DX1,DP2DX2,DY2DX3,DPHI2DX4,DP3DX5 DOUBLE PRECISION DY2DX1,DP1DX2,DY1DX3,DPHI1DX4 DOUBLE PRECISION DXIDX5,DWDX6,DPHI3DX7 DOUBLE PRECISION DOMEGADX67,DY3PHI3DOMEGA DOUBLE PRECISION CONVERT C DOUBLE PRECISION JACOBIAN,SINGULARITY,COEF DOUBLE PRECISION FEXP C INTEGER IPJ,IYJ INTEGER NPJ,NYJ DOUBLE PRECISION PJGRID(1:30),YJGRID(1:30) COMMON /GRIDS/ PJGRID,YJGRID,NPJ,NYJ C INTEGER PROCESS,NPERM,I,J C DOUBLE PRECISION PI PARAMETER ( PI = 3.141592654 D0) C DOUBLE PRECISION TEMPARRAY(1:30,1:30) DOUBLE PRECISION XSECTARRAY(1:30,1:30),ERRARRAY(1:30,1:30) COMMON /ARRAYS/ TEMPARRAY,XSECTARRAY,ERRARRAY C C Arrays for RESIDUE, initialized by RESIDUEINIT C INTEGER JA(5,5,5),JB(5,5,5),JC(5,5,5),JD(5,5,5) INTEGER SUBSCRIPTA(5,5,5),SUBSCRIPTB(5,5,5) INTEGER SUBSCRIPTC(5,5,5),SUBSCRIPTD(5,5,5) INTEGER PERMA(5,5,5,5),PERMB(5,5,5,5) INTEGER PERMC(5,5,5,5),PERMD(5,5,5,5) COMMON /RESIDUEHELP/ JA,JB,JC,JD, > SUBSCRIPTA,SUBSCRIPTB,SUBSCRIPTC,SUBSCRIPTD, > PERMA,PERMB,PERMC,PERMD C C Permutations for each process. The separate definitions for C each process, with an equivalence statement, are made in order C to keep under 20 continuation lines, for reasons of program C protability. C INTEGER PERMS(30,5,4) INTEGER PERMSA(30,5),PERMSB(30,5),PERMSC(30,5),PERMSD(30,5) EQUIVALENCE (PERMS(1,1,1),PERMSA(1,1)) EQUIVALENCE (PERMS(1,1,2),PERMSB(1,1)) EQUIVALENCE (PERMS(1,1,3),PERMSC(1,1)) EQUIVALENCE (PERMS(1,1,4),PERMSD(1,1)) C C The permutaions will be translated into these arrays for later use. C INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4) COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR C C -- SET UP SUM OVER PERMUTATIONS FOR EACH PROCESS IN ORDER: C C -- FIRST PROCESS A: USE (1,3) <-> (2,4) AND C -- ALSO USE C INVARIANCE (1,2) <-> (3,4) C -- TO MAKE PI(1)1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2, >2,2,2,2,2,2,3,3,4,4,5,5,3,3,4,4,5,5,3,3,4,4,5,5,3,3,4,4,5,5, >3,3,4,4,5,5,2,2,2,2,2,2,4,5,3,5,3,4,4,5,3,5,3,4,4,5,3,5,3,4, >4,5,3,5,3,4,4,5,3,5,3,4,2,2,2,2,2,2,5,4,5,3,4,3,5,4,5,3,4,3, >5,4,5,3,4,3,5,4,5,3,4,3,5,4,5,3,4,3,2,2,2,2,2,2,1,1,1,1,1,1/ C C -- FOR PROCESS B LET PI(3) < PI(4) ALSO USING EXTRA SYMMETRY C DATA PERMSB/ >1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >2,2,2,3,3,4,4,5,5,3,4,5,3,4,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >3,3,4,2,2,2,2,2,2,4,3,3,4,3,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >4,5,5,4,5,3,5,3,4,5,5,4,5,5,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >5,4,3,5,4,5,3,4,3,2,2,2,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ C C -- FOR PROCESS C USE 345 SYMMETRY TO MAKE PI(3) < PI(4) < PI(5) C AND C INVARIANCE TO MAKE PI(1)1,1,1,1,2,2,2,3,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >2,3,4,5,3,4,5,4,5,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >3,2,2,2,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >4,4,3,3,4,3,3,2,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >5,5,5,4,5,5,4,5,4,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ C C -- FOR PROCESS D SYMMETRIC IN EVERYTHING C DATA PERMSD/ >1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, >5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0/ C C Initializations. C NPERMS(1) = 30 NPERMS(2) = 15 NPERMS(3) = 10 NPERMS(4) = 1 C C Initialize logic arrays for RESIDUE C CALL RESIDUEINIT C C Initialize results arrays C DO I = 1,10 DO J = 1,10 XSECTARRAY(I,J) = 0.0D0 ERRARRAY(I,J) = 0.0D0 TEMPARRAY(I,J) = 0.0D0 ENDDO ENDDO C C Initialize permutation arrays C DO PROCESS = 1,4 DO NPERM = 1,NPERMS(PROCESS) DO I = 1,5 J = PERMS(NPERM,I,PROCESS) PERMINVARR(J,NPERM,PROCESS) = I PERMARR(I,NPERM,PROCESS) = J ENDDO ENDDO ENDDO C C Set parameters for choosing kinematic variables C TINY = 1.0 D-10 C C Calculation starts here. C C Start loop over ICALL1, ICALL2. C DO ICALL2 = 1,NRENO DO ICALL1 = 1,1000 C C ----------------------- Do 2-to-2 terms ---------------------------- C IF (S2TO2) THEN C X1 = RANDOM(1.0) X2 = RANDOM(1.0) X3 = RANDOM(1.0) X4 = RANDOM(1.0) C Y1 = YSCALE * DLOG( (X1 + TINY) /(1.0D0 + TINY - X1) ) DY1DX1 = YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1)) C P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS) DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS) > /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2 C Y2 = YSCALE * DLOG( (X3 + TINY) /(1.0D0 + TINY - X3) ) DY2DX3 = YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3)) C PHI2 = 2.0D0 * PI * X4 - PI DPHI2DX4 = 2.0D0 * PI C C JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4 C CALL INTEGRATE2TO2(Y1,Y2,PHI2,P2,JACOBIAN) C ENDIF C C -------------------------- Do Term A ---------------------------- C IF (SAFINITE) THEN C X1 = RANDOM(1.0) X2 = RANDOM(1.0) X3 = RANDOM(1.0) X4 = RANDOM(1.0) X5 = RANDOM(1.0) X6 = RANDOM(1.0) X7 = RANDOM(1.0) C Y1 = YSCALE * DLOG( (X1 + TINY) /(1.0D0 + TINY - X1) ) DY1DX1 = YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1)) C P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS) DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS) > /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2 C Y2 = YSCALE * DLOG( (X3 + TINY) /(1.0D0 + TINY - X3) ) DY2DX3 = YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3)) C PHI2 = 2.0D0 * PI * X4 - PI DPHI2DX4 = 2.0D0 * PI C XI = X5**2 DXIDX5 = 2.0D0 * X5 C W = RTS * X6**2 /(1.0D0 - X6) DWDX6 = RTS*X6*(2.0D0 - X6) /(1.0D0 - X6)**2 C PHI3 = 2.0D0 * PI * X7 - PI DPHI3DX7 = 2.0D0 * PI C JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4*DXIDX5*DWDX6*DPHI3DX7 C SINGULARITY = XI * W C C Proceed as long as SINGULARITY is not too small: C IF (SINGULARITY.GT.1.0D-6) THEN C COEF = JACOBIAN /SINGULARITY C CALL INTEGRATEA(Y1,P2,Y2,PHI2,XI, W, PHI3, COEF) IF (W.LT.RTS) THEN CALL INTEGRATEA(Y1,P2,Y2,PHI2,XI, 0.0D0,PHI3,-COEF) ENDIF IF (XI.LT.(1.0D0 - ( FEXP(Y1) + FEXP(Y2) ) * P2 / RTS)) THEN CALL INTEGRATEA(Y1,P2,Y2,PHI2,0.0D0,W, PHI3,-COEF) IF (W.LT.RTS) THEN CALL INTEGRATEA(Y1,P2,Y2,PHI2,0.0D0,0.0D0,PHI3, COEF) ENDIF ENDIF C C End IF for SINGULARITY not too small. C ENDIF C C End IF for SAFINITE switch on C ENDIF C C -------------------------- Do Term B ---------------------------- C IF (SBFINITE) THEN C X1 = RANDOM(1.0) X2 = RANDOM(1.0) X3 = RANDOM(1.0) X4 = RANDOM(1.0) X5 = RANDOM(1.0) X6 = RANDOM(1.0) X7 = RANDOM(1.0) C Y1 = YSCALE * DLOG( (X1 + TINY) /(1.0D0 + TINY - X1) ) DY1DX1 = YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1)) C P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS) DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS) > /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2 C Y2 = YSCALE * DLOG( (X3 + TINY) /(1.0D0 + TINY - X3) ) DY2DX3 = YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3)) C PHI2 = 2.0D0 * PI * X4 - PI DPHI2DX4 = 2.0D0 * PI C XI = X5**2 DXIDX5 = 2.0D0 * X5 C W = RTS * X6**2 /(1.0D0 - X6) DWDX6 = RTS*X6*(2.0D0 - X6) /(1.0D0 - X6)**2 C PHI3 = 2.0D0 * PI * X7 - PI DPHI3DX7 = 2.0D0 * PI C JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4*DXIDX5*DWDX6*DPHI3DX7 C SINGULARITY = XI * W C C Proceed as long as SINGULARITY is not too small: C IF (SINGULARITY.GT.1.0D-6) THEN C COEF = JACOBIAN /SINGULARITY C CALL INTEGRATEB(Y1,P2,Y2,PHI2,XI, W, PHI3, COEF) IF (W.LT.RTS) THEN CALL INTEGRATEB(Y1,P2,Y2,PHI2,XI, 0.0D0,PHI3,-COEF) ENDIF IF (XI.LT.(1.0D0 - ( FEXP(-Y1) + FEXP(-Y2) ) * P2 / RTS)) THEN CALL INTEGRATEB(Y1,P2,Y2,PHI2,0.0D0,W, PHI3,-COEF) IF (W.LT.RTS) THEN CALL INTEGRATEB(Y1,P2,Y2,PHI2,0.0D0,0.0D0,PHI3, COEF) ENDIF ENDIF C C End IF for SINGULARITY not too small. C ENDIF C C End IF for SBFINITE switch on. C ENDIF C C -------------------------- Do Term 1 ---------------------------- C IF (S1FINITE) THEN C X1 = RANDOM(1.0) X2 = RANDOM(1.0) X3 = RANDOM(1.0) X4 = RANDOM(1.0) X5 = RANDOM(1.0) X6 = RANDOM(1.0) X7 = RANDOM(1.0) C Y2 = YSCALE * DLOG( (X1 + TINY) /(1.0D0 + TINY - X1) ) DY2DX1 = YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1)) C P1 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS) DP1DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS) > /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2 C Y1 = YSCALE * DLOG( (X3 + TINY) /(1.0D0 + TINY - X3) ) DY1DX3 = YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3)) C PHI1 = 2.0D0 * PI * X4 - PI DPHI1DX4 = 2.0D0 * PI C P3 = P1 * X5**2 DP3DX5 = 2.0D0 * P1 * X5 C COS0 = 1.0D0 - 2.0D0 * X6**2 PHI0 = 2.0D0 * PI * X7 - PI DOMEGADX67 = 8.0D0 * PI * X6 C SIN0 = DSQRT(1.0D0 - COS0**2) X0 = COS0 Y0 = SIN0 * DCOS(PHI0) Z0 = SIN0 * DSIN(PHI0) C Y3 = Y1 + 2.0D0 * DLOG( (1.0D0 + Z0)/(1.0D0 - Z0) ) PHI3 = CONVERT( PHI1 + DATAN2(Y0,X0)) C DY3PHI3DOMEGA = 4.0D0/(1.0D0 - Z0**2) C JACOBIAN = DY2DX1*DP1DX2*DY1DX3*DPHI1DX4*DP3DX5 > * DOMEGADX67 * DY3PHI3DOMEGA C SINGULARITY = P3 * ( DCOSH(Y1-Y3) - DCOS(PHI1-PHI3) ) C C Proceed as long as SINGULARITY is not too small: C IF (SINGULARITY.GT.1.0D-6) THEN C COEF = JACOBIAN /SINGULARITY C CALL INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF) CALL INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y1,PHI1,-COEF) IF (P3.LT.0.5D0*P1) THEN CALL INTEGRATE1(Y2,P1,Y1,PHI1,0.0D0,Y3,PHI3,-COEF) CALL INTEGRATE1(Y2,P1,Y1,PHI1,0.0D0,Y1,PHI1,COEF) ENDIF C C End IF for SINGULARITY not too small. C ENDIF C C End IF for S1FINITE switch on. C ENDIF C C -------------------------- Do Term 2 ---------------------------- C IF (S2FINITE) THEN C X1 = RANDOM(1.0) X2 = RANDOM(1.0) X3 = RANDOM(1.0) X4 = RANDOM(1.0) X5 = RANDOM(1.0) X6 = RANDOM(1.0) X7 = RANDOM(1.0) C Y1 = YSCALE * DLOG( (X1 + TINY) /(1.0D0 + TINY - X1) ) DY1DX1 = YSCALE*(1.0D0+2.0D0*TINY)/((X1+TINY)*(1.0D0+TINY-X1)) C P2 = PSCALE * X2 /(1.0D0 - X2 + 2.0D0*PSCALE/RTS) DP2DX2 = PSCALE * (1.0D0 + 2.0D0*PSCALE/RTS) > /(1.0D0 - X2 + 2.0D0*PSCALE/RTS)**2 C Y2 = YSCALE * DLOG( (X3 + TINY) /(1.0D0 + TINY - X3) ) DY2DX3 = YSCALE*(1.0D0+2.0D0*TINY)/((X3+TINY)*(1.0D0+TINY-X3)) C PHI2 = 2.0D0 * PI * X4 - PI DPHI2DX4 = 2.0D0 * PI C P3 = P2 * X5**2 DP3DX5 = 2.0D0 * P2 * X5 C COS0 = 1.0D0 - 2.0D0 * X6**2 PHI0 = 2.0D0 * PI * X7 - PI DOMEGADX67 = 8.0D0 * PI * X6 C SIN0 = DSQRT(1.0D0 - COS0**2) X0 = COS0 Y0 = SIN0 * DCOS(PHI0) Z0 = SIN0 * DSIN(PHI0) C Y3 = Y2 + 2.0D0 * DLOG( (1.0D0 + Z0)/(1.0D0 - Z0) ) PHI3 = CONVERT( PHI2 + DATAN2(Y0,X0)) C DY3PHI3DOMEGA = 4.0D0/(1.0D0 - Z0**2) C JACOBIAN = DY1DX1*DP2DX2*DY2DX3*DPHI2DX4*DP3DX5 > * DOMEGADX67 * DY3PHI3DOMEGA C SINGULARITY = P3 * ( DCOSH(Y2-Y3) - DCOS(PHI2-PHI3) ) C C Proceed as long as SINGULARITY is not too small: C IF (SINGULARITY.GT.1.0D-6) THEN C COEF = JACOBIAN /SINGULARITY C CALL INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF) CALL INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y2,PHI2,-COEF) IF (P3.LT.0.5D0*P2) THEN CALL INTEGRATE2(Y1,P2,Y2,PHI2,0.0D0,Y3,PHI3,-COEF) CALL INTEGRATE2(Y1,P2,Y2,PHI2,0.0D0,Y2,PHI2,COEF) ENDIF C C End IF for SINGULARITY not too small. C ENDIF C C End IF for S2FINITE switch on. C ENDIF C C End loop over ICALL1 = 1,1000 C ENDDO C C Update arrays C DO IPJ = 1,NPJ DO IYJ = 1,NYJ XSECTARRAY(IPJ,IYJ) = XSECTARRAY(IPJ,IYJ) > + (TEMPARRAY(IPJ,IYJ)/1000)/NRENO ERRARRAY(IPJ,IYJ) = ERRARRAY(IPJ,IYJ) > + (TEMPARRAY(IPJ,IYJ)/1000)**2/NRENO TEMPARRAY(IPJ,IYJ) = 0.0D0 ENDDO ENDDO C C End loop over ICALL2 = 1,NRENO C ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C C Routines for 2-to-2 integrals C C*********************************************************************** C*********************************************************************** C SUBROUTINE INTEGRATE2TO2(Y1,Y2,PHI2,P2,JACOBIAN) C DOUBLE PRECISION Y1,Y2,PHI2,P2,JACOBIAN C C This function calculates the integrand multiplying C DY1 DP2 DY2 DPHI2 C in the calculation of the 2 -> 2 cross section C DOUBLE PRECISION INTEGRAND2TO2 C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C Jet parameters LOGICAL OK DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE C Switches LOGICAL SBORN,SVIRTUAL LOGICAL SACOLL,SASOFT LOGICAL SBCOLL,SBSOFT LOGICAL S1COLL,S1SOFT LOGICAL S2COLL,S2SOFT COMMON /SW2TO2/ SBORN,SVIRTUAL, > SACOLL,SASOFT,SBCOLL,SBSOFT, > S1COLL,S1SOFT,S2COLL,S2SOFT LOGICAL STYPE(4) COMMON /SWTYPE/ STYPE C Switch for p-pbar collisions (if true) or p-p collisions (if false) LOGICAL PPBAR COMMON /BEAMS/ PPBAR C Common block for mu values for 2-to-2 subroutines DOUBLE PRECISION MUUV,MUCO COMMON /MUVALS/ MUUV,MUCO C Loop controls,permutations,parton labels INTEGER PERMS(12,4,4),NPERMS(4),PERMINV(4),PERM(4) INTEGER PROCESS,NPERM,FLAVORS1,FLAVORS2,A1,A2 INTEGER N1,N2,I,J C Parton distributions DOUBLE PRECISION PRTNA(-6:6),PRTNB(-6:6) DOUBLE PRECISION EFFPRTNA(-6:6),EFFPRTNB(-6:6) C Useful parameters DOUBLE PRECISION PI,V,N,LOG2 PARAMETER ( PI = 3.141592654 D0) PARAMETER ( V = 8.0 D0) PARAMETER ( N = 3.0 D0) PARAMETER ( LOG2 = 0.693147181 D0) C Units for error reports: INTEGER NIN,NOUT,NWRT COMMON /IOUNIT/ NIN,NOUT,NWRT C Variables for collinear integration DOUBLE PRECISION EFFPARTON C Parton flavors INTEGER AA,AB C Kinematical variables DOUBLE PRECISION SHAT,UHAT,THAT DOUBLE PRECISION XA,XB DOUBLE PRECISION PIN(4,4),K(4,4) DOUBLE PRECISION FEXP,EY1,EY2 C Ellis and Sexton scale parameter DOUBLE PRECISION QES COMMON /KEITHSQ/ QES C QCD functions DOUBLE PRECISION D,D0,ALPHAS DOUBLE PRECISION L,LEFFA,LEFFB DOUBLE PRECISION PSI4,PSITILDE4,SIGN DOUBLE PRECISION PSI6NS,PSITILDE6NS DOUBLE PRECISION GAMMAPRIME DOUBLE PRECISION PARTON,WNUM,COLOR,GAMMA DOUBLE PRECISION JAB,JA1,JA2,JBA,JB1,JB2,J1A,J1B,J12,J2A,J2B,J21 DOUBLE PRECISION PSICAB,PSICA1,PSICA2,PSICB1,PSICB2,PSIC12 DOUBLE PRECISION PSICTILDE DOUBLE PRECISION LI2 C Precalculated logs DOUBLE PRECISION LOGS,LOGT,LOGU,LOGMUSQ,LOG16P2SQ DOUBLE PRECISION LOGXA,LOGXB C C Each column is one of the permutations we need. C DATA PERMS/1,1,1,1,1,1,2,2,2,2,3,3, > 2,2,3,3,4,4,3,3,4,4,4,4, > 3,4,2,4,2,3,1,4,1,3,1,2, > 4,3,4,2,3,2,4,1,3,1,2,1, C > 1,1,1,2,2,3,0,0,0,0,0,0, > 2,3,4,3,4,4,0,0,0,0,0,0, > 3,2,2,1,1,1,0,0,0,0,0,0, > 4,4,3,4,3,2,0,0,0,0,0,0, C > 1,1,1,2,2,2,3,3,3,4,4,4, > 2,3,4,1,3,4,1,2,4,1,2,3, > 3,2,2,3,1,1,2,1,1,2,1,1, > 4,4,3,4,4,3,4,4,2,3,3,2, C > 1,0,0,0,0,0,0,0,0,0,0,0, > 2,0,0,0,0,0,0,0,0,0,0,0, > 3,0,0,0,0,0,0,0,0,0,0,0, > 4,0,0,0,0,0,0,0,0,0,0,0/ C DATA NPERMS/12,6,12,1/ C C -- Initialize computation ------------ C C CALL JETDEFBORN(Y1,P2,Y2,PHI2,OK,JETVAR1,JETVAR2,JETVAR3,SVALUE) IF(.NOT.OK) THEN RETURN ENDIF CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUCO) C EY1 = FEXP(Y1) EY2 = FEXP(Y2) XA = (EY1 + EY2)* P2/RTS XB = (1.0D0/EY1 + 1.0D0/EY2)* P2/RTS SHAT = P2**2 * ( 2.0D0 + EY2/EY1 + EY1/EY2 ) THAT = - P2**2 * (1.0D0 + EY2/EY1) UHAT = - P2**2 * (1.0D0 + EY1/EY2) C IF( (XA.GT.1.0D0).OR.(XB.GT.1.0D0))THEN RETURN ENDIF C QES = RTS/10.0D0 C C Momentum matrix (incoming) C PIN(1,1) = 0.0D0 PIN(1,2) = SHAT /2.0D0 PIN(1,3) = THAT /2.0D0 PIN(1,4) = UHAT /2.0D0 PIN(2,1) = SHAT /2.0D0 PIN(2,2) = 0.0D0 PIN(2,3) = UHAT /2.0D0 PIN(2,4) = THAT /2.0D0 PIN(3,1) = THAT /2.0D0 PIN(3,2) = UHAT /2.0D0 PIN(3,3) = 0.0D0 PIN(3,4) = SHAT /2.0D0 PIN(4,1) = UHAT /2.0D0 PIN(4,2) = THAT /2.0D0 PIN(4,3) = SHAT /2.0D0 PIN(4,4) = 0.0D0 C D0 = ALPHAS(MUUV)**2 * P2 /2.0D0 /RTS**4 D = D0 * ALPHAS(MUUV)/2.0D0/PI C LOGS = DLOG(SHAT /QES**2) LOGT = DLOG(-THAT /QES**2) LOGU = DLOG(-UHAT /QES**2) LOGMUSQ = DLOG(MUCO**2 /QES**2) LOGXA = DLOG( (XA/(1.0D0 - XA))**2 ) LOGXB = DLOG( (XB/(1.0D0 - XB))**2 ) LOG16P2SQ = DLOG(16.0D0*P2**2/QES**2) C C The tilde J_MN soft integrals C J21 = 0.25D0 * (2.0D0 * LOG2 - LOGS)**2 > - 0.25D0 * (2.0D0 * LOGS - LOGT - LOGU)**2 > + (LOGS - LOGT)*(LOGS - LOGU) C J2A = 0.25D0 * (2.0D0 * LOG2 - LOGU)**2 > - 0.25D0 * (LOGS - LOGT)**2 > - LOG2 * DLOG(2.0D0 + UHAT/SHAT) C J2B = 0.25D0 * (2.0D0 * LOG2 - LOGT)**2 > - 0.25D0 * (LOGS - LOGU)**2 > - LOG2 * DLOG(2.0D0 + THAT/SHAT) C J12 = J21 J1A = J2B J1B = J2A C JAB = 0.25D0 * ( LOGXA - LOGS )**2 + PI**2 / 12.0D0 C JA1 = 0.25D0 * (LOGXA - LOGT)**2 > - 0.25D0 * DLOG(0.5D0 + 0.25D0 * THAT / SHAT)**2 > + LOG2**2 + PI**2 / 12.0D0 > - 0.5D0 * LI2( (SHAT + THAT)/(2.0D0*SHAT + THAT) ) C JA2 = 0.25D0 * (LOGXA - LOGU)**2 > - 0.25D0 * DLOG(0.5D0 + 0.25D0 * UHAT / SHAT)**2 > + LOG2**2 + PI**2 / 12.0D0 > - 0.5D0 * LI2( (SHAT + UHAT)/(2.0D0*SHAT + UHAT) ) C C JBA = 0.25D0 * ( LOGXB - LOGS )**2 + PI**2 / 12.0D0 C JB1 = JA2 - 0.25D0 * (LOGXA - LOGU)**2 > + 0.25D0 * (LOGXB - LOGU)**2 C JB2 = JA1 - 0.25D0 * (LOGXA - LOGT)**2 > + 0.25D0 * (LOGXB - LOGT)**2 C C ---- Calculate parton distributions. For p-pbar collisions, exchange C partons -> antipartons for hadron B. C DO AA = -NFL,NFL PRTNA(AA) = PARTON(AA,XA,MUCO)/(XA*WNUM(AA)) EFFPRTNA(AA) = EFFPARTON(AA,XA,MUCO)/(XA*WNUM(AA)) ENDDO C IF (PPBAR) THEN DO AB = -NFL,NFL PRTNB(AB) = PARTON(-AB,XB,MUCO)/(XB*WNUM(AB)) EFFPRTNB(AB) = EFFPARTON(-AB,XB,MUCO)/(XB*WNUM(AB)) ENDDO ELSE DO AB = -NFL,NFL PRTNB(AB) = PARTON(AB,XB,MUCO)/(XB*WNUM(AB)) EFFPRTNB(AB) = EFFPARTON(AB,XB,MUCO)/(XB*WNUM(AB)) ENDDO ENDIF C C ----- Initialize sum: C INTEGRAND2TO2 = 0.0D0 C C ----- Sum over processes A,B,C,D (PROCESS = 1,2,3,4) ------- C DO PROCESS = 1,4 IF (STYPE(PROCESS)) THEN C C ----- Sum over independent permutations the partons C DO NPERM = 1,NPERMS(PROCESS) C C Calculate inverse permutation C DO I = 1,4 J = PERMS(NPERM,I,PROCESS) PERMINV(J) = I PERM(I) = J ENDDO C C Set K matrix as permutation of P matrix: C K(I) = P(J) with J = PERM(I). C DO I = 1,4 DO J = 1,4 K(I,J) = PIN(PERMS(NPERM,I,PROCESS), > PERMS(NPERM,J,PROCESS)) ENDDO ENDDO C C Here we calculate functions that don't depend on the choice C of quark flavors. Note that these are the 'tilde' functions, C so we need a crossing sign when we use them. C PSI4 = PSITILDE4(PROCESS,K) PSI6NS = PSITILDE6NS(PROCESS,K) PSICAB = PSICTILDE(PROCESS,PERMINV(1),PERMINV(2),K) PSICA1 = PSICTILDE(PROCESS,PERMINV(1),PERMINV(3),K) PSICA2 = PSICTILDE(PROCESS,PERMINV(1),PERMINV(4),K) PSICB1 = PSICTILDE(PROCESS,PERMINV(2),PERMINV(3),K) PSICB2 = PSICTILDE(PROCESS,PERMINV(2),PERMINV(4),K) PSIC12 = PSICTILDE(PROCESS,PERMINV(3),PERMINV(4),K) C C Find how many flavors for incoming partons to sum over. C IF (PROCESS.EQ.1) THEN FLAVORS1 = NFL FLAVORS2 = NFL ELSE IF (PROCESS.EQ.2) THEN FLAVORS1 = NFL FLAVORS2 = 1 ELSE IF (PROCESS.EQ.3) THEN FLAVORS1 = NFL FLAVORS2 = 1 ELSE IF (PROCESS.EQ.4) THEN FLAVORS1 = 1 FLAVORS2 = 1 ENDIF C C ----- Sum over particular parton flavors to calculate luminosities C L = 0.0D0 LEFFA = 0.0D0 LEFFB = 0.0D0 C DO N1 = 1,FLAVORS1 DO N2 = 1,FLAVORS2 IF (.NOT.((PROCESS.EQ.1).AND.(N1.EQ.N2))) THEN C C Find flavors for incoming partons: a(J) = a_0(I) with J = PERM(I) C and set outgoing flavors to 0 (gluon) or 101 (generic C quark or antiquark). C IF (PROCESS.EQ.1) THEN IF (PERM(1).EQ.1) THEN AA = N1 ELSE IF (PERM(2).EQ.1) THEN AA = N2 ELSE IF (PERM(3).EQ.1) THEN AA = - N1 ELSE IF (PERM(4).EQ.1) THEN AA = - N2 END IF IF (PERM(1).EQ.2) THEN AB = N1 ELSE IF (PERM(2).EQ.2) THEN AB = N2 ELSE IF (PERM(3).EQ.2) THEN AB = - N1 ELSE IF (PERM(4).EQ.2) THEN AB = - N2 END IF A1 = 101 A2 = 101 ELSE IF (PROCESS.EQ.2) THEN IF (PERM(1).EQ.1) THEN AA = N1 ELSE IF (PERM(2).EQ.1) THEN AA = N1 ELSE IF (PERM(3).EQ.1) THEN AA = - N1 ELSE IF (PERM(4).EQ.1) THEN AA = - N1 END IF IF (PERM(1).EQ.2) THEN AB = N1 ELSE IF (PERM(2).EQ.2) THEN AB = N1 ELSE IF (PERM(3).EQ.2) THEN AB = - N1 ELSE IF (PERM(4).EQ.2) THEN AB = - N1 END IF A1 = 101 A2 = 101 ELSE IF (PROCESS.EQ.3) THEN IF (PERM(1).EQ.1) THEN AA = N1 ELSE IF (PERM(2).EQ.1) THEN AA = - N1 ELSE IF (PERM(3).EQ.1) THEN AA = 0 ELSE IF (PERM(4).EQ.1) THEN AA = 0 END IF IF (PERM(1).EQ.2) THEN AB = N1 ELSE IF (PERM(2).EQ.2) THEN AB = - N1 ELSE IF (PERM(3).EQ.2) THEN AB = 0 ELSE IF (PERM(4).EQ.2) THEN AB = 0 END IF IF (PERM(1).EQ.3) THEN A1 = 101 ELSE IF (PERM(2).EQ.3) THEN A1 = 101 ELSE IF (PERM(3).EQ.3) THEN A1 = 0 ELSE IF (PERM(4).EQ.3) THEN A1 = 0 END IF IF (PERM(1).EQ.4) THEN A2 = 101 ELSE IF (PERM(2).EQ.4) THEN A2 = 101 ELSE IF (PERM(3).EQ.4) THEN A2 = 0 ELSE IF (PERM(4).EQ.4) THEN A2 = 0 END IF ELSE IF (PROCESS.EQ.4) THEN AA = 0 AB = 0 A1 = 0 A2 = 0 ENDIF C C Calculate luminosity factors C L = L + PRTNA(AA)*PRTNB(AB) LEFFA = LEFFA + EFFPRTNA(AA)*PRTNB(AB) LEFFB = LEFFB + PRTNA(AA)*EFFPRTNB(AB) C C Here we close the loops for flavor sums, having calculated the C luminosity factors. In the following parts of the calculation, C the indices AA,AB,A1,A2 are used, but the calculation needs only C to remember if Ax.EQ.0 for gluon or Ax.NE.0 for quark or antiquark. C As we exit these loops, the Ax are correctly set for this purpose. C C ----- Close IF that omits the case N1=N2 in process A. ENDIF C ----- Close DO for sum over flavors N2. ENDDO C ----- Close DO for sum over flavors N1. ENDDO C C Get crossing sign C IF (((AA.EQ.0).AND.(AB.NE.0)).OR. > ((AB.EQ.0).AND.(AA.NE.0)) ) THEN SIGN = - 1.0D0 ELSE SIGN = 1.0D0 ENDIF C C Calculate contribution from Born graph C IF (SBORN) THEN INTEGRAND2TO2 = INTEGRAND2TO2 + D0 * L * SIGN * PSI4 ENDIF C C Calculate contribution from virtual graphs C IF (SVIRTUAL) THEN INTEGRAND2TO2 = INTEGRAND2TO2 + D * L * SIGN * PSI6NS ENDIF C C Calculate contribution from term A C IF(SACOLL) THEN C INTEGRAND2TO2 = INTEGRAND2TO2 > + D * SIGN * PSI4 > *( - L * LOGMUSQ > *( GAMMA(AA) - COLOR(AA) *LOGXA ) > + LEFFA ) C ENDIF C IF(SASOFT) THEN INTEGRAND2TO2 = INTEGRAND2TO2 > + D * L * SIGN > *( PSICAB * JAB + PSICA1 * JA1 + PSICA2 * JA2) C ENDIF C C Calculate contribution from term B C IF(SBCOLL) THEN C INTEGRAND2TO2 = INTEGRAND2TO2 > + D * SIGN * PSI4 > *( - L * LOGMUSQ > *( GAMMA(AB) - COLOR(AB) * LOGXB ) > + LEFFB ) C ENDIF C IF(SBSOFT) THEN INTEGRAND2TO2 = INTEGRAND2TO2 > + D * L * SIGN > *( PSICAB * JBA + PSICB1 * JB1 + PSICB2 * JB2) C ENDIF C C C Calculate contribution from term 1 C C IF(S1COLL) THEN C INTEGRAND2TO2 = INTEGRAND2TO2 > + D * L * SIGN * PSI4 > *( - LOG16P2SQ > *( GAMMA(A1) - 2.0D0 * COLOR(A1) * LOG2 ) > + GAMMAPRIME(A1) - 2.0D0 * COLOR(A1) * LOG2**2 ) C ENDIF C IF(S1SOFT) THEN INTEGRAND2TO2 = INTEGRAND2TO2 > + D * L * SIGN > *( PSICA1 * J1A + PSICB1 * J1B + PSIC12 * J12) C ENDIF C C C Calculate contribution from term 2 C IF(S2COLL) THEN C INTEGRAND2TO2 = INTEGRAND2TO2 > + D * L * SIGN * PSI4 > *( - LOG16P2SQ > *( GAMMA(A2) - 2.0D0 * COLOR(A2) * LOG2 ) > + GAMMAPRIME(A2) - 2.0D0 * COLOR(A2) * LOG2**2 ) C ENDIF C IF(S2SOFT) THEN INTEGRAND2TO2 = INTEGRAND2TO2 > + D * L * SIGN > *( PSICA2 * J2A + PSICB2 * J2B + PSIC12 * J21) C ENDIF C C ---- Done! C C ----- Close DO for sum over permutations. ENDDO C ----- Close IF that omits process if switch is off. ENDIF C ----- Close DO for sum over processes A,B,C,D. ENDDO C INTEGRAND2TO2 = INTEGRAND2TO2*JACOBIAN*SVALUE CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND2TO2) RETURN END C C*********************************************************************** C -------------------------------------------- C Functions for Born and Virtual Contributions C -------------------------------------------- C*********************************************************************** C DOUBLE PRECISION FUNCTION PSITILDE4(PROCESS,K) C INTEGER PROCESS DOUBLE PRECISION K(4,4) C C This function gives the psitilde^(4) functions from E.S. C DOUBLE PRECISION S,T,U DOUBLE PRECISION N,V PARAMETER ( V = 8.0 D0) PARAMETER ( N = 3.0 D0) C S = 2.0D0 * K(1,2) T = 2.0D0 * K(1,3) U = 2.0D0 * K(1,4) C IF (PROCESS.EQ.1) THEN PSITILDE4 = 2.0D0 * V * (S**2 + U**2) / T**2 ELSEIF (PROCESS.EQ.2) THEN PSITILDE4 = 2.0D0 * V * (S**2 + U**2) / T**2 > + 2.0D0 * V * (S**2 + T**2) / U**2 > - 4.0D0 * V / N * S**2 / U / T ELSEIF (PROCESS.EQ.3) THEN PSITILDE4 = 2.0D0 * V / N > *( V / U / T - 2.0D0 * N**2 / S**2 ) > *( T**2 + U**2 ) ELSEIF (PROCESS.EQ.4) THEN PSITILDE4 = 4.0D0 * V * N**2 > *( S**2 + T**2 + U**2 ) > /S**2 /T**2 /U**2 > *( S**4 + T**4 + U**4 ) ELSE PSITILDE4 = 0.0D0 ENDIF RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION PSITILDE6NS(PROCESS,K) C INTEGER PROCESS DOUBLE PRECISION K(4,4) C C This function gives the psitilde^(6)_NS functions calculated C by subtracting the singular pieces from from E.S. C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R DOUBLE PRECISION MUUV,MUCO COMMON /MUVALS/ MUUV,MUCO DOUBLE PRECISION QES COMMON /KEITHSQ/ QES DOUBLE PRECISION PI,N,V PARAMETER ( PI = 3.141592654 D0) PARAMETER ( V = 8.0 D0) PARAMETER ( N = 3.0 D0) DOUBLE PRECISION S,T,U DOUBLE PRECISION THETA DOUBLE PRECISION LS,LT,LU,LMU,L2S,L2T,L2U C S = 2.0D0 * K(1,2) T = 2.0D0 * K(1,3) U = 2.0D0 * K(1,4) C C Here are the Ellis and Sexton log functions C LS = DLOG(DABS(S/QES**2)) LT = DLOG(DABS(T/QES**2)) LU = DLOG(DABS(U/QES**2)) LMU = DLOG(DABS(MUUV**2/QES**2)) L2S = LS**2 - PI**2 * THETA(S.GT.0.0D0) L2T = LT**2 - PI**2 * THETA(T.GT.0.0D0) L2U = LU**2 - PI**2 * THETA(U.GT.0.0D0) C C Find what type, then calculate C IF (PROCESS.EQ.1) THEN PSITILDE6NS = > V/(9.0D0*N*T**2) > *(- 9*PI**2*(N**2-4 )*(S**2 - U**2) > + 2*(72 + 13*N**2 - 10*N*NFLAVOR + 9*N**2*PI**2) > *(S**2 + U**2) > + 6*N*(11*N - 2*NFLAVOR)*(S**2 + U**2)*LMU > - 36*T*U*LS > + 3*(- 6*S**2 - 30*U**2 + 4*N*NFLAVOR*(S**2 + U**2) > - N**2*(7*S**2 + 3*T**2 + U**2))*LT > - 18*S*T*(N**2 -2)*LU > - 36*(3*S**2 + U**2)*LS*LT > - 18*(N**2 -2)*(S**2 + 3*U**2)*LT*LU > + 18*(S**2 - U**2)*L2S > + 9*(N**2*(S**2 + 3*U**2) + 2*(3*S**2 - U**2))*L2T > - 9*(N**2 - 2)*(S**2 - U**2)*L2U) ELSEIF (PROCESS.EQ.2) THEN C PSITILDE6NS = > 40*N*NFLAVOR*S**2*T*U - 4*N**2*S**2*T*U*(13 + 9*PI**2) > - 9*T*U*(32*S**2 + PI**2*(S**2 + 2*T*U)) > + 36*N*(S**2*(4 + PI**2) * (T**2 + U**2) > + (4 - PI**2) * (T**4 + U**4)) > + N**3*(S**2*(26 + 9*PI**2)*(T**2 + U**2) > + (26 + 27*PI**2)*(T**4 + U**4)) > - 20*N**2*NFLAVOR*(T**4 + U**4 + S**2*(T**2 + U**2)) C PSITILDE6NS = PSITILDE6NS > + 6*N*(11*N - 2*NFLAVOR) > *(-2*S**2*T*U + N*(T**4 + U**4 + S**2*(T**2 + U**2))) *LMU > - 36*N*T*U*(T**2 + U**2)*LS > - 6*U*(- 2*N**2*S**2*T + 2*N*NFLAVOR*S**2*T > - 2*N**2*NFLAVOR*U*(S**2 + U**2) > + 6*S*T*(T + 2*U) + N**3*(-3*T**3 + 2*T**2*U > + 7*T*U**2 + 4*U**3) > + 3*N*(2*T**3 + 3*T**2*U + 2*T*U**2 + 6*U**3)) *LT > - 6*T*(- 2*N**2*S**2*U + 2*N*NFLAVOR*S**2*U > - 2*N**2*NFLAVOR*T*(S**2 + T**2) > + 6*S*U*(2*T + U) > + N**3*(4*T**3 + 7*T**2*U + 2*T*U**2 - 3*U**3) > + 3*N*(6*T**3 + 2*T**2*U + 3*T*U**2 + 2*U**3)) * LU C PSITILDE6NS = PSITILDE6NS > - 36*U*(- S**2*T - N**2*S**2*T + N*U*(3*S**2 + U**2)) *LS*LT > - 36*T*(- S**2*U - N**2*S**2*U + N*T*(3*S**2 + T**2)) *LS*LU > - 18*( 2*N*(-2 + N**2)*(T**2 - T*U + U**2) > *(2*T**2 + 3*T*U + 2*U**2) > + T*U*(3*T**2 + 4*T*U + 3*U**2)) * LT*LU > + 18*N*T*U*(S**2 + T**2 + U**2) * L2S > + 9*U*(- 2*N**2*S**2*T - 4*N*(S**3 - S*T**2 - T**3) > - T*(3*T**2 + 8*T*U + 3*U**2) > - 2*N**3*(T**3 - T*U**2 - 2*U**3)) * L2T > + 9*T*(- 2*N**2*S**2*U - 4*N*(S**3 - S*U**2 - U**3) > - U*(3*T**2 + 8*T*U + 3*U**2) > + 2*N**3*(2*T**3 + T**2*U - U**3)) * L2U C PSITILDE6NS = V/(9*N**2*T**2*U**2) * PSITILDE6NS C ELSEIF (PROCESS.EQ.3) THEN C PSITILDE6NS = > 3*S**2*( -7*(T**2 + U**2) + PI**2*(3*S**2 -2*T*U) ) > + 3*N**2*( 14*(T**4 + U**4) + T*U*(13*S**2 + 4*T*U ) > + PI**2*(T**3*U + 2*T**2*U**2 + T*U**3) ) > - 3*N**4*( 7*(T**4 + U**4)+ T*U*(S**2 + 10*T*U) > - PI**2*(T - U)**2*(T**2 + U**2) ) > + 2*N*(11*N - 2*NFLAVOR)*(-S**2 + N**2*(T**2 + U**2)) > *(T**2 + U**2) * LMU > - 3*T*U*( 4*T**2 + 8*T*U + 4*U**2 > + N**2*(-1 + N**2)*(T**2 - 10*T*U + U**2) ) * LS > + 3*S*U*(S + N**2*U)*(2*T + 3*U + N**2*(2*T - 3*U)) * LT > + 3*S*T*(S + N**2*T)*(3*T + 2*U + N**2*(-3*T + 2*U)) * LU C PSITILDE6NS = PSITILDE6NS > - 6*( N**4*U*(U - T)*(2*T**2 + T*U + U**2) > + S*(S - N**2*T)*(2*T**2 + 2*T*U + U**2) ) * LS*LT > - 6*( N**4*T*(T - U)*(T**2 + T*U + 2*U**2) > + S*(S - N**2*U)*(T**2 + 2*T*U + 2*U**2) ) * LS*LU > + 12*N**2*S**2*(T**2 + U**2) * LT*LU > + 3*( 2*S**4 - 2*N**4*T*U*(T**2 + U**2) > + N**2*(T**2 + T*U + 2*U**2)*(2*T**2 + T*U + U**2) ) *L2S > - 3*S*(- N**4*U*(2*T**2 - T*U + U**2) > + (N**2*T - S)*(2*T**2 + 2*T*U + U**2) ) *L2T > - 3*S*(- N**4*T*(2*U**2 - T*U + T**2 ) > + (N**2*U - S)*( 2*U**2 + 2*T*U + T**2)) *L2U C PSITILDE6NS = V/(3.0D0*N**2*S**2*T*U) * PSITILDE6NS C ELSEIF (PROCESS.EQ.4) THEN C PSITILDE6NS = > 2*N**2*NFLAVOR*( - (66 + 27*PI**2)*S**2*T**2*U**2 > + 40*(S**6 + T**6 + U**6)) > + 2*N**3*( 6*(125 - 27*PI**2)*S**2*T**2*U**2 > - 4*(67 - 9*PI**2)*(S**6 + T**6 + U**6)) > + 6*N**2*(11*N - 2*NFLAVOR)*(S**2 + T**2 + U**2)**3 *LMU > + 6*N**2*T*U*(S**2 + T**2 + U**2) > *( NFLAVOR *(5*T**2 + 2*T*U + 5*U**2) > - 2*N*(7*T**2 - 8*T*U + 7*U**2) ) *LS C PSITILDE6NS = PSITILDE6NS > + 6*N**2*S*U*(S**2 + T**2 + U**2) > *( NFLAVOR *(5*S**2 + 2*S*U + 5*U**2) > - 2*N*(7*S**2 - 8*S*U + 7*U**2) ) *LT > + 6*N**2*S*T*(S**2 + T**2 + U**2) > *( NFLAVOR* (5*S**2 + 2*S*T + 5*T**2) > - 2*N*(7*S**2 - 8*S*T + 7*T**2) ) *LU > - 36*N**2*U**2*( NFLAVOR*S*T*(S**2 + T**2) > + 2*N*(2*S**4 + 2*S**3*T + 3*S**2*T**2 + 2*S*T**3 > + 2*T**4))*LS*LT C PSITILDE6NS = PSITILDE6NS > - 36*N**2*S**2*( NFLAVOR*T*U*(T**2 + U**2) > + 2*N*(2*T**4 + 2*T**3*U + 3*T**2*U**2 + 2*T*U**3 > + 2*U**4))*LT*LU > - 36*N**2*T**2*( NFLAVOR*S*U*(S**2 + U**2) > + 2*N*(2*S**4 + 2*S**3*U + 3*S**2*U**2 + 2*S*U**3 > + 2*U**4))*LS*LU > + 18*N**2*S**2*T*U*(4*N*(T**2 + U**2) > - NFLAVOR*(T**2 + 3*T*U + U**2))*L2S > + 18*N**2*S*T**2*U*(4*N*(S**2 + U**2) > - NFLAVOR*(S**2 + 3*S*U + U**2))*L2T > + 18*N**2*S*T*U**2*(4*N*(S**2 + T**2) > - NFLAVOR*(S**2 + 3*S*T + T**2))*L2U C PSITILDE6NS = V/(9.0D0*S**2*T**2*U**2) * PSITILDE6NS C ELSE PSITILDE6NS = 0.0D0 ENDIF RETURN END C C*********************************************************************** C ------------------------------------- C Functions for Collinear Contributions C ------------------------------------- C*********************************************************************** C DOUBLE PRECISION FUNCTION EFFPARTON(A,X,SCALE) C INTEGER A DOUBLE PRECISION X,SCALE C DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL C DOUBLE PRECISION RANDOM,Z,DZDX1 INTEGER APRIME DOUBLE PRECISION FCOLL DOUBLE PRECISION PARTON,COLOR DOUBLE PRECISION ALTARELLI,ALTRLIPRIME DOUBLE PRECISION FACTOR,SUBTRACTION C C Integrates FCOLL(z) from X to 1, Monte Carlo style, C if called many times. C Z = X + (1.0D0 - X) * RANDOM(1.0) DZDX1 = 1.0D0 - X C------ C Calculate FCOLL(Z), the integrand for the 2 -> 2 coll pieces. C FACTOR = 2.0D0 * DLOG(RTS*X*(1.0D0 - Z)/SCALE/Z) SUBTRACTION = FACTOR * 2.0D0 * COLOR(A) > * PARTON(A,X,SCALE) /Z/(1.0D0 - Z) C FCOLL = - SUBTRACTION C IF (A.NE.0) THEN FCOLL = > FCOLL > + (FACTOR * ALTARELLI(A,A,Z) - ALTRLIPRIME(A,A,Z)) > * PARTON(A,X/Z,SCALE) /Z > + (FACTOR * ALTARELLI(A,0,Z)- ALTRLIPRIME(A,0,Z)) > * PARTON(0,X/Z,SCALE) /Z ELSE DO APRIME = - NFL,+NFL FCOLL = > FCOLL > + (FACTOR * ALTARELLI(A,APRIME,Z)- ALTRLIPRIME(A,APRIME,Z)) > * PARTON(APRIME,X/Z,SCALE) /Z C ENDDO ENDIF C------ EFFPARTON = FCOLL * DZDX1 RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION ALTARELLI(AOUT,AIN,Z) C INTEGER AOUT,AIN DOUBLE PRECISION Z DOUBLE PRECISION V,N,TWON,CF PARAMETER( V = 8.0D0 ) PARAMETER( N = 3.0D0 ) PARAMETER( TWON = 6.0D0 ) PARAMETER( CF = 4.0D0/3.0D0 ) C C The Altarelli-Parisi kernel P_{AOUT/AIN}(z) C IF ( (Z.GT.1.0D0).OR.(Z.LT.0.0D0) ) THEN ALTARELLI = 0.0D0 RETURN ENDIF C IF ( (AOUT.EQ.0).AND.(AIN.EQ.0) ) THEN ALTARELLI = TWON > *( (1.0D0 - Z) /Z + Z /(1.0D0 - Z) + Z*(1.0D0 - Z) ) ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.AOUT) ) THEN ALTARELLI = CF *(1.0D0 + Z**2) /(1.0D0 - Z) ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.0) ) THEN ALTARELLI = 0.5D0 * ( Z**2 + (1.0D0 - Z)**2 ) ELSE IF ( (AOUT.EQ.0).AND.(AIN.NE.0) ) THEN ALTARELLI = CF *( 1.0D0 + (1.0D0 - Z)**2 ) /Z ELSE ALTARELLI = 0.0D0 ENDIF RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION ALTRLIPRIME(AOUT,AIN,Z) C INTEGER AOUT,AIN DOUBLE PRECISION Z DOUBLE PRECISION V,N,TWON,CF PARAMETER( V = 8.0D0 ) PARAMETER( N = 3.0D0 ) PARAMETER( TWON = 6.0D0 ) PARAMETER( CF = 4.0D0/3.0D0) C C The derivative of the Altarelli-Parisi kernel C P_{AOUT/AIN}(z,\epsilon) w.r.t. \epsilon (\epsilon = 0) C IF ( (Z.GT.1.0D0).OR.(Z.LT.0.0D0) ) THEN ALTRLIPRIME = 0.0D0 RETURN ENDIF C IF ( (AOUT.EQ.0).AND.(AIN.EQ.0) ) THEN ALTRLIPRIME = 0.0D0 ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.AOUT) ) THEN ALTRLIPRIME = - CF * (1.0D0 - Z) ELSE IF ( (AOUT.NE.0).AND.(AIN.EQ.0) ) THEN ALTRLIPRIME = 0.5D0 * ( Z**2 + (1.0D0 - Z)**2 - 1.0D0 ) ELSE IF ( (AOUT.EQ.0).AND.(AIN.NE.0) ) THEN ALTRLIPRIME = - CF * Z ELSE ALTRLIPRIME = 0.0D0 ENDIF RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION GAMMA(NPARTON) C INTEGER NPARTON C C This function gives the A-P delta function coefficient for C parton NPARTON: 2 for quarks and antiquarks, C 11/2 - Nfl/3 for gluons C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL C IF (NPARTON.EQ.0) THEN GAMMA = 11.0D0/2.0D0 - NFLAVOR/3.0D0 ELSE GAMMA = 2.0D0 ENDIF RETURN END C C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION GAMMAPRIME(NPARTON) C INTEGER NPARTON C C This function gives the remainder from the collinear integration C for the coll-1 and coll-2 contributions for parton NPARTON C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL C DOUBLE PRECISION COLOR DOUBLE PRECISION PI,LOG2,V,N PARAMETER ( PI = 3.141592654 D0) PARAMETER ( LOG2 = 0.693147181 D0) PARAMETER ( V = 8.0 D0) PARAMETER ( N = 3.0 D0) C IF (NPARTON.EQ.0) THEN GAMMAPRIME = 67.0D0*N/9.0D0 - 23.0D0/18.0D0 * NFLAVOR > -2.0D0 * PI**2 * COLOR(NPARTON) /3.0D0 ELSE GAMMAPRIME = 13.0D0*V/(4.0D0*N) > -2.0D0 * PI**2 * COLOR(NPARTON) /3.0D0 ENDIF RETURN END C C*********************************************************************** C -------------------------------- C Functions for Soft Contributions C -------------------------------- C*********************************************************************** C DOUBLE PRECISION FUNCTION PSICTILDE(PROCESS,I1,I2,K) C INTEGER PROCESS INTEGER I1,I2 DOUBLE PRECISION K(4,4) C C This function gives the lambdatilde(a1,a2,s,t,u) functions C INTEGER II1,II2 DOUBLE PRECISION S,T,U DOUBLE PRECISION HAT,HAU,HB,HCT,HCU,HCS,HDS,HDT,HDU DOUBLE PRECISION N,V,CF PARAMETER( N = 3.0D0 ) PARAMETER( V = 8.0D0 ) PARAMETER( CF = 4.0D0 / 3.0D0 ) C S = 2.0D0 * K(1,2) T = 2.0D0 * K(1,3) U = 2.0D0 * K(1,4) C C Make it a symmetric function of I1,I2 C IF (I1.LT.I2) THEN II1 = I1 II2 = I2 ELSE II1 = I2 II2 = I1 ENDIF C IF (PROCESS.EQ.1) THEN HAT = 2.0D0 * (S**2 + U**2) /T**2 IF ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN PSICTILDE = 4.0D0 * CF * HAT ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN PSICTILDE = - 2.0D0 * CF * HAT ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN PSICTILDE = 2.0D0 * CF * (N**2 - 2.0D0) * HAT ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN PSICTILDE = 4.0D0 * CF * HAT ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN PSICTILDE = - 2.0D0 * CF* HAT ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN PSICTILDE = 2.0D0 * CF * (N**2 - 2.0D0) * HAT ENDIF C ELSEIF (PROCESS.EQ.2) THEN HAT = 2.0D0 * (S**2 + U**2) /T**2 HAU = 2.0D0 * (S**2 + T**2) /U**2 HB = 2.0D0 * S**2 /T /U IF ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN PSICTILDE = 4.0D0 * CF /N > *( - HB + N*(HAT + HAU) - N**2 * HB ) ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN PSICTILDE = 4.0D0 * CF /N > *( HB - N*(HAU + 0.5D0*HAT) + 0.5D0 * N**3 * HAU ) ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN PSICTILDE = 4.0D0 * CF /N > *( HB - N*(HAT + 0.5D0*HAU) + 0.5D0 * N**3 * HAT ) ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN PSICTILDE = 4.0D0 * CF /N > *( - HB + N*(HAT + HAU) - N**2 * HB ) ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN PSICTILDE = 4.0D0 * CF /N > *( HB - N*(HAU + 0.5D0*HAT) + 0.5D0 * N**3 * HAU ) ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN PSICTILDE = 4.0D0 * CF /N > *( HB - N*(HAT + 0.5D0*HAU) + 0.5D0 * N**3 * HAT ) ENDIF C ELSEIF (PROCESS.EQ.3) THEN HCT = 2.0D0 * (T**2 + U**2) /S**2 * T/U HCU = 2.0D0 * (T**2 + U**2) /S**2 * U/T HCS = 2.0D0 * (T**2 + U**2) /T /U IF ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN PSICTILDE = V * ( - HCT - HCU + HCS + HCS /N**2 ) ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN PSICTILDE = V * ( N**2 * HCU - HCS ) ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN PSICTILDE = V * ( N**2 * HCT - HCS ) ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN PSICTILDE = V * N**2 * ( HCT + HCU ) ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN PSICTILDE = V * ( N**2 * HCU - HCS ) ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN PSICTILDE = V * ( N**2 * HCT - HCS ) ENDIF C ELSEIF (PROCESS.EQ.4) THEN HDS = 2.0D0 * (T**2 + U**2) /S**2 /T**2 /U**2 > * (S**4 + T**4 + U**4) HDT = 2.0D0 * (U**2 + S**2) /S**2 /T**2 /U**2 > * (S**4 + T**4 + U**4) HDU = 2.0D0 * (S**2 + T**2) /S**2 /T**2 /U**2 > * (S**4 + T**4 + U**4) IF ( (II1.EQ.1).AND.(II2.EQ.2) ) THEN PSICTILDE = 2.0D0 * V * N**3 * HDS ELSEIF ( (II1.EQ.1).AND.(II2.EQ.3) ) THEN PSICTILDE = 2.0D0 * V * N**3 * HDT ELSEIF ( (II1.EQ.1).AND.(II2.EQ.4) ) THEN PSICTILDE = 2.0D0 * V * N**3 * HDU ELSEIF ( (II1.EQ.3).AND.(II2.EQ.4) ) THEN PSICTILDE = 2.0D0 * V * N**3 * HDS ELSEIF ( (II1.EQ.2).AND.(II2.EQ.4) ) THEN PSICTILDE = 2.0D0 * V * N**3 * HDT ELSEIF ( (II1.EQ.2).AND.(II2.EQ.3) ) THEN PSICTILDE = 2.0D0 * V * N**3 * HDU ENDIF C ELSE PSICTILDE = 0.0D0 ENDIF RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION LI2(X) C DOUBLE PRECISION X C C Let ERR be the target fractional error C DOUBLE PRECISION XN INTEGER N DOUBLE PRECISION ERR PARAMETER(ERR = 1.0D-6) C DOUBLE PRECISION PI PARAMETER ( PI = 3.141592654 D0) C XN = 1.0D0 LI2 = 0.0D0 C DO N=1,100 XN = XN * X IF (XN.LT.ERR) THEN RETURN ENDIF LI2 = LI2 + XN/N**2 ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C C End of functions for 1- and 2-dimensional integrals C C*********************************************************************** C*********************************************************************** C C SUBROUTINE INTEGRATEA(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF) C DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,COEF C C Integrates Term A. First compute Y3 and P3. (For Y3, set C Y3 = ln(\infty) = 30 if W is zero or very near to zero.) C We sum over terms in the jet-defining function, labelled by CASE. C For each term, JETDEF returns OK = .TRUE. if the theta functions C for that term are satisfied and then returns the values of PJ and C YJ determined by the delta functions for that term. We then set C MUUV and MUCO accordingly and call the BINIT subroutine, which C multiplies by gaussian smearing for each PJ,YJ grid point and C adds the resulting integrand to the integral for that grid point. C DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C DOUBLE PRECISION Y3,P3 DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE DOUBLE PRECISION MUUV,MUCO DOUBLE PRECISION FA,INTEGRAND LOGICAL OK INTEGER CASE C INTEGER NJETDEF COMMON /DEFJET/ NJETDEF C IF(W.GT.1.0D-12)THEN Y3 = DLOG(RTS/W) ELSE Y3 = 30 ENDIF C P3 = XI*W C DO CASE = 1,NJETDEF CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3, > OK,JETVAR1,JETVAR2,JETVAR3,SVALUE) IF(OK) THEN CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUCO) INTEGRAND = COEF * SVALUE > * FA(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO) CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND) ENDIF ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C SUBROUTINE INTEGRATEB(Y1,P2,Y2,PHI2,XI,W,PHI3,COEF) C DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,COEF C C Integrates Term A. First compute Y3 and P3. (For Y3, set C Y3 = -ln(\infty) = -30 if W is zero or very near to zero.) C We sum over terms in the jet-defining function, labelled by CASE. C For each term, JETDEF returns OK = .TRUE. if the theta functions C for that term are satisfied and then returns the values of PJ and C YJ determined by the delta functions for that term. We then set C MUUV and MUCO accordingly and call the BINIT subroutine, which C multiplies by gaussian smearing for each PJ,YJ grid point and C adds the resulting integrand to the integral for that grid point. C DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C DOUBLE PRECISION Y3,P3 DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE DOUBLE PRECISION MUUV,MUCO DOUBLE PRECISION FB,INTEGRAND LOGICAL OK INTEGER CASE C INTEGER NJETDEF COMMON /DEFJET/ NJETDEF C IF(W.GT.1.0D-12)THEN Y3 = - DLOG(RTS/W) ELSE Y3 = - 30 ENDIF C P3 = XI*W C DO CASE = 1,NJETDEF CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3, > OK,JETVAR1,JETVAR2,JETVAR3,SVALUE) IF(OK) THEN CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUCO) INTEGRAND = COEF * SVALUE > * FB(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO) CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND) ENDIF ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C SUBROUTINE INTEGRATE1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF) C DOUBLE PRECISION Y2,P1,Y1,PHI1,P3,Y3,PHI3,COEF C C Integrates Term 1. First we compute new variables P2 and PHI2. C We sum over terms in the jet-defining function, labelled by CASE. C For each term, JETDEF returns OK = .TRUE. if the theta functions C for that term are satisfied and then returns the values of PJ and C YJ determined by the delta functions for that term. We then set C MUUV and MUCO accordingly and call the BINIT subroutine, which C multiplies by gaussian smearing for each PJ,YJ grid point and C adds the resulting integrand to the integral for that grid point. C DOUBLE PRECISION P2X,P2Y,P2,PHI2 DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE DOUBLE PRECISION MUUV,MUCO DOUBLE PRECISION F1,INTEGRAND LOGICAL OK INTEGER CASE C INTEGER NJETDEF COMMON /DEFJET/ NJETDEF C P2X = - P1 * DCOS(PHI1) - P3 * DCOS(PHI3) P2Y = - P1 * DSIN(PHI1) - P3 * DSIN(PHI3) P2 = DSQRT(P2X**2 + P2Y**2) PHI2 = DATAN2(P2Y,P2X) C DO CASE = 1,NJETDEF CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3, > OK,JETVAR1,JETVAR2,JETVAR3,SVALUE) IF(OK) THEN CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUCO) INTEGRAND = COEF * SVALUE > * F1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUCO) CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND) ENDIF ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C SUBROUTINE INTEGRATE2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF) C DOUBLE PRECISION Y1,P2,Y2,PHI2,P3,Y3,PHI3,COEF C C Integrates Term 2. C We sum over terms in the jet-defining function, labelled by CASE. C For each term, JETDEF returns OK = .TRUE. if the theta functions C for that term are satisfied and then returns the values of PJ and C YJ determined by the delta functions for that term. We then set C MUUV and MUCO accordingly and call the BINIT subroutine, which C multiplies by gaussian smearing for each PJ,YJ grid point and C adds the resulting integrand to the integral for that grid point. C DOUBLE PRECISION JETVAR1,JETVAR2,JETVAR3,SVALUE DOUBLE PRECISION MUUV,MUCO DOUBLE PRECISION F2,INTEGRAND LOGICAL OK INTEGER CASE C INTEGER NJETDEF COMMON /DEFJET/ NJETDEF C DO CASE = 1,NJETDEF CALL JETDEF(CASE,Y1,P2,Y2,PHI2,P3,Y3,PHI3, > OK,JETVAR1,JETVAR2,JETVAR3,SVALUE) IF(OK) THEN CALL SETMU(JETVAR1,JETVAR2,JETVAR3,MUUV,MUCO) INTEGRAND = COEF * SVALUE > * F2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUCO) CALL BINIT(JETVAR1,JETVAR2,JETVAR3,INTEGRAND) ENDIF ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION > FA(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO) C DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C LOGICAL STYPE(4) COMMON /SWTYPE/ STYPE C INTEGER PERM(5),PERMINV(5) INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4) COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR C DOUBLE PRECISION PIN(5,5),K(5,5) DOUBLE PRECISION PAB,PA1,PA2,PA3MOD, > PB1,PB2,PB3MOD,P12,P13MOD,P23MOD DOUBLE PRECISION P1,XA,XB DOUBLE PRECISION COS23,EY1,EY2,P3,FEXP DOUBLE PRECISION ALPHAS,FACOEF,FASUM,LUMINOSITY DOUBLE PRECISION LMNSTY(7) INTEGER I,J,PROCESS,NPERM C DOUBLE PRECISION SIGN,CROSSINGSIGN DOUBLE PRECISION PSIAB,PSIA1,PSIA2,RESIDUE C DOUBLE PRECISION PISQUARE PARAMETER (PISQUARE = 9.869604404D0) C C COS23 = DCOS(PHI2-PHI3) EY1 = FEXP(Y1) EY2 = FEXP(Y2) P3 = XI*W P1 = DSQRT( P2**2 + P3**2 + 2.0D0*P2*P3*COS23 ) XA = (1/RTS)*( P1*EY1 + P2*EY2 + XI*RTS ) XB = (1/RTS)*( P1/EY1 + P2/EY2 + XI * W**2 /RTS ) C C Test for particle identification conditions and XA,XB: C IF (P3.GT.P1) THEN FA = 0.0D0 RETURN ELSE IF (P3.GT.P2) THEN FA = 0.0D0 RETURN ELSE IF (XA.GT.1.0D0) THEN FA = 0.0D0 RETURN ELSE IF (XB.GT.1.0D0) THEN FA = 0.0D0 RETURN ENDIF C C Construct the dot products of the P_J vectors, J = A,B,1,2,3. C The 'MOD' variables are the dot products divided by XI. C This division is done algebraicly because XI may be zero. C PAB = XA*XB*RTS**2 /2.0D0 PA1 = XA*RTS*P1 /EY1 /2.0D0 PA2 = XA*RTS*P2 /EY2 /2.0D0 PB1 = XB*RTS*P1 *EY1 /2.0D0 PB2 = XB*RTS*P2 *EY2 /2.0D0 P12 = P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P2**2 + P2*P3*COS23 PA3MOD = XA * W**2 /2.0D0 PB3MOD = XB * RTS**2 /2.0D0 P13MOD = 0.5D0 * P1 * RTS /EY1 + 0.5D0 * P1 * W**2 * EY1 /RTS > + XI * W**2 + P2 * W * COS23 P23MOD = 0.5D0 * P2 * RTS /EY2 + 0.5D0 * P2 * W**2 * EY2 /RTS > - P2 * W * COS23 C C Construct the PIN matrix (incoming momenta) C PIN(1,2) = PAB PIN(1,3) = -PA1 PIN(1,4) = -PA2 PIN(1,5) = -XI*PA3MOD PIN(2,3) = -PB1 PIN(2,4) = -PB2 PIN(2,5) = -XI*PB3MOD PIN(3,4) = P12 PIN(3,5) = XI*P13MOD PIN(4,5) = XI*P23MOD C DO I = 1,5 DO J = 1,I-1 PIN(I,J) = PIN(J,I) ENDDO ENDDO C C Calculate coefficient C FACOEF = ALPHAS(MUUV)**3 * P2 /(4.0D0 * PISQUARE * XA * RTS**4) C C Generate LMNSTY array for the given XA,XB C CALL PARTONSIN(LMNSTY,XA,XB,MUCO) C C Initialize FA C FA = 0.0D0 C C Main loops over processes and permutations C DO PROCESS = 1,4 C IF ( .NOT.STYPE(PROCESS) ) GOTO 300 C DO NPERM = 1,NPERMS(PROCESS) C C Initialize for each process and permutation C DO I = 1,5 PERM(I) = PERMARR(I,NPERM,PROCESS) PERMINV(I) = PERMINVARR(I,NPERM,PROCESS) ENDDO C DO I = 1,5 DO J = 1,5 K(I,J) = PIN(PERM(I),PERM(J)) ENDDO ENDDO C C Complete calculation C SIGN = CROSSINGSIGN(PROCESS,PERMINV) C PSIAB = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(2),K) PSIA1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(3),K) PSIA2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(4),K) C FASUM = PAB/(PA3MOD + PB3MOD) * PSIAB > + PA1/(PA3MOD + P13MOD) * PSIA1 > + PA2/(PA3MOD + P23MOD) * PSIA2 C FA = FA + FACOEF * FASUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV) C C Close loop over NPERM ENDDO C Close loop over PROCESS 300 ENDDO C RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION > FB(Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO) C DOUBLE PRECISION Y1,P2,Y2,PHI2,XI,W,PHI3,MUUV,MUCO C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C LOGICAL STYPE(4) COMMON /SWTYPE/ STYPE C INTEGER PERM(5),PERMINV(5) INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4) COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR C DOUBLE PRECISION PIN(5,5),K(5,5) DOUBLE PRECISION PAB,PA1,PA2,PA3MOD, > PB1,PB2,PB3MOD,P12,P13MOD,P23MOD DOUBLE PRECISION P1,XA,XB DOUBLE PRECISION COS23,EY1,EY2,P3,FEXP DOUBLE PRECISION ALPHAS,FBCOEF,FBSUM,LUMINOSITY DOUBLE PRECISION LMNSTY(7) INTEGER I,J,PROCESS,NPERM C DOUBLE PRECISION SIGN,CROSSINGSIGN DOUBLE PRECISION PSIAB,PSIB1,PSIB2,RESIDUE C DOUBLE PRECISION PISQUARE PARAMETER (PISQUARE = 9.869604404D0) C C COS23 = DCOS(PHI2-PHI3) EY1 = FEXP(Y1) EY2 = FEXP(Y2) P3 = XI*W P1 = DSQRT( P2**2 + P3**2 + 2.0D0*P2*P3*COS23 ) XA = (1/RTS)*( P1*EY1 + P2*EY2 + XI * W**2 / RTS) XB = (1/RTS)*( P1/EY1 + P2/EY2 + XI*RTS ) C C Test for particle identification conditions and XA,XB: C IF (P3.GT.P1) THEN FB = 0.0D0 RETURN ELSE IF (P3.GT.P2) THEN FB = 0.0D0 RETURN ELSE IF (XA.GT.1.0D0) THEN FB = 0.0D0 RETURN ELSE IF (XB.GT.1.0D0) THEN FB = 0.0D0 RETURN ENDIF C C Construct the dot products of the P_J vectors, J = A,B,1,2,3. C The 'MOD' variables are the dot products divided by XI. C This division is done algebraicly because XI may be zero. C PAB = XA*XB*RTS**2 /2.0D0 PA1 = XA*RTS*P1 /EY1 /2.0D0 PA2 = XA*RTS*P2 /EY2 /2.0D0 PB1 = XB*RTS*P1* EY1 /2.0D0 PB2 = XB*RTS*P2* EY2 /2.0D0 P12 = P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P2**2 + P2*P3*COS23 PA3MOD = XA * RTS**2 /2.0D0 PB3MOD = XB * W**2 /2.0D0 P13MOD = 0.5D0 * P1 * RTS * EY1 + 0.5D0 * P1* W**2 / EY1 / RTS > + XI * W**2 + P2 * W * COS23 P23MOD = 0.5D0 * P2 * RTS * EY2 + 0.5D0 * P2 * W**2 / EY2 / RTS > - P2 * W * COS23 C C Construct the PIN matrix (incoming momenta): C PIN(1,2) = PAB PIN(1,3) = -PA1 PIN(1,4) = -PA2 PIN(1,5) = -XI*PA3MOD PIN(2,3) = -PB1 PIN(2,4) = -PB2 PIN(2,5) = -XI*PB3MOD PIN(3,4) = P12 PIN(3,5) = XI*P13MOD PIN(4,5) = XI*P23MOD C DO I = 1,5 DO J = 1,I-1 PIN(I,J) = PIN(J,I) ENDDO ENDDO C C Calculate coefficient for factor HB C FBCOEF = ALPHAS(MUUV)**3 * P2 /(4.0D0 * PISQUARE * XB * RTS**4) C C Generate LMNSTY array for the given XA,XB C CALL PARTONSIN(LMNSTY,XA,XB,MUCO) C C Initialize FB C FB = 0.0D0 C C Main loops over processes and permutations C DO PROCESS = 1,4 C IF ( .NOT.STYPE(PROCESS) ) GOTO 300 C DO NPERM = 1,NPERMS(PROCESS) C C Initialize for each process and permutation C DO I = 1,5 PERM(I) = PERMARR(I,NPERM,PROCESS) PERMINV(I) = PERMINVARR(I,NPERM,PROCESS) ENDDO C DO I = 1,5 DO J = 1,5 K(I,J) = PIN(PERM(I),PERM(J)) ENDDO ENDDO C C Complete calculation: C SIGN = CROSSINGSIGN(PROCESS,PERMINV) C PSIAB = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(2),K) PSIB1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(3),K) PSIB2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(4),K) C FBSUM = PAB/(PA3MOD + PB3MOD) * PSIAB > + PB1/(PB3MOD + P13MOD) * PSIB1 > + PB2/(PB3MOD + P23MOD) * PSIB2 C FB = FB + FBCOEF * FBSUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV) C C Close loop over NPERM ENDDO C Close loop over PROCESS 300 ENDDO C RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION > F1(Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUCO) C DOUBLE PRECISION Y2,P1,Y1,PHI1,P3,Y3,PHI3,MUUV,MUCO C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C LOGICAL STYPE(4) COMMON /SWTYPE/ STYPE C INTEGER PERM(5),PERMINV(5) INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4) COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR C DOUBLE PRECISION PIN(5,5),K(5,5) DOUBLE PRECISION PAB,PA1,PA2,PA3MOD, > PB1,PB2,PB3MOD,P12,P13MOD,P23MOD DOUBLE PRECISION P2,XA,XB DOUBLE PRECISION COS13,EY1,EY2,EY3,FEXP DOUBLE PRECISION ALPHAS,F1COEF,F1SUM,LUMINOSITY DOUBLE PRECISION LMNSTY(7) INTEGER I,J,PROCESS,NPERM C DOUBLE PRECISION SIGN,CROSSINGSIGN DOUBLE PRECISION PSIA1,PSIB1,PSI12,RESIDUE C DOUBLE PRECISION PISQUARE PARAMETER (PISQUARE = 9.869604404D0) C COS13 = DCOS(PHI1-PHI3) EY1 = FEXP(Y1) EY2 = FEXP(Y2) EY3 = FEXP(Y3) P2 = DSQRT( P1**2 + P3**2 + 2.0D0*P1*P3*COS13 ) XA = (1/RTS)*( P1*EY1 + P2*EY2 + P3*EY3 ) XB = (1/RTS)*( P1/EY1 + P2/EY2 + P3/EY3 ) C C Test for particle identification conditions and XA,XB: C IF (P3.GT.P1) THEN F1 = 0.0D0 RETURN ELSE IF (P3.GT.P2) THEN F1 = 0.0D0 RETURN ELSE IF (XA.GT.1.0D0) THEN F1 = 0.0D0 RETURN ELSE IF (XB.GT.1.0D0) THEN F1 = 0.0D0 RETURN ENDIF C C Construct the dot products of the P_J vectors, J = A,B,1,2,3. C The 'MOD' variables are the dot products divided by P3. C This division is done algebraicly because P3 may be zero. C PAB = XA*XB*(RTS)**2 /2.0D0 PA1 = XA*RTS*P1 /EY1 /2.0D0 PA2 = XA*RTS*P2 /EY2 /2.0D0 PA3MOD = XA*RTS /EY3 /2.0D0 PB1 = XB*RTS*P1* EY1 /2.0D0 PB2 = XB*RTS*P2* EY2 /2.0D0 PB3MOD = XB*RTS * EY3 /2.0D0 P12 = P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P1**2 + P1*P3*COS13 P13MOD = P1 *0.5D0*(EY1/EY3 + EY3/EY1) - P1 *COS13 P23MOD = P2 *0.5D0*(EY2/EY3 + EY3/EY2) + P3 + P1 *COS13 C C Construct the PIN matrix (incoming momenta): C PIN(1,2) = PAB PIN(1,3) = -PA1 PIN(1,4) = -PA2 PIN(1,5) = -P3*PA3MOD PIN(2,3) = -PB1 PIN(2,4) = -PB2 PIN(2,5) = -P3*PB3MOD PIN(3,4) = P12 PIN(3,5) = P3*P13MOD PIN(4,5) = P3*P23MOD C DO I = 1,5 DO J = 1,I-1 PIN(I,J) = PIN(J,I) ENDDO ENDDO C C Calculate coefficient for factor H1 C F1COEF = ALPHAS(MUUV)**3 /(8.0D0 * PISQUARE * RTS**4) C C Generate LMNSTY array for the given XA,XB C CALL PARTONSIN(LMNSTY,XA,XB,MUCO) C C Initialize psi variables C F1 = 0.0D0 C C Main loops over processes and permutations C DO PROCESS = 1,4 C IF ( .NOT.STYPE(PROCESS) ) GOTO 300 C DO NPERM = 1,NPERMS(PROCESS) C C Initialize for each process and permutation C DO I = 1,5 PERM(I) = PERMARR(I,NPERM,PROCESS) PERMINV(I) = PERMINVARR(I,NPERM,PROCESS) ENDDO C DO I = 1,5 DO J = 1,5 K(I,J) = PIN(PERM(I),PERM(J)) ENDDO ENDDO C C Complete calculation: C SIGN = CROSSINGSIGN(PROCESS,PERMINV) C PSIA1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(3),K) PSIB1 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(3),K) PSI12 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(3),PERMINV(4),K) C F1SUM = PA1/(PA3MOD + P13MOD) * PSIA1 > + PB1/(PB3MOD + P13MOD) * PSIB1 > + P12/(P13MOD + P23MOD) * PSI12 C F1 = F1 + F1COEF * F1SUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV) C C Close loop over NPERM ENDDO C Close loop over PROCESS 300 ENDDO C RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION > F2(Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUCO) C DOUBLE PRECISION Y1,P2,Y2,PHI2,P3,Y3,PHI3,MUUV,MUCO C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL DOUBLE PRECISION RTS,R COMMON /PHYSICS/ RTS,R C LOGICAL STYPE(4) COMMON /SWTYPE/ STYPE C INTEGER PERM(5),PERMINV(5) INTEGER NPERMS(4),PERMARR(5,30,4),PERMINVARR(5,30,4) COMMON/PERMS5/NPERMS,PERMARR,PERMINVARR C DOUBLE PRECISION PIN(5,5),K(5,5) DOUBLE PRECISION PAB,PA1,PA2,PA3MOD, > PB1,PB2,PB3MOD,P12,P13MOD,P23MOD DOUBLE PRECISION P1,XA,XB DOUBLE PRECISION COS23,EY1,EY2,EY3,FEXP DOUBLE PRECISION ALPHAS,F2COEF,F2SUM,LUMINOSITY DOUBLE PRECISION LMNSTY(7) INTEGER I,J,PROCESS,NPERM C DOUBLE PRECISION SIGN,CROSSINGSIGN DOUBLE PRECISION PSIA2,PSIB2,PSI12,RESIDUE C DOUBLE PRECISION PISQUARE PARAMETER (PISQUARE = 9.869604404D0) C COS23 = DCOS(PHI2-PHI3) EY1 = FEXP(Y1) EY2 = FEXP(Y2) EY3 = FEXP(Y3) P1 = DSQRT( P2**2 + P3**2 + 2.0D0*P2*P3*COS23 ) XA = (1/RTS)*( P1*EY1 + P2*EY2 + P3*EY3 ) XB = (1/RTS)*( P1/EY1 + P2/EY2 + P3/EY3 ) C C Test for particle identification conditions and XA,XB: C IF (P3.GT.P1) THEN F2 = 0.0D0 RETURN ELSE IF (P3.GT.P2) THEN F2 = 0.0D0 RETURN ELSE IF (XA.GT.1.0D0) THEN F2 = 0.0D0 RETURN ELSE IF (XB.GT.1.0D0) THEN F2 = 0.0D0 RETURN ENDIF C C Construct the dot products of the P_J vectors, J = A,B,1,2,3. C The 'MOD' variables are the dot products divided by P3. C This division is done algebraicly because P3 may be zero. C PAB = XA*XB*(RTS)**2 /2.0D0 PA1 = XA*RTS*P1 /EY1 /2.0D0 PA2 = XA*RTS*P2 /EY2 /2.0D0 PA3MOD = XA*RTS /EY3 /2.0D0 PB1 = XB*RTS*P1* EY1 /2.0D0 PB2 = XB*RTS*P2* EY2 /2.0D0 PB3MOD = XB*RTS * EY3 /2.0D0 P12 = P1*P2*0.5D0*(EY1/EY2 + EY2/EY1) + P2**2 + P2*P3*COS23 P13MOD = P1 *0.5D0*(EY1/EY3 + EY3/EY1) + P3 + P2 *COS23 P23MOD = P2 *0.5D0*(EY2/EY3 + EY3/EY2) - P2 *COS23 C C Construct the PIN matrix (incoming momenta): C PIN(1,2) = PAB PIN(1,3) = -PA1 PIN(1,4) = -PA2 PIN(1,5) = -P3*PA3MOD PIN(2,3) = -PB1 PIN(2,4) = -PB2 PIN(2,5) = -P3*PB3MOD PIN(3,4) = P12 PIN(3,5) = P3*P13MOD PIN(4,5) = P3*P23MOD C DO I = 1,5 DO J = 1,I-1 PIN(I,J) = PIN(J,I) ENDDO ENDDO C C Calculate coefficient for factor H2 C F2COEF = ALPHAS(MUUV)**3 /(8.0D0 * PISQUARE * RTS**4) C C Generate LMNSTY array for the given XA,XB C CALL PARTONSIN(LMNSTY,XA,XB,MUCO) C C Initialize F2 C F2 = 0.0D0 C C Main loops over processes and permutations C DO PROCESS = 1,4 C IF ( .NOT.STYPE(PROCESS) ) GOTO 300 C DO NPERM = 1,NPERMS(PROCESS) C C Initialize for each process and permutation C DO I = 1,5 PERM(I) = PERMARR(I,NPERM,PROCESS) PERMINV(I) = PERMINVARR(I,NPERM,PROCESS) ENDDO C DO I = 1,5 DO J = 1,5 K(I,J) = PIN(PERM(I),PERM(J)) ENDDO ENDDO C C Complete calculation: C SIGN = CROSSINGSIGN(PROCESS,PERMINV) C PSIA2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(1),PERMINV(4),K) PSIB2 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(2),PERMINV(4),K) PSI12 = SIGN * RESIDUE(PROCESS,PERMINV(5),PERMINV(3),PERMINV(4),K) C F2SUM = PA2/(PA3MOD + P23MOD) * PSIA2 > + PB2/(PB3MOD + P23MOD) * PSIB2 > + P12/(P13MOD + P23MOD) * PSI12 C F2 = F2 + F2COEF * F2SUM * LUMINOSITY(LMNSTY,PROCESS,PERMINV) C C Close loop over NPERM ENDDO C Close loop over PROCESS 300 ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION CROSSINGSIGN(PROCESS,PERMINV) C INTEGER PROCESS,PERMINV(5) C C This function calculated the crossing sign needed in the C the four dimensional integrals. SIGN is -1 if one of C particles A and B is a gluon while the other is not. Otherwise C SIGN = +1. C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL C INTEGER A,B C C Start with a plus sign, then change it if necessary: C CROSSINGSIGN = 1.0D0 A = PERMINV(1) B = PERMINV(2) C C----------------------------- C IF ((PROCESS.EQ.1).OR.(PROCESS.EQ.2)) THEN C IF ((A.EQ.5).OR.(B.EQ.5)) THEN CROSSINGSIGN = - CROSSINGSIGN ENDIF C ELSE IF (PROCESS.EQ.3) THEN C IF ((A.EQ.1).OR.(A.EQ.2)) THEN CROSSINGSIGN = - CROSSINGSIGN ENDIF IF ((B.EQ.1).OR.(B.EQ.2)) THEN CROSSINGSIGN = - CROSSINGSIGN ENDIF C ENDIF C C------------------------ C RETURN END C C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION CONVERT(PHI) DOUBLE PRECISION PHI C C Gets PHI + 2*pi*N such that the result lies between - pi and pi. C Here IDNINT( x ) is the nearest integer to x. C DOUBLE PRECISION TWOPI PARAMETER ( TWOPI = 6.283185308 D0) C CONVERT = PHI - TWOPI * IDNINT( PHI/TWOPI ) C RETURN END C C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION RESIDUE(PROCESS,J,N,M,KIN) C INTEGER PROCESS,J,N,M DOUBLE PRECISION KIN(5,5) C C This function calculates the function tilde-Psi^j_{nm}, C which is the coefficient of (k_n*k_m) / (k_n*k_j)(k_m*k_j)) C in the expansion of tilde-\Psi(K) for the PROCESS. C C In the computation,only certain standard values of j,n,m are used. C Then tilde-\Psi for other values of j,n,m are obtained by C using symmetries. The relation is C C Psi^j_{n m}(KIN) = Psi^jj_{nn mm}(K) C C with j = PERM(jj), n = PERM(nn), m = PERM(nm) C KIN(I,J) = K(PERM(I),PERM(J)). C The proper value of JJ for process A is given in JA(J,N,M). C The proper value of NN,MM for process A is in SUBSCRIPTA(J,N,M). C as a two digit integer SUBSCRIPTA = 10*min(NN,MM) + max(NN,MM). C (Recall that Psi^j_{n m} = Psi^j_{m n} ). C The proper value of PERM(I) for process A is in PERMA(I,J,N,M). C Simarly for processes B,C,D. C These values are put into COMMON /RESIDUEHELP/ by the C subroutine RESIDUEINIT at the start of the program. C INTEGER I INTEGER JJ,SUBSCRIPT,PERM(5) C INTEGER JA(5,5,5),JB(5,5,5),JC(5,5,5),JD(5,5,5) INTEGER SUBSCRIPTA(5,5,5),SUBSCRIPTB(5,5,5) INTEGER SUBSCRIPTC(5,5,5),SUBSCRIPTD(5,5,5) INTEGER PERMA(5,5,5,5),PERMB(5,5,5,5) INTEGER PERMC(5,5,5,5),PERMD(5,5,5,5) COMMON /RESIDUEHELP/ JA,JB,JC,JD, > SUBSCRIPTA,SUBSCRIPTB,SUBSCRIPTC,SUBSCRIPTD, > PERMA,PERMB,PERMC,PERMD C C Variables for the calculation of all processes C DOUBLE PRECISION K12,K13,K14,K15,K23,K24,K25,K34,K35,K45 DOUBLE PRECISION NCOLOR,V,CF,TWOCF,FOURCF PARAMETER( NCOLOR = 3.0D0) PARAMETER( V = 8.0D0) PARAMETER( CF = V/2.0D0/NCOLOR) PARAMETER(TWOCF = 2.0D0 * CF) PARAMETER(FOURCF = 4.0D0 * CF) DOUBLE PRECISION NCF,TWONCF,FOURNCF2,FOURCF2 PARAMETER(NCF = NCOLOR * CF) PARAMETER(TWONCF = 2.0D0 * NCOLOR * CF) PARAMETER(FOURNCF2 = 4.0D0 * NCOLOR * CF**2) PARAMETER(FOURCF2 = 4.0D0 * CF**2) C C Variables for process A: C DOUBLE PRECISION COEF512,COEF513,COEF514,FCTOR PARAMETER(COEF512 = 4.0D0 * CF) PARAMETER(COEF513 = - 2.0D0 * CF) PARAMETER(COEF514 = 2.0D0 * CF * ( 2.0D0 * NCOLOR * CF - 1.0D0 )) DOUBLE PRECISION Q1,Q2,Q3 C C Variables for process B: C DOUBLE PRECISION FCTOR1,FCTOR2,FCTOR3,COMBO,COMBO2 DOUBLE PRECISION H1,H2,H3,H4,H5,H6,H7 DOUBLE PRECISION SFACTOR,TFACTOR,UFACTOR,STUFACTOR DOUBLE PRECISION ONELESS2NCF,TWOPLUS2NCF PARAMETER(ONELESS2NCF = 1.0D0 - 2.0D0*NCOLOR*CF) PARAMETER(TWOPLUS2NCF = 2.0D0 + 2.0D0*NCOLOR*CF) C C Variables for process C: C DOUBLE PRECISION VOVERN2,ONEPLUSN2,N2,N2OVER3,N2COMBO PARAMETER(VOVERN2 = V/NCOLOR**2) PARAMETER(ONEPLUSN2 = 1.0D0 + NCOLOR**2) PARAMETER(N2 = NCOLOR**2) PARAMETER(N2OVER3 = NCOLOR**2/3.0D0) PARAMETER(N2COMBO = (1.0D0 + NCOLOR**2)/NCOLOR**2) DOUBLE PRECISION H123,H124 DOUBLE PRECISION U123,U124,V124 DOUBLE PRECISION W123,W124,W125 DOUBLE PRECISION S1314,S1525,S1535,S1545 DOUBLE PRECISION S2324,S2334,S2434,S2535,S2545,S3545 DOUBLE PRECISION N2OVERK12 DOUBLE PRECISION X123,X124,X125,Y123,Y124,Y125 DOUBLE PRECISION L2345,L3435,L3445,L3545 C C Variables for process D: C DOUBLE PRECISION SUM4,DENOMS,TWOVN3 PARAMETER(TWOVN3 = 2.0D0 * V * NCOLOR**3) C C-------------- C C Program must change this! C RESIDUE = 1.0D15 C IF (PROCESS.EQ.1) THEN C JJ = JA(J,N,M) SUBSCRIPT = SUBSCRIPTA(J,N,M) C DO I = 1,5 PERM(I) = PERMA(I,J,N,M) ENDDO C IF ( JJ.EQ.5 ) THEN C IF ( SUBSCRIPT.EQ.12 ) THEN C C The function tilde-Psi^5_{12} for process A: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- FCTOR = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24) RESIDUE = COEF512 * FCTOR C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.13 ) THEN C C The function tilde-Psi^5_{13} for process A: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- FCTOR = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24) RESIDUE = COEF513 * FCTOR C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.14 ) THEN C C The function tilde-Psi^5_{14} for process A: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- FCTOR = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24) RESIDUE = COEF514 * FCTOR C-- RETURN C C End SUBSCRIPT test, continue JJ test C ENDIF ELSEIF ( JJ.EQ.1 ) THEN C IF ( SUBSCRIPT.EQ.35 ) THEN C C The function tilde-Psi^1_{35} for process A: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- Q1 = K14*K23 - K12*K34 Q2 = K13*K24 Q3 = K12**2 + K14**2 + K23**2 + K34**2 RESIDUE = FOURCF * Q3 /(K24 * K25 * K35**2 * K45) > *( - K14*K23*(K13 + K24) > - K13*K24*(K14 + K23) > + (K12 + K34)*(Q1 + Q2)/2.0D0 > + NCF * ( K14*(K12*K13 + K24*K34) > + K23*(K12*K24 + K13*K34) > - (K14 + K23)*(Q1 - Q2) )) C-- RETURN C ELSEIF((SUBSCRIPT.EQ.23).OR.(SUBSCRIPT.EQ.24).OR. > (SUBSCRIPT.EQ.25).OR.(SUBSCRIPT.EQ.34).OR. > (SUBSCRIPT.EQ.45) ) THEN C C The function tilde-Psi^1_{23} for process A: C The function tilde-Psi^1_{24} for process A: C The function tilde-Psi^1_{25} for process A: C The function tilde-Psi^1_{34} for process A: C The function tilde-Psi^1_{45} for process A: C-- RESIDUE = 0.0D0 C-- RETURN C C End SUBSCRIPT test and JJ test C ENDIF ENDIF C C-------------- C ELSE IF (PROCESS.EQ.2) THEN C JJ = JB(J,N,M) SUBSCRIPT = SUBSCRIPTB(J,N,M) C DO I = 1,5 PERM(I) = PERMB(I,J,N,M) ENDDO C IF ( JJ.EQ.5 ) THEN C IF ( SUBSCRIPT.EQ.12 ) THEN C C The function tilde-Psi^5_{12} for process B: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- FCTOR1 = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24) FCTOR2 = (K12**2 + K13**2 + K24**2 + K34**2)/(K14*K23) FCTOR3 = (K12**2 + K34**2) > *(K12*K34 - K14*K23 - K13*K24)/(K13*K14*K23*K24) COMBO2 = FOURCF * (FCTOR1 + FCTOR2 - FCTOR3 /NCOLOR) RESIDUE = COMBO2 - FOURCF2 * FCTOR3 C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.13 ) THEN C C The function tilde-Psi^5_{13} for process B: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- FCTOR1 = (K12**2 + K14**2 + K23**2 + K34**2)/(K13*K24) FCTOR2 = (K12**2 + K13**2 + K24**2 + K34**2)/(K14*K23) FCTOR3 = (K12**2 + K34**2) > *(K12*K34 - K14*K23 - K13*K24)/(K13*K14*K23*K24) COMBO = TWOCF * (FCTOR1 + FCTOR2 - FCTOR3 /NCOLOR) RESIDUE = - COMBO + FOURNCF2 * FCTOR2 C-- RETURN C C End SUBSCRIPT test, continue JJ test C ENDIF ELSEIF ( JJ.EQ.1 ) THEN IF ( SUBSCRIPT.EQ.23 ) THEN C C The function tilde-Psi^1_{23} for process B: C-- RESIDUE = 0.0D0 C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.34 ) THEN C C The function tilde-Psi^1_{34} for process B: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- H2 = K24*K35 H3 = K23*K45 H4 = K35*K45 H5 = K25*K34 H1 = 2.0D0*H5 - H2 - H3 SFACTOR = K12**2 + K34**2 TFACTOR = K13**2 + K24**2 UFACTOR = K14**2 + K23**2 STUFACTOR = K12*K34 - K14*K23 - K13*K24 C RESIDUE = TWOCF/(H5*H4) > *( ( 2.0D0*H4 - H1 + TWONCF*(H4-H5) )/(K23*K24*NCOLOR) > *SFACTOR*STUFACTOR > + (K13/K23)*(SFACTOR + TFACTOR)*(H1 + TWONCF*H2) > + (K14/K24)*(SFACTOR + UFACTOR)*(H1 + TWONCF*H3) ) C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.25 ) THEN C C The function tilde-Psi^1_{25} for process B: C-- RESIDUE = 0.0D0 C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.35 ) THEN C C The function tilde-Psi^1_{35} for process B: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- H6 = K24*K45 H7 = K24*K25 SFACTOR = K12**2 + K34**2 UFACTOR = K14**2 + K23**2 STUFACTOR = K12*K34 - K14*K23 - K13*K24 C RESIDUE = ( TWOCF/K35) > *((SFACTOR + UFACTOR) > *( (2.0D0*K12 - K13*K25/K35)/H7 - K14*ONELESS2NCF/H6 ) > + SFACTOR*STUFACTOR/(K23*NCOLOR)*(1.0D0/H6 + TWOPLUS2NCF/H7)) C-- RETURN C C End SUBSCRIPT test and JJ test C ENDIF ENDIF C C-------------- C ELSE IF (PROCESS.EQ.3) THEN C JJ = JC(J,N,M) SUBSCRIPT = SUBSCRIPTC(J,N,M) C DO I = 1,5 PERM(I) = PERMC(I,J,N,M) ENDDO C IF ( JJ.EQ.5 ) THEN C IF ( SUBSCRIPT.EQ.12 ) THEN C C The function tilde-Psi^5_{12} for process C: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- H123 = K13**2 + K23**2 H124 = K14**2 + K24**2 C RESIDUE = VOVERN2 > *( ONEPLUSN2*(H123/(K14*K24) + H124/(K13*K23)) > - (N2/K34/K12) > *( H123*(K13/K14 + K23/K24)+ H124*(K24/K23 + K14/K13))) C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.13 ) THEN C C The function tilde-Psi^5_{13} for process C: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- H123 = K13**2 + K23**2 H124 = K14**2 + K24**2 U123 = K13**2 + K23**2 U124 = K14**2 + K24**2 V124 = K14*K24 C RESIDUE = V/K13 > *( - U124/K23 - K13*U123/V124 > + N2/(K12*K24*K34) > *( V124*H124 + K13*K23*H123) ) C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.34 ) THEN C C The function tilde-Psi^5_{34} for process C: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- W123 = K13**2 + K23**2 W124 = K14**2 + K24**2 W125 = K15**2 + K25**2 S1314 = K13*K14 S1525 = K15*K25 S1535 = K15*K35 S1545 = K15*K45 S2324 = K23*K24 S2334 = K23*K34 S2434 = K24*K34 S2535 = K25*K35 S2545 = K25*K45 S3545 = K35*K45 N2OVERK12 = N2/K12 C RESIDUE = V/K34 *( - W125 > *( (S1535/K23 + S1545/K24 )/S1314 > + (S2535/S2324 + S3545/S2434)/K13 > + (S2545/S2324 + S3545/S2334)/K14 > - S3545*ONEPLUSN2 /(S1314*S2324*N2OVERK12) > - N2OVERK12 > *( (S1525/K23 + S1545/K34)/K14 > + (S1525/K24 + S1535/K34)/K13 > + S2535/S2334 + S2545/S2434 ) ) > + N2OVERK12 > *( W123*(K13/K14 + K23/K24) > + W124*(K14/K13 + K24/K23) ) ) C-- RETURN C C End SUBSCRIPT test, continue JJ test C ENDIF ELSEIF ( JJ.EQ.1 ) THEN C IF ( SUBSCRIPT.EQ.23 ) THEN C C The function tilde-Psi^1_{23} for process C: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- X123 = K13**2 + K23**2 X124 = K14**2 + K24**2 X125 = K15**2 + K25**2 Y123 = K23 * X123 Y124 = K24 * X124 Y125 = K25 * X125 L2345 = K23 * K45 L3435 = K34 * K35 L3445 = K34 * K45 L3545 = K35 * K45 C RESIDUE = V/K23 > *( N2*(K25*K34 + K24*K35)/(L3435*K45*K24*K25) > *( K15*Y125 + K14*Y124) > - K12/K24/K25/L2345*( Y125*(L2345/K34 + K25) > + Y124*(L2345/K35 + K24) ) > + K13*N2OVER3 > *( (Y124/L3545 + Y125/L3445)/K23 > + (Y125/L3435 + Y123/L3545)/K24 > + (Y123/L3445 + Y124/L3435)/K25 )) C-- RETURN C ELSEIF ( SUBSCRIPT.EQ.34 ) THEN C C The function tilde-Psi^1_{34} for process C: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- RESIDUE = V*(K15**2 + K25**2)/K34/K23/K24 > *( K12*N2COMBO - K15*(K23/K35 + K24/K45)) C-- RETURN C C End SUBSCRIPT test and JJ test C ENDIF ENDIF C-------------- C ELSE IF (PROCESS.EQ.4) THEN C JJ = JD(J,N,M) SUBSCRIPT = SUBSCRIPTD(J,N,M) C DO I = 1,5 PERM(I) = PERMD(I,J,N,M) ENDDO C IF ( JJ.EQ.5 ) THEN C IF ( SUBSCRIPT.EQ.34 ) THEN C C The function tilde-Psi^5_{34} for process D: C CALL PERMUTE(KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C-- SUM4 = K12**4 + K13**4 + K14**4 + K15**4 > + K23**4 + K24**4 + K25**4 > + K34**4 + K35**4 > + K45**4 C DENOMS = (1.0D0 /K13/K24 + 1.0D0 /K23/K14)/K12/K34 C RESIDUE = TWOVN3 * SUM4 * DENOMS C-- RETURN C C End SUBSCRIPT test and JJ test C ENDIF ENDIF C C-------------- C C End PROCESS test C ENDIF C RETURN END C C*********************************************************************** C*********************************************************************** C SUBROUTINE PERMUTE > (KIN,K12,K13,K14,K15,K23,K24,K25,K34,K35,K45,PERM) C DOUBLE PRECISION KIN(5,5) DOUBLE PRECISION K12,K13,K14,K15,K23,K24,K25,K34,K35,K45 INTEGER PERM(5) C C Produces K(I,J) = KIN(PERM(I),PERM(J)) C K12 = KIN( PERM(1),PERM(2) ) K13 = KIN( PERM(1),PERM(3) ) K14 = KIN( PERM(1),PERM(4) ) K15 = KIN( PERM(1),PERM(5) ) K23 = KIN( PERM(2),PERM(3) ) K24 = KIN( PERM(2),PERM(4) ) K25 = KIN( PERM(2),PERM(5) ) K34 = KIN( PERM(3),PERM(4) ) K35 = KIN( PERM(3),PERM(5) ) K45 = KIN( PERM(4),PERM(5) ) C RETURN END C C*********************************************************************** C*********************************************************************** C SUBROUTINE RESIDUEINIT C C C This is to initialize arrays for the function RESIDUE, C which calculates the function tilde-Psi^j_{nm}, C which is the coefficient of (k_n*k_m) / (k_n*k_j)(k_m*k_j)) C in the expansion of tilde-\Psi(K) for the PROCESS. C C In the computation,only certain standard values of j,n,m are used. C Then tilde-\Psi for other values of j,n,m are obtained by C using symmetries. The relation is C C Psi^j_{n m}(KIN) = Psi^jj_{nn mm}(K) C C with jj = PERM(j), nn = PERM(n), mm = PERM(n) C K(I,J) = KIN(PERM(I),PERM(J)). C The proper value of JJ for process A is given in JA(J,N,M). C The proper value of NN,MM for process A is in SUBSCRIPTA(J,N,M) C as a two digit integer SUBSCRIPTA = 10*min(NN,MM) + max(NN,MM). C (Recall that Psi^j_{n m} = Psi^j_{m n} ). C The proper value of PERM(I) for process A is in PERMA(I,J,N,M). C Simarly for processes B,C,D. C These values are put into COMMON /RESIDUEHELP/ by the C subroutine RESIDUEINIT at the start of the program. C C INTEGER JA(5,5,5),JB(5,5,5),JC(5,5,5),JD(5,5,5) INTEGER SUBSCRIPTA(5,5,5),SUBSCRIPTB(5,5,5) INTEGER SUBSCRIPTC(5,5,5),SUBSCRIPTD(5,5,5) INTEGER PERMA(5,5,5,5),PERMB(5,5,5,5) INTEGER PERMC(5,5,5,5),PERMD(5,5,5,5) COMMON /RESIDUEHELP/ JA,JB,JC,JD, > SUBSCRIPTA,SUBSCRIPTB,SUBSCRIPTC,SUBSCRIPTD, > PERMA,PERMB,PERMC,PERMD C INTEGER PROCESS,J,N,M LOGICAL MATCH INTEGER I,L,LL INTEGER A(4,5),B(8,5),C(12,5) DATA A/1,2,3,4, > 2,1,4,3, > 3,4,1,2, > 4,3,2,1, > 5,5,5,5/ DATA B/1,1,2,2,3,3,4,4, > 2,2,1,1,4,4,3,3, > 3,4,3,4,1,2,1,2, > 4,3,4,3,2,1,2,1, > 5,5,5,5,5,5,5,5/ DATA C/1,1,1,1,1,1,2,2,2,2,2,2, > 2,2,2,2,2,2,1,1,1,1,1,1, > 3,3,4,4,5,5,3,3,4,4,5,5, > 4,5,3,5,3,4,4,5,3,5,3,4, > 5,4,5,3,4,3,5,4,5,3,4,3/ C C-------------- C DO J = 1,5 DO N = 1,5 DO M = 1,5 C IF((N.NE.J).AND.(M.NE.J).AND.(M.NE.N) ) THEN C PROCESS=1 C DO L = 1,4 C IF (MATCH(J,N,M,A(L,5),A(L,1),A(L,2))) THEN C JA(J,N,M) = 5 SUBSCRIPTA(J,N,M) = 12 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,5),A(L,1),A(L,3))) THEN C JA(J,N,M) = 5 SUBSCRIPTA(J,N,M) = 13 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,5),A(L,1),A(L,4))) THEN C JA(J,N,M) = 5 SUBSCRIPTA(J,N,M) = 14 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,1),A(L,2),A(L,3))) THEN C JA(J,N,M) = 1 SUBSCRIPTA(J,N,M) = 23 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,1),A(L,2),A(L,4))) THEN C JA(J,N,M) = 1 SUBSCRIPTA(J,N,M) = 24 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,1),A(L,3),A(L,4))) THEN C JA(J,N,M) = 1 SUBSCRIPTA(J,N,M) = 34 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,1),A(L,2),A(L,5))) THEN C JA(J,N,M) = 1 SUBSCRIPTA(J,N,M) = 25 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,1),A(L,3),A(L,5))) THEN C JA(J,N,M) = 1 SUBSCRIPTA(J,N,M) = 35 GOTO 101 C ELSE IF (MATCH(J,N,M,A(L,1),A(L,4),A(L,5))) THEN C JA(J,N,M) = 1 SUBSCRIPTA(J,N,M) = 45 GOTO 101 C C If none of the tests was successful, try a different L C ENDIF ENDDO C C This gets executed when test is successful for this L C 101 DO I = 1,5 PERMA(I,J,N,M) = A(L,I) ENDDO C C C-------------- C PROCESS=2 C DO L = 1,8 C IF (MATCH(J,N,M,B(L,5),B(L,1),B(L,2))) THEN C JB(J,N,M) = 5 SUBSCRIPTB(J,N,M) = 12 GOTO 102 C ELSE IF (MATCH(J,N,M,B(L,5),B(L,1),B(L,3))) THEN C JB(J,N,M) = 5 SUBSCRIPTB(J,N,M) = 13 GOTO 102 C ELSE IF (MATCH(J,N,M,B(L,1),B(L,2),B(L,3))) THEN C JB(J,N,M) = 1 SUBSCRIPTB(J,N,M) = 23 GOTO 102 C ELSE IF (MATCH(J,N,M,B(L,1),B(L,3),B(L,4))) THEN C JB(J,N,M) = 1 SUBSCRIPTB(J,N,M) = 34 GOTO 102 C ELSE IF (MATCH(J,N,M,B(L,1),B(L,2),B(L,5))) THEN C JB(J,N,M) = 1 SUBSCRIPTB(J,N,M) = 25 GOTO 102 C ELSE IF (MATCH(J,N,M,B(L,1),B(L,3),B(L,5))) THEN C JB(J,N,M) = 1 SUBSCRIPTB(J,N,M) = 35 GOTO 102 C C If none of the tests was successful, try a different L C ENDIF ENDDO C C This gets executed when test is successful for this L C 102 DO I = 1,5 PERMB(I,J,N,M) = B(L,I) ENDDO C C-------------- C PROCESS=3 C DO L = 1,12 C IF (MATCH(J,N,M,C(L,5),C(L,1),C(L,2))) THEN C JC(J,N,M) = 5 SUBSCRIPTC(J,N,M) = 12 GOTO 103 C ELSE IF (MATCH(J,N,M,C(L,5),C(L,1),C(L,3))) THEN C JC(J,N,M) = 5 SUBSCRIPTC(J,N,M) = 13 GOTO 103 C ELSE IF (MATCH(J,N,M,C(L,5),C(L,3),C(L,4))) THEN C JC(J,N,M) = 5 SUBSCRIPTC(J,N,M) = 34 GOTO 103 C ELSE IF (MATCH(J,N,M,C(L,1),C(L,2),C(L,3))) THEN C JC(J,N,M) = 1 SUBSCRIPTC(J,N,M) = 23 GOTO 103 C ELSE IF (MATCH(J,N,M,C(L,1),C(L,3),C(L,4))) THEN C JC(J,N,M) = 1 SUBSCRIPTC(J,N,M) = 34 GOTO 103 C C If none of the tests was successful, try a different L C ENDIF ENDDO C C This gets executed when test is successful for this L C 103 DO I = 1,5 PERMC(I,J,N,M) = C(L,I) ENDDO C C-------------- C PROCESS=4 C JD(J,N,M) = 5 SUBSCRIPTD(J,N,M) = 34 C PERMD(5,J,N,M) = J PERMD(4,J,N,M) = M PERMD(3,J,N,M) = N C C Now assign PERMD(2)=LL, PERMD(1)=L where L and LL are chosen C from the remaining numbers (i.e. not J,N,M) with L < LL. C When done, get out of the Do loops with a GOTO. C DO L = 1,4 IF ((L.NE.J).AND.(L.NE.N).AND.(L.NE.M)) THEN PERMD(1,J,N,M) = L DO LL = L+1,5 IF ((LL.NE.J).AND.(LL.NE.N).AND.(LL.NE.M)) THEN PERMD(2,J,N,M) = LL GOTO 702 ENDIF ENDDO ENDIF ENDDO C C C End of loops over J,N,M and ENDIF for (N.NE.J),(M.NE.J),(M.NE.N) C C 702 ENDIF ENDDO ENDDO ENDDO C RETURN END C C*********************************************************************** C LOGICAL FUNCTION MATCH(J,N,M,JJ,NN,MM) C INTEGER J,N,M,JJ,NN,MM C IF (J.EQ.JJ) THEN IF ( (N.EQ.NN).AND.(M.EQ.MM) ) THEN MATCH = .TRUE. ELSE IF ( (N.EQ.MM).AND.(M.EQ.NN) ) THEN MATCH = .TRUE. ELSE MATCH = .FALSE. ENDIF ELSE MATCH = .FALSE. ENDIF C RETURN END C C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION LUMINOSITY(LMNSTY,PROCESS,PERMINV) C DOUBLE PRECISION LMNSTY(7) INTEGER PROCESS,PERMINV(5) C C This function calculated the luminosity function needed in the C the four dimensional integrals. LUMINOSITY is the sum over C quark flavors of L(\vec a_A, \vec a_B, x_A, x_B), where the C sum runs over flavors a_A and a_B of the initial state particles C and over flavors a_1, a_2, a_3 of the final state particles, C consistent with PROCESS and PERMINV. For each case, the charge C conjugate luninosity is included also. When the sum over final C state particles is nontrivial, it produces a factor NFLAVOR -1 C (for process A) or 2*NFLAVOR (for process C). C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL C INTEGER A,B C C Initialize: C A = PERMINV(1) B = PERMINV(2) C C----------------------------- C IF (PROCESS.EQ.1) THEN C C Possible cases are (A,B) = (1,2),(1,3),(1,4),(1,5),(5,1). C IF (A.EQ.1) THEN C IF (B.EQ.2) THEN C C Luminosity for (PROCESS;A,B) = (1;1,2) i.e. (A,B) = (q,Q): C We use WNUM(101) to get the counting factor for a quark or antiquark. C LUMINOSITY = LMNSTY(1) C ELSE IF (B.EQ.3) THEN C C Luminosity for (PROCESS;A,B) = (1;1,3) i.e. (A,B) = (q,qbar): C LUMINOSITY = LMNSTY(4) C LUMINOSITY = (NFLAVOR - 1.0D0) * LUMINOSITY C ELSE IF (B.EQ.4) THEN C C Luminosity for (PROCESS;A,B) = (1;1,4) i.e. (A,B) = (q,Qbar): C LUMINOSITY = LMNSTY(2) C ELSE IF (B.EQ.5) THEN C C Luminosity for (PROCESS;A,B) = (1;1,5) i.e. (A,B) = (q,g): C LUMINOSITY = LMNSTY(5) C LUMINOSITY = (NFLAVOR - 1.0D0) * LUMINOSITY C ENDIF C ELSE IF ((A.EQ.5).AND.(B.EQ.1)) THEN C C Luminosity for (PROCESS;A,B) = (1;5,1) i.e. (A,B) = (g,q): C LUMINOSITY = LMNSTY(6) C LUMINOSITY = (NFLAVOR - 1.0D0) * LUMINOSITY C ENDIF C C------------------ C ELSE IF (PROCESS.EQ.2) THEN C C Possible cases are (A,B) = (1,2),(1,3),(1,5),(5,1). C IF (A.EQ.1) THEN C IF (B.EQ.2) THEN C C Luminosity for (PROCESS;A,B) = (2;1,2) i.e. (A,B) = (q,q): C LUMINOSITY = LMNSTY(3) C ELSE IF (B.EQ.3) THEN C C Luminosity for (PROCESS;A,B) = (2;1,3) i.e. (A,B) = (q,qbar): C LUMINOSITY = LMNSTY(4) C ELSE IF (B.EQ.5) THEN C C Luminosity for (PROCESS;A,B) = (2;1,5) i.e. (A,B) = (q,g): C LUMINOSITY = LMNSTY(5) C ENDIF C ELSE IF ((A.EQ.5).AND.(B.EQ.1)) THEN C C Luminosity for (PROCESS;A,B) = (2;5,1) i.e. (A,B) = (g,q): C LUMINOSITY = LMNSTY(6) C ENDIF C C------------------ C ELSE IF (PROCESS.EQ.3) THEN C C Possible cases are (A,B) = (1,2),(1,3),(3,1),(3,4). C IF (A.EQ.1) THEN C IF (B.EQ.2) THEN C C Luminosity for (PROCESS;A,B) = (3;1,2) i.e. (A,B) = (q,qbar): C LUMINOSITY = LMNSTY(4) C ELSE IF (B.EQ.3) THEN C C Luminosity for (PROCESS;A,B) = (3;1,3) i.e. (A,B) = (q,g): C LUMINOSITY = LMNSTY(5) C ENDIF C ELSE IF(A.EQ.3) THEN C IF (B.EQ.1) THEN C C Luminosity for (PROCESS;A,B) = (3;3,1) i.e. (A,B) = (g,q): C LUMINOSITY = LMNSTY(6) C ELSE IF (B.EQ.4) THEN C C Luminosity for (PROCESS;A,B) = (3;3,4) i.e. (A,B) = (g,g): C LUMINOSITY = LMNSTY(7) C LUMINOSITY = 2.0D0 * NFLAVOR * LUMINOSITY C ENDIF C ENDIF C C------------------------ C ELSE IF (PROCESS.EQ.4) THEN C C Only case is (A,B) = (g,g). C LUMINOSITY = LMNSTY(7) C C------------------------ C ENDIF C RETURN END C C*********************************************************************** C*********************************************************************** C SUBROUTINE PARTONSIN(LMNSTY,XA,XB,MUCO) C DOUBLE PRECISION LMNSTY(7) DOUBLE PRECISION XA,XB,MUCO C C This function computes the product, LMNSTY(KIND), of a parton C distribution function for A and a parton distribution function C for B, summed over flavors according to KIND: C A B C 1 = quark + quark (different flavors) C antiquark + antiquark (different flavors) C C 2 = quark + antiquark (different flavors) C antiquark + quark (different flavors) C C 3 = quark + quark (same flavor) C antiquark + antiquark (same flavor) C C 4 = quark + antiquark (same flavor) C antiquark + quark (same flavor) C C 5 = quark + gluon C antiquark + gluon C C 6 = gluon + quark C gluon + antiquark C C 7 = gluon + gluon C C INTEGER I,J,KIND DOUBLE PRECISION PARTON,PRTNS,WNUM C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL C Switch for p-pbar collisions (if true) or p-p collisions (if false) LOGICAL PPBAR COMMON /BEAMS/ PPBAR C DOUBLE PRECISION PRTNA(-6:+6),PRTNB(-6:+6) C C For p-pbar collisions, replace partons -> antipartons for hadron B. C IF (PPBAR) THEN DO I = -NFL,+NFL PRTNA(I) = PARTON( I,XA,MUCO)/XA/WNUM(I) PRTNB(I) = PARTON(-I,XB,MUCO)/XB/WNUM(I) ENDDO ELSE DO I = -NFL,+NFL PRTNA(I) = PARTON( I,XA,MUCO)/XA/WNUM(I) PRTNB(I) = PARTON( I,XB,MUCO)/XB/WNUM(I) ENDDO ENDIF C DO KIND = 1,7 C PRTNS = 0.0D0 C IF (KIND.EQ.1) THEN C i.e. (A,B) = (q,Q): C DO I = 1,NFL DO J = 1,NFL IF (J.NE.I) THEN PRTNS = PRTNS > + PRTNA( I) > *PRTNB( J) ENDIF ENDDO ENDDO C DO I = -NFL,-1 DO J = -NFL,-1 IF (J.NE.I) THEN PRTNS = PRTNS > + PRTNA( I) > *PRTNB( J) ENDIF ENDDO ENDDO C ELSE IF (KIND.EQ.2) THEN C i.e. (A,B) = (q,Qbar): C DO I = 1,NFL DO J = 1,NFL IF (J.NE.I) THEN PRTNS = PRTNS > + PRTNA( I) > *PRTNB(-J) ENDIF ENDDO ENDDO C DO I = -NFL,-1 DO J = -NFL,-1 IF (J.NE.I) THEN PRTNS = PRTNS > + PRTNA( I) > *PRTNB(-J) ENDIF ENDDO ENDDO C ELSE IF (KIND.EQ.3) THEN C i.e. (A,B) = (q,q): C DO I = 1,NFL PRTNS = PRTNS > + PRTNA( I) > *PRTNB( I) ENDDO C DO I = -NFL,-1 PRTNS = PRTNS > + PRTNA( I) > *PRTNB( I) ENDDO C ELSE IF (KIND.EQ.4) THEN C i.e. (A,B) = (q,qbar): C DO I = 1,NFL PRTNS = PRTNS > + PRTNA( I) > *PRTNB(-I) ENDDO C DO I = -NFL,-1 PRTNS = PRTNS > + PRTNA( I) > *PRTNB(-I) ENDDO C ELSE IF (KIND.EQ.5) THEN C i.e. (A,B) = (q,g): C DO I = 1,NFL PRTNS = PRTNS > + PRTNA( I) > *PRTNB( 0) ENDDO C DO I = -NFL,-1 PRTNS = PRTNS > + PRTNA( I) > *PRTNB( 0) ENDDO C ELSE IF (KIND.EQ.6) THEN C i.e. (A,B) = (g,q): C DO I = 1,NFL PRTNS = PRTNS > + PRTNA( 0) > *PRTNB( I) ENDDO C DO I = -NFL,-1 PRTNS = PRTNS > + PRTNA( 0) > *PRTNB( I) ENDDO C ELSE IF (KIND.EQ.7) THEN C i.e. (A,B) = (g,g): C PRTNS = PRTNA( 0) > *PRTNB( 0) C ENDIF C LMNSTY(KIND) = PRTNS C C End loop over KIND ENDDO C RETURN END C C*********************************************************************** C*********************************************************************** C END OF SUBROUTINES FOR 4 DIMENSIONAL INTEGRATIONS C*********************************************************************** C C*********************************************************************** C MISCELLANEOUS LIBRARY ROUTINES FOR JET CALCULATION C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION ALPHAS(Q) C DOUBLE PRECISION Q C INTEGER NFL DOUBLE PRECISION NFLAVOR,LAMBDA COMMON /QCDPARS/ NFLAVOR,LAMBDA,NFL C DOUBLE PRECISION PI PARAMETER ( PI = 3.141592654 D0) DOUBLE PRECISION FACTOR,LOG C FACTOR = 33 - 2.0D0*NFLAVOR LOG = DLOG( Q**2 / LAMBDA**2) ALPHAS = 12.0D0 * PI /FACTOR /LOG > - 72.0D0 * PI * (153.0D0 - 19.0D0*NFLAVOR) * DLOG(LOG) > /FACTOR**3 /LOG**2 C RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION WNUM(NPARTON) C INTEGER NPARTON C C This function gives the spin and color weighting for C parton NPARTON: 6 for quarks and antiquarks, 16 for gluons C IF (NPARTON.EQ.0) THEN WNUM = 16.0D0 ELSE WNUM = 6.0D0 ENDIF RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION COLOR(NPARTON) C INTEGER NPARTON C C This function gives the color charge for C parton NPARTON: 4/3 for quarks and antiquarks, 3 for gluons C DOUBLE PRECISION CA,CF PARAMETER( CA = 3.0D0 ) PARAMETER( CF = 4.0D0/3.0D0) C IF (NPARTON.EQ.0) THEN COLOR = CA ELSE COLOR = CF ENDIF RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION THETA(BOOL) C LOGICAL BOOL IF (BOOL) THEN THETA = 1.0D0 ELSE THETA = 0.0D0 ENDIF RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION FEXP(Y) C DOUBLE PRECISION Y C THE PURPOSE OF THIS FUNCTION IS JUST TO PREVENT A PROGRAM CRASH INTEGER NIN,NOUT,NWRT COMMON /IOUNIT/ NIN,NOUT,NWRT C IF (Y.GT.35.0D0) THEN FEXP=DEXP(35.0D0) WRITE(NWRT,101)Y ELSE IF (Y.LT.-35.0D0) THEN FEXP=DEXP(-35.0D0) WRITE(NWRT,101)Y ELSE FEXP=DEXP(Y) ENDIF RETURN 101 FORMAT(' Y = ',G10.2,'IN FEXP(Y); NO HARM IF THIS IS RARE.') END C C************************************************************************ C*********************************************************************** C Fitting subroutines for unsmearing calculation with FITJET C*********************************************************************** C*********************************************************************** C DOUBLE PRECISION FUNCTION UNSMEAR(M,J,DELTA) INTEGER M,J DOUBLE PRECISION DELTA C C Calculates the 'unsmearing coefficients' C C UNSMEAR(M,J,DELTA) = (-1 DELTA^2 /2)^(J-M) (2J)! / (2M)! /(J-M)! C C where DELTA is the smearing length. C DOUBLE PRECISION FACTORIAL C UNSMEAR = ( - DELTA**2 / 2.0D0 )**(J-M) > * FACTORIAL(2*J) / FACTORIAL(2*M) / FACTORIAL(J-M) C RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION DELTA(M,J) INTEGER M,J C C Kroneker delta C IF(M.EQ.J) THEN DELTA = 1.0D0 ELSE DELTA = 0.0D0 ENDIF C RETURN END C C*********************************************************************** C DOUBLE PRECISION FUNCTION FACTORIAL(N) INTEGER N C INTEGER J C FACTORIAL = 1.0D0 DO J = 1,N FACTORIAL = FACTORIAL * J ENDDO RETURN END C C*********************************************************************** C SUBROUTINE FITVAL(FUN,VAL,DEV,NPT,NPAR,PARAMS,CHISQ) INTEGER MAXPAR,MAXPT PARAMETER (MAXPAR=100) PARAMETER (MAXPT=300) C Input INTEGER NPT,NPAR DOUBLE PRECISION VAL(MAXPT),DEV(MAXPT) DOUBLE PRECISION FUN(MAXPAR,MAXPT) C Output DOUBLE PRECISION PARAMS(MAXPAR),CHISQ C C C Fit a function f(x) = sum(f_j(x)*params(j), j=1,npar) C at a set of points x(k) with k=1,npt by minimizing C the chisq = sum( ( f(x(k))-VAL(k) )^2 / DEV(k)^2, k=1,npt), C where VAL(k) is the value to fit at each point x(k), and DEV(k) C is its standard deviation. C C ARGUMENTS: C C FUN(J,K) : values of user given function f_j(x(k)) C VAL(K) : the set of y values. C DEV(K) : the standard deviations. C NPT : number of x values. C NPAR : number of parameters. C PARAMS(J): the parameters of the fit. -output C CHISQ : chisquare of the fit. - output C DOUBLE PRECISION GRAD(MAXPAR),HES(MAXPAR,MAXPAR) DOUBLE PRECISION TMP INTEGER J,K,J1,J2 C IF((NPAR.GT.MAXPAR).OR.(NPT.GT.MAXPT)) THEN WRITE(6,*) 'Error: NPT or NPAR too big in subroutine FITVAL' WRITE(6,*) 'NPT: ',NPT,' NPAR: ',NPAR STOP ENDIF C DO J = 1,NPAR GRAD(J) = 0.0D0 DO K=1,NPT GRAD(J) = GRAD(J) - 2.0D0*VAL(K)*FUN(J,K)/DEV(K)**2 ENDDO ENDDO C DO J1=1,NPAR DO J2=J1,NPAR HES(J1,J2) = 0.0D0 DO K=1,NPT HES(J1,J2) = HES(J1,J2) > + 2.0D0*FUN(J1,K)*FUN(J2,K)/DEV(K)**2 ENDDO HES(J2,J1)=HES(J1,J2) ENDDO ENDDO C C MMIN computes PARAMS according to C PARAMS(J) = - Sum HES^{-1}(J,L) * GRAD(L) C CALL MMIN(GRAD,HES,PARAMS,NPAR) C C Check that GRAD(J) = - Sum HES(J,L) * PARAMS(L): C C DO J1 = 1,NPAR C TMP = 0.0D0 C DO J2 = 1,NPAR C TMP = TMP - HES(J1,J2) * PARAMS(J2) C ENDDO C WRITE(6,*)TMP,GRAD(J1) C ENDDO C CHISQ = 0.0D0 DO K=1,NPT TMP = 0.0D0 DO J=1,NPAR TMP = TMP + PARAMS(J)*FUN(J,K) ENDDO CHISQ = CHISQ + (TMP-VAL(K))**2/DEV(K)**2 ENDDO CHISQ = CHISQ/NPAR C RETURN END C C------------------------------------------------------------------ C------------------------------------------------------------------ C SUBROUTINE MMIN(FUNGRD,FUNHES,DELVEC,N) C INTEGER N,MAXPAR PARAMETER(MAXPAR=100) DOUBLE PRECISION FUNGRD(MAXPAR),DELVEC(MAXPAR) DOUBLE PRECISION FUNHES(MAXPAR,MAXPAR) C DOUBLE PRECISION XMAT(MAXPAR,MAXPAR),CVEC(MAXPAR) INTEGER J,K,IMX,L,NTEST DOUBLE PRECISION TMP,XMX IF(N.GT.MAXPAR) THEN WRITE(6,*) 'Error: N>MAXPAR in subroutine MMIN' STOP ENDIF DO J=1,N CVEC(J) = - FUNGRD(J) DO K=1,N XMAT(K,J) = FUNHES(J,K) ENDDO ENDDO C C start elimination C loop over rows C DO J=1,N XMX = ABS(XMAT(1,J)) IMX = 1 C C Loop over columns: find maximum coef. for that row C DO K=1,N TMP = ABS(XMAT(K,J)) IF(TMP.GT.XMX) THEN XMX = TMP IMX = K ENDIF ENDDO IF(XMX.EQ.0) GOTO 110 C C Add row J to all others so to cancel column IMX C DO L=1,N IF(L.NE.J) THEN CVEC(L) = CVEC(L)-CVEC(J)*XMAT(IMX,L)/XMAT(IMX,J) DO K=1,N IF(K.NE.IMX) THEN XMAT(K,L) = # XMAT(K,L)-XMAT(K,J)*XMAT(IMX,L)/XMAT(IMX,J) ENDIF ENDDO XMAT(IMX,L) = 0 ENDIF ENDDO 110 ENDDO C DO J=1,N C WRITE(6,'(10(1X,D10.4))')(XMAT(K,J),K=1,N),CVEC(J) C ENDDO C C Now only one term per row is different from zero. C DO K=1,N NTEST = 0 DELVEC(K) = 0 DO J=1,N IF(XMAT(K,J).NE.0) THEN DELVEC(K) = CVEC(J)/XMAT(K,J) GOTO 111 ENDIF ENDDO 111 ENDDO RETURN END C C*********************************************************************** C SUBROUTINE RANDOMINIT(IRAN) INTEGER IRAN C C Initialize the shift-register random number generator from a random C seed IRAN, 0 <= IRAN < 259200, with the help of a portable "quick C and dirty" generator C INTEGER LBIT DOUBLE PRECISION FAC PARAMETER(LBIT=31,FAC= 2.0D0**(-LBIT)) C Here LBIT is the number of bits used in the shift register generator INTEGER J,IR(250) DOUBLE PRECISION RR(250) COMMON/RANDO/ RR,J,IR SAVE /RANDO/ INTEGER I,LB,JRAN,IFAC1,ISK,IDI,I1S C Configuration of the shift register generator INTEGER IM,IA,IC DATA IM,IA,IC /259200,7141,54773/ C C IF(IRAN.LT.0.OR.IRAN.GE.IM)STOP > 'RINI: IRAN OUT OF RANGE' JRAN=IRAN C Warm up the auxiliary generator a little DO I=1,10 JRAN=MOD(JRAN*IA+IC,IM) ENDDO IFAC1=((2**(LBIT-1)-1)*2+1)/IM DO I=2,250 JRAN=MOD(JRAN*IA+IC,IM) IR(I)=IFAC1*JRAN ENDDO C Guarantee LBIT linearly independent (over the field of 2 el.) C elements in IR(I): IR(1)=1 I=1 ISK=250/LBIT IDI=1 I1S=1 DO LB=1,LBIT-1 I=I+ISK IDI=IDI*2 I1S=I1S+IDI IR(I)=IR(I).OR.IDI IR(I)=IR(I).AND.I1S ENDDO CALL NEWRAN END C C*********************************************************************** C SUBROUTINE NEWRAN C C Fills IR(I),RR(I) with 250 new random numbers, resets J=0. C Increment J before use! C INTEGER LBIT DOUBLE PRECISION FAC PARAMETER(LBIT=31,FAC= 2.0D0**(-LBIT)) INTEGER J,IR(250) DOUBLE PRECISION RR(250) COMMON/RANDO/ RR,J,IR SAVE /RANDO/ C INTEGER N C DO N=1,103 IR(N)=IR(N+147).XOR.IR(N) RR(N)=FAC*(DBLE(IR(N)) + 0.5D0) ENDDO C DO N=104,250 IR(N)=IR(N-103).XOR.IR(N) RR(N)=FAC*(DBLE(IR(N)) + 0.5D0) ENDDO C J=0 END C C*********************************************************************** C DOUBLE PRECISION FUNCTION RANDOM(DUMMY) REAL DUMMY C C Random number between 2**(-32) and 1 - 2**(-32): C INTEGER J,IR(250) DOUBLE PRECISION RR(250) COMMON/RANDO/ RR,J,IR SAVE /RANDO/ C IF(J.GE.250)CALL NEWRAN J=J+1 RANDOM = RR(J) END C C*********************************************************************** C END OF LIBRARY ROUTINES FOR JET CALCULATION C***********************************************************************