Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap1_198111 - decus/20-0009/nvevnt.for
There is 1 other file named nvevnt.for in the archive. Click here to see a list.
      SUBROUTINE EVENT(NERR)
C      ******************    COMMON COMMON   ***************************    0030
      COMMON    MAP(2000),PARS(1000),MISC(27),KLIST(500),MTABLE(2)
      DIMENSION ZMAP(2000)
      DIMENSION  HTABLE(7, 100),  KTABLE(7, 100), OTABLE(7, 50)
      DIMENSION  RTABLE(9, 20, 2),  ITABLE(6, 20), LTABLE(9, 20, 2)
      DIMENSION  JTABLE(7, 50), RBANK(10), DCTGT(3), DCTRA(3)
      DIMENSION  TITLE(48), NC(48), DIRINC(3)
      DIMENSION PARA(1000),NPARA(1000),SNAME(1000),NAME(1000),TABLE(100)
      DIMENSION  TBLMS (30), ADIS (40), ZEROS (10) , AD(50)
      DIMENSION    NTABLE(100), HEAD(11), NBRNCH(10)
      DIMENSION OBANK(90), IBANK(90), OB(90),VM(3), VN(3)
      DIMENSION PPLANE(3),FMASS(8),PTBL(8),ETBL(8),DCSO(3,8),DCSINC(3)
      EQUIVALENCE  (KTABLE,HTABLE,MAP)
      EQUIVALENCE  (OTABLE,JTABLE,MAP(701)), (RTABLE,LTABLE,MAP(1051)),
     1             (ITABLE,MAP(1411)), (VAL,IVAL,MAP(1531)),
     2             (WGT,MAP(1631))
      EQUIVALENCE ( MAP (1731), NC(1) ),( MAP (1870), TITLE(1) )
      EQUIVALENCE  (NCFLAG,MAP(1869)), (WEIGHT,MAP(1978)),
     1             (NTAPE,MAP(1988)), (EINC,MAP(1998)),
     2             (PINC,MAP(1999)), (BINC,MAP(2000))
      EQUIVALENCE  (IPS,MAP(1973)), (WPS,MAP(1974))
      EQUIVALENCE  (MAP(1975), GMINP), (MAP(1976), GMAXP)
      EQUIVALENCE  (MAP(1977), GSCALE), (WSCALE,MAP(1972))
      EQUIVALENCE (MTOT,MAP(1987))                                      7/13/68
      EQUIVALENCE (PARA,NPARA,PARS),(SNAME,NAME,MAP(1))
      EQUIVALENCE (ZMAP,MAP),(TABLE,PARS(101),NTABLE)
      EQUIVALENCE ( TBLMS (1), PARA (1) ),( NTOT, NPARA (99) )
      EQUIVALENCE ( AD   (1), PARA (40) ),( ZEROS (1), PARA (40) )
      EQUIVALENCE (RSCALE, PARA(98)), (DIRINC(1), PARA(95))
      EQUIVALENCE (ADIS(1),AD(11)), (OBANK(1),OB(1))
      EQUIVALENCE (OBANK(1), MAP(1779)),(OBANK(1), IBANK(1))
      EQUIVALENCE (OBANK(1), PPLANE(1)), (OBANK(4), FMASS(1))
      EQUIVALENCE (OBANK(12), PTBL(1)), (OBANK(20), ETBL(1) )
      EQUIVALENCE (OBANK(28), DCSINC(1)),(OBANK(31), DCSO(1))
      EQUIVALENCE  (OBANK(55), VM(1)),  (OBANK(58), VN(1))
      EQUIVALENCE  ( OBANK(61), DCTGT(1)),  (OBANK(64), RBANK(1 ))
      EQUIVALENCE   (OBANK(74), DCTRA(1)) , (BENUC, PARA(94))
      EQUIVALENCE  (PI, MISC), (RADIAN, MISC(2)), (NIT, MISC(3)),
     1             (NOT, MISC(4)), (HEAD, MISC(5)), (NBRNCH, MISC(16)),
     2             (NPAGE, MISC(26)), (NORD, MISC(27))
C      ***************    END  COMMON  COMMON  *************************
      NTRY = 0
 10   NTRY = NTRY+1
C     TRY 50 TIMES BEFORE GIVING UP --- THE FAILURE MAY HAVE BEEN
C          THE RESULT OF AN UNPHYSICAL THROW FOR A RESONANCE
      IF (NTRY-50) 12, 12, 500
 12   CONTINUE
      NERR = 0                                                              0610
      WPS = 1.0                                                             0620
      IPS = IPS                                                             0630
C     FILL UP RESONANCE BANK                                                0640
      DO 50  K = 20, 29                                                     0650
      KK = K + 10                                                           0660
      LL = K - 19                                                           0670
      CALL RANRES (TBLMS(K), TBLMS(KK), RBANK(LL))                          0680
 50   CONTINUE                                                              0690
      DO 150 K = 1,20                                                       0700
      NPROD = ITABLE(1,K)                                                   0710
      IF (NPROD -1) 150,150,75                                              0720
 75   TMASS = RTABLE(9,K,1)                                                 0730
      IF (TMASS) 800, 76, 76
C     A NEGATIVE MASS
 800  NERR = 23
      GO TO 10
 76   CALL FERMI (RTABLE(9, K, 2), BENUC, TP, DCTGT(1))                     0740
      DO 80 KK = 1, NPROD                                                   0750
      CALL GETMS(LTABLE(KK,K,1), FMASS(KK))                                 0760
      IF(FMASS(KK))800,80,80                                                0770
   80 CONTINUE                                                              0800
C     TEST FOR VERTEX A                                                     0810
      IF (K-1) 90,85,90                                                     0820
 85   DO 86 KK = 1,3                                                        0830
 86   DCSINC(KK) = DIRINC(KK)                                               0840
      BMASS = BINC                                                          0850
      EB = EINC                                                             0860
      IF (BMASS)  892, 87, 892
 87   IF (GMINP)  892, 892, 89
C     INVOKE BREMMSTRAHLUNG SPECTRUM                                        0910
 89   CALL RANDOM (RAN)
      PINC = GMINP * EXP( ABS( GSCALE * RAN) )                              0930
      EINC = PINC                                                           0940
      EB = PINC                                                             0950
C    FILL OTABLE(I,1) WITH BEAM QUANTITIES
 892  DO 894 KK=1,3
 894  OTABLE(KK,1) = DCSINC(KK)
      OTABLE(4,1) = PINC
      OTABLE(5,1) = EINC
      GO TO 100                                                             0960
 90   II = ITABLE(3,K)                                                      0970
      KK = ITABLE(4,K)                                                      0980
      CALL GETMS (LTABLE(KK, II, 1), BMASS)                                 0990
      JJ = ITABLE(2,II) + KK                                                1000
      DO 92 KK = 1,3                                                        1010
 92   DCSINC(KK) = OTABLE(KK,JJ)                                            1020
      EB = OTABLE(5,JJ)                                                     1030
 100  IF (BMASS) 800, 1002, 1002
 1002 IF (EB) 108, 1004, 1004
 1004 CONTINUE
C     GET ADIS ENTRY                                                        1050
      LL = ITABLE(5, K) + 1                                                 1060
C     SET PPLANE                                                            1070
      IF (ITABLE(6,K)) 103,103,105
 103  DO 104 KK = 1,3                                                       1120
 104  PPLANE(KK) = 0.0                                                      1130
      GO TO 107                                                             1140
 105  MM = ITABLE(6,K)                                                      1150
       DO 106 KK= 1,3                                                       1160
 106  PPLANE(KK) = OTABLE(KK, MM)                                           1170
 107  CONTINUE                                                              1180
      CALL PSEGEN (TMASS, TP, DCTGT(1), BMASS, EB, DCSINC(1), NPROD, FMA    1200
     1SS(1),      PTBL(1), ETBL(1), DCSO(1,1), AD(LL), PPLANE(1), WATE)     1210
      IF (ETBL(1) + 2.0) 108,108,109                                        1220
C     NOT ENOUGH ENERGY
 108  NERR = 21                                                             1230
      GO TO 10
 109  IF (ETBL(1) )       110,110,115                                       1250
C     REPEATED ZERO MOMENTA (OR DIVIDE CHECKS)
 110  NERR = 22                                                             1260
      GO TO 500                                                             1270
 115  GO TO (118, 116), IPS                                                 1280
 116  WPS = WPS * WATE                                                      1290
 118  CONTINUE                                                              1300
      DO 120 NN = 1, NPROD                                                  1310
      II  =  ITABLE(2, K) + NN                                              1320
      OTABLE(5, II) = ETBL(NN)                                              1330
      OTABLE(4,II) = PTBL(NN)                                               1340
      DO 120 KK = 1,3                                                       1350
 120  OTABLE(KK,II) = DCSO(KK,NN)                                           1360
 150  CONTINUE                                                              1370
 500  CONTINUE                                                              1380
      WPS = WPS*WSCALE                                                      1390
      RETURN                                                                1400
      END                                                                   1410
      SUBROUTINE PSEGEN ( TMASS, TP, DCTGT, BMASS, EB, DCSIN, NPROD,        0010
     1     FMASS, PTBL, ETBL, DCSO, A, PPLANE, W )                          0040
C     VERSION OF 1/69--COMPATIBLE WITH INVARIANT RDECAY                 1/69
C     ALSO RETURNS WEIGHT TO CONVERT TO NON-INVARIANT PHASE SPASE       1/69
      DIMENSION   DCSINC(3)    ,FMASS(8)     ,PTBL(8), DCSIN(3)             0050
      DIMENSION   ETBL(8)      ,DCSO  (3,8)  ,R(3,8) , DCTGT(3)             0060
      DIMENSION   P(8,4)       ,X(8)                                    2/69
      DIMENSION   PP(3,8)      ,ROT(3,3)                                1/69
      DIMENSION A(10)  , PPLANE(3)                                          0090
      CALL DVCHK (K000FX)                                                   0140
      NDVCHK = -1                                                           0150
C     CALCULATE ENERGY IN CENTER OF MASS                                    0190
      PB =  SQRT( EB**2 - BMASS**2)                                         0200
      IF (TP)  11, 3, 11                                                    0210
 3    DO 4 K = 1, 3                                                         0220
 4    DCSINC(K) = DCSIN(K)                                                  0230
      DEN = PB                                                              0240
      TE = TMASS                                                            0250
      IF (TMASS)  5, 5, 10                                                  0260
    5 OMEGA = BMASS                                                         0270
      GO TO 12                                                              0280
   10 OMEGA = SQRT( TMASS ** 2 + BMASS ** 2 + 2.0 * TMASS * EB )            0290
      GO TO 12                                                              0300
 11   DEN = 0                                                               0310
      DO 121 K = 1, 3                                                       0320
      DCSINC(K) = DCSIN(K) * PB + DCTGT(K) * TP                             0330
 121  DEN = DEN + DCSINC(K)**2                                              0340
      TE = SQRT( TMASS**2 + TP**2 )                                         0350
      OMEGA = SQRT(( EB + TE ) **2 - DEN )                                  0360
      DEN = SQRT(DEN)                                                       0370
      DO 131 K = 1, 3                                                       0380
 131  DCSINC(K) = DCSINC(K) / DEN                                           0390
 12   ESUM = OMEGA                                                          0400
      DO 13  K = 1, NPROD                                                   0410
      ESUM = ESUM - FMASS(K)                                                0420
      I = NPROD-K+1                                                     2/69
      X(I) = FMASS(K)                                                   2/69
 13   CONTINUE                                                              0430
      BETA = DEN / ( TE + EB)                                               0440
      GAMMA = (EB + TE ) / OMEGA                                            0450
      IF (ESUM)  14, 15, 15                                                 0460
C     NOT ENOUGH ENERGY * EXIT AND SET ETBL(1) = -10.0                      0470
 14   ETBL(1) = -10.0                                                       0480
      GO TO  401                                                            0490
 15   NDVCHK = NDVCHK + 1                                                   0500
      IF (NDVCHK - 50)  17, 16, 16                                          0510
C     EXIT IF DIVIDE CHECK 50 TIMES AND SET ETBL(1) = -1.0                  0520
 16   ETBL(1) = -1.0                                                        0530
      GO TO 401                                                             0540
 17   CONTINUE                                                              0550
      CALL RDECAY (NPROD, OMEGA,   X,   P, FLAG, A, TMASS, BMASS)       2/69
      IF (FLAG)  14, 18, 16                                             1/69
 18   CONTINUE                                                              0580
C     CALCULATE NON-INVARIANT WEIGHT                                    1/69
      W = 1.0                                                               0600
      DO 100  I = 1, NPROD                                                  0610
      W = W*P(I,4)                                                      1/69
  100 CONTINUE                                                              0640
      IF (NPROD - 2)         401, 2001, 201                                 0650
 2001      CONTINUE                                                         0660
      DO       2002     K = 1, 3                                            0670
      IF (PPLANE(K) )     2003, 2002, 2003                                  0680
 2002    CONTINUE                                                           0690
      GO TO 201                                                             0700
 2003      CONTINUE                                                         0710
C     CALCULATE ROTATION MATRIX FOR TRANSVERSE POLARIZATION                 0720
      ALFC = PPLANE(2) * DCSINC(3) - PPLANE(3) * DCSINC(2)                  0730
      BETC = PPLANE(3) * DCSINC(1) - PPLANE(1) * DCSINC(3)                  0740
      GAMC = PPLANE(1) * DCSINC(2) - PPLANE(2) * DCSINC(1)                  0750
C     NORMALOIZE NORMAL VECTOR                                              0760
      DEN = SQRT( ALFC**2 + BETC**2 + GAMC**2 )                             0770
      IF (DEN) 2008, 2008, 2009                                             0780
C     SET RANDOM NORMAL                                                     0790
 2008 CONTINUE                                                              0800
      CALL RANDOM ( ALFC )                                              6/68
      CALL RANDOM ( RSN )                                               6/68
      BETC  = SIGN( SQRT( 1.0 - ALFC **2), RSN )                            0820
      GAMC  = 0.0                                                           0830
      GO TO 2007                                                            0840
 2009 ALFC = ALFC / DEN                                                     0850
      BETC = BETC / DEN                                                     0860
      GAMC = GAMC / DEN                                                     0870
 2007      CONTINUE                                                         0880
C     SET UP X-AXIS ALONG BEAM DIRECTION AND Z-AXIS ALONG POLARIZATION D    0890
      COST = GAMC                                                           0900
      SINT = SQRT( 1.0 - GAMC**2)                                           0910
      IF (SINT)  2012, 2011, 2012                                           0920
 2011 COSP = 1.0                                                            0930
      SINP = 0.0                                                            0940
      GO TO 2014                                                            0950
 2012 COSP = ALFC/SINT                                                      0960
      SINP = BETC / SINT                                                    0970
 2014 CONTINUE                                                              0980
      ALFCP = DCSINC(1) * COST * COSP + DCSINC(2) * COST * SINP             0990
     1    - DCSINC(3) * SINT                                                1000
      BETCP = - DCSINC(1) * SINP + DCSINC(2) * COSP                         1010
      GAMCP = DCSINC(1) * SINT * COSP + DCSINC(2) * SINT * SINP             1020
     1    + DCSINC(3) * COST                                                1030
C     LORENTZ TRANSFORMATION  ALONG THE X AXIS                              1040
      DO 210  K = 1, NPROD                                                  1050
      I = NPROD-K+1                                                     2/69
      PP(1, K) = GAMMA * ( P(I,1) + BETA * P(I,4) )                     2/69
      PP(2, K) = P(I,2)                                                 2/69
      PP(3, K) = P(I,3)                                                 2/69
      ETBL(K) = GAMMA * ( P(I,4) + BETA * P(I,1) )                      2/69
      PTBL(K) = SQRT( ETBL(K) ** 2 - FMASS(K) ** 2)                         1100
  210 CONTINUE                                                              1110
      ROT (1,1) =-BETCP * SINP + ALFCP * COST * COSP                        1120
      ROT (1, 2) = - ALFCP * SINP - BETCP * COST * COSP                     1130
      ROT(1, 3) = SINT * COSP                                               1140
      ROT (2, 1) = + BETCP * COSP + ALFCP * COST * SINP                     1150
      ROT (2, 2) = ALFCP * COSP - BETCP * COST * SINP                       1160
      ROT (2,3) = SINT * SINP                                               1170
      ROT (3,1) =-ALFCP * SINT                                              1180
      ROT (3,2) = BETCP * SINT                                              1190
      ROT (3,3) = COST                                                      1200
      GO TO 220                                                             1210
 201       CONTINUE                                                         1220
C     LONGITUDINAL POLARIZATION * PROD ANGULAR DIST * NO ASSYMMETRY         1230
C     LORENTZ TRANSFORM ALONG THE Z AXIS                                    1240
      DO 200  K = 1, NPROD                                                  1250
      I = NPROD-K+1                                                     2/69
      PP(1, K) = P(I,1)                                                 2/69
      PP(2, K) = P(I,2)                                                 2/69
      PP(3, K) = GAMMA * ( P(I,3) + BETA * P(I,4) )                     2/69
      ETBL(K) = GAMMA * ( P(I,4) + BETA * P(I,3) )                      2/69
      PTBL(K) = SQRT( ETBL(K) ** 2 - FMASS(K) ** 2)                         1300
  200 CONTINUE                                                              1310
      COST = DCSINC(3)                                                      1320
      SINT = SQRT( 1.0 - COST ** 2)                                         1330
      IF (SINT)    203, 202, 203                                            1340
 202  COSP = 1.0                                                            1350
      SINP = 0.0                                                            1360
      GO TO 204                                                             1370
 203  COSP = DCSINC(1) / SINT                                               1380
      SINP = DCSINC(2) / SINT                                               1390
 204  CONTINUE                                                              1400
C     NEXT ROTATE Z AXIS TO COINCIDE WITH DIRINC                            1410
C     SET UP ROTATION MATRIX   ROT                                          1420
      ROT(1,1) = COST * COSP                                                1430
      ROT(1,2) = -SINP                                                      1440
      ROT(1,3) = SINT * COSP                                                1450
      ROT(2,1) = COST * SINP                                                1460
      ROT(2,2) = COSP                                                       1470
      ROT(2,3) = SINT * SINP                                                1480
      ROT(3,1) = -SINT                                                      1490
      ROT(3,2) = 0.0                                                        1500
      ROT(3,3) = COST                                                       1510
C     CALCULATE MOMENTUM COMPONENTS                                         1520
 220  CONTINUE                                                              1530
      DO 300  L = 1, NPROD                                                  1540
      DO 300  M = 1, 3                                                      1550
      R(M, L) = 0.0                                                         1560
      DO 300  N = 1, 3                                                      1570
      R(M, L) = R(M, L) + ROT(M, N) * PP(N, L)                              1580
  300 CONTINUE                                                              1590
C     CALCULATE DIRECTIONS FOR PRODUCT PARTICLES                            1600
      DO 400  I = 1, NPROD                                                  1610
      IF (PTBL(I))  312, 15, 312                                            1640
  312 CONTINUE                                                              1650
      DO 400 K = 1, 3                                                       1690
      DCSO(K,I) = R(K,I) / PTBL(I)                                          1700
  400 CONTINUE                                                              1710
      CALL DVCHK (K000FX)                                                   1712
       GO TO(15,401,15),K000FX                                              1714
 401  CONTINUE                                                              1720
      RETURN                                                                1730
      END                                                                   1740
      SUBROUTINE RDECAY (N, E, X, P, FLAG, A, TMASS, BMASS)
C     THIS VERSION OF 1/69 GENERATES INVARIANT PHASE SPACE EVENTS
C     IT IS A MODIFICATION OF THE J.H.U. ROUTINES PSEVNT AND ISODK
C     REFERENCE-- D.T. GILLESPIE, DISSERTATION(1968), APPENDICES C AND D
C                 JOHNS HOPKINS UNIVERSITY (HIGH ENERGY PHYSICS GROUP)
C   *    *    *    *    *    *    *    *    *    *    *    *    *    *
      DIMENSION P(8,4), Q(8,4), X(8), Z(8)                              1/69
      DIMENSION ZL(8), POWER(7), A(10)                                  1/69
      DIMENSION NSAVE(5), XSAVE(8,5), ESAVE(5), WMSAVE(5)               1/69
      DOUBLE PRECISION P10(3),BETA(3),TLOR(4,4),GAMMA,DELTA,PJ3
      PS(A,B,C) = ((A+B+C)*(A-B-C)*(A+B-C)*(A-B+C))/(A*A)
C   *    *    *    *    *    *    *    *    *    *    *    *    *    *
      IF (N .GT. 0)  GO TO 10                                           1/69
      IMAX = 0                                                          1/69
      DO 4 I = 1,5                                                      1/69
    4 NSAVE(I) = 0                                                      1/69
      GO TO 400                                                         1/69
C
C     INITIALIZE.
   10 FLAG = 0.                                                         1/69
      NP1 = N + 1
      NM1 = N - 1
      NM2 = N - 2
C     DETERMINE DIST TYPE--ITYP = 0 (UNIFORM), 1 (E**AT), 2 (UNIFORM    1/69
C     WITHIN LIMITS), 3 (COSINE SERIES)                                 1/69
      ITYP = 0                                                          1/69
      IF (A(1) .NE. 0)  ITYP = 1                                        1/69
      IF (A(10) .NE. 0)  ITYP = ITYP+2                                  1/69
      DO 14 I = 2,9                                                     1/69
   14 IF (A(I) .NE. 0)  ITYP = 3                                        1/69
      DO 20 I = 1,3                                                     1/69
   20 Q(N,I) = 0.                                                       1/69
      Q(N,4) = E                                                        1/69
      Z(N) = E
      Z(1) = X(1)
      ZL(1) = X(1)
      DO 22 I=2,N
   22 ZL(I) = ZL(I-1) + X(I)
      T = Z(N) - ZL(N)                                                  1/69
      IF (T .GT. 0.) GO TO 26                                           1/69
      FLAG = -1.                                                        1/69
      GO TO 400                                                         1/69
   26 IF (NM2 .EQ. 0) GO TO 300                                         1/69
      DO 28 I=2,NM1
      FIM1 = I - 1
   28 POWER(I) = 1./FIM1
C     FIND WMAXS.
C     FIRST SEE IF IT IS BEING SAVED                                    1/69
      DO 86 I = 1,IMAX                                                  1/69
      IF (NSAVE(I) .NE. N)  GO TO 86                                    1/69
      DO 82 J = 1,N                                                     1/69
      IF (XSAVE(J,I) .NE. X(J))  GO TO 86                               1/69
   82 CONTINUE                                                          1/69
      IF (ESAVE(I) .LT. E)  GO TO 84                                    1/69
      WMAXS = WMSAVE(I)                                                 1/69
      GO TO 200                                                         1/69
   84 ISAV = I                                                          1/69
      GO TO 100                                                         1/69
   86 CONTINUE                                                          1/69
      ISAV = IMAX+1                                                     1/69
      IF (ISAV .LE. 5) IMAX = ISAV                                      1/69
C     CALCULATE WMAXS                                                   1/69
  100 DO 104 J=1,NM2
      I = N - J
  104 Z(I) = ZL(I) + (Z(I+1) - ZL(I+1)) * .5**POWER(I)
      WMAXS = 1.                                                        2/69
      DO 106 K=2,N                                                      2/69
  106 WMAXS = WMAXS*PS(Z(K),Z(K-1),X(K))
      DZ = T/3.                                                         1/69
  108 DZ = DZ/5.
      IF (DZ .LT. T/10000.) GO TO 122                                   1/69
  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 = 1.                                                           2/69
      DO 113 K=2,N                                                      2/69
  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 = 1.                                                           2/69
      DO 117 K=2,N                                                      2/69
  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
C     SAVE WMAXS                                                        1/69
      IF (ISAV .GT. 5)  GO TO 200                                       1/69
      NSAVE(ISAV) = N                                                   1/69
      DO 124 J = 1,N                                                    1/69
  124 XSAVE(J,ISAV) = X(J)                                              1/69
      ESAVE(ISAV) = E                                                   1/69
      WMSAVE(ISAV) = WMAXS                                              1/69
C     END INITIALIZATION
C
C     GENERATE EVENTS.  FIRST GENERATE INTERMEDIATE MASSES.
  200 DO 206 J=1,NM2                                                    1/69
      IF (J.NE.1 .OR. ITYP.NE.1)  GO TO 204                             1/69
      CALL EXPDIS (A, Z(N-1), E, TMASS, BMASS, X(N), ZL(N-1), N)        2/69
      GO TO 206                                                         1/69
  204 I = N - J                                                         1/69
      CALL RANDOM (R)                                                   1/69
      Z(I) = ZL(I) + (Z(I+1) - ZL(I+1)) * ABS(R)**POWER(I)              1/69
  206 CONTINUE                                                          1/69
      WS = 1.                                                           2/69
      DO 208 K=2,N                                                      2/69
  208 WS = WS*PS(Z(K),Z(K-1),X(K))
  210 IF (WS .LE. WMAXS) GO TO 214
C     WEIGHTING FAULT
      FLAG = 1.                                                         1/69
      GO TO 400                                                         1/69
  214 CALL RANDOM (R)                                                   1/69
      IF ((WS/WMAXS) .LT. R**2) GO TO 200
C     HAVE INTERMEDIATE MASSES.  NOW GENERATE VERTEX DECAYS.
  300 DO 303 J=1,NM1                                                    1/69
      I = NP1 - J
      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)))
      IF (J.NE.1 .OR. ITYP.EQ.0)  GO TO 301                             1/69
      CALL RANDIS (ITYP, A, CPOL, E, TMASS, BMASS, X(N), Z(N-1))        1/69
      GO TO 302                                                         1/69
  301 CALL RANDOM (CPOL)                                                1/69
  302 CALL RANDOM (R)                                                   1/69
      AZM = 3.141593*R                                                  1/69
      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 303 K = 1,4                                                    1/69
      PJ3 = TLOR(K,1)*P10(1) + TLOR(K,2)*P10(2) + TLOR(K,3)*P10(3)      1/69
      Q(IM1,K) = +PJ3 + TLOR(K,4)*E10                                   1/69
  303 P(I,K) = -PJ3 + TLOR(K,4)*E20                                     1/69
      DO 304 I=1,4
  304 P(1,I) = Q(1,I)                                                   1/69
  400 RETURN                                                            1/69
      END