*
* $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 17/09/90  23.22.39  by  Paul Avery
*CMZ :  1.00/00 04/06/90  18.55.47  by  Jorge L. Rodriguez
*-- Author :

      SUBROUTINE QED032(QPP,QM,PH)
      IMPLICIT DOUBLE PRECISION(A-H,K-Z)

      DIMENSION QPP(4),QM(4),PH(4),AA(4)
      COMMON / GETAL / BARN,PI,ALFA,BETA,DELTA,SIG0,AME,AMT,EME,S
      COMMON / GRENS / KMAX,A ,CMAX,SOFT,TRSHLD,K0
      COMMON / OUT11 / K,Z,ZM,ZP,SZ,FG,SFG,CFG
      COMMON / INV1 / EPLMN,PYPLMN,PZPLMN,BOOST1,BOOST2,BOOST3,IPLUS
      EXTERNAL Y3KQ,Y3ZQ
      COMMON/TRY5/ITRY,IEV
      REAL RANP
      EXTERNAL RANP
      IEV=IEV+1
      CMXX=1.D0-CMAX
      DMXX=CMAX
111   ITRY=ITRY+1
      IF(RANP(0).GT.TRSHLD) GOTO 911
C----SOFT REGIME------
      K=0.D0
      Z=1.D0
      ZM=0.D0
      ZP=2.D0
      SZ=0.D0
      FG=0.D0
      CFG=1.D0
      SFG=0.D0
      FP=1.D0
      FM=1.D0
      S1=4.D0
      W1=2.D0
      GOTO 910
C----HARD REGIME------
911   CALL EVGENH(Y3KQ,0.D0,DLOG(KMAX/K0),QKQK)
      K=K0*DEXP(QKQK)
      EPSLN=AME/2.
      CALL EVGENH(Y3ZQ,1./(EPSLN+2.),1./EPSLN,ZQZQ)
      ZM=1./ZQZQ-EPSLN
      Z=1.D0-ZM
      ZP=2.D0-ZM
      IF(RANP(0).LT.0.5D0) GOTO 311
            ZP=ZM
            Z=ZP-1.D0
            ZM=2.-ZP
311   SZ=DSQRT(ZM*ZP)
      FG=2.D0*PI*RANP(0)
      SFG=DSIN(FG)
      CFG=DCOS(FG)
C CALCULATE K-DEPENDENT INVARIANTS IN LAB FRAME
      DELDEL=AME/(EME+1.D0)      !=EME-1
      ZETA=Z/EME
      IF(ZETA.LT.0.) GOTO 312
            OMZETA=(DELDEL+ZM)/EME
            OPZETA=1.D0+ZETA
            GOTO 313
312         OMZETA=1.D0-ZETA
            OPZETA=(DELDEL+ZP)/EME
313   KM=1.D0-K
      S1=4.D0*KM
      W1=DSQRT(S1)
C....FP AND FM ARE NEEDED ONLY AS RELATIVE WEIGHTS IN Y3C,Y3F
      FP=-2.D0*AME*S1*(OMZETA/OPZETA)+(S1+2.D0*OPZETA)**2
      FM=-2.D0*AME*S1*(OPZETA/OMZETA)+(S1+2.D0*OMZETA)**2
C
C PICK AXIS ALONG WHICH TO MEASURE ANGLE OF FINAL LEPTON
C (EITHER R(+) OR R(-) )
C.......IPLUS=1 SELECT R(+),  IPLUS=-1 SELECT R(-)
C
910   PH(4)=K
      PH(3)=K*Z
      PH(1)=K*SZ*SFG
      PH(2)=K*SZ*CFG
      IPLUS=1
      IF(RANP(0)*(FP+FM).GT.FP) IPLUS=-1
C---ENERGY/MOMENTUM VECTOR OF E+ (OR E-)
      IF(K.NE.0.D0) GOTO 411
            EPLMN=EME
            PZPLMN=IPLUS

            GOTO 412
411   EPLMN=W1/2.+K/W1*(EME+IPLUS*Z)
      PYPLMN=EME+K*IPLUS*Z/(2.*EME-K+W1)
      PZPLMN=IPLUS + PYPLMN*K/W1*Z
      PYPLMN=        PYPLMN*K/W1*SZ
      PPPLMN=DSQRT(PYPLMN**2+PZPLMN**2)
      PYPLMN=PYPLMN/PPPLMN
      PZPLMN=PZPLMN/PPPLMN
      BOOST1=(2.*EME-K)/W1
      BOOST2=-K*SZ/W1
      BOOST3=-K*Z/W1
C
C HAVE NOW THE BOOST FACTORS AND ROTATION STUFF TOGETHER
C
412   IF(AMT.EQ.0.D0) RETURN
C
C
C NOW PICK TETA AND FAI FOR THE FINAL LEPTONS IN THE C.M. FRAME
C     OF THE FINAL STATE.
      A=4.D0*AMT/S1
      A1=1.D0+A
      A2=1.D0-A
CCCCCCCALL EVGENH(Y3CC,-1.D0,1.D0,C)
      C=RANP(0)*2.D0-1.D0
      IF(RANP(0)*(A1+A2/3.D0).GT.A1) C=DSIGN(DABS(C)**0.333333D0,C)
CCCCCCCCCCCCCCCCCCC
      SC=DSQRT(1.D0-C*C)
      FI=2.D0*PI*RANP(0)
      QQQ=W1/2.*DSQRT(A2)
      QPP(1)=QQQ*SC*DCOS(FI)
      QPP(2)=QQQ*SC*DSIN(FI)
      QPP(3)=QQQ*C
      QPP(4)=W1/2.
C...APPLY ROTATION+BOOST+FIROTATION
      CALL QED03T(QPP)
      IF(QPP(3)**2*CMXX.GT.(QPP(1)**2+QPP(2)**2)*DMXX) GOTO 111
      QM(1)=-QPP(1)-PH(1)
      QM(2)=-QPP(2)-PH(2)
      QM(3)=-QPP(3)-PH(3)
      QM(4)=2.D0-QPP(4)-PH(4)
      IF(QM(3)**2*CMXX.GT.(QM(1)**2+QM(2)**2)*DMXX) GOTO 111
C
C
C
      RETURN
      END