*
* $Id$
*
* $Log$
* Revision 1.1  2000/06/19 20:00:34  eugenio
* Initial revision
*
* Revision 1.1.1.1  1994/10/08  02:21:28  zfiles
* first version of qqlib in CVS
*
*
#include "sys/CLEO_machine.h"
#include "pilot.h"
*CMZ :  1.01/00 03/11/90  16.54.51  by  Paul Avery
*CMZ :  1.00/00 04/06/90  18.55.47  by  Jorge L. Rodriguez
*-- Author :

      SUBROUTINE QED031(ARG1,ARG2,ARG3,IARG4,ARG5)
C  ---------------------------------------------------------------------------
C   MAIN ROUTINE FOR RADIATIVE TAU PRODUCTION.  SPECIAL FORM
C   FOR THE HARD BREMSTRAHLUNG INTEGRAL CONTAINING THE
C   CROSS SECTION FOR THE TURN-ON OF A HEAVY LEPTON PAIR
C   IS USED HERE.  (SEE B+K PAPER ON INITIAL STATE RADIATION
C   APPENDIX B).
C  ---------------------------------------------------------------------------
C beam energy and mass to be given in gev
      IMPLICIT DOUBLE PRECISION(A-H,K-Z)

      COMMON/QED000/IMODDD,IOUT
      COMMON / GETAL / BARN,PI,ALFA,BETA,DELTA,SIG0,AME,AMT,EME,S
      COMMON / GRENS / KMAX,A ,CMAX,SOFT,TRSHLD,K0
      COMMON/TRY5/ITRY,IEV
C constants
      ENERGY=ARG1
      AMAX  =ARG2
      MASS  =ARG3
      INCL  =IARG4
      KMAXX =ARG5
      K0   =1.D-2
      ITRY=0
      IEV =0
      BARN=(1.9732858D0)**2*(1.D+05)
      PI=3.1415927D0
      ME=0.5110034D-03
      MM=105.65948D-03
      MT=1.784D0
      ALFA=1./137.035982D0
      CMAX=DCOS(AMAX*PI/180.D0)
C DIMNESIONLESS QUANTITIES AND S
      S=4.D0*ENERGY**2
      LE=DLOG(S/ME**2)
      LM=DLOG(S/MM**2)
      FAKTOR=2.D0*ALFA/PI
      BETA=FAKTOR*(LE-1.D0)
      AME=(ME/ENERGY)**2
      EME=DSQRT(1.D0+AME)
C..NOTE...AT=STANDARD TAU MASS, OF NO CONCERN TO OUTSIDE ROUTINES
      AT=MT/ENERGY
C....AMT=THE INPUTTED FINAL LEPTON MASS, NEEDED EVERYWHERE
C...SIGNAL AMT=0.D0 IF JUST INITIAL STATE RADIATION WANTED (BREVNT)
      AMT=(MASS/ENERGY)**2
C ELECTRON VERTEX CORRECTION AND VACUUM POLARIZATION
      DVTX=FAKTOR*(0.75D0*LE+PI**2/6.D0-1.D0)
      DELEC=FAKTOR*(LE/3.D0-5.D0/9.D0)
C MUON VACUUM POLARIZATION
      DMUON=FAKTOR*(LM/3.D0-5.D0/9.D0)
C TAU VACUUM POLARIZATION
      A2=1.D0-AT*AT
      A=DSQRT(A2)
      LOGB=DLOG(AT**2/(1.D0+A)**2)
      DTAU=-FAKTOR*(8.D0/9.D0-A2/3.D0+A*(0.5D0-A2/6.D0)*LOGB)
C HADRONIC VACUUM POLARIZATION
C-----DELTH=-2.D0*PIH(S*1.D+6)----IS TOO SLOW
      DELTH=-2.D0*PIHINT(S)
C FINITE PART ANALYTIC CORRECTION
      DELTA=DVTX+DELEC+DMUON+DTAU+DELTH
      SOFT=1.D0+DELTA+BETA*DLOG(K0)
C LOWEST ORDER CROSSECTION
      SIG0=4.D0*ALFA**2*PI/3.D0/S*BARN
C----------INCLUDE MASSEFFECTS FROM OUTGOING LEPTON-------
      SIG0=SIG0*DSQRT(1.D0-AMT)*(1.D0+0.5D0*AMT)
C---APPLY ANGULAR CUT
      SIG00=SIG0 / (4./3.) * CMAX*(1.+CMAX**2/3.)
C TOTAL CROSSECTION AND CORRECTION
      KMAX=1.D0-AMT
C---THE FOLLOWING FOR BRSTRT PURPOSES
C---FIXUP 5.JUNE.85   0.5 GEV INTO 0.5 MEV
      IF(AMT.EQ.0.D0) KMAX=1.D0-0.0005/ENERGY
C------FIXUP 5.JUNE.85------
      IF (KMAXX.GT.0.D0.AND.KMAXX.LT.1.D0) KMAX=DMIN1(KMAX,KMAXX)
      SIGT=Y3K(KMAX)
      DELTOT=SIGT/SIG00-1.D0
      TRSHLD=SOFT*SIG0/SIGT
C PRINT OUT RESULTS
      WRITE (IOUT,2) ENERGY,AMAX,SIG00,KMAX,DVTX,DELEC,DMUON,DTAU,DELTH
C*****TYPE         2,ENERGY,AMAX,SIG00,KMAX,DVTX,DELEC,DMUON,DTAU,DELTH
    2 FORMAT(1H ,/
     Z      ,20X,'BEAM ENERGY                        =',F10.5,' GEV',/
     Z      ,20X,'ANGLE CUT AT                       =',F10.5,
     Z' DEGREES',/
     Z      ,20X,'LOWEST ORDER CROSSECTION           =',F10.5,' NB',/
     Z      ,20X,'UPPER LIMIT ON K/E                 =',F10.5,/
     Z      ,20X,'VERTEX CORRECTIONS                 =',2PF10.5,' %',/
     Z      ,20X,'ELECTRON VACUUM POLARIZATION       =',2PF10.5,' %',/
     Z      ,20X,'MUON   VACUUM POLARIZATION         =',2PF10.5,' %',/
     Z      ,20X,'TAU    VACUUM POLARIZATION         =',2PF10.5,' %',/
     Z      ,20X,'HADRON VACUUM POLARIZATION         =',2PF10.5,' %')
      IF(MASS.GT.0.)WRITE(IOUT,233) MASS
C*****IF(MASS.GT.0.)TYPE 233,       MASS
  233 FORMAT(20X,'HEAVY PARTICLE MASS                =',F10.5,' GEV')
C
C HEAVY PARTICLE VACUUM POLARIZATION
C
      IF(INCL.EQ.0)GOTO 43
      A2=1.D0-AMT
      A=DSQRT(A2)
      LOGB=DLOG(AMT/(1.D0+A)**2)
      DHEAVY=-FAKTOR*(8.D0/9.D0-A2/3.D0+A*(0.5D0-A2/6.D0)*LOGB)
      DELTA=DELTA+DHEAVY
      WRITE (IOUT,3) DHEAVY
C*****TYPE         3,DHEAVY
    3 FORMAT(20X,'HEAVY PARTICLE VACUUM POLARIZATION =',2PF10.5,' %')
   43 CONTINUE
C
      WRITE (IOUT,4) DELTA,SIGT,DELTOT
C*****TYPE         4,DELTA,SIGT,DELTOT
    4 FORMAT(20X,'FINITE ANALYTIC CORRECTION         =',2PF10.5,' %',/
     Z      ,20X,'CORRECTED CROSSECTION              =',0PF10.5,' NB',/
     Z      ,20X,'TOTAL CORRECTION                   =',2PF10.5,' %')
      RETURN
C
C
      ENTRY QED034
      RETURN
C
C
      ENTRY QED033
      WRITE (IOUT,3331) ITRY,IEV
C*****  TYPE         3331,ITRY,IEV
3331  FORMAT(//,' TRIED',I10,' EVENTS',/,
     *    ' ACCEPTED',I10,' EVENTS AFTER ANGLE CUT')
      CROSS=SIGT*(IEV+0.D0)/(ITRY+0.D0)
      DCROSS=CROSS/SIG00-1.D0
      WRITE (IOUT,3332) CROSS,DCROSS
C*****TYPE         3332,CROSS,DCROSS
3332  FORMAT(20X,'FINAL CROSSECTION     =',0PF10.5,' NB',/,
     *         20X,'TOTAL CORRECTION      =',2PF10.5,' %')
      RETURN
      END