C----------------------------------------------------------------------
C
C                        ##       ##
C                       ###     #####
C                      ######     #####
C                        ######    #####
C                         /######   /###
C                        /   ######/
C                      ####    ######
C                        ####    ######
C                        ####       ####
C                        ###       ##
C                       ##        #
C
C
C----------------------------------------------------------------------
C
      PROGRAM FIT
C----------------------------------------------------------------------
C!  -
C!
C!   Authors   :- Manel Martinez        JUN-1992
C!                Frederic Teubert      FEB-1995
C!
C!   Outputs: FOR006
C!
C!   Libraries required: CERNLIB
C!
C!   Description
C!   ===========
C!
C!      * Z0 lineshape + AFB fit
C!
C!======================================================================
C  Book arrays
      IMPLICIT DOUBLE PRECISION(A-H,M,O-Z)
      REAL*8 RT,JT,RFB,JFB,RT0,JT0,RFB0,JFB0
      CHARACTER*10,CHNAM0,CHNAM(16),PARNAM(5)
C Fitted parameters
      COMMON /UNKNOWN/ MZ,GZ,RT(4),JT(4),RFB(3),JFB(3)
      COMMON /INITIAL/ RT0(4),RFB0(3),JT0(4),JFB0(3)
C Communication with FCN
      COMMON /FCNXSEC/ NTOT,NECM(20),IYEAR(300),IECMS(300),
     .    X(300),SS(300),
     .    SEF(300),SSUB(300),ACC(300),XNTCOR(300),ITYPE(20)
      COMMON /FCNENER/ ECMS(300),DECMS(300),EABS(20),EALFA(20),EREF(20),
     .  CENER(20,20),IENEF
      COMMON /FCNLUMI/ XLUMI(300),NBHA(300),SYLUM(20),SYCORR(20,20),
     &  SLPTP(20)
cft      COMMON /FCNLUMI/ XLUMI(300),NBHA(300),SYLUM(20),SYCORR(20,20)
      COMMON /CHI/ V1(300,300),STHER(300),CHISD(300)
      COMMON /OUTLVL/ IOUT
      COMMON /TAUPOL/ IVAT,VATM,SVAT,IVAE,VAEM,SVAE
C Information for combined fits
      COMMON /DECIDE/ NL,NA,NP,IFIT,IEE,UNIV,IASY
C Effective sinus for the S.M. modelled parts (very small)
      COMMON /INTER/ S2TW,MMIN1,MMINV(20),XMINV(20)
C Refit with a different parametrization
      COMMON /COVAR/ ISTAT
      COMMON /MIZA/IMIZA
C
      DIMENSION SP(16),SSP(16),SMIN(16),SMAX(16)
      DIMENSION JTYPE(20)
      DIMENSION SX(300),BSC(300)
      DIMENSION ERRVEC(50)
      DIMENSION R(3),AFB0(3)
      DIMENSION EMAT(300),CMAT(9,9)
      EXTERNAL FCSMT
C
      DATA PARNAM/'HADRONS   ',
     . 'ELECTRONS ','MUONS     ','TAUS      ','COM. LEPT.'/
      DATA CHNAM/
     . 'MZ-91     ','GZ        ','RT_h    ',
     . 'JT_h      ','RT_e      ','JT_e      ',
     . 'RT_m      ','JT_m      ','RT_t      ',
     . 'JT_t      ','RFB_e      ','JFB_e      ',
     . 'RFB_m      ','JFB_m      ','RFB_t      ',
     . 'JFB_t   '/
      DATA SMIN/0.D0,2.20D0,2.0D0,-5.D0,-0.5D0,-5.D0,-0.5D0,-5.D0,-0.
     &  5D0,-5.D0,-0.1D0,-0.5D0,-0.1D0,-0.5D0,-0.1D0,-0.5D0/
      DATA SMAX/0.5D0,2.80D0,4.0D0,5.D0,0.5D0,5.D0,0.5D0,5.D0,0.5D0,
     &  5.D0,0.1D0,1.5D0,0.1D0,1.5D0,0.1D0,1.5D0/
C***********************************************************************
C  Setting initial conditions for the fit:
C***********************************************************************
C
C=======================================================================
C  READ from control cards
C  # Parameters fitted, Output level
      READ(5,*) ITHEOR,ISTAT,IOUT,IEE
      READ(5,*) IMIZA
      READ(5,*) NP,UNIV
C  # Fit with F-B asymmetry also ?
      READ(5,*) IASY
C  READ initial values for the fit
      READ(5,*) MZ,GZ
      READ(5,*) (RT0(J1),JT0(J1),J1=1,4)
      READ(5,*) (RFB0(J1),JFB0(J1),J1=1,3)
CFT      READ(5,*) S2TW,MZ,GZ,S0,(R(J1),AFB0(J1),J1=1,3)
C=======================================================================
C READ from energy-luminosity data file
C=======================================================================
      OPEN(9,STATUS='OLD',READONLY)
C READ from energy-luminosity data file
      READ(9,*) NY,NPE,IENEF
C READ luminosity systematics information
      DO I = 1,NY
        READ(9,*) SLPTP(I)
      ENDDO
      DO I = 1,NY
        READ(9,*) SYLUM(I),(SYCORR(I,J),J=1,NY)
      ENDDO
C  READ LEP ENERGY systematics information
      IF(NY.GT.4)THEN
        NYM = 4
      ELSE
        NYM = NY
      ENDIF
C
      DO I = 1,NYM
        READ(9,*) EABS(I),EALFA(I),EREF(I)
      ENDDO
      IF(NY.GE.5)THEN
        DO I = 1,NPE
          READ(9,*) (CENER(I,J),J=1,NPE)
        ENDDO
      ENDIF
C
      READ(9,*) NE
      DO I = 1,NE
        READ(9,*) IYEAR(I),IDUMMY,ECMS(I),DECMS(I),XLUMI(I)
     .           ,NBHA(I),BSC(I)
      ENDDO
      CLOSE(9)
C=======================================================================
C READ from cross-section data file
C=======================================================================
      OPEN(7,STATUS='OLD',READONLY)
      READ(7,*) NL
      NTOT = 0
      I = 0
      DO J1 = 1,NL
C ITYPE = 0  --> hadrons
C       = 2  --> muons
C       = 3  --> taus
        READ(7,*) ITYPE(J1),NECM(J1),MMINV(J1),XMINV(J1)
C Choose the right name
        JTYPE(J1) = ITYPE(J1) + 1
        IF(ITYPE(J1).EQ.123) JTYPE(J1) = 5
C
        NTOT = NTOT+NECM(J1)
        DO K1 = 1,NECM(J1)
          I = I+1
C
        IF(ITHEOR.EQ.1)THEN
            READ(7,*) IECMS(I),X(I),SS(I)
            ISTAT = 0
        ELSE
            READ(7,*) IECMS(I),X(I),
     .                 SS(I),SEFR,SSUB(I),ACC(I),XNTCOR(I)
            SEF(I) = X(I)*SEFR
CFT Consider statistical errors of 4 LEP experiments...
cft            IF (IECMS(I).LT.35) THEN
cft              SS(I) = SS(I) / 2.D0
cft            ENDIF
        ENDIF
C
C BSC is the beam spread correction
        X(I) = X(I)/(1.D0 + BSC(IECMS(I)))
        ENDDO
      ENDDO
      CLOSE(7)
C=======================================================================
C READ from F-B asymmetry data file
C=======================================================================
      OPEN(8,STATUS='OLD',READONLY)
      NA = 0
      IF(IASY.EQ.1)THEN
C # of asymmetry courbes
        READ(8,*) NA
C Constraint on (gv/ga)_tau from P_tau
        READ(8,*) IVAT,VATM,SVAT
C Constraint on (gv/ga)_e from P_tau_fb
        READ(8,*) IVAE,VAEM,SVAE
        DO J1 = NL+1,NL+NA
          READ(8,*) ITYPE(J1),NECM(J1),MMINV(J1),XMINV(J1)
          JTYPE(J1) = ITYPE(J1) + 1
          IF(ITYPE(J1).EQ.10) JTYPE(J1) = 2
          IF(ITYPE(J1).EQ.123) JTYPE(J1) = 5
          NTOT = NTOT+NECM(J1)
          DO K1 = 1,NECM(J1)
            I = I+1
            READ(8,*) IECMS(I),X(I),SS(I),SEF(I)
          ENDDO
        ENDDO
      ENDIF
      CLOSE(8)
C
C=======================================================================
C  Writing the relevant input
C=======================================================================
C
      WRITE(6,318)MZ,GZ,RT0(1),JT0(1),RT0(2),JT0(2),RFB0(1),JFB0(1)
 318  FORMAT(/' ***** INPUT PARAMETERS:'/
     . ' INITIAL MZ          =',G12.5/
     . '         GZ          =',G12.5/
     . '         RT_h        =',G12.5/
     . '         JT_h        =',G12.5,'  '/
     . '         RT_l        =',G12.5,'  '/
     . '         JT_l        =',G12.5,'  '/
     . '         RFB_l       =',G12.5,'  '/
     . '         JFB_l       =',G12.5,'  '//)
C
CFT      NP = 4
      NDF = NTOT+IVAT+IVAE-NP
      WRITE(6,319)NTOT,IVAT+IVAE,NP,NDF
 319  FORMAT(/' ***** FIT INFORMATION:'/
     . ' Number of points      =',I8/
     . ' Number of constraints =',I8/
     . ' Number of fitting par.=',I8/
     . '           ===>   NDF  =',I8//)
C
C
      IF(IEE.NE.0)WRITE(6,150)(I,EABS(I),EALFA(I),EREF(I),I=1,NY)
      IF(IEE.NE.0) THEN
      IF(NY.GE.5)THEN
cft        DO I = 1,NPE
        DO I = 1,10
          K = IENEF + (I - 1)
          EFTEN = 1.D3*SQRT(CENER(I,I))
          WRITE(6,155)K,EFTEN,(1.D6*CENER(I,J),J=1,10)
        ENDDO
      ENDIF
      ENDIF
 150  FORMAT(/'  Global LEP energy uncertainties:'/
     .      14(I4,' Abs.scale=',G12.3,
     .           ' Non-linear.=',G12.3,' E_ref=',G12.5/))

      WRITE(6,160)
 155    FORMAT(I4,' Ener.(MeV)=',G12.3,
     .            ' Corr.(MeV**2)=',10F8.2/)
 160  FORMAT(/'  Global X-SECT uncertainties:'/)
      DO I=1,NY
        WRITE(6,161)I,SYLUM(I),(SYCORR(I,J),J=1,NY)
 161    FORMAT(I4,'   Lumi.=',G12.3,
     .            '   Corr.=',11F8.2/)
      ENDDO
      IF(IVAT.NE.0)WRITE(6,170)VATM,SVAT
 170  FORMAT(/'  Additionnal constraints: TAU POLARIZATION:'/
     .      '   gv(tau)/ga(tau) =',G12.3,'  +-  ',G12.3)
      IF(IVAE.NE.0)WRITE(6,171)VAEM,SVAE
 171  FORMAT(/'  Additionnal constraints: FORWARD-BACKWARD TAU POL.:'/
     .      '   gv(e)/ga(e)     =',G12.3,'  +-  ',G12.3)
C
C  Lineshape
      I = 0
      DO J1=1,NL
        IF(ITHEOR.EQ.1)THEN
            WRITE(6,320)PARNAM(JTYPE(J1)),NECM(J1),MMINV(J1),XMINV(J1)
320         FORMAT(/' LINESHAPE: ',A10,' #POINTS =',I4,
     .     '  MMIN =',G9.2,' GEV','  XMIN =',G9.2/
     .      '  N  ECM     +-      SIGMA        +-')
        ELSE
            WRITE(6,321)PARNAM(JTYPE(J1)),NECM(J1),MMINV(J1),XMINV(J1)
321         FORMAT(/' ',A10,' #POINTS =',I4,'  MMIN =',G12.5,'  GEV',
     .        '  XMIN =',G12.5/
     .   '  N  ECM     +-        LUMI     NBHA    SIGMA +-(STA)',
     .   ' +-(TOT) SIG_CORR MUL.SYS ADD.SYS    ACC.       SUB.')
        ENDIF
        DO K1=1,NECM(J1)
          I = I + 1
          SX(I) =  DSQRT( SS(I)**2 + SSUB(I)**2 + SEF(I)**2 +
     .          X(I)**2*SYLUM(IYEAR(I))**2  )
          IECM = IECMS(I)
          DECM = ECMS(IECM)*DECMS(IECM)
          IF(IEE.EQ.0)DECM = 0.D0
C
          IF(ITHEOR.EQ.1)THEN
            WRITE(6,330)IECM,ECMS(IECM),DECM,X(I),SX(I)
330         FORMAT(' ',I2,1X,F10.3,1X,F7.4,1X,F8.4,1X,F8.4)
          ELSE
            WRITE(6,331)IECM,ECMS(IECM),DECM,
     .        XLUMI(IECM),NBHA(IECM),
     .        X(I)*(1.+BSC(IECM)),SS(I),SX(I),X(I),SEF(I),SSUB(I),
     .        ACC(I),XNTCOR(I)
331         FORMAT(' ',I2,1X,F10.3,1X,F7.4,1X,F9.2,1X,I8,1X,F7.3,1X,
     .        F7.3,1X,F7.3,1X,F7.3,1X,
     .        F7.3,1X,F7.3,1X,F8.3,1X,F10.3)
          ENDIF
C
        ENDDO
      ENDDO
C
C  A_FB
      DO J1=NL+1,NL+NA
        WRITE(6,420)PARNAM(JTYPE(J1)),NECM(J1),MMINV(J1),XMINV(J1)
420     FORMAT(/' F-B ASYMMETRY: ',A10,' #POINTS =',I4,
     .     '  MMIN =',G9.2,' GEV','  XMIN =',G9.2/
     .     '  N  ECM     +-        AFB        +-     MUL.SYS')
        DO K1=1,NECM(J1)
          I = I + 1
          IECM = IECMS(I)
          DECM = ECMS(IECM)*DECMS(IECM)
          IF(IEE.EQ.0)DECM = 0.D0
          WRITE(6,430)IECM,ECMS(IECM),DECM,X(I),SS(I),SEF(I)
430       FORMAT(' ',I2,1X,F10.3,1X,F7.4,1X,F8.4,1X,F8.4,1X,F8.4)
        ENDDO
      ENDDO
C
C
C****************************************************************
C FITTING PROCEDURE
C****************************************************************
C
C INITIALIZE parameters
      SP(1) = MZ-91.d0
      SP(2) = GZ
      SP(3) = RT0(1)
      SP(4) = JT0(1)
      SP(5) = RT0(2)
      SP(6) = JT0(2)
      SP(7) = RT0(3)
      SP(8) = JT0(3)
      SP(9) = RT0(4)
      SP(10) = JT0(4)
      SP(11) = RFB0(1)
      SP(12) = JFB0(1)
      SP(13) = RFB0(2)
      SP(14) = JFB0(2)
      SP(15) = RFB0(3)
      SP(16) = JFB0(3)
C
C MINUIT minimization
      DO J2=0,1
CFT      DO J2=1,1
        CALL MNINIT(5,6,8)
        CALL MNEXCM(FCSMT,'SET PRINT',1.D0,1,IERR)
        IF(J2.EQ.0)THEN
C Fit # 0 uses data errors
          ISTAT= 0
        ELSE
          ISTAT= 1
C Use input parameters from prefit.
CFT          CALL CONV(0,SP)
          SP(1) = MZ-91.d0
          SP(2) = GZ
          SP(3) = RT(1)
          SP(4) = JT(1)
          SP(5) = RT(2)
          SP(6) = JT(2)
          SP(7) = RT(3)
          SP(8) = JT(3)
          SP(9) = RT(4)
          SP(10) = JT(4)
          SP(11) = RFB(1)
          SP(12) = JFB(1)
          SP(13) = RFB(2)
          SP(14) = JFB(2)
          SP(15) = RFB(3)
          SP(16) = JFB(3)
        ENDIF
        WRITE(6,444)J2,ISTAT
444     FORMAT(//' =============== FIT # ',I4,' ISTAT =',I2,
     .' ==============='//)
        CALL MNPARM(1,CHNAM(1),SP(1),1.D-5,SMIN(1),SMAX(1),IERR)
        CALL MNPARM(2,CHNAM(2),SP(2),1.D-5,SMIN(2),SMAX(2),IERR)
        CALL MNPARM(3,CHNAM(3),SP(3),1.D-5,SMIN(3),SMAX(3),IERR)
        IF ((NL.EQ.1).AND.(NP.EQ.4)) THEN
          CALL MNPARM(4,CHNAM(4),SP(4),1.D-5,SMIN(4),SMAX(4),IERR)
        ENDIF
        IF ((NL.GT.1).AND.(NP.EQ.4)) THEN
          CALL MNPARM(4,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
        ENDIF
        IF ((UNIV.EQ.1).AND.(NP.EQ.6).AND.(IMIZA.EQ.0)) THEN
          CALL MNPARM(4,CHNAM(4),SP(4),1.D-5,SMIN(4),SMAX(4),IERR)
          CALL MNPARM(5,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(6,CHNAM(6),SP(6),1.D-5,SMIN(6),SMAX(6),IERR)
        ENDIF
        IF ((UNIV.EQ.1).AND.(NP.EQ.6).AND.(IMIZA.EQ.2)) THEN
          CALL MNPARM(4,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(5,CHNAM(11),SP(11),1.D-5,SMIN(11),SMAX(11),IERR)
          CALL MNPARM(6,CHNAM(12),SP(12),1.D-5,SMIN(12),SMAX(12),IERR)
        ENDIF
        IF ((UNIV.EQ.0).AND.(NP.EQ.6)) THEN
          CALL MNPARM(4,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(5,CHNAM(7),SP(7),1.D-5,SMIN(7),SMAX(7),IERR)
          CALL MNPARM(6,CHNAM(9),SP(9),1.D-5,SMIN(9),SMAX(9),IERR)
        ENDIF
        IF ((UNIV.EQ.0).AND.(NP.EQ.10)) THEN
          CALL MNPARM(4,CHNAM(4),SP(4),1.D-5,SMIN(4),SMAX(4),IERR)
          CALL MNPARM(5,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(6,CHNAM(6),SP(6),1.D-5,SMIN(6),SMAX(6),IERR)
          CALL MNPARM(7,CHNAM(7),SP(7),1.D-5,SMIN(7),SMAX(7),IERR)
          CALL MNPARM(8,CHNAM(8),SP(8),1.D-5,SMIN(8),SMAX(8),IERR)
          CALL MNPARM(9,CHNAM(9),SP(9),1.D-5,SMIN(9),SMAX(9),IERR)
          CALL MNPARM(10,CHNAM(10),SP(10),1.D-5,SMIN(10),SMAX(10),IERR)
        ENDIF
        IF ((IASY.EQ.1).AND.(UNIV.EQ.1).AND.(NP.EQ.5)) THEN
          CALL MNPARM(4,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(5,CHNAM(11),SP(11),1.D-5,SMIN(11),SMAX(11),IERR)
        ENDIF
        IF ((IASY.EQ.1).AND.(UNIV.EQ.1).AND.(NP.EQ.8)) THEN
          CALL MNPARM(4,CHNAM(4),SP(4),1.D-5,SMIN(4),SMAX(4),IERR)
          CALL MNPARM(5,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(6,CHNAM(6),SP(6),1.D-5,SMIN(6),SMAX(6),IERR)
          CALL MNPARM(7,CHNAM(11),SP(11),1.D-5,SMIN(11),SMAX(11),IERR)
          CALL MNPARM(8,CHNAM(12),SP(12),1.D-5,SMIN(12),SMAX(12),IERR)
        ENDIF
        IF ((IASY.EQ.1).AND.(UNIV.EQ.0).AND.(NP.EQ.9)) THEN
          CALL MNPARM(4,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(5,CHNAM(7),SP(7),1.D-5,SMIN(7),SMAX(7),IERR)
          CALL MNPARM(6,CHNAM(9),SP(9),1.D-5,SMIN(9),SMAX(9),IERR)
          CALL MNPARM(7,CHNAM(11),SP(11),1.D-5,SMIN(11),SMAX(11),IERR)
          CALL MNPARM(8,CHNAM(13),SP(13),1.D-5,SMIN(13),SMAX(13),IERR)
          CALL MNPARM(9,CHNAM(15),SP(15),1.D-5,SMIN(15),SMAX(15),IERR)
        ENDIF
        IF ((IASY.EQ.1).AND.(UNIV.EQ.0).AND.(NP.EQ.16)) THEN
          CALL MNPARM(4,CHNAM(4),SP(4),1.D-5,SMIN(4),SMAX(4),IERR)
          CALL MNPARM(5,CHNAM(5),SP(5),1.D-5,SMIN(5),SMAX(5),IERR)
          CALL MNPARM(6,CHNAM(6),SP(6),1.D-5,SMIN(6),SMAX(6),IERR)
          CALL MNPARM(7,CHNAM(7),SP(7),1.D-5,SMIN(7),SMAX(7),IERR)
          CALL MNPARM(8,CHNAM(8),SP(8),1.D-5,SMIN(8),SMAX(8),IERR)
          CALL MNPARM(9,CHNAM(9),SP(9),1.D-5,SMIN(9),SMAX(9),IERR)
          CALL MNPARM(10,CHNAM(10),SP(10),1.D-5,SMIN(10),SMAX(10),IERR)
          CALL MNPARM(11,CHNAM(11),SP(11),1.D-5,SMIN(11),SMAX(11),IERR)
          CALL MNPARM(12,CHNAM(12),SP(12),1.D-5,SMIN(12),SMAX(12),IERR)
          CALL MNPARM(13,CHNAM(13),SP(13),1.D-5,SMIN(13),SMAX(13),IERR)
          CALL MNPARM(14,CHNAM(14),SP(14),1.D-5,SMIN(14),SMAX(14),IERR)
          CALL MNPARM(15,CHNAM(15),SP(15),1.D-5,SMIN(15),SMAX(15),IERR)
          CALL MNPARM(16,CHNAM(16),SP(16),1.D-5,SMIN(16),SMAX(16),IERR)
        ENDIF
        CALL MNEXCM(FCSMT,'MIGRAD',0.D0,0,IERR)
C       CALL MNEXCM(FCSMT,'MINOS',0.D0,0,IERR)
C-----------------------------------------
C Report the Diagonal CHIS after fit
C-----------------------------------------
        IF(J2.EQ.1)THEN
C===========================
          I = 0
C Lineshape
          CHISTL = 0.D0
          DO J1 = 1,NL
            WRITE(6,200) J1,NECM(J1)
 200        FORMAT(/'  LINESHAPE#',I5,' #POINTS =',I6/
     .   '  N  ECM     SIGMA   +-(TOT) STHEOR  CHIS    ST.DEV.')
            CHISSUM = 0.D0
            DO K1 = 1,NECM(J1)
              I = I + 1
              IECM = IECMS(I)
              CHISSUM = CHISSUM + CHISD(I)
              WRITE(6,210) IECM,ECMS(IECM),X(I),1.D0/DSQRT(V1(I,I)),
     .          STHER(I),
     .          CHISD(I),DSIGN( DSQRT(CHISD(I)) , (X(I)-STHER(I)) )
 210          FORMAT(' ',I2,1X,F10.3,1X,F7.3,1X,F7.3,
     .           1X,F7.3,1X,F7.3,1X,F7.3)
            ENDDO
          WRITE(6,220)CHISSUM
 220      FORMAT('  Total DIAGONAL CHIS for this lineshape =',G12.5)
          CHISTL = CHISTL + CHISSUM
          ENDDO
          IL = I
          WRITE(6,230)IL,CHISTL
 230      FORMAT(/'  Total DIAGONAL CHIS for ',I6,
     .                ' lineshape points =',G12.5)
C A_fb
          CHISTA = 0.D0
          DO J1 = NL+1,NL+NA
            WRITE(6,201) J1,NECM(J1)
 201        FORMAT(/'  A_FB #',I5,' #POINTS =',I6/
     .   '  N  ECM       AFB   +-(TOT) STHEOR  CHIS    ST.DEV.')
            CHISSUM = 0.D0
            DO K1 = 1,NECM(J1)
              I = I + 1
              IECM = IECMS(I)
              CHISSUM = CHISSUM + CHISD(I)
              WRITE(6,210) IECM,ECMS(IECM),X(I),1.D0/DSQRT(V1(I,I)),
     .          STHER(I),
     .          CHISD(I),DSIGN( DSQRT(CHISD(I)) , (X(I)-STHER(I)) )
 211          FORMAT(' ',I2,1X,F6.3,1X,F7.3,1X,F7.3,
     .           1X,F7.3,1X,F7.3,1X,F7.3)
            ENDDO
          WRITE(6,221)CHISSUM
 221      FORMAT('  Total DIAGONAL CHIS for this A_fb =',G12.5)
          CHISTA = CHISTA + CHISSUM
          ENDDO
          IA = I-IL
          WRITE(6,231)IA,CHISTA
 231      FORMAT(/'  Total DIAGONAL CHIS for ',I6,
     .                ' A_fb points =',G12.5)
          WRITE(6,232)I,CHISTL+CHISTA
 232      FORMAT(//'  Total DIAGONAL CHIS for ',I6,
     .                ' points =',G12.5)
C
C
C-----------------------------------------
C Write output file for refitting
C-----------------------------------------
          CALL MNSTAT(CHIS,FEDM,ERRDEF,NPARI,NPARX,ISTA)
          WRITE(10,*)NP,CHIS,NDF
C FITTED VALUES
CFT          IF(NP.NE.5)THEN
            DO J1=1,NP
              CALL MNPOUT(J1,CHNAM0,SP(J1),SSP(J1),BND1,BND2,IVARBL)
              WRITE(10,*)SP(J1),SSP(J1)
            ENDDO
CFT          ELSE
CFT            DO J1=1,4
CFT              CALL MNPOUT(J1,CHNAM0,SP(J1),SSP(J1),BND1,BND2,IVARBL)
CFT              WRITE(10,*)SP(J1),SSP(J1)
CFT            ENDDO
CFT            CALL MNPOUT(7,CHNAM0,SP(7),SSP(7),BND1,BND2,IVARBL)
CFT            WRITE(10,*)SP(7),SSP(7)
CFT          ENDIF
C CORRELATION MATRIX
          CALL MNEMAT(EMAT,NP)
          DO I1 = 1,NP
          DO I2 = 1,NP
            CMAT(I1,I2) = EMAT((I1-1)*NP+I2)/
     .         SQRT(EMAT((I1-1)*NP+I1)*EMAT((I2-1)*NP+I2))
          ENDDO
          ENDDO
          DO J1 = 1,NP
            WRITE(10,*) (CMAT(I1,J1),I1=1,NP)
          ENDDO
C===========================
        ENDIF
C
C
      ENDDO
      STOP
      END
C
C
C
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
      SUBROUTINE FCSMT(NPAR,GRAD,CHIS,SP,IFLAG)
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
      REAL*8 RT,JT,RFB,JFB,RT0,JT0,RFB0,JFB0
C
      COMMON /UNKNOWN/ MZ,GZ,RT(4),JT(4),RFB(3),JFB(3)
      COMMON /INITIAL/ RT0(4),RFB0(3),JT0(4),JFB0(3)
C
      COMMON /FCNXSEC/ NTOT,NECM(20),IYEAR(300),IECMS(300),
     .    X(300),SS(300),
     .    SEF(300),SSUB(300),ACC(300),XNTCOR(300),ITYPE(20)
      COMMON /FCNENER/ ECMS(300),DECMS(300),EABS(20),EALFA(20),EREF(20),
     .  CENER(20,20),IENEF
      COMMON /FCNLUMI/ XLUMI(300),NBHA(300),SYLUM(20),SYCORR(20,20),
     &  SLPTP(20)
cft      COMMON /FCNLUMI/ XLUMI(300),NBHA(300),SYLUM(20),SYCORR(20,20)
      COMMON /TAUPOL/ IVAT,VATM,SVAT,IVAE,VAEM,SVAE
C
      COMMON /INTER/ S2TW,MMIN1,MMINV(20),XMINV(20)
      COMMON /OUTLVL/ IOUT
      COMMON /DECIDE/ NL,NA,NP,IFIT,IEE,UNIV,IASY
      COMMON /FINFER/ IFLAV
C
      COMMON /COVAR/ ISTAT
      COMMON /MIZA/IMIZA
      COMMON /CHI/ V1(300,300),STHER(300),CHISD(300)
C
      DIMENSION AUX(300),DSDE(300)
      DIMENSION GRAD(*),SP(*)
      DATA ISTAT0/-1/
C Filling the input parameters
CFT      CALL CONV(1,SP)
      MZ = 91.D0 + SP(1)
      GZ = SP(2)
      RT(1) = SP(3)
      IF ((IASY.EQ.0).AND.(NL.EQ.1).AND.(NP.EQ.4)) THEN
        JT(1) = SP(4)
      ENDIF
      IF ((IASY.EQ.0).AND.(NL.GT.1).AND.(NP.EQ.4)) THEN
        JT(1) = JT0(1)
        JT(2) = JT0(2)
        JT(3) = JT0(3)
        JT(4) = JT0(4)
        RT(2) = SP(4)
        RT(3) = RT(2)
        RT(4) = RT(2)
      ENDIF
      IF ((IASY.EQ.0).AND.(UNIV.EQ.1).AND.(NP.EQ.6)) THEN
        JT(1) = SP(4)
        RT(2) = SP(5)
        JT(2) = SP(6)
        RT(3) = RT(2)
        RT(4) = RT(2)
        JT(3) = JT(2)
        JT(4) = JT(2)
      ENDIF
      IF ((IASY.EQ.0).AND.(UNIV.EQ.0).AND.(NP.EQ.6)) THEN
        JT(1) = JT0(1)
        JT(2) = JT0(2)
        JT(3) = JT0(3)
        JT(4) = JT0(4)
        RT(2) = SP(4)
        RT(3) = SP(5)
        RT(4) = SP(6)
      ENDIF
      IF ((IASY.EQ.0).AND.(UNIV.EQ.0).AND.(NP.EQ.10)) THEN
        JT(1) = SP(4)
        JT(2) = SP(6)
        JT(3) = SP(8)
        JT(4) = SP(10)
        RT(2) = SP(5)
        RT(3) = SP(7)
        RT(4) = SP(9)
      ENDIF
      IF ((IASY.EQ.1).AND.(UNIV.EQ.1).AND.(NP.EQ.5).AND.(IMIZA.EQ.0))
     &  THEN
        JT(1) = JT0(1)
        JT(2) = JT0(2)
        JT(3) = JT0(3)
        JT(4) = JT0(4)
        JFB(1) = JFB0(1)
        JFB(2) = JFB0(2)
        JFB(3) = JFB0(3)
        RT(2) = SP(4)
        RT(3) = RT(2)
        RT(4) = RT(2)
        RFB(1) = SP(5)
        RFB(2) = RFB(1)
        RFB(3) = RFB(1)
      ENDIF
C This is a special fit to compare with the MIZA output, where the leptonic
C interference is not fixed but related with the other parameters.
      IF ((IASY.EQ.1).AND.(UNIV.EQ.1).AND.(NP.EQ.5).AND.(IMIZA.EQ.1))
     &  THEN
        ALMZR = 1.064027D0
        RT(2) = SP(4)
        RT(3) = RT(2)
        RT(4) = RT(2)
        RFB(1) = SP(5)
        RFB(2) = RFB(1)
        RFB(3) = RFB(1)
        JT(1) = JT0(1)
        JT(2) = ALMZR * 2.D0/3.D0 * RFB(1) / DSQRT(RT(2))
        JT(3) = JT(2)
        JT(4) = JT(2)
        JFB(1) = ALMZR * 2.D0 * DSQRT(RT(2))
        JFB(2) = JFB(1)
        JFB(3) = JFB(1)
      ENDIF
      IF ((IASY.EQ.1).AND.(UNIV.EQ.1).AND.(NP.EQ.6).AND.(IMIZA.EQ.2))
     &  THEN
        RT(2) = SP(4)
        RT(3) = RT(2)
        RT(4) = RT(2)
        RFB(1) = SP(5)
        RFB(2) = RFB(1)
        RFB(3) = RFB(1)
        JT(1) = JT0(1)
        JT(2) = JT0(2)
        JT(3) = JT(2)
        JT(4) = JT(2)
        JFB(1) = SP(6)
        JFB(2) = JFB(1)
        JFB(3) = JFB(1)
      ENDIF
      IF ((IASY.EQ.1).AND.(UNIV.EQ.1).AND.(NP.EQ.8)) THEN
        JT(1) = SP(4)
        JT(2) = SP(6)
        JT(3) = JT(2)
        JT(4) = JT(2)
        JFB(1) = SP(8)
        JFB(2) = JFB(1)
        JFB(3) = JFB(1)
        RT(2) = SP(5)
        RT(3) = RT(2)
        RT(4) = RT(2)
        RFB(1) = SP(7)
        RFB(2) = RFB(1)
        RFB(3) = RFB(1)
      ENDIF
      IF ((IASY.EQ.1).AND.(UNIV.EQ.0).AND.(NP.EQ.9)) THEN
        JT(1) = JT0(1)
        JT(2) = JT0(2)
        JT(3) = JT0(3)
        JT(4) = JT0(4)
        JFB(1) = JFB0(1)
        JFB(2) = JFB0(2)
        JFB(3) = JFB0(3)
        RT(2) = SP(4)
        RT(3) = SP(5)
        RT(4) = SP(6)
        RFB(1) = SP(7)
        RFB(2) = SP(8)
        RFB(3) = SP(9)
      ENDIF
      IF ((IASY.EQ.1).AND.(UNIV.EQ.0).AND.(NP.EQ.16)) THEN
        JT(1) = SP(4)
        JT(2) = SP(6)
        JT(3) = SP(8)
        JT(4) = SP(10)
        JFB(1) = SP(12)
        JFB(2) = SP(14)
        JFB(3) = SP(16)
        RT(2) = SP(5)
        RT(3) = SP(7)
        RT(4) = SP(9)
        RFB(1) = SP(11)
        RFB(2) = SP(13)
        RFB(3) = SP(15)
      ENDIF
C-------------------------------
C Compute theoretical expectation
C-------------------------------
      I = 0
C-------------------
C Lineshape
C-------------------
      DO J1 = 1,NL
C Lineshape type
         IFLAV = ITYPE(J1)
C
         DO J2 = 1,NECM(J1)
           I = I + 1
           ECM = ECMS(IECMS(I))
C Minimum invariant mass cut: take the largest one
           MMIN2 = MMINV(J1)
           MMIN3 = XMINV(J1) * ECM/2.D0
           MMIN1 = MAX(MMIN2,MMIN3)
C Only for 1993 data the minimum invariant mass cut for electrons is the tau
C mass
           IF (ITYPE(J1).EQ.1) THEN
             IF (IECMS(I).GE.28.AND.IECMS(I).LE.31) MMIN1=1.777
           ENDIF
C Only for >140 GeV data the minimum invariant mass cut is variable.
           IF (IECMS(I).GE.38) THEN
             IF (ITYPE(J1).EQ.0) MMIN1=0.90D0*ECM/2.D0
             IF (ITYPE(J1).EQ.2) MMIN1=0.90D0*ECM/2.D0
             IF (ITYPE(J1).EQ.3) MMIN1=0.90D0*ECM/2.D0
           ENDIF
c
C
           IF(IFLAV.EQ.123)THEN
C Common leptons
cft             STHER(I) = 3.D0 * MIZA(ECM)
           ELSE
C Muon-like
cft             do km = 1,15
cft             mmin1 = 0.d0
cft             ecm = 20.000 + (km - 1)*10
CFT Consider TOPAZ data deconvoluted...
cft             IF (IECMS(I).EQ.38) THEN
cft               XEC2 = ECM**2
cft               STHER(I) = CSX(XEC2)
cft             ELSE
cft               STHER(I) = SMATR(ECM)
cft             ENDIF
cft             write (18,*) ecm,sther(i)
cft             enddo
               STHER(I) = SMATR(ECM)
           ENDIF
         ENDDO
      ENDDO
C-------------------
C A_FB
C-------------------
      DO J1 = NL+1,NL+NA
C Lineshape type
         IFLAV = ITYPE(J1)
C
         DO J2 = 1,NECM(J1)
           I = I + 1
           ECM = ECMS(IECMS(I))
C Minimum invariant mass cut: take the largest one
           MMIN2 = MMINV(J1)
           MMIN3 = XMINV(J1) * ECM/2.D0
           MMIN1 = MAX(MMIN2,MMIN3)
C Only for 1993 data the minimum invariant mass cut for electrons is the tau
C mass
           IF (ITYPE(J1).EQ.1) THEN
             IF (IECMS(I).GE.28.AND.IECMS(I).LE.31) MMIN1=1.777
           ENDIF
C Only for >140 GeV data the minimum invariant mass cut is variable.
           IF (IECMS(I).GE.38) THEN
             IF (ITYPE(J1).EQ.2) MMIN1=0.90D0*ECM/2.D0
             IF (ITYPE(J1).EQ.3) MMIN1=0.90D0*ECM/2.D0
           ENDIF

C
           STHER(I) = FBSMATR(ECM)
         ENDDO
      ENDDO
C----------------------------------------------------------------
C Covariance matrix calculation
C----------------------------------------------------------------
      IF(ISTAT.NE.ISTAT0)THEN
      ISTAT0 = ISTAT
C*********************************************
C
C==== Derivative for energy uncertainty propagation
        IF(IEE.NE.0)THEN
           DATA EPS/1.D-5/
           I = 0
C Lineshape
           DO J1 = 1,NL
           IFLAV = ITYPE(J1)
             DO J2 = 1,NECM(J1)
               I = I + 1
               ECM = ECMS(IECMS(I))
C Minimum invariant mass cut: take the largest one
               MMIN2 = MMINV(J1)
               MMIN3 = XMINV(J1) * ECM/2.D0
               MMIN1 = MAX(MMIN2,MMIN3)
               MIZAPE = smatr(ECM+EPS/2.D0)
               MIZAME = smatr(ECM-EPS/2.D0)
               DSDE(I) = ( MIZAPE - MIZAME )/EPS
               IF(IFLAV.EQ.123)DSDE(I) = 3.D0*DSDE(I)
             ENDDO
           ENDDO
C A_FB
           DO J1 = NL+1,NL+NA
           IFLAV = ITYPE(J1)
             DO J2 = 1,NECM(J1)
               I = I + 1
               ECM = ECMS(IECMS(I))
C Minimum invariant mass cut: take the largest one
               MMIN2 = MMINV(J1)
               MMIN3 = XMINV(J1) * ECM/2.D0
               MMIN1 = MAX(MMIN2,MMIN3)
               FBMPE = FBSMATR(ECM+EPS/2.D0)
               FBMME = FBSMATR(ECM-EPS/2.D0)
               DSDE(I) = ( FBMPE - FBMME )/EPS
             ENDDO
           ENDDO
        ENDIF
C==================================================
C
      DATA EPSE/0.020/
      I = 0
C
      DO J1=1,NL+NA
C
        DO K1=1,NECM(J1)
        I = I + 1
        J = 0
C
          DO J2=1,NL+NA
C
            DO K2=1,NECM(J2)
            J = J + 1
            IECM = IECMS(I)
            JECM = IECMS(J)
            IY = IYEAR(IECM)
            JY = IYEAR(JECM)
            V1(I,J) = 0.D0
C
C=======================================================
C    X-sect. uncertainty
C=======================================================
            IF(J1.LE.NL.AND.J2.LE.NL)THEN
C Decide how to compute the covariance matrix
              IF(ISTAT.EQ.0)THEN
C Use data
                 XIXJ = X(I)*X(J)
                 SXI = SS(I)**2
                 XNBI = 0.D0
                 XNBJ = 0.D0
              ELSE
C Use theory
                 XIXJ = STHER(I)*STHER(J)
                 SXI = STHER(I)/(ACC(I)*XLUMI(IECM))
                 XNBI = 1.D0/NBHA(IECM)
                 XNBJ = 1.D0/NBHA(JECM)
              ENDIF
C
C Compute diagonal terms
              IF(J.EQ.I) THEN
                 V1(I,J) = SXI +  SSUB(I)**2 + SEF(I)**2 +
     .             XIXJ * ( XNBI + SYLUM(IY)**2 ) +
     .             XNTCOR(I) * SXI**2 /XIXJ
              ELSE
C
C Compute off-diagonal terms
C    Luminosity systematic
C       Correlation SYCORR between different years
                   SYL2 = SYLUM(IY)*SYLUM(JY)*SYCORR(IY,JY)
                V1(I,J) = XIXJ*SYL2
C       Take care of the point to point error in the same year.
                IF (IY.EQ.JY) THEN
                  V1(I,J) = V1(I,J) - XIXJ*SLPTP(IY)*SLPTP(IY)
                ENDIF
C
CFT take 140 GeV run , as completely uncorrelated with the rest...
                IF ((IY.EQ.8).and.(JY.NE.8)) GOTO 884
                IF ((IY.NE.8).and.(JY.EQ.8)) GOTO 884
                IF(J1.EQ.J2) THEN
C    Global multiplicative systematic (for each specie)
                   V1(I,J) = V1(I,J) + SEF(I)*SEF(J)
C    Subtraction systematic (for each specie)
C    Specific modification for the correlation in the gamma-gamma
C      background subtraction in hadrons and taus that are uncorrelated
C      for different periods
                   IF(ITYPE(J1).NE.0.OR.SSUB(I).EQ.SSUB(J)) THEN
                      V1(I,J) = V1(I,J) + SSUB(I)*SSUB(J)
                   ENDIF
C                   IF(ITYPE(J1).NE.0.OR.ITYPE(J1).NE.3.OR.SSUB(I).EQ.
C     &               SSUB(J)) THEN
C                      V1(I,J) = V1(I,J) + SSUB(I)*SSUB(J)
C                   ENDIF
                ENDIF
C    Common error due to number of bhabha events for normalization
C    (equal for points at the same ECMS)
                IF(IECM.EQ.JECM)THEN
                   V1(I,J) = V1(I,J) + XIXJ*MIN(XNBI,XNBJ)
                ENDIF
              ENDIF
            ELSEIF(J1.GT.NL.AND.J2.GT.NL)THEN
C
C=======================================================
C    A_FB uncertainty
C=======================================================
CFT take 140 GeV run , as completely uncorrelated...
                IF ((IY.EQ.8).and.(JY.NE.8)) GOTO 884
                IF ((IY.NE.8).and.(JY.EQ.8)) GOTO 884
              IF(I.EQ.J)THEN
                V1(I,I) = SS(I)**2 + SEF(I)**2
              ELSE
                IF(J1.EQ.J2) THEN
C    Global multiplicative systematic (for each specie)
                   V1(I,J) = SEF(I)*SEF(J)
C    Subtraction systematic (for each specie)
                   V1(I,J) = V1(I,J) + SSUB(I)*SSUB(J)
                ENDIF
              ENDIF
            ENDIF
C
C=======================================================
C    Energy uncertainty (X-sect and A_FB)
C=======================================================
  884       CONTINUE
              IF(IEE.NE.0)THEN
              IF(IY.LT.5.AND.JY.LT.5)THEN
                ECMI = ECMS(IECM)
                ECMJ = ECMS(JECM)
                  IF(IECM.EQ.JECM)THEN
C    Points at the same energy
                    DECMI = DECMS(IECM)
                    EVX1 = ECMI**2 * ( DECMI**2 + EABS(IY)**2
     .               + EALFA(IY)**2 * ( (EREF(IY)-ECMI)/ECMI )**2 )
                    V1(I,J) = V1(I,J) + DSDE(I)*DSDE(J)*EVX1
                  ELSE
C    Points at different energy
C       Absolute energy scale
C          No correlation between energies in different years
                    IF(IY.EQ.JY)THEN
                      EVX1 = ECMI*ECMJ*EABS(IY)**2
                    ELSE
                      EVX1 = 0.D0
                    ENDIF
C          Exception: full correlation between the two absolute energy
C           uncertainties in 1992
                    IF((IY.EQ.3.AND.JY.EQ.4).OR.(IY.EQ.4.AND.JY.EQ.3))
     .                  THEN
                      EVX1 = ECMI*ECMJ*EABS(IY)*EABS(JY)
                    ENDIF
C       Non-linearity
C          Full correlation between energies in different years
                    EVX1 = EVX1 +
     .                ECMI * EALFA(IY) * ( (EREF(IY)-ECMI)/ECMI ) *
     .                ECMJ * EALFA(JY) * ( (EREF(JY)-ECMJ)/ECMJ )
                    V1(I,J) = V1(I,J) + DSDE(I)*DSDE(J)*EVX1
C       Point-to-point
C          Full correlation between energies in the same "point"
C          (only data in period 89-92)
                    IF(IY.LE.4.AND.JY.LE.4)THEN
                      IF(DABS(ECMI-ECMJ).LT.0.25D0)THEN
                        EPTP = 3.D-5
                        V1(I,J) = V1(I,J) +
     .                      DSDE(I)*DSDE(J)*ECMI*ECMJ*EPTP**2
                      ENDIF
                    ENDIF
                  ENDIF
                ELSEIF(IY.GE.5.AND.JY.GE.5)THEN
                  V1(I,J) = V1(I,J) +
     .                DSDE(I)*DSDE(J)*CENER(IECM-IENEF,JECM-IENEF)
                ENDIF
              ENDIF
C=======================================================
C
            ENDDO
          ENDDO
        ENDDO
      ENDDO
C-----------------------------
C Invert the covariance matrix
C-----------------------------
C
C  Write down (if required)
      IF(IOUT.GE.2)THEN
        WRITE(6,886)
886     FORMAT(//'   COVARIANCE MATRIX:')
        DO J=1,NTOT
          WRITE(6,*) J
          WRITE(6,887)( V1(I,J),I=1,NTOT )
887       FORMAT(8(1X,G9.3))
        ENDDO
      ENDIF
C
      CALL DINV(NTOT,V1,300,AUX,IFAIL)
      IF(IFAIL.NE.0) WRITE(6,*)IFAIL
C*********************************************
      ENDIF
C
C-----------------------------
C Chi squared calculation
C-----------------------------
      CHIS = 0.D0
C Diagonal terms
      DO I = 1,NTOT
         CHISD(I) = (X(I)-STHER(I))**2 * V1(I,I)
         CHIS = CHIS + CHISD(I)
      ENDDO
C Off-diagonal terms
      DO I = 1,NTOT-1
         DO J = I+1,NTOT
             CHIS = CHIS +
     .    2.D0 * (X(I)-STHER(I)) * V1(I,J) * (X(J)-STHER(J))
         ENDDO
      ENDDO
C-----------------------------
C Additional constraints
C-----------------------------
C  (gv/ga)_tau from tau polarization
CFT      IF(IVAT.EQ.1) CHIS = CHIS + ( (VATM - GV(3)/GA(3))/SVAT )**2
C  (gv/ga)_e from F-B tau polarization
CFT      IF(IVAE.EQ.1) CHIS = CHIS + ( (VAEM - GV(1)/GA(1))/SVAE )**2
C
      RETURN
      END
C
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C             SMATR block
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C
      DOUBLE PRECISION FUNCTION SMATR(ECM)
*******************************************************************
*                                                                 *
*                                     Manel Martinez              *
*                                     Frederic Teubert            *
*                                     Feb. 1995                   *
*                                                                 *
*******************************************************************
*     exponentiation as in Nucl. Phys. B297(1988)429.             *
*     The differences between the various Born approximations are *
*     explained in Phys.Lett. B203(1988)177 by F.A. Berends et al.*
*******************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,K-Z)
      INTEGER*4 NL,NA,NP
      REAL*8 RT,JT,RFB,JFB,RT0,JT0,RFB0,JFB0
      EXTERNAL DSIGVS,DSIGH
      COMMON /UNKNOWN/ MZ,GZ,RT(4),JT(4),RFB(3),JFB(3)
      COMMON /INITIAL/ RT0(4),RFB0(3),JT0(4),JFB0(3)
      COMMON /FINFER/ IFLAV
      COMMON /DECIDE/ NL,NA,NP,IFIT,IEE,UNIV,IASY
      COMMON /SIG/ SIG0
      COMMON /SIGMAS/ SGA
      COMMON /BEAM/S, LNS, BB, AP
      COMMON /INTER/ S2TW,MMIN1,MMINV(20),XMINV(20)
      COMMON /WHICH/ISEL
      COMMON /FINAL/ DQED,CQED,CQCD
      PARAMETER (Z2=1.6449340668482264D0,Z3=1.2020569031595943D0,
     .           PI=3.1415926535897932D0,
     .           ALFINV=137.036D0,
     .           GEVNB=0.38937966D6,
     .           M=0.0005110034D0,
     .           MMU=0.10565839D0,
     .           MTAU=1.7841D0,
     .           MPI=0.1349734D0)
C
C
C Initialize constant parameters
      DATA INIT/0/
      IF(INIT.EQ.0) THEN
        INIT = 1
        AP    = 1.D0/(ALFINV*PI)
        LN2 = LOG(2.D0)
C
C
C  Parameters of CERN INTEGRATION SUBROUTINE:
C EPSREL is the requested accuracy.
C A and B are the lower and upper integration boundary.
        EPSREL=1.D-06
        A=0.D0
        B=1.0D0
C QCD and QED correction factors
C QCD
      ALFAS = 0.D0
      IF (ALFAS .EQ. 0.D0)THEN
        DQCD = 0.045D0
      ELSEIF(ALFAS.GT.0)THEN
        DQCD = ALFAS/PI + (ALFAS/PI)**2*( 1.98D0 - 5D0*.115D0 )
      END IF
        CQCD = 1.D0 + DQCD
C QED
      DQED = 3.D0/(4.D0*PI*ALFINV)
        CQED = 1.D0 + DQED
C
C Compute relevant factors for the interference size
        T3 = 1.D0/2
        QE = -1.D0
        QU = 2.D0/3
        QD = -1.D0/3
C Charged leptons
        SGAL = 4.D0*PI / 3.D0 * CQED
cft        CALL COUPLIGS(-T3,QE,S2TW,GLP,GLM)
cft        SITL = GEVNB * 4.D0*PI / 6.D0 * CQED
cft     .     * ( GLP + GLM )**2
C Hadrons
        SGAH = SGAL/CQED * CQCD *(11./3 +  35./27. * DQED)
C    Up-type quarks couplings
cft        CALL COUPLIGS(T3,QU,S2TW,GUP,GUM)
C    Down-type quarks couplings
cft        CALL COUPLIGS(-T3,QD,S2TW,GDP,GDM)
cft        SITH = GEVNB * 4.D0*PI / 6.D0 * CQCD
cft     .     * ( GLP + GLM )
cft     . * 3.D0* ( 2.D0*QE*QU*(1.+ 4./9.*DQED)*( GUP + GUM )
cft     .         + 3.D0*QE*QD*(1.+ 1./9.*DQED)*( GDP + GDM ) )
      ENDIF
C Decide among different processes
      IF(IFLAV.EQ.0)THEN
C Hadrons
        MMIN0 = MPI
        SGA = SGAH
      ELSEIF(IFLAV.EQ.1.OR.IFLAV.EQ.123)THEN
C Electrons
        MMIN0 = MMU
        SGA = SGAL
      ELSEIF(IFLAV.EQ.2)THEN
C Muons
        MMIN0 = MMU
        SGA = SGAL
      ELSEIF(IFLAV.EQ.3)THEN
C Taus
        MMIN0 = MTAU
CMMR        MMIN0 = MMU
        SGA = SGAL
      ENDIF

C
       EB = ECM/2.
       S = 4.*EB**2

C
C Asking for lower integration bound of s'
       IF(MMIN1.EQ.0.D0)THEN
C Take the internal one
        MMIN = MMIN0
       ELSE
C Take the external one
        MMIN = MMIN1
       ENDIF
C
       SMIN  = (2.*MMIN)**2
       B   = 1.D0-SMIN/S
       LNS = 2.D0*LOG(ECM/M)
       BB = LNS -1.D0
       BETA  = 2.*AP*BB
       BBETA =  EXP(BETA*LOG(B))
C
       SIG0 = CSX(S)
       DVS  = AP*(1.5D0*LNS+2.D0*Z2-2.D0)+
     .        AP**2*( (9.D0/8.D0-2.D0*Z2)*LNS**2+
     .        (-45.D0/16.D0+11.D0*Z2/2.D0+3.D0*Z3)*LNS
     .        -6.D0*Z2**2/5.D0-9.D0*Z3/2.D0-6.D0*Z2*LN2+
     .          3.D0*Z2/8.D0+ 513.D0/108.D0)
C  Integration of (regularized) differential cross sections:
C  CALLS TO  *** C E R N ***SUBROUTINE:
       SIGVS = DGAUSS(DSIGVS,A,B,EPSREL)
C
       SIGH  = DGAUSS(DSIGH ,A,B,EPSREL)
C
       SMATR = (BBETA*SIG0+SIGVS)*(1.D0+DVS)+SIGH
C
C Initial-final state interference correction
       XMIN = DSQRT(SMIN/S)
       S2TW = 0.231
       IF(IFLAV.EQ.0)THEN
            DIF = DIFH(MZ,GZ,S2TW,XMIN,ECM)
       ELSE
            DIF = DIFL(MZ,GZ,S2TW,XMIN,ECM)
       ENDIF
       SMATR = SMATR * (1.D0 + DIF )
C
      RETURN
      END
C
C
C
C-----------------------------------------------------------------------
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      FUNCTION DIFL(XMZ,XGZ,S2TW,XMIN,ECM)
C=======================================================================
C
C    INITIAL/FINAL STATE INTERFERENCE FOR LEPTONS
C
C                      M.MARTINEZ AND R.MIQUEL   CERN/1990
C
C=======================================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16 Z
      PARAMETER (ALPHA=1.D0/137.036D0, PI=3.141592654D0,
     .            QE=-1.D0, T3E=-0.5D0)
C
      XMAX = 1.D0 - XMIN
      IPART = 1
C----------------
      S = ECM**2
      Z = DCMPLX(XMZ**2/S,-XGZ/XMZ)
      XI2ZGZ = DREAL( Z*(Z+1.)*CDLOG((XMAX+Z-1.)/Z) + (Z-1.)*(1.-XMAX) )
     .    - DLOG(CDABS(Z)) - 2.*DLOG(XMAX)
      XI2ZGG = DREAL( (1.-DCONJG(Z))*
     .             ( Z*(Z+1.)*CDLOG((XMAX+Z-1.)/Z) +
     .              (Z-1.)*(1.-XMAX)-CDLOG(Z)-2.*DLOG(XMAX) ) )
      XI2GGZ = DREAL(1.-Z)*(-2.*DLOG(XMAX)+XMAX+0.5)
C----------------
C  COUPLINGS
C----------------
      C2TW = 1.D0 - S2TW
      BEL = (T3E-S2TW*QE)/(DSQRT(S2TW)*DSQRT(C2TW))
      BER = (   -S2TW*QE)/(DSQRT(S2TW)*DSQRT(C2TW))
      BEL2=BEL*BEL
      BER2=BER*BER
C
      IF     (IPART.EQ.1) THEN
        T3F = -.5D0
        QF  = -1.D0
      ELSEIF (IPART.EQ.2) THEN
        T3F = .5D0
        QF  = 2.D0/3.D0
      ELSEIF (IPART.EQ.3) THEN
        T3F = -.5D0
        QF  = -1.D0/3.D0
      ENDIF
C
      BFL = (T3F-S2TW*QF)/(DSQRT(S2TW)*DSQRT(C2TW))
      BFR = (   -S2TW*QF)/(DSQRT(S2TW)*DSQRT(C2TW))
      BFL2=BFL*BFL
      BFR2=BFR*BFR
C----------------
      DELZGZ = 3.*ALPHA/PI*QE*QF*(BEL2-BER2)/(BEL2+BER2)
     .             *(BFL2-BFR2)/(BFL2+BFR2)*XI2ZGZ
C
      DELZGG = 3.*ALPHA/PI*QE**2*QF**2*(BEL-BER)/(BEL2+BER2)
     .             *(BFL-BFR)/(BFL2+BFR2)*XI2ZGG
C
      DELGGZ = 3.*ALPHA/PI*QE**2*QF**2*(BEL-BER)/(BEL2+BER2)
     .             *(BFL-BFR)/(BFL2+BFR2)*XI2GGZ
C----------------
      DIFL = DELZGZ + DELZGG + DELGGZ
      END
C
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C
      FUNCTION DIFH(XMZ,XGZ,S2TW,XMIN,ECM)
C=======================================================================
C
C    INITIAL/FINAL STATE INTERFERENCE FOR HADRONS
C
C                      M.MARTINEZ AND R.MIQUEL   CERN/1990
C
C=======================================================================
      IMPLICIT REAL*8 (A-H,O-Z)
      COMPLEX*16 Z
      PARAMETER (ALPHA=1.D0/137.036D0, PI=3.141592654D0,
     .   QE=-1.D0, T3E=-0.5D0)
      PARAMETER (GZU=0.297, GZD=0.383, GZB=0.379)
C
      XMAX = 1.D0 - XMIN
C----------------
      S = ECM**2
      Z = DCMPLX(XMZ**2/S,-XGZ/XMZ)
      XI2ZGZ = DREAL( Z*(Z+1.)*CDLOG((XMAX+Z-1.)/Z) + (Z-1.)*(1.-XMAX) )
     .    - DLOG(CDABS(Z)) - 2.*DLOG(XMAX)
      XI2ZGG = DREAL( (1.-DCONJG(Z))*
     .             ( Z*(Z+1.)*CDLOG((XMAX+Z-1.)/Z) +
     .              (Z-1.)*(1.-XMAX)-CDLOG(Z)-2.*DLOG(XMAX) ) )
      XI2GGZ = DREAL(1.-Z)*(-2.*DLOG(XMAX)+XMAX+0.5)
C----------------
C COUPLINGS
C----------------
      C2TW = 1 - S2TW
      BEL = (T3E-S2TW*QE)/(DSQRT(S2TW)*DSQRT(C2TW))
      BER = (   -S2TW*QE)/(DSQRT(S2TW)*DSQRT(C2TW))
      BEL2=BEL*BEL
      BER2=BER*BER
C------------
C UP QUARKS
C------------
      T3F = .5
      QF = 2./3.
      BFL = (T3F-S2TW*QF)/(DSQRT(S2TW)*DSQRT(C2TW))
      BFR = (   -S2TW*QF)/(DSQRT(S2TW)*DSQRT(C2TW))
      BFL2=BFL*BFL
      BFR2=BFR*BFR
C
      DUZGZ = 3.*ALPHA/PI*QE*QF*(BEL2-BER2)/(BEL2+BER2)
     .                          *(BFL2-BFR2)/(BFL2+BFR2)*XI2ZGZ
      DUZGG = 3.*ALPHA/PI*QE**2*QF**2*(BEL-BER)/(BEL2+BER2)
     .                                *(BFL-BFR)/(BFL2+BFR2)*XI2ZGG
      DUGGZ = 3.*ALPHA/PI*QE**2*QF**2*(BEL-BER)/(BEL2+BER2)
     .                                *(BFL-BFR)/(BFL2+BFR2)*XI2GGZ
C------------
C DOWN QUARKS
C------------
      T3F = -.5
      QF = -1./3.
      BFL = (T3F-S2TW*QF)/(DSQRT(S2TW)*DSQRT(C2TW))
      BFR = (   -S2TW*QF)/(DSQRT(S2TW)*DSQRT(C2TW))
      BFL2=BFL*BFL
      BFR2=BFR*BFR
      DDZGZ = 3.*ALPHA/PI*QE*QF*(BEL2-BER2)/(BEL2+BER2)
     .                          *(BFL2-BFR2)/(BFL2+BFR2)*XI2ZGZ
      DDZGG = 3.*ALPHA/PI*QE**2*QF**2*(BEL-BER)/(BEL2+BER2)
     .                                *(BFL-BFR)/(BFL2+BFR2)*XI2ZGG
      DDGGZ = 3.*ALPHA/PI*QE**2*QF**2*(BEL-BER)/(BEL2+BER2)
     .                                *(BFL-BFR)/(BFL2+BFR2)*XI2GGZ
C------------
      DLZGZH = (2.*GZU*DUZGZ + (2.*GZD+GZB)*DDZGZ)/
     .                    (2.*GZU + 2.*GZD+GZB)
      DLZGGH = (2.*GZU*DUZGG + (2.*GZD+GZB)*DDZGG)/
     .                    (2.*GZU + 2.*GZD+GZB)
      DLGGZH = (2.*GZU*DUGGZ + (2.*GZD+GZB)*DDGGZ)/
     .                    (2.*GZU + 2.*GZD+GZB)
C
      DIFH = DLZGZH + DLZGGH + DLGGZH
      END
C
C
C

      SUBROUTINE COUPLIGS(T3,Q,S2TW,GP,GM)
*********************************************************************
* couplings as a function of the effective sinus
*********************************************************************
      IMPLICIT DOUBLE PRECISION (A-Z)
C + - helicity couplings squared
      DEN = DSQRT(S2TW*(1.-S2TW))
      GP = (T3-Q*S2TW)/DEN
      GM = (-Q*S2TW)/DEN
      RETURN
      END
C
C
C
      DOUBLE PRECISION FUNCTION CSX(SS)
      IMPLICIT DOUBLE PRECISION (A-H,M,O-Z)
      REAL*8 RT,JT,RFB,JFB,RT0,JT0,RFB0,JFB0
      REAL*8 RTTAU,JTTAU,JTTOT
      COMMON /SIGMAS/ SGA
      COMMON /FINFER/ IFLAV
      COMMON /UNKNOWN/ MZ,GZ,RT(4),JT(4),RFB(3),JFB(3)
      COMMON /INITIAL/ RT0(4),RFB0(3),JT0(4),JFB0(3)
      COMMON /FINAL/ DQED,CQED,CQCD
      PARAMETER (Z2=1.6449340668482264D0,Z3=1.2020569031595943D0,
     .           PI=3.1415926535897932D0,
     .           ALFINV=137.036D0,
     .           GEVNB=0.38937966D6,
     .           M=0.0005110034D0,
     .           MMU=0.10565839D0,
     .           MTAU=1.7841D0,
     .           MPI=0.1349734D0)
      PARAMETER (RH0 = 0.D0, RH1 = 0.D0)
C
C
      ALFA1 = 1.D0/ALF(SS)
CFT      ALFMZ = 1.D0/ALF(MZ**2)
      ALFA2 = (ALFA1)**2
C Hadrons
      IF (IFLAV.EQ.0) THEN
CFT        JTTOT = JT(1) * (ALFA1/ALFMZ)
        CSX = ALFA2 * SGA / SS +
     .        4.D0*PI/(3.D0*ALFINV**2)* (
     .        (SS*RT(1)+(SS - MZ**2)*JT(1))/((SS-MZ**2)**2 + (GZ*MZ)**2)
     .        + rh0 + rh1*(1 - MZ**2/SS) ) 
      ENDIF
C Electrons
      IF (IFLAV.EQ.1) THEN
        CSX = ALFA2 * SGA / SS +
     .        4.D0*PI/(3.D0*ALFINV**2)*(SS*RT(2)+(SS - MZ**2)*JT(2))
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2)
ccc                                         * CQED
      ENDIF
C Muons
      IF (IFLAV.EQ.2) THEN
        CSX = ALFA2 * SGA / SS +
     .        4.D0*PI/(3.D0*ALFINV**2)*(SS*RT(3)+(SS - MZ**2)*JT(3))
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2) 
ccc                                         * CQED
      ENDIF
C Taus
      IF (IFLAV.EQ.3) THEN
CFT take into account mass effect for taus.
        RTTAU = RT(4)*0.998D0
        JTTAU = JT(4)*0.998D0
        CSX = ALFA2 * SGA / SS +
     .        4.D0*PI/(3.D0*ALFINV**2)*(SS*RTTAU+(SS - MZ**2)*JTTAU)
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2) 
ccc                                         * CQED
      ENDIF
      CSX = GEVNB * CSX
      RETURN
      END
C
C
C
      DOUBLE PRECISION FUNCTION ALF(S)
*********************************************************************
* Calculation of running ALPHA QED
*********************************************************************
      IMPLICIT REAL*8(A-Z)
C
      PARAMETER (S0=92.D0**2,A=0.0280D0,B=0.002980D0,C=0.006307D0)
      PARAMETER (PI=3.1415926536D0)
      PARAMETER (ALINV0=137.0359895D0,AL=1.D0/PI/ALINV0)
      PARAMETER (ME=0.000511D0,MMU=0.10566D0,MTAU=1.777D0)
      DATA INIT/0/
      IF(INIT.EQ.0)THEN
        INIT = 1
        LS0 = DLOG(S0)
        PIG0 = AL/3.D0*(-2.D0*(DLOG(ME)+DLOG(MMU)+DLOG(MTAU))-5.D0)+A
      ENDIF
C
      DS = DABS(S)
      LQ2 = DLOG(DS)
      PIG = PIG0 + AL*LQ2 + B*((LQ2-LS0) + C*(DS/S0-1.D0))
      ALF = (1.D0-PIG)*ALINV0
      END
C
C
C
      DOUBLE PRECISION FUNCTION DSIGVS(X)
*********************************************************************
* Calculation of integral over the soft-photon csx dsig(s')/ds':
* beta*(1-s'/s)**(beta-1).  The almost-pole at s'=s has been split off.
* This calculation sums over all orders of alfa ("exponentiation")
*********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,K-Z)
      COMMON/BEAM/S, LNS, BB, AP
      COMMON/SIG/ SIG0
        BETA = 2.D0*AP*BB
        FBETA  = BETA*DEXP(BETA*LOG(X))/X
        SX  = S*(1.D0-X)
        DSIGVS = FBETA*(CSX(SX)-SIG0)
      END
C
C
C
      DOUBLE PRECISION FUNCTION DSIGH(X)
*********************************************************************
* Calculation of the non-IR part of dsig/ds'.                       *
* 1st+2nd order calculation.                                        *
*********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,K-Z)
      COMMON/BEAM/S, LNS, BB, AP
      COMMON/SIG/ SIG0
      PARAMETER (Z2=1.6449340668482264D0,Z3=1.2020569031595943D0)
       Z   = 1.D0-X
       SX  = S*Z
       LNZ = LOG(Z)
       LN1Z= LOG(X)
       LNZ2= LNZ**2
       LNV2= LN1Z**2
       LNVZ= LN1Z*LNZ
       IF(X.LE. 0.5D0) THEN
         U =-LNZ
       ELSE
         U =-LN1Z
       ENDIF
       LI2 = ((((((-1.D0/10886400.D0)*U**2+1.D0/211680.D0)*U**2-
     .         1.D0/3600)*U**2+1.D0/36.D0)*U-1.D0/4.D0)*U+1.D0)*U
       S12 = ((((((-1.D0/24192000.D0)*U**2+1.D0/483840D0)*U**2
     .  -1.D0/8640.D0)*U**2+1.D0/96.D0)*U-1.D0/12.D0)*U+1.D0/4.D0)*U**2
       LI3 = ((((((((((11.D0/435456000D0)*U+11.D0/127008000D0)*U-
     .         1.D0/752640.D0)*U-19.D0/5556600.D0)*U+
     .         7.D0/86400.D0)*U+7.D0/54000.D0)*U-5.D0/576.D0)*U+
     .         17.D0/216.D0)*U-3.D0/8.D0)*U+1.D0)*U
       IF(X.GT. 0.5D0) THEN
          SSS= S12
          S12=-LI3+LNZ*LI2+LN1Z*LNZ**2/2.D0+Z3
          LI3=-SSS-LN1Z*LI2-LNZ*LN1Z**2/2.D0+LN1Z*Z2+Z3
          LI2=-LI2+Z2-LNVZ
       ENDIF
       RZ  = 2./X-2.+X
      C2= RZ*(-LNZ)+(1.D0+Z)*(0.5*LNZ-2.*LN1Z)-5./2.-0.5*Z
      C1= RZ*(+LNZ*LN1Z+3.5D0*LNZ+LI2-LNZ**2/2.D0)+
     .       3.D0*(1.D0-Z)+0.5D0-LNZ+
     .     (1.D0+Z)*(0.25*LNZ**2+4.D0*LN1Z-2.D0*Z2+3.5D0)
      DX  = (-1.D0-2.D0*LNZ/X-LNZ2/X**2)/3.D0
      C0= RZ*(LNZ2*LN1Z/2.D0+LI2*LNZ/2.D0-LNZ**3/6.D0+Z2*LNZ-
     .            1.5D0*LI2-1.5D0*LNVZ-17.D0*LNZ/6.D0+LNZ2/6.D0
     .           +LNZ2/X+DX)
     .   +(1.D0+Z)*(1.5D0*LI3-2.D0*S12-LI2*LN1Z-Z2+LNZ2/6.D0-1.D0/6.D0+
     .            LNVZ/2.D0-LNV2/4.D0+5.D0*LNZ/2.D0+3.D0*LN1Z/2.D0+
     .            2.D0*Z2-2.D0-2.D0*LN1Z)
     .   +   Z  *(16.D0*Z2/3.D0-25.D0*LI2/6.D0-13.D0*LNZ2/12.D0-
     .    5.D0*LNZ/6.D0-4.D0*LNVZ+1.5D0*LNV2-15.D0*LN1Z/6.D0-2.D0/3.D0)
      CSXSX = CSX(SX)
        DSIGH1 = -AP*BB*(1.D0+Z)*CSXSX
        DSIGH2 =  AP**2*(C2*LNS**2+C1*LNS+C0)*CSXSX
      DSIGH  = DSIGH1+DSIGH2
      END
C
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C             FBSMATR block
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C
      DOUBLE PRECISION FUNCTION FBSMATR(ECM)
*******************************************************************
*                 Authors:   Manel MARTINEZ                       *
*                          Frederic TEUBERT                       *
*                                  Feb 1995                       *
*******************************************************************
*     exponentiation as in Nucl. Phys. B297(1988)429.             *
*     The differences between the various Born approximations are *
*     explained in Phys.Lett. B203(1988)177 by F.A. Berends et al.*
*******************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,L,M,O-Z)
      EXTERNAL DSIGVSA,DSIGHA
      COMMON /UNKNOWN/ MZ,GZ,RT(4),JT(4),RFB(3),JFB(3)
      COMMON /INITIAL/ RT0(4),RFB0(3),JT0(4),JFB0(3)
      COMMON /RAD/ IRAD,SMIN
      COMMON /FINFER/ IFLAV
      COMMON /INTER/ S2TW,MMIN1,MMINV(20),XMINV(20)
      COMMON /SIG/ SIG0
      COMMON /BEAM/S, LNS, BB, AP
      COMMON /WHICH/ ISEL
      COMMON /BORN/ FBSMATRB
      PARAMETER (Z2=1.6449340668482264D0,Z3=1.2020569031595943D0,
     .           PI=3.1415926535897932D0,
     .           ALFINV=137.036D0,
     .           M=0.0005110034D0,
     .           MMU=0.10565839D0,
     .           MTAU=1.7841D0,
     .           MPI=0.1349734D0)
C
C
       DATA INIT/0/,IFLAV0/0/
       IF(INIT.EQ.0) THEN
         INIT = 1
         AP    = 1.D0/(ALFINV*PI)
         LN2 = DLOG(2.D0)
C
C  Parameters of CERN INTEGRATION SUBROUTINE:
C EPSREL is the requested accuracy.
C A and B are the lower and upper integration boundary.
         EPSREL=1.D-04
         A=0.D0
         B=1.0D0
C QED correction factors
C QED
        DQED = 3.D0/(4.D0*PI*ALFINV)
        CQED = 1.D0 + DQED
C Charged leptons
        SGAL = 4.D0*PI / 3.D0 * CQED
       ENDIF
C Decide among different processes
      IF(IFLAV.EQ.1)THEN
C Electrons
        MMIN0 = MMU
        SGA = SGAL
      ELSEIF(IFLAV.EQ.2)THEN
C Muons
        MMIN0 = MMU
        SGA = SGAL
      ELSEIF(IFLAV.EQ.3)THEN
C Taus
        MMIN0 = MTAU
CMMR        MMIN0 = MMU
        SGA = SGAL
      ENDIF
C
C Asking for lower integration bound of s'
       IF(MMIN1.EQ.0.D0)THEN
C Take the internal one
        MMIN = MMIN0
       ELSE
C Take the external one
        MMIN = MMIN1
       ENDIF
C
C
        SMIN  = (2.*MMIN)**2
C
       EB = ECM/2.D0
       S = 4.D0*EB**2
       IRAD = 1
       IF(IRAD.GE.1)THEN
         B = 1.D0-SMIN/S
         LNS = 2.D0*DLOG(ECM/M)
         BB = LNS - 1.D0
         BETA  = 2.D0*AP*BB
         BBETA =  EXP(BETA*DLOG(B))
*
         DVS  = AP*(1.5D0*LNS+2.D0*Z2-2.D0)+
     .        AP**2*( (9.D0/8.D0-2.D0*Z2)*LNS**2+
     .        (-45.D0/16.D0+11.D0*Z2/2.D0+3.D0*Z3)*LNS
     .        -6.D0*Z2**2/5.D0-9.D0*Z3/2.D0-6.D0*Z2*LN2-
     .          3.D0*Z2/8.D0+ 513.D0/108.D0)
       ENDIF
C Compute numerator
       ISEL = 1
       SIG0FB = CSXA(S)
       SIG0 = SIG0FB
       IF(IRAD.GE.1)THEN
         SIGVS = DGAUSS(DSIGVSA,A,B,EPSREL)
         SIGH  = DGAUSS(DSIGHA ,A,B,EPSREL)
         SIGFB = (BBETA*SIG0+SIGVS)*(1.D0+DVS)+SIGH
       ELSE
         SIGFB = SIG0FB
       ENDIF
C Compute denominator
       ISEL = 2
       SIG0TOT = CSXA(S)
       SIG0 = SIG0TOT
       IF(IRAD.GE.1)THEN
         SIGVS = DGAUSS(DSIGVSA,A,B,EPSREL)
         SIGH  = DGAUSS(DSIGHA ,A,B,EPSREL)
         SIGTOT = (BBETA*SIG0+SIGVS)*(1.D0+DVS)+SIGH
       ELSE
         SIGTOT = SIG0TOT
       ENDIF
C Theoretical asymmetry
       FBSMATR = SIGFB/SIGTOT
       FBSMATRB = SIG0FB/SIG0TOT
      RETURN
      END
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      DOUBLE PRECISION FUNCTION CSXA(SS)
* one of the simplest approximations for the lowest-order csx
      IMPLICIT DOUBLE PRECISION (A-H,L,M,O-Z)
      REAL*8 RT,JT,RFB,JFB,RT0,JT0,RFB0,JFB0
      REAL*8 RFBTAU,JFBTAU
      COMMON /SIGMAS/ SGA
      COMMON /FINFER/ IFLAV
      COMMON /UNKNOWN/ MZ,GZ,RT(4),JT(4),RFB(3),JFB(3)
      COMMON /INITIAL/ RT0(4),RFB0(3),JT0(4),JFB0(3)
      COMMON /RAD/ IRAD,SMIN
      COMMON /BEAM/S, LNS, BB, AP
      COMMON /WHICH/ISEL
      PARAMETER (Z2=1.6449340668482264D0,Z3=1.2020569031595943D0,
     .           PI=3.1415926535897932D0,
     .           ALFINV=137.036D0,
     .           GEVNB=0.38937966D6,
     .           M=0.0005110034D0,
     .           MMU=0.10565839D0,
     .           MTAU=1.7841D0,
     .           MPI=0.1349734D0)
C
C
      ALFA2 = 1.D0/ALF(SS)**2
      IF (ISEL.EQ.1) THEN
C Electrons
      IF (IFLAV.EQ.1) THEN
        CSXA = 4.D0*PI/(3.D0*ALFINV**2)*(SS*RFB(1)+(SS - MZ**2)*JFB(1))
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2)
      ENDIF
C Muons
      IF (IFLAV.EQ.2) THEN
        CSXA = 4.D0*PI/(3.D0*ALFINV**2)*(SS*RFB(2)+(SS - MZ**2)*JFB(2))
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2)
      ENDIF
C Taus
      IF (IFLAV.EQ.3) THEN
CFT take into account mass effect for taus.
        RFBTAU = RFB(3)*0.9992D0
        JFBTAU = JFB(3)*0.9992D0
        CSXA = 4.D0*PI/(3.D0*ALFINV**2)*(SS*RFBTAU+(SS - MZ**2)*JFBTAU)
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2)
      ENDIF
      IF(IRAD.EQ.2)THEN
        Z = SS/S
        CSXA = CSXA * 4.D0*Z/(1.D0+Z)**2
      ENDIF
      CSXA = 3.D0/4.D0 * CSXA
      ELSEIF (ISEL.EQ.2) THEN
C Electrons
      IF (IFLAV.EQ.1) THEN
        CSXA = ALFA2 * SGA / SS +
     .        4.D0*PI/(3.D0*ALFINV**2)*(SS*RT(2)+(SS - MZ**2)*JT(2))
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2)
      ENDIF
C Muons
      IF (IFLAV.EQ.2) THEN
        CSXA = ALFA2 * SGA / SS +
     .        4.D0*PI/(3.D0*ALFINV**2)*(SS*RT(3)+(SS - MZ**2)*JT(3))
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2)
      ENDIF
C Taus
      IF (IFLAV.EQ.3) THEN
CFT take into account mass effect for taus.
        RTTAU = RT(4)*0.998D0
        CSXA = ALFA2 * SGA / SS +
     .        4.D0*PI/(3.D0*ALFINV**2)*(SS*RTTAU+(SS - MZ**2)*JT(4))
     .        /((SS-MZ**2)**2 + (GZ*MZ)**2)
      ENDIF
      ENDIF
      CSXA = GEVNB * CSXA
      RETURN
      END

C
CFT      IF(ISEL.EQ.1)THEN
CFT        CSXA =
CFT     .  ( 2.D0*GAE*GAF* ( (SS-MZ**2) + SS*GZ/MZ*DELFTI ) * CGMU +
CFT     .     4.D0*GVE*GVF*GAE*GAF*SS *CGMU**2 ) * BW * QED * 3.D0/4.D0
CFT        IF(IRAD.EQ.2)THEN
CFT          Z = SS/S
CFT          CSXA = CSXA * 4.D0*Z/(1.D0+Z)**2
CFT        ENDIF
CFT      ELSEIF(ISEL.EQ.2)THEN
CFT        CSXA = ( ( 2.D0*GVE*GVF*(SS-MZ**2) *CGMU +
CFT     .    (GVE**2+GAE**2)*(GVF**2+GAF**2)*SS *CGMU**2 ) * BW
CFT     .            + 1.D0 ) * QED
CFT      ENDIF
CFT      END
C
C
C
      DOUBLE PRECISION FUNCTION DSIGVSA(X)
*********************************************************************
* Calculation of integral over the soft-photon csx dsig(s')/ds':
* beta*(1-s'/s)**(beta-1).  The almost-pole at s'=s has been split off.
* This calculation sums over all orders of alfa ("exponentiation")
*********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,L-Z)
      COMMON/BEAM/S, LNS, BB, AP
      COMMON/SIG/ SIG0
        BETA = 2.D0*AP*BB
        FBETA  = BETA*DEXP(BETA*LOG(X))/X
        SX  = S*(1.D0-X)
        DSIGVSA = FBETA*(CSXA(SX)-SIG0)
      END
C
C
C
      DOUBLE PRECISION FUNCTION DSIGHA(X)
*********************************************************************
* Calculation of the non-IR part of dsig/ds'.                       *
* 1st+2nd order calculation.                                        *
*********************************************************************
      IMPLICIT DOUBLE PRECISION (A-H,L-Z)
      COMMON/BEAM/S, LNS, BB, AP
      COMMON/SIG/ SIG0
      COMMON /RAD/ IRAD,SMIN
      COMMON /WHICH/ISEL
      PARAMETER (Z2=1.6449340668482264D0,Z3=1.2020569031595943D0)
       Z   = 1.D0-X
       SX  = S*Z
       LNZ = LOG(Z)
       LN1Z= LOG(X)
       LNZ2= LNZ**2
       LNV2= LN1Z**2
       LNVZ= LN1Z*LNZ
       IF(X.LE. 0.5D0) THEN
         U =-LNZ
       ELSE
         U =-LN1Z
       ENDIF
       LI2 = ((((((-1.D0/10886400.D0)*U**2+1.D0/211680.D0)*U**2-
     .         1.D0/3600)*U**2+1.D0/36.D0)*U-1.D0/4.D0)*U+1.D0)*U
       S12 = ((((((-1.D0/24192000.D0)*U**2+1.D0/483840D0)*U**2
     .  -1.D0/8640.D0)*U**2+1.D0/96.D0)*U-1.D0/12.D0)*U+1.D0/4.D0)*U**2
       LI3 = ((((((((((11.D0/435456000D0)*U+11.D0/127008000D0)*U-
     .         1.D0/752640.D0)*U-19.D0/5556600.D0)*U+
     .         7.D0/86400.D0)*U+7.D0/54000.D0)*U-5.D0/576.D0)*U+
     .         17.D0/216.D0)*U-3.D0/8.D0)*U+1.D0)*U
       IF(X.GT. 0.5D0) THEN
          SSS= S12
          S12=-LI3+LNZ*LI2+LN1Z*LNZ**2/2.D0+Z3
          LI3=-SSS-LN1Z*LI2-LNZ*LN1Z**2/2.D0+LN1Z*Z2+Z3
          LI2=-LI2+Z2-LNVZ
       ENDIF
       RZ  = 2./X-2.+X
      C2= RZ*(-LNZ)+(1.D0+Z)*(0.5*LNZ-2.*LN1Z)-5./2.-0.5*Z
      C1= RZ*(+LNZ*LN1Z+3.5D0*LNZ+LI2-LNZ**2/2.D0)+
     .       3.D0*(1.D0-Z)+0.5D0-LNZ+
     .     (1.D0+Z)*(0.25*LNZ**2+4.D0*LN1Z-2.D0*Z2+3.5D0)
      DX  = (-1.D0-2.D0*LNZ/X-LNZ2/X**2)/3.D0
      C0= RZ*(LNZ2*LN1Z/2.D0+LI2*LNZ/2.D0-LNZ**3/6.D0+Z2*LNZ-
     .            1.5D0*LI2-1.5D0*LNVZ-17.D0*LNZ/6.D0+LNZ2/6.D0
     .           +LNZ2/X+DX)
     .   +(1.D0+Z)*(1.5D0*LI3-2.D0*S12-LI2*LN1Z-Z2+LNZ2/6.D0-1.D0/6.D0+
     .            LNVZ/2.D0-LNV2/4.D0+5.D0*LNZ/2.D0+3.D0*LN1Z/2.D0+
     .            2.D0*Z2-2.D0-2.D0*LN1Z)
     .   +   Z  *(16.D0*Z2/3.D0-25.D0*LI2/6.D0-13.D0*LNZ2/12.D0-
     .    5.D0*LNZ/6.D0-4.D0*LNVZ+1.5D0*LNV2-15.D0*LN1Z/6.D0-2.D0/3.D0)
      CSXSX = CSXA(SX)
      DSIGH1 = -AP*BB*(1.D0+Z)*CSXSX
      IF(ISEL.EQ.1.AND.IRAD.EQ.2)THEN
       DSIGH1 = DSIGH1
     .         -AP*(1.D0+Z**2)/(1.D0-Z)*DLOG(4.D0*Z/(1.D0+Z)**2)*CSXSX
      ENDIF
      DSIGH2 =  AP**2*(C2*LNS**2+C1*LNS+C0)*CSXSX
      DSIGHA  = DSIGH1+DSIGH2
      END

