      COMMON/PSIN/XINC,XTAR,PINC,PTAR,NPTCLS,XFIN
      COMMON/PSOUT/PFIN
      DIMENSION PINC(4),PTAR(4),XFIN(20),PFIN(20,4)
C        APRIL 1975
C                  PHASE SPACE EVENTS
C
C
      COMMON /LUNS/ ICRD,ILPT,IKBD,ITTY
	call TIMRB
C        T1=SECNDS(0.)
        write (itty,*) 'BENCHMARK -- LUSTY/10'
        CALL TIMRB
       XINC=0.0
      XTAR=0.938
      PINC(1)=0.0
      PINC(2)=0.0
      PINC(3)=50.0
      PTAR(1)=0.0
      PTAR(2)=0.0
      PTAR(3)=0.0
C
C      NEV=6812
        NEV=681
      NPTCLS=10
      XFIN(1)=0.1396
      XFIN(2)=0.1396
      XFIN(3)=0.1396
      XFIN(4)=0.1396
      XFIN(5)=0.1396
      XFIN(6)=0.1396
      XFIN(7)=0.1396
      XFIN(8)=0.1396
      XFIN(9)=0.1396
      XFIN(10)=0.938
C
      INIT=1234
      CALL PSEVNT(INIT)
C
      DO 200 I=1,1
      DO 100 J=1,NEV
      CALL PSEVNT(0)
100   CONTINUE
200   CONTINUE
C
        CALL TIMRE
c        WRITE(6,999) TIMU
999     FORMAT(5X,E10.4)
      STOP
      END
      SUBROUTINE PSEVNT(INIT)
C           *  GENERATES PURE PHASE SPACE MONTE CARLO EVENTS  *
C           *     WITH FROM 2 TO 20 FINAL STATE PARTICLES     *
      COMMON/PSIN/XINC,XTAR,PINC,PTAR,NPTCLS,XFIN
      DIMENSION PINC(4),PTAR(4),XFIN(20),PFIN(20,4)
      COMMON/PSOUT/PFIN
      COMMON/DTG/Q,X,Z
      DIMENSION Q(20,4),X(20),Z(20)
      DIMENSION ZL(20),POWER(19)
      PS(A,B,C) = ((A+B+C)*(A-B-C)*(A+B-C)*(A-B+C))/(A*A)
      RS(A,B,C) = ((A+B+C)*(A-B-C)*(A+B-C)*(A-B+C))
    1 IF (INIT .EQ. 0) GO TO 200
C
C     INITIALIZE.
    2 N = NPTCLS
      NP1 = N + 1
      NM1 = N - 1
      NM2 = N - 2
      DO 4 I=1,N
    4 X(I) = XFIN(I)
      PINC(4) = SQRT(PINC(1)**2 + PINC(2)**2 + PINC(3)**2 + XINC**2)
      PTAR(4) = SQRT(PTAR(1)**2 + PTAR(2)**2 + PTAR(3)**2 + XTAR**2)
      DO 5 I=1,4
    5 Q(N,I) = PINC(I) + PTAR(I)
      E = SQRT(XINC**2 + XTAR**2 + 2.*PINC(4)*PTAR(4))
C      WRITE(6,6)
    6 FORMAT (1H1,10X,51H*  INITIALIZING PARAMETERS FOR SUBROUTINE PSEVN
     1T  *//)
C      WRITE(6,8) XINC
    8 FORMAT (1H0,18X,24HINCIDENT PARTICLE MASS =,F8.5,4H GEV)
C      WRITE(6,10) XTAR
   10 FORMAT (1H0,18X,24H  TARGET PARTICLE MASS =,F8.5,4H GEV)
C      WRITE(6,12) (PINC(J), J=1,3)
   12 FORMAT (1H0,3X,32HINCIDENT PARTICLE (PX,PY,PZ) = (,F8.4,1H,,F8.4,1
     1H,,F8.4,7H) GEV/C)
C      WRITE(6,13) (PTAR(J), J=1,3)
   13 FORMAT (1H0,3X,32H  TARGET PARTICLE (PX,PY,PZ) = (,F8.4,1H,,F8.4,1
     1H,,F8.4,7H) GEV/C)
C      WRITE(6,14) E
   14 FORMAT (1H0/19X,21HTOTAL C.O.M. ENERGY =,F9.5,4H GEV)
C      WRITE(6,16)
   16 FORMAT (1H0/21X,33HFINAL STATE PARTICLE MASSES (GEV)/)
C      WRITE(6,18) (I,X(I), I=1,N)
   18 FORMAT (31X,I2,2X,F8.5)
C      WRITE(6,20) INIT
   20 FORMAT (1H0/30X,7HINIT = ,I6///)
      Z(N) = E
      Z(1) = X(1)
      ZL(1) = X(1)
      DO 22 I=2,N
   22 ZL(I) = ZL(I-1) + X(I)
      AA= RRN(INIT)
      AA= Z(N) - ZL(N)
      A=AA
      IF (AA.GT. 0.) GO TO 26
C      WRITE(6,24)
   24 FORMAT (71H0*****  REACTION NOT ENERGETICALLY POSSIBLE.  PROGRAM T
     1ERMINATED. *****/1H0,18X,36HCHECK ABOVE PARAMETERS FOR MISTAKES.)
      STOP
   26 IF (NM2 .EQ. 0) GO TO 124
      DO 28 I=2,NM1
      FIM1 = I - 1
   28 POWER(I) = 1./FIM1
C     FIND WMAXS.
  100 DO 104 J=1,NM2
      I = N - J
  104 Z(I) = ZL(I) + (Z(I+1) - ZL(I+1)) * .5**POWER(I)
      WMAXS = RS(Z(N),Z(NM1),X(N))
      DO 106 KK=1,NM2
      K = N - KK
  106 WMAXS = WMAXS*PS(Z(K),Z(K-1),X(K))
      DZ =AA/3.
  108 DZ = DZ/5.
      IF (DZ .LT. A/10000.) GO TO 122
  110 KDZ = 0
      DO 120 J=1,NM2
      I = N - J
  112 Z(I) = Z(I) + DZ
      IF (Z(I) .GE. (Z(I+1)-X(I))) GO TO 115
      WS = RS(Z(N),Z(NM1),X(N))
      DO 113 KK=1,NM2
      K = N - KK
  113 WS = WS*PS(Z(K),Z(K-1),X(K))
      IF (WS .LE. WMAXS) GO TO 115
      WMAXS = WS
      KDZ = 1
      GO TO 112
  115 Z(I) = Z(I) - DZ
  116 Z(I) = Z(I) - DZ
      IF (Z(I) .LE. (Z(I-1)+X(I))) GO TO 119
      WS = RS(Z(N),Z(NM1),X(N))
      DO 117 KK=1,NM2
      K = N - KK
  117 WS = WS*PS(Z(K),Z(K-1),X(K))
      IF (WS .LE. WMAXS) GO TO 119
      WMAXS = WS
      KDZ = 1
      GO TO 116
  119 Z(I) = Z(I) + DZ
  120 CONTINUE
      IF (KDZ .EQ. 1) GO TO 110
      GO TO 108
  122 WMAXS = 1.001*WMAXS
  124 CONTINUE
C     WRITE(6,126)
  126 FORMAT (34X,5H* * *//////)
      RETURN
C     END INITIALIZATION
C
C     GENERATE EVENTS.  FIRST GENERATE INTERMEDIATE MASSES.
  200 IF (N - 3) 300,203,204
  203 Z(2) = ZL(2) +AA*RRN(0)
      WS = RS(Z(3),Z(2),X(3))*PS(Z(2),Z(1),X(2))
      GO TO 210
  204 DO 206 J=1,NM2
      I = N - J
      R = RRN(0)
      Z(I) = ZL(I) + (Z(I+1) - ZL(I+1)) * R**POWER(I)
C     CALL ZMOD(I,Z)
206   CONTINUE
      WS = RS(Z(N),Z(NM1),X(N))
      DO 208 KK=1,NM2
      K = N - KK
  208 WS = WS*PS(Z(K),Z(K-1),X(K))
C     CALL WSMOD(WS,WMAXS)
210   CONTINUE
      IF (WS .LE. WMAXS) GO TO 214
C      WRITE(6,212)
  212 FORMAT (72H1*****  WEIGHTING FAULT IN SUBROUTINE PSEVNT.  PROGRAM
     1TERMINATED. *****////)
      STOP
  214 R = RRN(0)
C     NOTE---   (WT .LT. RRN) IF AND ONLY IF (WT**2 .LT. RRN**2).
      IF ((WS/WMAXS) .LT. R**2) GO TO 200
C     HAVE INTERMEDIATE MASSES.  NOW GENERATE VERTEX DECAYS.
  300 DO 302 J=1,NM1
      I = NP1 - J
  302 CALL ISODK(I)
      DO 304 I=1,4
  304 PFIN(1,I) = Q(1,I)
      RETURN
      END
      SUBROUTINE ISODK(I)
      COMMON /GJL/ITSTA,DELSQ
      COMMON/PSIN/XINC,XTAR,PINC,PTAR,NPTCLS,XFIN
      COMMON/PSOUT/P
      COMMON/DTG/Q,X,Z
      DIMENSION PINC(4),PTAR(4),XFIN(20)
      DIMENSION P(20,4),Q(20,4),X(20),Z(20)
      DOUBLE PRECISION P10(3),BETA(3),TLOR(4,4),GAMMA,DELTA,PJ3
      IM1 = I - 1
      U02 = Z(I)/2.
      V = (Z(IM1) - X(I))*(Z(IM1) + X(I))/(2.*Z(I))
      E10 = U02 + V
      E20 = U02 - V
      P120 = SQRT((E10 - Z(IM1))*(E10 + Z(IM1)))
      CPOL = -1.0 + 2.0*RRN(0)
      AZM = 6.2831853*RRN(0)
C     ***************************************************
C     SPECIAL FOR FIRST VERTEX   MAIN ROUTINE FIXES 4-MOM TRANSFER
C     IF (I.EQ.NPTCLS.AND.ITSTA.EQ.1) CALL XXCPOL(CPOL)
C     ***************************************************
      P120XY = P120*SQRT((1.0-CPOL)*(1.0+CPOL))
      P10(1) = P120XY*COS(AZM)
      P10(2) = P120XY*SIN(AZM)
      P10(3) = P120*CPOL
      BETA(1) = Q(I,1)/Q(I,4)
      BETA(2) = Q(I,2)/Q(I,4)
      BETA(3) = Q(I,3)/Q(I,4)
      GAMMA = Q(I,4)/Z(I)
      DELTA = (GAMMA*GAMMA)/(1.0+GAMMA)
      TLOR(1,1) = 1.0 + DELTA*BETA(1)*BETA(1)
      TLOR(2,2) = 1.0 + DELTA*BETA(2)*BETA(2)
      TLOR(3,3) = 1.0 + DELTA*BETA(3)*BETA(3)
      TLOR(1,2) = DELTA*BETA(1)*BETA(2)
      TLOR(1,3) = DELTA*BETA(1)*BETA(3)
      TLOR(2,3) = DELTA*BETA(2)*BETA(3)
      TLOR(1,4) = GAMMA*BETA(1)
      TLOR(2,4) = GAMMA*BETA(2)
      TLOR(3,4) = GAMMA*BETA(3)
      TLOR(4,4) = GAMMA
      TLOR(2,1) = TLOR(1,2)
      TLOR(3,1) = TLOR(1,3)
      TLOR(4,1) = TLOR(1,4)
      TLOR(3,2) = TLOR(2,3)
      TLOR(4,2) = TLOR(2,4)
      TLOR(4,3) = TLOR(3,4)
      DO 5 J=1,4
      PJ3 = TLOR(J,1)*P10(1) + TLOR(J,2)*P10(2) + TLOR(J,3)*P10(3)
      Q(IM1,J) = +PJ3 + TLOR(J,4)*E10
    5 P(I,J) = -PJ3 + TLOR(J,4)*E20
      RETURN
      END
      FUNCTION RRN(K)
C
      IF (K.GT.0) GO TO 100
C
10     RRN=RANDR(START)
       RRN=ABS(RRN)
      RETURN
C
100   START=FLOAT(K)
      START=1./START
      RRN=RANDR(START)
      GO TO 10
      END
      FUNCTION RANDR(STARTR)
      DATA JCALL,ISTART,JSTART,KSTART,LSTART/0,317,4631,1387,6033/
      DATA JCALL/0/
      IF(JCALL.NE.0) GO TO 1
      JCALL=1
      CALL SETUP(ISTART,JSTART)
    1 CALL SHUFEL(ISTART,JSTART,KSTART,LSTART,RAN)
      RANDR=2.*RAN-1.
      RETURN
      END
      SUBROUTINE SETUP(ISTART,JSTART)
      COMMON/RAND/LOC(128),MODUL,MODULL,MDIV,FMODUL
      MODUL=16384
      MDIV=64
      MODULL=8192
      FMODUL=MODUL
      DO 1 I=1,128
      ISTART=ISTART+JSTART
      ISTART=ISTART-(ISTART/MODUL)*MODUL
      JSTART=ISTART+JSTART
      JSTART=JSTART-(JSTART/MODUL)*MODUL
    1 LOC(I)=JSTART
C      WRITE(6,101)(LOC(I),I=1,128)
  101 FORMAT(1H0,16H SET UP COMPLETE,(/,1X,8I10))
      RETURN
      END
      SUBROUTINE SHUFEL(ISTART,JSTART,KSTART,LSTART,RAN)
      COMMON/RAND/LOC(128),MODUL,MODULL,MDIV,FMODUL
      KT=KSTART+LSTART
      KSTART=LSTART
      LSTART=KT-(KT/MODULL)*MODULL
      INDEX=(LSTART/MDIV)+1
      RAN=LOC(INDEX)
      RAN=RAN/FMODUL + 0.000001
      ISTART=ISTART+JSTART
      ISTART=ISTART-(ISTART/MODUL)*MODUL
      JSTART=ISTART+JSTART
      JSTART=JSTART-(JSTART/MODUL)*MODUL
      LOC(INDEX)=JSTART
      RETURN
      END
