C REMEMBER TO LINK WITH LATEST VERSION OF LIBEXCA.._AND_ LUDATA74.O
C SAMPLE CODE:
C      EVERY EVENT/CHANGE OF COUPLINGS
C      SET FREE COUPLINGS (DZ,XG,XZ,YG,YZ) IN AN DBLE-ARRAY WITH 5 ELEMENTS
C      CALL SETCON.. TO GET CONSTRAINS (IE. TO CONVERT TO ALPHA.. SCHEME)
C      CALL QEXMTX TO GET XSEC WITH (NON-)ZERO COUPLINGS, MW AND ALPHA_QED
C _IF_ YOU WANT TO INITIALIZE EXCALIBUR YOUR-SELF AND NOT USE VALUES
C GIVEN HERE, THEN USE: CALL INIEXCA....(SEE BELOW)
C     
C QEXMTX IS AN "ALPHA"-ROUTINE -- USE EXMTX FOR STAND-ALONE ROUTINE
      SUBROUTINE QEXMTX(IDIACL,IPARIN,ALQED,M_W,DTGC,XSEC)
C IPARIN: ALPHA TRACK NUMBERS OF FOUR FERMIONS
C IDIACL=IEXDIA=KREL (SEE INIEXCA)
C DTGC: DBLE-ARRAY WITH 5 ELEMENTS CONTAINING COUPLINGS (DZ,XG,XZ,YG,YZ)
      IMPLICIT NONE
      
      INCLUDE '/aleph/phy/qdecl.inc'
      INCLUDE '/aleph/phy/qcde.inc'
      
      INTEGER IDIACL,IPARIN(4)
      DOUBLE PRECISION DTGC(5),XSEC,DTEMP(5),PEVT(4,4),ALQED,M_W
      
      LOGICAL FIRST
      DATA FIRST /.TRUE./
      INTEGER IFSC,IDCL,I,IPAR(4),IPTRAN
      DOUBLE PRECISION S
* Original common
      common/area1/sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +             pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      double precision sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +   pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      
      SAVE FIRST
      
      INCLUDE '/aleph/phy/qmacro.inc'
      
      IDCL = IDIACL
      DTEMP(1)=M_W
      DTEMP(2)=ALQED
      IF (FIRST.AND.WMI.LT.0.0D0) THEN
        CALL INIEXCA(KUPRNT,IDCL,1,1,1,DTEMP(2),0.12D0,91.1888D0
     &    ,DTEMP(1),2.4974D0)
      ELSE
        CALL EXSET(DTEMP(1),DTEMP(2))
      ENDIF
      FIRST=.FALSE.
      IFSC = 0
      DO I=1,3,2
        IF (KTPCOD(IPARIN(I)).EQ.3.OR.KTPCOD(IPARIN(I)).EQ.6.OR
     &    .KTPCOD(IPARIN(I)).EQ.34.OR.KTPCOD(IPARIN(I)).EQ.58.OR
     &    .KTPCOD(IPARIN(I)).EQ.60.OR.KTPCOD(IPARIN(I)).EQ.62.OR
     &    .(KTPCOD(IPARIN(I)).GE.101.AND.KTPCOD(IPARIN(I)).LE.106)) THEN
          IFSC = IFSC+IPTRAN(KTPCOD(IPARIN(I)))*1000000/(100**(I-1))
     &      +IPTRAN(KTPCOD(IPARIN(I+1)))*1000000/(100**I)
        ELSE
          IFSC = IFSC+IPTRAN(KTPCOD(IPARIN(I+1)))*1000000/(100**(I-1))
     &      +IPTRAN(KTPCOD(IPARIN(I)))*1000000/(100**I)
        ENDIF
      ENDDO
      
      DO I=1,4
        IPAR(I)=IPARIN(I)
      ENDDO      
      CALL EXCCNV(IFSC,IPAR)
      S = DBLE(QMSQ4(IPAR(1),IPAR(2),IPAR(3),IPAR(4)))
      DO I=1,4
        PEVT(1,I) = DBLE(QX(IPAR(I)))
        PEVT(2,I) = DBLE(QY(IPAR(I)))
        PEVT(3,I) = DBLE(QZ(IPAR(I)))
        PEVT(4,I) = DBLE(QM(IPAR(I)))
      ENDDO
      
      CALL JBHSETPRO(1,IDCL,IFSC,DTGC)
      CALL JBHSQMTX(S,XSEC,IFSC,PEVT)
      
      RETURN
      END
      
      SUBROUTINE EXMTX(IDIACL,IPATY,PMOM,ALQED,M_W,DTGC,XSEC)
C IDIACL=IEXDIA=KREL (SEE INIEXCA)
C IPATY: ARRAY GIVING THE FERMION TYPE (ALPHA-code in Appendix F in the
C        alpha manual)
C PMOM: FOUR FERMION MOMENTA (px,py,pz,mass)
C DTGC: DBLE-ARRAY WITH 5 ELEMENTS CONTAINING COUPLINGS (DZ,XG,XZ,YG,YZ)
      IMPLICIT NONE
      
      INTEGER IDIACL,IPATY(4)
      DOUBLE PRECISION DTGC(5),XSEC,XSECSM,DTEMP(5),PEVT(4,4),PMOM(4,4)
     &  ,ALQED,M_W
      
      LOGICAL FIRST
      DATA FIRST /.TRUE./
      INTEGER IFSC,IDCL,I,IPAR(4),IPTRAN
      DOUBLE PRECISION S
* Original common
      common/area1/sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +             pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      double precision sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +   pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      
      SAVE FIRST
      
      IDCL = IDIACL
      DTEMP(1)=M_W
      DTEMP(2)=ALQED
      IF (FIRST.AND.WMI.LT.0.0D0) THEN
        CALL INIEXCA(6,IDCL,1,1,1,DTEMP(2),0.12D0,91.1888D0
     &    ,DTEMP(1),2.4974D0)
      ELSE
        CALL EXSET(DTEMP(1),DTEMP(2))
      ENDIF
      FIRST=.FALSE.
      IFSC = 0
      DO I=1,3,2
        IF (IPATY(I).EQ.3.OR.IPATY(I).EQ.6.OR
     &    .IPATY(I).EQ.34.OR.IPATY(I).EQ.58.OR
     &    .IPATY(I).EQ.60.OR.IPATY(I).EQ.62.OR
     &    .(IPATY(I).GE.101.AND.IPATY(I).LE.106)) THEN
          IFSC = IFSC+IPTRAN(IPATY(I))*1000000/(100**(I-1))
     &      +IPTRAN(IPATY(I+1))*1000000/(100**I)
        ELSE
          IFSC = IFSC+IPTRAN(IPATY(I+1))*1000000/(100**(I-1))
     &      +IPTRAN(IPATY(I))*1000000/(100**I)
        ENDIF
      ENDDO
      DO I=1,5
        DTEMP(I)=0.0D0
      ENDDO
      DO I=1,4
        IPAR(I)=I
        DTEMP(1)=DTEMP(1)+PMOM(1,I)
        DTEMP(2)=DTEMP(2)+PMOM(2,I)
        DTEMP(3)=DTEMP(3)+PMOM(3,I)
        DTEMP(4)=DTEMP(4)+SQRT(PMOM(1,I)*PMOM(1,I)+PMOM(2,I)*PMOM(2,I)
     &    +PMOM(3,I)*PMOM(3,I)+PMOM(4,I)*PMOM(4,I))
      ENDDO      
      CALL EXCCNV(IFSC,IPAR)
      S = MAX(0.0D0,DTEMP(4)*DTEMP(4)-DTEMP(1)*DTEMP(1)-DTEMP(2)
     &  *DTEMP(2)-DTEMP(3)*DTEMP(3))
      DO I=1,4
        PEVT(1,I) = PMOM(1,IPAR(I))
        PEVT(2,I) = PMOM(2,IPAR(I))
        PEVT(3,I) = PMOM(3,IPAR(I))
        PEVT(4,I) = PMOM(4,IPAR(I))
      ENDDO
      
      CALL JBHSETPRO(1,IDCL,IFSC,DTGC)
      CALL JBHSQMTX(S,XSEC,IFSC,PEVT)
      
      RETURN
      END
      
      INTEGER FUNCTION IPTRAN(IAPCOD)
      IMPLICIT NONE
      
c      INCLUDE 'excini.inc'
      INTEGER LMCWRT,IPROC,NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD,IEXUWT
     &  ,IEXCOU,IEVT,INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG,IEXQDW,IEXNDB
     &  ,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT,IEXINT,IEXFCR,IEXUCR,IEXUWP
      DOUBLE PRECISION ECM,GMU,ALPHAR,DANOMC,CUTS,SDVRT,EXCFAC,EXALSZ
     &  ,EXUMWT,EXUWEI
      COMMON /EXINII/ LMCWRT,IPROC(300),NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD
     &  ,IEXUWT,IEXCOU,IEVT(152),INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG
     &  ,IEXQDW,IEXNDB,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT(153)
     &  ,IEXINT,IEXFCR,IEXUCR,IEXUWP(152)
      COMMON /EXINID/ ECM,GMU,ALPHAR,DANOMC(14),CUTS(26,153) 
     &  ,EXCFAC,EXALSZ,EXUMWT(152),EXUWEI(152)
      
      INTEGER IAPCOD
      
      IF (IAPCOD.EQ.2.OR.IAPCOD.EQ.3) THEN
        IPTRAN=11
      ELSEIF (IAPCOD.EQ.5.OR.IAPCOD.EQ.6) THEN
        IPTRAN=13
      ELSEIF (IAPCOD.EQ.33.OR.IAPCOD.EQ.34) THEN
        IPTRAN=15
      ELSEIF (IAPCOD.EQ.58.OR.IAPCOD.EQ.59) THEN
        IPTRAN=12
      ELSEIF (IAPCOD.EQ.60.OR.IAPCOD.EQ.61) THEN
        IPTRAN=14
      ELSEIF (IAPCOD.EQ.62.OR.IAPCOD.EQ.63) THEN
        IPTRAN=16
      ELSEIF (IAPCOD.EQ.101.OR.IAPCOD.EQ.107) THEN
        IPTRAN=2
      ELSEIF (IAPCOD.EQ.102.OR.IAPCOD.EQ.108) THEN
        IPTRAN=1
      ELSEIF (IAPCOD.EQ.103.OR.IAPCOD.EQ.109) THEN
        IPTRAN=3
      ELSEIF (IAPCOD.EQ.104.OR.IAPCOD.EQ.110) THEN
        IPTRAN=4
      ELSEIF (IAPCOD.EQ.105.OR.IAPCOD.EQ.111) THEN
        IPTRAN=5
      ELSEIF (IAPCOD.EQ.106.OR.IAPCOD.EQ.102) THEN
        IPTRAN=6
      ELSE
        WRITE(LMCWRT,*)
     &    'Fatal error: Bad FERMION id... Program stopped'
        STOP        
      ENDIF
 999  RETURN
      END
          
      SUBROUTINE EXCCNV(IFSC,IPAR)
      IMPLICIT NONE
      
      INTEGER IFSC,IPAR(4)
      
      INTEGER I,H,ID1,ID2,ID3,ID4
* Final state codes -- using PDG numbering scheme (only positive numbers)
      INTEGER IFSCOD(151)
      DATA IFSCOD 
* Leptonic
     &  /11111212,11121413,11121615,12111314,12111516,13131414,15151616
     &  ,13141615,15161413,11111111,11111313,11111515,13131313,15151515
     &  ,13131515,11111414,11111616,12121313,12121515,16161313,14141515
     &  ,12121212,12121414,12121616,14141414,16161616,14141616
* Semi-leptonic
     &  ,11120201,11120403,12110102,12110304,13140201,13140403,14130102
     &  ,14130304,15160201,15160403,16150102,16150304,11110202,11110404
     &  ,11110101,11110303,11110505,13130202,13130404,15150202,15150404
     &  ,13130101,13130303,13130505,15150101,15150303,15150505,12120202
     &  ,12120404,12120101,12120303,12120505,14140202,14140404,16160202
     &  ,16160404,14140101,14140303,14140505,16160101,16160303,16160505
* Hadronic
     &  ,02020101,04040303,02010304,01020403,02020202,04040404,01010101
     &  ,03030303,05050505,02020404,02020303,02020505,04040101,04040505
     &  ,01010303,01010505,03030505
* New mixing processes - first semi-leptonic...
     &  ,11120203,11120205,11120401,11120405,12110302,12110502,12110104
     &  ,12110504,13140203,13140205,13140401,13140405,14130302,14130502
     &  ,14130104,14130504,15160203,15160205,15160401,15160405,16150302
     &  ,16150502,16150104,16150504,12120204,12120402,12120103,12120105
     &  ,12120305,12120301,12120501,12120503
* Hadronic
     &  ,02010302,02010502,02030102,02050102,02050302,02030502,04030104
     &  ,04030504,04010304,04050304,04050104,04010504,02030304,02050304
     &  ,02010104,02010504,02030104,02030504,02050104,02050504,03020403
     &  ,05020403,01020401,01020405,03020401,03020405,05020401,05020405
* QQgg
     &  ,02022121,04042121,01012121,03032121,05052121/
*
* Translate specified channels into EXCALIBUR codes and remove redundant
* codes 
*
      IF (IFSC.LE.0) GOTO 60
      DO H=1,151
        IF (IFSC.EQ.IFSCOD(H)) GOTO 60
        IF (IFSC.NE.IFSCOD(H)) THEN
          ID1=INT(IFSCOD(H)/1000000)
          ID2=MOD(INT(IFSCOD(H)/10000),100)
          ID3=MOD(INT(IFSCOD(H)/100),100)
          ID4=MOD(IFSCOD(H),100)
          IF (ID1.EQ.ID2.AND.ID3.EQ.ID4.AND.MOD(MAX(ID1,ID3),2).EQ.0
     &      .AND.MAX(ID1,ID3).EQ.MIN(ID1,ID3)+1) THEN
            IF ((IFSC.EQ.(ID1*1000000+ID3*10000+ID4*100+ID2)).OR
     &        .(IFSC.EQ.(ID3*1000000+ID4*10000+ID1*100+ID2)).OR
     &        .(IFSC.EQ.(ID3*1000000+ID1*10000+ID2*100+ID4)))
     &        THEN
              IF (IFSC.EQ.(ID1*1000000+ID3*10000+ID4*100+ID2))
     &          THEN
                I=IPAR(2)
                IPAR(2)=IPAR(4)
                IPAR(4)=I
                I=IPAR(3)
                IPAR(3)=IPAR(4)
                IPAR(4)=I
              ELSEIF (IFSC.EQ.(ID3*1000000+ID4*10000+ID1*100+ID2
     &            )) THEN
                I=IPAR(1)
                IPAR(1)=IPAR(3)
                IPAR(3)=I
                I=IPAR(2)
                IPAR(2)=IPAR(4)
                IPAR(4)=I
              ELSEIF (IFSC.EQ.(ID3*1000000+ID1*10000+ID2*100+ID4
     &            )) THEN
                I=IPAR(1)
                IPAR(1)=IPAR(2)
                IPAR(2)=I
                I=IPAR(2)
                IPAR(2)=IPAR(3)
                IPAR(3)=I
              ENDIF
              IFSC=IFSCOD(H)
            ENDIF
          ELSEIF (((ID1.EQ.ID2.AND.ID3.EQ.ID4).OR.(MOD(MAX(ID1,ID2)
     &        ,2).EQ.0.AND.MAX(ID1,ID2).EQ.MIN(ID1,ID2)+1.OR
     &        .MOD(MAX(ID3,ID4),2).EQ.0.AND.MAX(ID3,ID4).EQ.MIN(ID3
     &        ,ID4)+1)).AND.ID1.NE.ID3.AND.IFSC.EQ.(ID3*1000000
     &        +ID4*10000+ID1*100+ID2)) THEN
            I=IPAR(1)
            IPAR(1)=IPAR(3)
            IPAR(3)=I
            I=IPAR(2)
            IPAR(2)=IPAR(4)
            IPAR(4)=I
            IFSC=IFSCOD(H)
          ENDIF
        ENDIF
        IF (IFSC.EQ.IFSCOD(H)) GOTO 60
      ENDDO
 60   CONTINUE
 999  RETURN
      END
      
      SUBROUTINE JBHSETPRO(IFLAG,IDCL,IFSC,ALPHA)
      IMPLICIT NONE
      
      INTEGER IFLAG,IDCL,IFSC
      DOUBLE PRECISION ALPHA(5)
      
C      INCLUDE 'excini.inc'
      INTEGER LMCWRT,IPROC,NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD,IEXUWT
     &  ,IEXCOU,IEVT,INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG,IEXQDW,IEXNDB
     &  ,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT,IEXINT,IEXFCR,IEXUCR,IEXUWP
      DOUBLE PRECISION ECM,GMU,ALPHAR,DANOMC,CUTS,SDVRT,EXCFAC,EXALSZ
     &  ,EXUMWT,EXUWEI
      COMMON /EXINII/ LMCWRT,IPROC(300),NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD
     &  ,IEXUWT,IEXCOU,IEVT(152),INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG
     &  ,IEXQDW,IEXNDB,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT(153)
     &  ,IEXINT,IEXFCR,IEXUCR,IEXUWP(152)
      COMMON /EXINID/ ECM,GMU,ALPHAR,DANOMC(14),CUTS(26,153) 
     &  ,EXCFAC,EXALSZ,EXUMWT(152),EXUWEI(152)
      
      CHARACTER*2 PAR(3:6)
      INTEGER ISTAT,NCM(1:150),NNM(1:150),NCT,NOPRO,IPRC
      DOUBLE PRECISION DZ,XG,YG,XZ,YZ,ZZ
      COMMON/AREA11/DZ,XG,YG,XZ,YZ,ZZ

      DZ=ALPHA(1)
      XG=ALPHA(2)
      XZ=ALPHA(3)
      YG=ALPHA(4)
      YZ=ALPHA(5)
      ZZ=0.0D0
      IF (IFLAG.NE.1) GOTO 999
*
* Decode final state
*
      IPRC=IFSC
      IF (IEXMIX.GT.0) CALL EXCDFS(IFSC,IPRC)
      CALL EXCFID(IPRC,PAR,ISTAT)
      IF (ISTAT.NE.0) THEN
        WRITE(LMCWRT,*)
     &    'Fatal error: Bad process id... Program stopped'
        STOP
      ENDIF
      CALL SETPRO(NCM,NNM,NCT,IDCL,NOPRO,PAR)
      IF (NOPRO.EQ.1) THEN
        WRITE(LMCWRT,*)
     &    'Fatal error: Wrong process... Program stopped'
        STOP
      ENDIF
 999  RETURN
      END
      
      SUBROUTINE JBHSQMTX(S,SQMTX,IFSC,PEVT)
      IMPLICIT NONE
      
      INTEGER IFSC
      DOUBLE PRECISION S,SQMTX,PEVT(4,4)
      
* Original common
      common/area1/sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +             pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      double precision sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +   pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      common/momenta/Droots,xr1,xr2,pm(0:4,0:900)
      double precision Droots,xr1,xr2,pm
* Additional commons...
C      INCLUDE 'excini.inc'
      INTEGER LMCWRT,IPROC,NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD,IEXUWT
     &  ,IEXCOU,IEVT,INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG,IEXQDW,IEXNDB
     &  ,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT,IEXINT,IEXFCR,IEXUCR,IEXUWP
      DOUBLE PRECISION ECM,GMU,ALPHAR,DANOMC,CUTS,SDVRT,EXCFAC,EXALSZ
     &  ,EXUMWT,EXUWEI
      COMMON /EXINII/ LMCWRT,IPROC(300),NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD
     &  ,IEXUWT,IEXCOU,IEVT(152),INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG
     &  ,IEXQDW,IEXNDB,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT(153)
     &  ,IEXINT,IEXFCR,IEXUCR,IEXUWP(152)
      COMMON /EXINID/ ECM,GMU,ALPHAR,DANOMC(14),CUTS(26,153) 
     &  ,EXCFAC,EXALSZ,EXUMWT(152),EXUWEI(152)
*     
      DOUBLE PRECISION M1,M2
      COMMON /BMASS/ M1,M2
      INTEGER IEXCLM,IEXABD,IEXNAD,NL,WWDIA
      COMMON /EXCSNG/ IEXCLM,IEXABD(144),IEXNAD(8),NL,WWDIA
*
      INTEGER I,ID(4)
      DOUBLE PRECISION RES,DTEMP1,DTEMP2,DTEMP3,PI,XFAC
      
      INTEGER IMIX(6)
      DOUBLE PRECISION QCD(2)
      
      PARAMETER (PI=3.14159265358979324D0)
      
      IF (STH2.LE.0.0D0) THEN
        WRITE(LMCWRT,*)
     &    'Fatal error: EXCALIBUR is not initialised.. Program stopped'
        STOP
      ENDIF
      
      IMIX(1)=0
      IF (IEXMIX.GT.0) THEN
        IMIX(1)=1
        IMIX(2)=1
        IMIX(3)=INT(IFSC/1000000)
        IMIX(4)=MOD(INT(IFSC/10000),100)
        IMIX(5)=MOD(INT(IFSC/100),100)
        IMIX(6)=MOD(IFSC,100)
        DO I=3,6
          IF (IMIX(I).LT.7) THEN
            IMIX(I)=IMIX(I)+6
          ELSEIF (IMIX(I).GT.10.AND.IMIX(I).LT.17) THEN
            IMIX(I)=IMIX(I)-10
          ENDIF
        ENDDO
      ENDIF
      QCD(1)=1.0D0
      QCD(2)=1.0D0
      IF (IEXQCD.GT.0) CALL EXQCDC(QCD)
      IEXCLM=IEXCOU
*     
      SQMTX=0.0D0
      ID(1)=INT(IFSC/1000000)
      ID(2)=MOD(INT(IFSC/10000),100)
      ID(3)=MOD(INT(IFSC/100),100)
      ID(4)=MOD(IFSC,100)
      CALL JBHSETMOM(S,PEVT)
      IF ((ID(1).NE.ID(2).OR.(ID(1).EQ.ID(2).AND.ID(1).GE.11
     &  .AND.ID(3).EQ.ID(1)+1).OR.(ID(1).EQ.ID(2).AND.ID(1).LE
     &  .6.AND.ID(3).EQ.ID(1)-1)).AND.ID(3).NE.21) THEN
        M1=SQRT(PM(4,34))
        M2=SQRT(PM(4,56))
        IF (ID(1).EQ.ID(2)) THEN
          M1=SQRT(PM(4,36))
          M2=SQRT(PM(4,45))
        ENDIF
      ENDIF
*
* Normalization factor
*
      XFAC= PI*PI*16.0D0*ALPHA*ALPHA/STH2/STH2
      FACNOR= XFAC*XFAC*FCOL*PBFAC*STATFAC/2048.0D0/PI**8
      CALL EXCMTX(RES,DTEMP1,DTEMP2,DTEMP3,IMIX,QCD)
      SQMTX=FACNOR*RES/DROOTS/DROOTS
 999  RETURN
      END
      
      SUBROUTINE JBHSETMOM(S,PEVT)
      IMPLICIT NONE
      
      INTEGER I,J,K,L,L1,L2,L3
      DOUBLE PRECISION S,PEVT(4,4),DTEMP(0:4),DTEMP1(0:4)
      COMMON/MOMENTA/ROOTS,XR1,XR2,PM(0:4,0:900)
      DOUBLE PRECISION ROOTS,XR1,XR2,PM

      ROOTS=SQRT(S)
      XR1=1.0D0
      XR2=1.0D0
      PM(0,1)= ROOTS/2.D0
      PM(1,1)= PM(0,1)

      PM(0,2)= PM(0,1)
      PM(1,2)= -PM(0,1)

      PM(0,12)= PM(0,1)+PM(0,2)
      PM(1,12)= PM(1,1)+PM(1,2)
      PM(2,12)= PM(2,1)+PM(2,2)
      PM(3,12)= PM(3,1)+PM(3,2)
      PM(4,12)= 4.0D0*PM(0,1)*PM(0,1)
      DO J= 0,4
        PM(J,21)= PM(J,12)
        DTEMP(J)=0.0D0
      ENDDO
      
      DO 5 I= 3,6
        PM(1,I)=-PEVT(3,I-2)
        PM(2,I)=PEVT(2,I-2)
        PM(3,I)=PEVT(1,I-2)
        PM(4,I)=PEVT(4,I-2)
        PM(0,I)=SQRT(PEVT(1,I-2)*PEVT(1,I-2)+PEVT(2,I-2)*PEVT(2,I-2)
     &    +PEVT(3,I-2)*PEVT(3,I-2)+PEVT(4,I-2)*PEVT(4,I-2))
        DTEMP(0)=DTEMP(0)+PM(0,I)
        DTEMP(1)=DTEMP(1)+PM(1,I)
        DTEMP(2)=DTEMP(2)+PM(2,I)
        DTEMP(3)=DTEMP(3)+PM(3,I)
 5    CONTINUE
C NOW BOOST INTO CM-FRAME OF FOUR FERMION SYSTEM
      DO 6 I= 3,6
        CALL BOOST(1,DTEMP,DTEMP1,PM(0,I))
        PM(1,I)=DTEMP1(1)
        PM(2,I)=DTEMP1(2)
        PM(3,I)=DTEMP1(3)
        PM(4,I)=0.0D0 !ZERO MASS
        PM(0,I)=SQRT(DTEMP1(1)*DTEMP1(1)+DTEMP1(2)*DTEMP1(2)+DTEMP1(3)
     &    *DTEMP1(3))
 6    CONTINUE
      
      DO 20 L= 0,4
      DO 10 I= 1,6
        DO 10 J= 1,6
          IF (I.NE.J) THEN
            L1= I*10+J
            IF (L.LE.3) THEN
              PM(L,L1)= PM(L,I)+PM(L,J)
            ELSE
              PM(L,L1)= PM(0,L1)*PM(0,L1)-PM(1,L1)*PM(1,L1)-
     +                  PM(2,L1)*PM(2,L1)-PM(3,L1)*PM(3,L1)
            ENDIF
          ENDIF
  10  CONTINUE
      DO 20 I= 1,6
        DO 20 J= 1,6
          DO 20 K= 1,6
            IF (I.NE.J.AND.I.NE.K.AND.K.NE.J) THEN
              L2= I*100+J*10+K
              IF (L.LE.3) THEN
                PM(L,L2)= PM(L,I)+PM(L,J*10+K)
              ELSE
                PM(4,L2)= PM(0,L2)*PM(0,L2)-PM(1,L2)*PM(1,L2)-
     +            PM(2,L2)*PM(2,L2)-PM(3,L2)*PM(3,L2)
              ENDIF
              IF (I.LE.2) THEN
                L3= L2+600
                IF (L.LE.3) THEN
                  PM(L,L3)= PM(L,I)-PM(L,J*10+K)
                ELSE
                  PM(4,L3)= PM(0,L3)*PM(0,L3)-PM(1,L3)*PM(1,L3)-
     +              PM(2,L3)*PM(2,L3)-PM(3,L3)*PM(3,L3)
                ENDIF
              ENDIF
            ENDIF
 20       CONTINUE
 999  RETURN
      END
      
      SUBROUTINE EXSET(M_W,AL_QED)
      IMPLICIT NONE
      
      DOUBLE PRECISION M_W,AL_QED
      
* Original common
      common/area1/sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +             pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      double precision sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +   pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
* Additional commons...
C      INCLUDE 'excini.inc'
      INTEGER LMCWRT,IPROC,NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD,IEXUWT
     &  ,IEXCOU,IEVT,INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG,IEXQDW,IEXNDB
     &  ,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT,IEXINT,IEXFCR,IEXUCR,IEXUWP
      DOUBLE PRECISION ECM,GMU,ALPHAR,DANOMC,CUTS,SDVRT,EXCFAC,EXALSZ
     &  ,EXUMWT,EXUWEI
      COMMON /EXINII/ LMCWRT,IPROC(300),NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD
     &  ,IEXUWT,IEXCOU,IEVT(152),INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG
     &  ,IEXQDW,IEXNDB,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT(153)
     &  ,IEXINT,IEXFCR,IEXUCR,IEXUWP(152)
      COMMON /EXINID/ ECM,GMU,ALPHAR,DANOMC(14),CUTS(26,153) 
     &  ,EXCFAC,EXALSZ,EXUMWT(152),EXUWEI(152)
C      INCLUDE 'excgen.inc'
      INTEGER NCONP,NNCONP,IEXWNL,IEXWNG,IGEVT,IFEVT,IEXCCF,IEXTPL
     &  ,IEXOFP,IEXPHG,IEXPRO
      COMMON /EXGENI/ NCONP,NNCONP,IEXWNL(152),IEXWNG(152),IGEVT(152)
     &  ,IFEVT(152),IEXCCF,IEXTPL(4),IEXOFP(15),IEXPHG,IEXPRO
      DOUBLE PRECISION EXXSEC,EXXERR,EXXSRN,EXWMAX,EXCUMC,EXSUMW
     &  ,EXSWSQ,EXMWEI,CL,ZETA,ZETA1,OMZ1,SHCUT,SAFETY,EXGMWT
     &  ,GAEUL,EXS,EXSWS3,EXSWS4,W,EX1234,EXINT,EX1432,VERTEX,EXHELI
     &  ,EXCSIG
      COMMON /EXGEND/ EXXSEC(152),EXXERR(152),EXXSRN(152),EXWMAX(152)
     &  ,EXCUMC(152),EXSUMW(152),EXSWSQ(152),EXSWS3(152),EXSWS4(152)
     &  ,EXMWEI(152),EXS(152),CL,ZETA,ZETA1,OMZ1,SHCUT(153),EXGMWT(153)
     &  ,GAEUL,SAFETY,W,EX1234,EXINT,EX1432,VERTEX(4),EXHELI(16)
     &  ,EXCSIG(3,8)

      
      DOUBLE PRECISION QCD(2),STH2R,PI
      PARAMETER (PI=3.14159265358979324D0)
*
* Unit number for output messages: LMCWRT
*
* Processing options:
*           Contributing diagrams (IEXDIA=Kleiss def.)
*           Debug printout (in org. EXCALIBUR code) (IEXDBG=1)
*           Unweighted events (IEXUWT=1)
*           ISR (IEXISR=1)
*           FSR (IEXFSR=1 (PHOTOS) or IEXFSR=2 (FERMISV))
*           Hadronization (IEXHAD=1)
*           Coulomb singularity (IEXCOU=1)
*           CKM-mixing (IEXMIX=1)
*           QCD corrections to boson widths (IEXQCD=1 (W))
*           Specified boson widths are EL-WEAK widths (IEXEWD=1)
*           Use TAUOLA for tau-decays (MSTJ(28)=2)
*
* IEXDIA  = 0  all diagrams
*         = 1  WW diagram
*         = 2  ZZ diagrams
*         = 3  W+ e- nue_bar production
*         = 4  Z e+ e- production
*         = 5  Z nue nue_bar production
* STH2.EQ.0 Forces use of value consistent with EXCALIBUR-scheme ("GMU")
      STH2=0.0D0
      WMI=M_W
* WWI.LE.0 Forces use of value consistent with EXCALIBUR-scheme ("GMU")
      WWI=-1.0D0
      ALPHA=AL_QED
      ALS=0.0D0
      STH2=PI*ALPHA/(SQRT(2.0D0)*WMI**2*GMU)
      WWI=3.0D0*GMU*WMI**3/SQRT(8.0D0)/PI
      IEXEWD=1
*
* Report Alpha_s corrected W boson width option...
*
      IF (IEXQCD.GT.0.AND.MOD(IEXQCD,2).EQ.1.AND.IEXEWD.EQ.1) THEN
        CALL EXQCDC(QCD)
* JBH Update total width of the W
        WWI=WWI*(1.0D0+2.0D0*QCD(2)*QCD(2))/3.0D0
      ENDIF
      RETURN
      END
      
      SUBROUTINE INIEXCA(UNIT,KREL,ICOU,IQDW,IQCD
     &  ,AL_QED,ALSZ,M_Z,M_W,G_Z)
C sample call: call iniexca(kuprnt,0,1,1,1,1.0d0/128.07d0,0.12d0
C ,91.1888d0,80.25d0,2.4974d0)
      IMPLICIT NONE
      
      INTEGER UNIT,KREL,ICOU,IQDW,IQCD,IEWD,IMIX,IGG
      DOUBLE PRECISION ENERGY,AL_QED,ALPHAS,ALSZ,M_Z,M_W,G_Z,G_W,SW2
      
* Original common
      common/area1/sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +             pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
      double precision sth2,zm,wz,wm,ww,alpha,grap,statfac,
     +   pbfac,fcol,facnor,als,zmi,wzi,wmi,wwi
* Additional commons...
C      INCLUDE 'excini.inc'
      INTEGER LMCWRT,IPROC,NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD,IEXUWT
     &  ,IEXCOU,IEVT,INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG,IEXQDW,IEXNDB
     &  ,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT,IEXINT,IEXFCR,IEXUCR,IEXUWP
      DOUBLE PRECISION ECM,GMU,ALPHAR,DANOMC,CUTS,SDVRT,EXCFAC,EXALSZ
     &  ,EXUMWT,EXUWEI
      COMMON /EXINII/ LMCWRT,IPROC(300),NEVT,IEXDIA,IEXISR,IEXFSR,IEXHAD
     &  ,IEXUWT,IEXCOU,IEVT(152),INFAIL,IEXMIX,IEXPAF,IEXCRC,IEXNNG
     &  ,IEXQDW,IEXNDB,IEXDBG,IEXQCD,IEXEWD,IEXGG,IEXCUT(153)
     &  ,IEXINT,IEXFCR,IEXUCR,IEXUWP(152)
      COMMON /EXINID/ ECM,GMU,ALPHAR,DANOMC(14),CUTS(26,153) 
     &  ,EXCFAC,EXALSZ,EXUMWT(152),EXUWEI(152)
C      INCLUDE 'excgen.inc'
      INTEGER NCONP,NNCONP,IEXWNL,IEXWNG,IGEVT,IFEVT,IEXCCF,IEXTPL
     &  ,IEXOFP,IEXPHG,IEXPRO
      COMMON /EXGENI/ NCONP,NNCONP,IEXWNL(152),IEXWNG(152),IGEVT(152)
     &  ,IFEVT(152),IEXCCF,IEXTPL(4),IEXOFP(15),IEXPHG,IEXPRO
      DOUBLE PRECISION EXXSEC,EXXERR,EXXSRN,EXWMAX,EXCUMC,EXSUMW
     &  ,EXSWSQ,EXMWEI,CL,ZETA,ZETA1,OMZ1,SHCUT,SAFETY,EXGMWT
     &  ,GAEUL,EXS,EXSWS3,EXSWS4,W,EX1234,EXINT,EX1432,VERTEX,EXHELI
     &  ,EXCSIG
      COMMON /EXGEND/ EXXSEC(152),EXXERR(152),EXXSRN(152),EXWMAX(152)
     &  ,EXCUMC(152),EXSUMW(152),EXSWSQ(152),EXSWS3(152),EXSWS4(152)
     &  ,EXMWEI(152),EXS(152),CL,ZETA,ZETA1,OMZ1,SHCUT(153),EXGMWT(153)
     &  ,GAEUL,SAFETY,W,EX1234,EXINT,EX1432,VERTEX(4),EXHELI(16)
     &  ,EXCSIG(3,8)

      
      DOUBLE PRECISION QCD(2),STH2R,PI
      PARAMETER (PI=3.14159265358979324D0)
*
* Unit number for output messages: LMCWRT
*
* Processing options:
*           Contributing diagrams (IEXDIA=Kleiss def.)
*           Debug printout (in org. EXCALIBUR code) (IEXDBG=1)
*           Unweighted events (IEXUWT=1)
*           ISR (IEXISR=1)
*           FSR (IEXFSR=1 (PHOTOS) or IEXFSR=2 (FERMISV))
*           Hadronization (IEXHAD=1)
*           Coulomb singularity (IEXCOU=1)
*           CKM-mixing (IEXMIX=1)
*           QCD corrections to boson widths (IEXQCD=1 (W))
*           Specified boson widths are EL-WEAK widths (IEXEWD=1)
*           Use TAUOLA for tau-decays (MSTJ(28)=2)
*
* IEXDIA  = 0  all diagrams
*         = 1  WW diagram
*         = 2  ZZ diagrams
*         = 3  W+ e- nue_bar production
*         = 4  Z e+ e- production
*         = 5  Z nue nue_bar production
C set defaults
      CALL EXCDEF
C
      LMCWRT=UNIT
      IEXDIA=KREL
      IEXCOU=ICOU
      IEXMIX=1
* Q-dependent width in boson-propagators?
      IEXQDW=IQDW
      IEXQCD=IQCD
      IEXEWD=1
* Include QQgg? ---- Test option only for cross-section calculations...
      IEXGG=0
      GMU=1.16639D-5
      ALPHAR=1.0D0/137.0359895D0
      ECM=172.0D0 ! NOT USED ANYWAY...
      EXALSZ=ALSZ
* STH2.EQ.0 Forces use of value consistent with EXCALIBUR-scheme ("GMU")
      STH2=0.0D0
      ZMI=M_Z
      WMI=M_W
      WZI=G_Z
* WWI.LE.0 Forces use of value consistent with EXCALIBUR-scheme ("GMU")
      WWI=-1.0D0
      ALPHA=AL_QED
      ALS=0.0D0
*
* Report processing options..
*
      WRITE(LMCWRT,1080) IEXDIA,(IEXCOU.GT.0),(IEXQCD.GT.0),(IEXQDW.GT.0
     &  ),(IEXGG.GT.0)
 1080 FORMAT(/,1X,'Processing options:',/,
     +  '       Krel option?                        :',I2,/,
     +  '       Coulomb-singularity?                :',L2,/,
     +  '       QCD correction to W width?          :',L2,/,
     +  '       Q dependent width?                  :',L2,/,
     +  '       Include QQgg contribution?          :',L2)
*
* Consistency check of Standard model parameters needed here
*
      IF (STH2.LE.0.0D0) THEN
        STH2=PI*ALPHA/(SQRT(2.0D0)*WMI**2*GMU)
      ELSE
        STH2R=PI*ALPHA/(SQRT(2.0D0)*WMI**2*GMU)
        IF (ABS(STH2-STH2R).GT.0.0001D0) THEN
          WRITE(LMCWRT,1070) STH2,STH2R
 1070     FORMAT(1X,'WARNING: sin2w value supplied:',F8.5
     &      ,' inconsistent with',/,1X
     &      ,'         EXCALIBUR recommended scheme, gives:',F8.5)
        ENDIF
      ENDIF
      SW2=STH2
      IF (WWI.LE.0.0D0) THEN
        WWI=3.0D0*GMU*WMI**3/SQRT(8.0D0)/PI
        IEXEWD=1
      ENDIF
      WRITE(LMCWRT,1040) GMU,STH2,ZMI,WZI,WMI,WWI,1.0D0/ALPHA,ALS
 1040 FORMAT(1X,'Standard model input parameters:',/
     &  ,'       G_F      = ',G14.8/,'       s^2_thet = ',G14.8/
     &  ,'       Z-mass   = ',G14.8/,'       Z-width  = ',G14.8/
     &  ,'       W-mass   = ',G14.8/,'       W-width  = ',G14.8/
     &  ,'       1/alpha  = ',G14.8/,'       alpha_s  = ',G14.8)
*
* Report Alpha_s corrected W boson width option...
*
      IF (IEXQCD.GT.0.AND.MOD(IEXQCD,2).EQ.1.AND.IEXEWD.EQ.1) THEN
        WRITE(LMCWRT,1082) EXALSZ
        CALL EXQCDC(QCD)
* JBH Update total width of the W
        WWI=WWI*(1.0D0+2.0D0*QCD(2)*QCD(2))/3.0D0
        WRITE(LMCWRT,1087) (IEXQCD+1)/2,WWI
      ELSEIF (IEXQCD.GT.0.AND.MOD(IEXQCD,2).EQ.1.AND.IEXEWD.EQ.0) THEN
        WRITE(LMCWRT,1082) EXALSZ
        WRITE(LMCWRT,1086) (IEXQCD+1)/2
      ELSE
        IF (IEXQCD.GT.0) IEXQCD=0
        WRITE(LMCWRT,1081)
      ENDIF
 1081 FORMAT(/,4X,'QCD corrections to boson widths are NOT applied',/)
 1082 FORMAT(/,4X,'QCD corrections to boson widths are applied using:',/
     &  ,/,10X,'Alpha_s(Mz) = ',G12.5,/)
 1086 FORMAT(4X
     &  ,'The specified total width of the W-boson remain unchanged,',/
     &  ,4X
     &  ,'while the partial widths of the boson are recalculated using'
     &  ,/,4X,I1,'. order QCD corrections.',/)
 1087 FORMAT(4X
     &  ,'The total and partial widths of the W-boson are calculated',/
     &  ,4X,'from the specified electro-weak width using',/,4X
     &  ,I1,'. order QCD corrections. The total widths are: ',//,10X
     &  ,'Gamma_W = ',G14.8,/)
      RETURN
      END
      
      SUBROUTINE SETCON (TYPE,ALPHA,SW,CW)
      IMPLICIT NONE
      
C THESE CONSTRAINS COME FROM THE LIST ON PAGE 30+31 BY BILENKY ET AL IN NUC.
C PHYS B409 (22-68). THE NUMBER REFER TO THEIR APPEARANCE IN THE LIST.
C WITH THIS NUMBERING SCHEME ALPHA_W_PHI IS 3, ALPHA_B_PHI IS 4 AND
C ALPHA_W IS 11
      
      INTEGER TYPE
      DOUBLE PRECISION ALPHA(5),SW,CW
      
      IF (TYPE.EQ.3.OR.TYPE.EQ.6.OR.TYPE.EQ.10) ALPHA(1) =
     &  ALPHA(2)/SW/CW
      IF (TYPE.EQ.6) ALPHA(1) = ALPHA(1)/2.0D0
      IF ((TYPE.GE.2.AND.TYPE.LE.4).OR.TYPE.EQ.6.OR.(TYPE.GE.8
     &  .AND.TYPE.LE.10)) ALPHA(3) = -ALPHA(2)*SW/CW
      IF (TYPE.GE.9.AND.TYPE.LE.11) ALPHA(5) = ALPHA(4)*CW/SW
 999  RETURN
      END
