*
* $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.03/10 12/06/91  11.31.52  by  Peter C Kim
*CMZ :  1.01/00 17/09/90  23.22.39  by  Paul Avery
*CMZ :  1.00/00 04/06/90  18.55.46  by  Jorge L. Rodriguez
*-- Author :

      SUBROUTINE QED026
C  ----------------------------------------------------------------------------
C...FORMERLY CALLED FIXFOT + SOFOT
      IMPLICIT DOUBLE PRECISION(A-H,K,O-Z)
C...GENERATES SOFT PHOTON APPROXIMATION WITH EXPONENTIATION
C  ----------------------------------------------------------------------------

      COMMON / OUTPAR / K,C,CG,FIG,FI
      COMMON / CONSTQ / E,AE,AM,ALE,ALM,ALL,PI,PI2,ALFA,BETA
      COMMON /CONST1/ CMAX,K0,KMIN,KMAX,SIG0,SOFT1,TRSHLD,CHI,GV2,DELW
      REAL RANP
      EXTERNAL RANP
C----DEFINE DOUBLE PRECISION RANDOM FUNCTION
      DRAN(IRAN)=RANP(0)+0.D0
C
C
C
      AE2=AE**2
      AM2=AM**2
      K0BETA=K0**BETA
      C0=BETA*SOFT1/K0BETA
      D1=BETA*(1.D0-BETA)*(SIG0-SOFT1)/K0
      C1=(BETA*SIG0-D1*K0)/K0BETA
      WMAX=SIG0/SOFT1
    1 GOTO 311
  312 WEIGH=(C1+D1*K**(1.D0-BETA))/C0
      IF(WMAX*DRAN(0).GT.WEIGH) GOTO 1
      CG=Y
      FIG=FG
      IF(I.EQ.1) RETURN
      FIG=DATAN2(DSIN(FG) ,
     *  C*DCOS(FG)-CG*DSQRT((1.-C)*(1.+C)/(1.-CG)/(1.+CG))   )
      RETURN
C--------------------
  311 CONTINUE
CCCCCCSUBROUTINE SOFOT(AE2,AM2,C,K0,BETA,PI,K,Z,FG)
C
C GENERATES FOTON ANGLE ACCORDING TO:
C
C D/DOMEGA = 1/DOT(E+,K) + 1/DOT(E-,K) + 1/DOT(MU+,K) + 1/DOT(MU-,K)
C
C DONE BY A) CHOOSING IF RADIATING ALONG E+- OR MU+-
C         B) IMPOSING SIGN SYMMETRY AFTERWARDS
C         C) THIS LEAVES SIMPLY A 1/DOT(E,K) TERM, WHATEVER THE 'E' IS.
C
C THEN A WEIGHTING FUNCTION IS IMPOSED AT THE END.
C
      K=K0*DRAN(IRAN)**(1./BETA)
      EPSE=AE2/2.D0
      EPSM=AM2/2.D0
      WMIXT=DLOG(2.D0/EPSE) / DLOG(4.D0/EPSE/EPSM)
  313 FG=2.D0*PI*DRAN(IRAN)
      IF(DRAN(IRAN).GT.WMIXT) GOTO 411
      I=0
C---RADIATE OFF E+,E-
      OMZ=EPSE*((2.D0/EPSE+1.D0)**DRAN(IRAN) - 1.D0)
      OPZ=2.D0-OMZ
      IF(DRAN(IRAN).GT.0.5D0) GOTO 412
      OPZ=OMZ
      OMZ=2.D0-OPZ
  412 Z=(OPZ-OMZ)/2.D0
      Y=C*Z + DSQRT((1.D0-C)*(1.D0+C)*OPZ*OMZ)*DCOS(FG)
      OMY=1.D0-Y
      OPY=1.D0+Y
      GOTO 500
  411 I=1
C---RADIATE OFF MU+,MU-
      OMY=EPSM*((2.D0/EPSM+1.D0)**DRAN(IRAN) - 1.D0)
      OPY=2.D0-OMY
      IF(DRAN(IRAN).GT.0.5D0) GOTO 417
      OPY=OMY
      OMY=2.D0-OPY
  417 Y=(OPY-OMY)/2.D0
      Z=C*Y + DSQRT((1.D0-C)*(1.D0+C)*OPY*OMY)*DCOS(FG)
      OMZ=1.D0-Z
      OPZ=1.D0+Z
  500 YP=OPY+EPSM
      YM=OMY+EPSM
      ZP=OPZ+EPSE
      ZM=OMZ+EPSE
      EXACT=OMZ*OPZ*YP*YM/ZP/ZM
     *     +OMY*OPY*ZP*ZM/YP/YM
     *     +0.5D0*(1.D0-C)*(ZP*YP+ZM*YM)
     *     -0.5D0*(1.D0+C)*(ZP*YM+ZM*YP)
      PROXIM=ZP*ZM+YP*YM
      WEIGH=EXACT/PROXIM
C8) TEST FOR EXIT ASSUMING MAXIMUM WEIGHT SMALLER THAN 10
      IF(DRAN(IRAN)*10.D0.LE.WEIGH) GOTO 312
      GO TO 313
      END