GEANT321/gheisha/genxpt.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:04  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.39  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GENXPT(IPPP,NFL,AVERN)
  13 C
  14 C *** GENERATION OF X- AND PT- VALUES FOR ALL PRODUCED PARTICLES ***
  15 C *** NVE 02-MAY-1988 CERN GENEVA ***
  16 C
  17 C ORIGIN : H.FESEFELDT 11-OCT-1987
  18 C
  19 C A SIMPLE SINGLE VARIABLE DESCRIPTION E D3S/DP3= F(Q) WITH
  20 C Q**2 = (M*X)**2 + PT**2 IS USED. FINAL STATE KINEMATIC IS PRODUCED
  21 C BY AN FF-TYPE ITERATIVE CASCADE METHOD
  22 C
  23 #include "geant321/s_defcom.inc"
  24 #include "geant321/s_genio.inc"
  25 C
  26       REAL ATNOL,ZNOL
  27       REAL MASPAR,LAMB,NUCSUP
  28       DIMENSION MASPAR(8),BP(8),PTEX(8),C1PAR(5),G1PAR(5),TAVAI(2),
  29      $          SIDE(MXGKCU),IAVAI(2),BINL(20),DNDL(20),TWSUP(8),
  30      $          NUCSUP(6),PSUP(6),IPAX(100)
  31       DIMENSION RNDM(3)
  32       DATA MASPAR/0.75,0.70,0.65,0.60,0.50,0.40,0.75,0.20/
  33       DATA     BP/3.50,3.50,3.50,6.00,5.00,4.00,3.50,3.50/
  34       DATA   PTEX/1.70,1.70,1.50,1.70,1.40,1.20,1.70,1.20/
  35       DATA  C1PAR/0.6,0.6,0.35,0.15,0.10/
  36       DATA  G1PAR/2.6,2.6,1.80,1.30,1.20/
  37       DATA BINL/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.11,1.25
  38      $         ,1.43,1.67,2.0,2.5,3.33,5.00,10.00/
  39       DATA TWSUP/1.,1.,0.7,0.5,0.3,0.2,0.1,0.0/
  40       DATA NUCSUP/1.00,0.7,0.5,0.4,0.35,0.3/
  41       DATA   PSUP/3.,6.,20.,50.,100.,1000./
  42 C
  43 C**
  44       CALL HIGSEL(ISEL)
  45       IF(ISEL.EQ.1) THEN
  46          CALL HIGXPT(IPPP,NFL,AVERN)
  47          RETURN
  48       ENDIF
  49 C**
  50 C**  FOR ANNIHILATION INTERACTIONS INTRODUCE PROPER KINEMATICS
  51 C**
  52       CALL CORANH(NIHIL,NFL)
  53 C**
  54 C**
  55 C** CHECK FIRST MASS-INDICES
  56 C**
  57       EK=ENP(5)
  58       EN=ENP(6)
  59       P=ENP(7)
  60       S=ENP(8)
  61       RS=ENP(9)
  62       NT=0
  63       DO 1 I=1,100
  64       IF(IPA(I).EQ.0) GOTO 1
  65       NT=NT+1
  66       IPA(NT)=IPA(I)
  67     1 CONTINUE
  68       CALL VZERO(IPA(NT+1),MXGKCU-NT)
  69       CALL UCOPY(IPA(1),IPAX(1),100)
  70 C**
  71 C** FOR LOW MULTIPLICITY USE TWO-BODY RESONANCE MODEL OR SINGLE/DOUBLE
  72 C** DIFFRACTION MODEL (--> TWOCLU (--> TWOB (--> COSCAT)))
  73 C**
  74       CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
  75       IF(NIHIL.GT.0) GOTO 200
  76       IF(NT.GE.8) GOTO 200
  77       IF(EK.LT.1.) GOTO 60
  78       CALL GRNDM(RNDM,1)
  79       RAN=RNDM(1)
  80       IF(IPART.GE.10.AND.IPART.LE.13.AND.RAN.LT.0.5) GOTO 200
  81       CALL GRNDM(RNDM,1)
  82       RAN=RNDM(1)
  83       WSUP=TWSUP(NT)
  84       IF(RAN.GT.WSUP) GOTO 200
  85    60 CALL UCOPY(IPAX,IPA,100)
  86       CALL TWOCLU(IPPP,NFL,AVERN)
  87       GO TO 9999
  88 C**
  89 C** SET EFFECTIVE 4-MOMENTUM OF PRIMARY PARTICLE
  90 C**
  91   200 MX =MXGKPV-20
  92       MX1=MX+1
  93       MX2=MX+2
  94       MX3=MX+3
  95       MX4=MX+4
  96       MX5=MX+5
  97       MX6=MX+6
  98       MX7=MX+7
  99       MX8=MX+8
 100       MX9=MX+9
 101       PV( 1,MXGKPV-1)=P*PX
 102       PV( 2,MXGKPV-1)=P*PY
 103       PV( 3,MXGKPV-1)=P*PZ
 104       PV( 4,MXGKPV-1)=EN
 105       PV( 5,MXGKPV-1)=AMAS
 106       PV( 6,MXGKPV-1)=NCH
 107       PV( 7,MXGKPV-1)=TOF
 108       PV( 8,MXGKPV-1)=IPART
 109       PV( 9,MXGKPV-1)=0.
 110       PV(10,MXGKPV-1)=USERW
 111       IER(49)=IER(49)+1
 112 C**
 113 C** SOME RANDOMISATION OF ORDER OF FINAL STATE PARTICLES
 114 C**
 115       DO 201 I=3,NT
 116       CALL GRNDM(RNDM,1)
 117       IPX=IFIX(3.+RNDM(1)*(NT-2.))
 118       IF(IPX.GT.NT) IPX=NT
 119       IPA1=IPA(IPX)
 120       IPA(IPX)=IPA(I)
 121   201 IPA(I)  =IPA1
 122 C**
 123 C** DISTRIBUTE IN FORWARD AND BACKWARD HEMISPHERE IN CMS
 124 C**
 125       SIDE(1)= 1.
 126       SIDE(2)=-1.
 127       NTB=1
 128       TARG=0.
 129       IF(IPART.LT.10.OR.IPART.GT.13) GOTO 53
 130       CALL GRNDM(RNDM,1)
 131       IF(RNDM(1).LT.0.7) GOTO 53
 132       IPA1=IPA(1)
 133       IPA(1)=IPA(2)
 134       IPA(2)=IPA1
 135    53 LEAD=0
 136       IF(IPART.LT.10.OR.IPART.EQ.14.OR.IPART.EQ.16) GOTO 532
 137       IPA1=ABS(IPA(1))
 138       IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 531
 139       LEAD=IPA1
 140       GOTO 532
 141   531 IPA1=ABS(IPA(2))
 142       IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 532
 143       LEAD=IPA1
 144   532 DO 3 I=1,NT
 145       IF(I.LE.2) GOTO 54
 146       SIDE(I)= 1.
 147       CALL GRNDM(RNDM,1)
 148       IF(RNDM(1).LT.0.5) SIDE(I)=-1.
 149       IF(SIDE(I).LT.-0.5) NTB=NTB+1
 150    54 CONTINUE
 151     3 CONTINUE
 152       TB=2.*NTB
 153       CALL GRNDM(RNDM,1)
 154       IF(RS.LT.(2.0+RNDM(1))) TB=(2.*NTB+NT)/2.
 155       ZNOL = ZNO2
 156       IF(NFL .EQ. 1)ZNOL = ZNOL -1
 157       ATNOL = ATNO2 - 1
 158 C**
 159 C** ADD PARTICLES FROM INTRANUCLEAR CASCADE
 160 C**
 161       AFC=0.312+0.200*LOG(LOG(S))+S**1.5/6000.
 162       IF(AFC.GT.0.75) AFC=0.75
 163       XTARG=AFC*(ATNO2**0.33 -1.0)*TB
 164       IF(XTARG.LE.0.) XTARG=0.01
 165       CALL POISSO(XTARG,NTARG)
 166       NT2=NT+NTARG
 167       IF(NT2.LE.MX) GOTO 2
 168       NT2=MX
 169       NTARG=NT2-NT
 170     2 CONTINUE
 171       IF (NPRT(4)) WRITE(NEWBCD,3001) NTARG,NT
 172       NT1=NT+1
 173       IF(NTARG.EQ.0) GOTO 51
 174 C**
 175 C** CHECK NUMBER OF EXTRA NUCLEONS AND PIONS
 176 C**
 177       DO 881 IPX=1,6
 178       IF(P.LE.PSUP(IPX)) GOTO 882
 179   881 CONTINUE
 180       IPX=6
 181   882 DO 4 I=NT1,NT2
 182       IF(ATNOL .GT. 0.99)THEN
 183       CALL GRNDM(RNDM,1)
 184       RAN=RNDM(1)
 185       IF(RAN.LT.NUCSUP(IPX)) GOTO 52
 186       ENDIF
 187       CALL GRNDM(RNDM,1)
 188       IPA(I)=-(7+IFIX(RNDM(1)*3.0))
 189       GOTO 4
 190    52 IPA(I)=-16
 191       PNRAT=1.-ZNOL/ATNOL
 192       CALL GRNDM(RNDM,1)
 193       IF(RNDM(1).GT.PNRAT)THEN
 194       IPA(I)=-14
 195       ZNOL = ZNOL -1
 196       ENDIF
 197       TARG=TARG+1.
 198       ATNOL = ATNOL -1
 199     4 SIDE(I)=-2.
 200       NT=NT2
 201 C**
 202 C** CHOOSE MASSES AND CHARGES FOR ALL PARTICLES
 203 C**
 204    51 DO 5 I=1,NT
 205       IPA1=ABS(IPA(I))
 206       PV(5,I)=RMASS(IPA1)
 207       PV(6,I)=RCHARG(IPA1)
 208       PV(7,I)=1.
 209       IF(PV(5,I).LT.0.) PV(7,I)=-1.
 210       PV(5,I)=ABS(PV(5,I))
 211     5 CONTINUE
 212 C**
 213 C** CHECK AVAILABLE KINETIC ENERGY, IN THIS MODEL CONSERVATION OF
 214 C** KINETIC ENERGY IN FORWARD AND BACKWARD HEMISPHERE IS ASSUMED
 215 C**
 216     6 IF(NT.LE.1) GOTO 60
 217       TAVAI(1)=RS/2.
 218       TAVAI(2)=(TARG+1.)*RS/2.
 219       IAVAI(1)=0
 220       IAVAI(2)=0
 221       DO 7 I=1,NT
 222       L=1
 223       IF(SIDE(I).LT.0.) L=2
 224       IAVAI(L)=IAVAI(L)+1
 225       TAVAI(L)=TAVAI(L)-ABS(PV(5,I))
 226     7 CONTINUE
 227       NTH=NT
 228       IF(NTH.GT.10) NTH=10
 229       IF (NPRT(4))
 230      $ WRITE(NEWBCD,3002) TAVAI,IAVAI,(IPA(I),SIDE(I),I=1,NTH)
 231       IF(IAVAI(1).LE.0) GOTO 60
 232       IF(IAVAI(2).LE.0) GOTO 60
 233       IF(TAVAI(1).GT.0.) GOTO 11
 234       CALL GRNDM(RNDM,1)
 235       ISKIP=IFIX(RNDM(1)*(IAVAI(1)-1))+1
 236       IS=0
 237       DO 10  I=1,NT
 238       II=NT-I+1
 239       IF(SIDE(II).LT.0.) GOTO 10
 240       IS=IS+1
 241       IF(IS.NE.ISKIP) GOTO 10
 242       IF(II.EQ.NT) GOTO 9
 243       NT1=II+1
 244       NT2=NT
 245       DO 8 J=NT1,NT2
 246       IPA(J-1)=IPA(J)
 247       SIDE(J-1)=SIDE(J)
 248       DO 71 K=1,10
 249    71 PV(K,J-1)=PV(K,J)
 250     8 CONTINUE
 251       GOTO 9
 252    10 CONTINUE
 253     9 IPA(NT)=0
 254       SIDE(NT)=0.
 255       NT=NT-1
 256       GOTO 6
 257    11 IF(TAVAI(2).GT.0.) GOTO 15
 258       CALL GRNDM(RNDM,1)
 259       ISKIP=IFIX(RNDM(1)*(IAVAI(2)-1))+1
 260       IS=0
 261       DO 14  I=1,NT
 262       II=NT-I+1
 263       IF(SIDE(II).GT.0.) GOTO 14
 264       IS=IS+1
 265       IF(IS.NE.ISKIP) GOTO 14
 266       IF(SIDE(II).LT.-1.5) NTARG=NTARG-1
 267       IF(NTARG.LT.0) NTARG=0
 268       IF(II.EQ.NT) GOTO 13
 269       NT1=II+1
 270       NT2=NT
 271       DO 12 J=NT1,NT2
 272       IPA(J-1)=IPA(J)
 273       SIDE(J-1)=SIDE(J)
 274       DO 74 K=1,10
 275    74 PV(K,J-1)=PV(K,J)
 276    12 CONTINUE
 277       GOTO 13
 278    14 CONTINUE
 279    13 IPA(NT)=0
 280       SIDE(NT)=0.
 281       NT=NT-1
 282       GOTO 6
 283    15 IF(NT.LE.1) GOTO 60
 284       IF(NT.EQ.MX) GOTO 29
 285       NT1=NT+1
 286       NT2=MX
 287       DO 28 I=NT1,NT2
 288    28 IPA(I)=0
 289    29 CONTINUE
 290 C**
 291 C** NOW THE PREPARATION IS FINISHED.
 292 C** DEFINE INITIAL STATE VECTORS FOR LORENTZ TRANSFORMATIONS.
 293 C**
 294       PV( 1,MX1)=0.
 295       PV( 2,MX1)=0.
 296       PV( 3,MX1)=P
 297       PV( 4,MX1)=SQRT(P*P+AMAS*AMAS)
 298       PV( 5,MX1)=ABS(AMAS)
 299       PV( 1,MX2)=0.
 300       PV( 2,MX2)=0.
 301       PV( 3,MX2)=0.
 302       PV( 4,MX2)=MP
 303       PV( 5,MX2)=MP
 304       PV( 1,MX4)=0.
 305       PV( 2,MX4)=0.
 306       PV( 3,MX4)=0.
 307       PV( 4,MX4)=MP*(1.+TARG)
 308       PV( 5,MX4)=PV(4,MX4)
 309       PV( 1,MX8)=0.
 310       PV( 2,MX8)=0.
 311       PV( 3,MX8)=0.
 312       PV( 1,MX9)=1.
 313       PV( 2,MX9)=0.
 314       PV( 3,MX9)=0.
 315       CALL ADD(MX1,MX2,MX3)
 316       CALL ADD(MX4,MX1,MX4)
 317       CALL LOR(MX1,MX3,MX1)
 318       CALL LOR(MX2,MX3,MX2)
 319 C**
 320 C** MAIN LOOP FOR 4-MOMENTUM GENERATION , SEE PITHA-REPORT (AACHEN)
 321 C** FOR A DETAILED DESCRIPTION OF THE METHOD.
 322 C**
 323       CALL GRNDM(RNDM,1)
 324       PHI=RNDM(1)*TWPI
 325       EKIN1=0.
 326       EKIN2=0.
 327       DO 39 J=1,10
 328       PV(J,MX5)=0.
 329    39 PV(J,MX6)=0.
 330       NPG=0
 331       TARG1=0.
 332       DO 16 III=1,NT
 333       I=NT-III+1
 334       IPA1=ABS(IPA(I))
 335 C**
 336 C** COUNT NUMBER OF BACKWARD NUCLEONS
 337 C**
 338       IF(I.EQ.2) GOTO 301
 339       IF(SIDE(I).LT.-1.5.AND.IPA1.GE.14) GOTO 301
 340       GOTO 38
 341   301 NPG=NPG+1
 342       IF(NPG.GT.18) GOTO 38
 343       SIDE(I)=-3.
 344       TARG1=TARG1+1.
 345       GOTO 16
 346    38 J=3
 347       IF(IPA1.LT.14) J=2
 348       IF(IPA1.LT.10) J=1
 349       IF(I.LE.2) J=J+3
 350       IF(SIDE(I).LT.-1.5) J=7
 351       IF(J.EQ.7.AND.IPA1.GE.14) J=8
 352 C**
 353 C** SET PT - AND PHI VALUES, THEY ARE CHANGED SOMEWHAT IN THE ITERATION
 354 C** LOOP, SET MASS PARAMETER FOR LAMBDA FRAGMENTATION MODEL
 355 C**
 356       CALL GRNDM(RNDM,1)
 357       RAN=RNDM(1)
 358       BPP=BP(J)
 359       BPE=PTEX(J)
 360       PT2=-LOG(1.-RAN)/BPP
 361       ASPAR=MASPAR(J)
 362       PT2=PT2**BPE
 363       PT =SQRT(PT2)
 364       IF(PT.LT.0.001) PT=0.001
 365       PV(1,I)=PT*COS(PHI)
 366       PV(2,I)=PT*SIN(PHI)
 367       PV(10,I)=PT
 368       BINL(1)=0.
 369       RLMAX=1./PV(10,I)
 370       DO 73 J=2,20
 371    73 BINL(J)=RLMAX*(J-1)/19.
 372       ET=PV(4,MX1)
 373       IF(SIDE(I).LT.0.) THEN
 374          ET=PV(4,MX2)
 375       ENDIF
 376       DNDL(1)=0.
 377       NTRIAL=0
 378 C**
 379 C** START OF BIG ITERATION LOOP
 380 C**
 381    30 NTRIAL=NTRIAL+1
 382       IF(NTRIAL.GT. 2) GOTO 169
 383       DO 17 L=2,20
 384       DNDL(L)=0.
 385       X=(BINL(L)+BINL(L-1))/2.
 386       IF(PV(10,I).LT.0.001) PV(10,I)=0.001
 387       IF(X.GT.1./PV(10,I)) GOTO 17
 388       DX=BINL(L)-BINL(L-1)
 389       DNDL(L)=ASPAR/SQRT((1.+(ASPAR*X)**2)**3)
 390       DNDL(L)=ET*DNDL(L)/SQRT((X*PV(10,I)*ET)**2+PV(10,I)**2
 391      *                             +PV(5,I)**2)
 392       DNDL(L)=DNDL(L)*DX
 393    17 DNDL(L)=DNDL(L-1)+DNDL(L)
 394       NTRI=0
 395    31 CALL GRNDM(RNDM,1)
 396       RAN=RNDM(1)*DNDL(20)
 397       DO 18 L=2,20
 398       IF(RAN.LT.DNDL(L)) GOTO 19
 399    18 CONTINUE
 400 C**
 401 C** START OF SMALL ITERATION LOOP
 402 C**
 403    19 NTRI=NTRI+1
 404       CALL GRNDM(RNDM,1)
 405       RAN=RNDM(1)
 406       DX=BINL(L)-BINL(L-1)
 407       LAMB=BINL(L-1)+RAN*DX/2.
 408       X=PV(10,I)*LAMB
 409       IF(X.GT.1.) X=1.
 410       X=X*SIDE(I)/ABS(SIDE(I))
 411       PV(3,I)=X*ET
 412       PV(4,I)=PV(3,I)**2+PV(10,I)**2+PV(5,I)**2
 413       PV(4,I)=SQRT(PV(4,I))
 414       IF(SIDE(I).LT.0.) GOTO 165
 415       IF(I.GT.2) GOTO 20
 416       EKIN=TAVAI(1)-EKIN1
 417       CALL NORMAL(RAN)
 418       IF(EKIN.LT.0.) EKIN=0.04*ABS(RAN)
 419       PV(4,I)=ABS(PV(5,I))+EKIN
 420       RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
 421       PP=SQRT(RNVE)
 422       CALL LENGTX(I,PP1)
 423 C
 424       IF (PP1 .GE. 1.0E-6) GO TO 8000
 425       CALL GRNDM(RNDM,2)
 426       RTHNVE=PI*RNDM(1)
 427       PHINVE=TWPI*RNDM(2)
 428       PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE)
 429       PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE)
 430       PV(3,I)=PP*COS(RTHNVE)
 431       GO TO 8001
 432  8000 CONTINUE
 433       PV(1,I)=PV(1,I)*PP/PP1
 434       PV(2,I)=PV(2,I)*PP/PP1
 435       PV(3,I)=PV(3,I)*PP/PP1
 436  8001 CONTINUE
 437 C
 438       CALL ADD(MX5,I,MX5)
 439       GOTO 16
 440    20 EKIN=EKIN1+PV(4,I)-ABS(PV(5,I))
 441       IF(EKIN.LT.0.95*TAVAI(1)) GOTO 161
 442       IF(NTRI.GT. 5) GOTO 167
 443       PV(10,I)=PV(10,I)*0.9
 444       PV( 1,I)=PV( 1,I)*0.9
 445       PV( 2,I)=PV( 2,I)*0.9
 446       DNDL(20)=DNDL(20)*0.9
 447       IF((TAVAI(2)-ABS(PV(5,I))).LT.0.) GOTO 31
 448       SIDE(I)=-1.
 449       TAVAI(1)=TAVAI(1)+ABS(PV(5,I))
 450       TAVAI(2)=TAVAI(2)-ABS(PV(5,I))
 451       GOTO 31
 452   161 CALL ADD(MX5,I,MX5)
 453       EKIN1=EKIN1+PV(4,I)-ABS(PV(5,I))
 454       GOTO 163
 455   165 EKIN=EKIN2+PV(4,I)-ABS(PV(5,I))
 456       XXX=0.95+0.05*TARG/20.
 457       IF(XXX.GT.0.999) X=0.999
 458       IF(EKIN.LT.XXX*TAVAI(2)) GOTO 166
 459       IF(NTRI.GT. 5) GOTO 167
 460       PV(10,I)=PV(10,I)*0.9
 461       PV( 1,I)=PV( 1,I)*0.9
 462       PV( 2,I)=PV( 2,I)*0.9
 463       DNDL(20)=DNDL(20)*0.9
 464       IF((TAVAI(1)-ABS(PV(5,I))).LT.0.) GOTO 31
 465       SIDE(I)=+1.
 466       TAVAI(1)=TAVAI(1)-ABS(PV(5,I))
 467       TAVAI(2)=TAVAI(2)+ABS(PV(5,I))
 468       GOTO 31
 469   166 CALL ADD(MX6,I,MX6)
 470       EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
 471   163 CALL ADD(MX5,MX6,MX7)
 472       PV(3,MX7)=0.
 473       CALL ANG(MX7,MX9,COST,PHIS)
 474       IF(PV(2,MX7).LT.0.) PHIS=TWPI-PHIS
 475       CALL NORMAL(RAN)
 476       RAN=RAN*PI/12.
 477       PHI=PHIS+PI+RAN
 478       IF(PHI.GT.TWPI) PHI=PHI-TWPI
 479       IF(PHI.LT.0.) PHI=TWPI-PHI
 480       GOTO 16
 481 C**
 482 C** PARTICLE MOMENTUM ZERO, REDUCE KINETIC ENERGY OF ALL OTHER
 483 C**
 484   167 EKIN1=0.
 485       EKIN2=0.
 486       DO 162 J=1,10
 487       PV(J,MX5)=0.
 488   162 PV(J,MX6)=0.
 489       II=I+1
 490       DO 168 L=II,NT
 491       IF(ABS(IPA(L)).GE.14.AND.SIDE(L).LT.0.) GOTO 168
 492       PV(4,L)=PV(4,L)*0.95+0.05*ABS(PV(5,L))
 493       IF(PV(4,L).LT.ABS(PV(5,L))) PV(4,L)=ABS(PV(5,L))
 494       RNVE=ABS(PV(4,L)**2-PV(5,L)**2)
 495       PP=SQRT(RNVE)
 496       CALL LENGTX(L,PP1)
 497 C
 498       IF (PP1 .GE. 1.0E-6) GO TO 8002
 499       CALL GRNDM(RNDM,2)
 500       RTHNVE=PI*RNDM(1)
 501       PHINVE=TWPI*RNDM(2)
 502       PV(1,L)=PP*SIN(RTHNVE)*COS(PHINVE)
 503       PV(2,L)=PP*SIN(RTHNVE)*SIN(PHINVE)
 504       PV(3,L)=PP*COS(RTHNVE)
 505       GO TO 8003
 506  8002 CONTINUE
 507       PV(1,L)=PV(1,L)*PP/PP1
 508       PV(2,L)=PV(2,L)*PP/PP1
 509       PV(3,L)=PV(3,L)*PP/PP1
 510  8003 CONTINUE
 511 C
 512       PV(10,L)=SQRT(PV(1,L)**2+PV(2,L)**2)
 513       IF(SIDE(L).LT.0.) GOTO 164
 514       EKIN1=EKIN1+PV(4,L)-ABS(PV(5,L))
 515       CALL ADD(MX5,L,MX5)
 516       GOTO 168
 517   164 EKIN2=EKIN2+PV(4,L)-ABS(PV(5,L))
 518       CALL ADD(MX6,L,MX6)
 519   168 CONTINUE
 520 C *** NEXT STMT. CHANGED TO PREVENT FROM INFINITE LOOPING ***
 521 C*************      GOTO 38
 522       GO TO 30
 523 C**
 524 C** SKIP PARTICLE, IF NOT ENOUGH ENERGY
 525 C**
 526   169 IPA(I)=0
 527       DO 170 J=1,10
 528   170 PV(J,I)=0.
 529       GOTO 163
 530    16 CONTINUE
 531       NTRI=0
 532       II=0
 533       DO 320 I=1,NT
 534       IF(IPA(I).EQ.0) GOTO 320
 535       II=II+1
 536       IPA(II)=IPA(I)
 537       SIDE(II)=SIDE(I)
 538       DO 321 J=1,10
 539   321 PV(J,II)=PV(J,I)
 540   320 CONTINUE
 541       NT=II
 542 C**
 543 C** BACKWARD NUCLEONS PRODUCED WITH A CLUSTER MODEL
 544 C**
 545       CALL LOR(MX4,MX3,MX7)
 546       CALL SUB(MX7,MX5,MX7)
 547       CALL SUB(MX7,MX6,MX7)
 548       IF(TARG1.GT.1.5) GOTO 310
 549   322 I=2
 550       CALL NORMAL(RAN)
 551       EKIN=TAVAI(2)-EKIN2
 552       EKINM=RS/2.-MP
 553       IF(EKIN.GT.EKINM) EKIN=EKINM
 554       CALL NORMAL(RAN)
 555       IF(EKIN.LT.0.04) EKIN=0.04*ABS(RAN)
 556       PV(4,I)=ABS(PV(5,I))+EKIN
 557       RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
 558       PP=SQRT(RNVE)
 559       CALL LENGTX(MX7,PP1)
 560 C
 561       IF (PP1 .GE. 1.0E-6) GO TO 8004
 562       CALL GRNDM(RNDM,2)
 563       RTHNVE=PI*RNDM(1)
 564       PHINVE=TWPI*RNDM(2)
 565       PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE)
 566       PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE)
 567       PV(3,I)=PP*COS(RTHNVE)
 568       GO TO 8005
 569  8004 CONTINUE
 570       PV(1,I)=PV(1,MX7)*PP/PP1
 571       PV(2,I)=PV(2,MX7)*PP/PP1
 572       PV(3,I)=PV(3,MX7)*PP/PP1
 573  8005 CONTINUE
 574 C
 575       CALL ADD(MX6,I,MX6)
 576       GOTO 330
 577   310 ITARG1=IFIX(TARG1+0.1)
 578       IF(ITARG1.GT.5) ITARG1=5
 579       RMB0=0.
 580       NPG=0
 581       DO 311 I=1,NT
 582       IF(SIDE(I).GT.-2.5) GOTO 311
 583       NPG=NPG+1
 584       RMB0=RMB0+ABS(PV(5,I))
 585   311 CONTINUE
 586       IF(NPG.LT.2) GOTO 322
 587       CALL GRNDM(RNDM,1)
 588       RAN=RNDM(1)
 589       RMB=-LOG(1.-RAN)
 590       GPAR=G1PAR(ITARG1)
 591       CPAR=C1PAR(ITARG1)
 592       RMB=RMB0+RMB**CPAR/GPAR
 593       PV(5,MX7)=RMB
 594       IF(PV(5,MX7).GT.PV(4,MX7)) PV(5,MX7)=PV(4,MX7)
 595       RNVE=ABS(PV(4,MX7)**2-PV(5,MX7)**2)
 596       PP=SQRT(RNVE)
 597       CALL LENGTX(MX7,PP1)
 598 C
 599       IF (PP1 .GE. 1.0E-6) GO TO 8006
 600       CALL GRNDM(RNDM,2)
 601       RTHNVE=PI*RNDM(1)
 602       PHINVE=TWPI*RNDM(2)
 603       PV(1,MX7)=PP*SIN(RTHNVE)*COS(PHINVE)
 604       PV(2,MX7)=PP*SIN(RTHNVE)*SIN(PHINVE)
 605       PV(3,MX7)=PP*COS(RTHNVE)
 606       GO TO 8007
 607  8006 CONTINUE
 608       PV(1,MX7)=PV(1,MX7)*PP/PP1
 609       PV(2,MX7)=PV(2,MX7)*PP/PP1
 610       PV(3,MX7)=PV(3,MX7)*PP/PP1
 611  8007 CONTINUE
 612 C
 613       I=MX7
 614       IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 615       PV(1,MX7)=-PV(1,MX7)
 616       PV(2,MX7)=-PV(2,MX7)
 617       PV(3,MX7)=-PV(3,MX7)
 618       KGENEV=1
 619       TECM=PV(5,MX7)
 620       NPG=0
 621       DO 312 I=1,NT
 622       IF(SIDE(I).GT.-2.5)GOTO 312
 623       NPG=NPG+1
 624       AMASS(NPG)=ABS(PV(5,I))
 625   312 CONTINUE
 626       CALL PHASP
 627       NPG=0
 628       DO 314 I=1,NT
 629       IF(SIDE(I).GT.-2.5) GOTO 314
 630       NPG=NPG+1
 631       PV(1,I)=PCM(1,NPG)
 632       PV(2,I)=PCM(2,NPG)
 633       PV(3,I)=PCM(3,NPG)
 634       PV(4,I)=PCM(4,NPG)
 635       CALL LOR(I,MX7,I)
 636       CALL ADD(MX6,I,MX6)
 637   314 CONTINUE
 638   330 IF (NPRT(4))
 639      $ WRITE(NEWBCD,2002) NTRIAL,EKIN1,EKIN2,TAVAI(1),TAVAI(2)
 640   175 IF (.NOT.NPRT(4)) GOTO 36
 641       CALL ADD(MX5,MX6,MX7)
 642       EKIN1=PV(4,MX1)+PV(4,MX2)
 643       EKIN2=PV(4,MX5)+PV(4,MX6)
 644       WRITE(NEWBCD,2000) EKIN1,EKIN2
 645       I=MX1
 646       WRITE(NEWBCD,2001) I,(PV(J,I),J=1,4)
 647       I=MX2
 648       WRITE(NEWBCD,2001) I,(PV(J,I),J=1,4)
 649       I=MX5
 650       WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 651       I=MX6
 652       WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 653       DO 37 I=1,NT
 654    37 WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
 655 C**
 656 C** LORENTZ TRANSFORMATION IN LAB SYSTEM
 657 C**
 658    36 IF(NT.LE.2) GOTO 60
 659       TARG=0.
 660       DO 601 I=1,NT
 661       IF(PV(5,I).GT.0.5) TARG=TARG+1.
 662       CALL LOR(I,MX2,I)
 663   601 CONTINUE
 664       IF(ABS(AMAS) .GT. 0.5)TARG = TARG - 1.
 665       IF(NIHIL .GT. 0)TARG = TARG + 2
 666       IF(TARG.LT.0.5) TARG=1.
 667       IF(LEAD.EQ.0) GOTO 6085
 668       DO 6081 I=1,NT
 669       IF(ABS(IPA(I)).EQ.LEAD) GOTO 6085
 670  6081 CONTINUE
 671       I=1
 672       IF(LEAD.GE.14.AND.ABS(IPA(2)).GE.14) I=2
 673       IF(LEAD.LT.14.AND.ABS(IPA(2)).LT.14) I=2
 674       IPA(I)=LEAD
 675       EKIN=PV(4,I)-ABS(PV(5,I))
 676       PV(5,I)=RMASS(LEAD)
 677       PV(7,I)=1.
 678       IF(PV(5,I).LT.0.) PV(7,I)=-1.
 679       PV(5,I)=ABS(PV(5,I))
 680       PV(6,I)=RCHARG(LEAD)
 681       PV(4,I)=PV(5,I)+EKIN
 682       CALL LENGTX(I,PP)
 683       RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
 684       PP1=SQRT(RNVE)
 685       PV(1,I)=PP1*PV(1,I)/PP
 686       PV(2,I)=PP1*PV(2,I)/PP
 687       PV(3,I)=PP1*PV(3,I)/PP
 688  6085 KGENEV=1
 689       PV(1,MX4)=0.
 690       PV(2,MX4)=0.
 691       PV(3,MX4)=P
 692       PV(4,MX4)=SQRT(P*P+AMAS*AMAS)
 693       PV(5,MX4)=ABS(AMAS)
 694       EKIN0=PV(4,MX4)-PV(5,MX4)
 695       PV(1,MX5)=0.
 696       PV(2,MX5)=0.
 697       PV(3,MX5)=0.
 698       PV(4,MX5)=MP*TARG
 699       PV(5,MX5)=PV(4,MX5)
 700       EKIN=PV(4,MX4)+PV(4,MX5)
 701       I=MX4
 702       IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 703       I=MX5
 704       IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 705       CALL ADD(MX4,MX5,MX6)
 706       CALL LOR(MX4,MX6,MX4)
 707       CALL LOR(MX5,MX6,MX5)
 708       TECM=PV(4,MX4)+PV(4,MX5)
 709       NPG=NT
 710       PV(1,MX8)=0.
 711       PV(2,MX8)=0.
 712       PV(3,MX8)=0.
 713       PV(4,MX8)=0.
 714       PV(5,MX8)=0.
 715       EKIN1=0.
 716       DO 598 I=1,NPG
 717       IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
 718       CALL ADD(MX8,I,MX8)
 719       EKIN1=EKIN1+PV(4,I)-PV(5,I)
 720       EKIN=EKIN-PV(5,I)
 721       IF(I.GT.18) GOTO 598
 722       AMASS(I)=PV(5,I)
 723   598 CONTINUE
 724       IF(NPG.GT.18) GOTO 597
 725       CALL PHASP
 726       EKIN=0.
 727       DO 599 I=1,NPG
 728       PV(1,MX7)=PCM(1,I)
 729       PV(2,MX7)=PCM(2,I)
 730       PV(3,MX7)=PCM(3,I)
 731       PV(4,MX7)=PCM(4,I)
 732       PV(5,MX7)=AMASS(I)
 733       CALL LOR(MX7,MX5,MX7)
 734   599 EKIN=EKIN+PV(4,MX7)-PV(5,MX7)
 735       CALL ANG(MX8,MX4,COST,TETA)
 736       IF (NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1,EKIN
 737 C**
 738 C** MAKE SHURE, THAT  KINETIC ENERGIES ARE CORRECT.
 739 C** EKIN= KINETIC ENERGY THEORETICALLY
 740 C** EKIN1= KINETIC ENERGY SIMULATED
 741 C**
 742   597 IF(EKIN1.EQ.0.) GOTO 600
 743       PV(1,MX7)=0.
 744       PV(2,MX7)=0.
 745       PV(3,MX7)=0.
 746       PV(4,MX7)=0.
 747       PV(5,MX7)=0.
 748       WGT=EKIN/EKIN1
 749       EKIN1=0.
 750       DO 602 I=1,NT
 751       EKIN=PV(4,I)-PV(5,I)
 752       EKIN=EKIN*WGT
 753       PV(4,I)=EKIN+PV(5,I)
 754       RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
 755       PP=SQRT(RNVE)
 756       CALL LENGTX(I,PP1)
 757 C
 758       IF (PP1 .GE. 1.0E-6) GO TO 8008
 759       CALL GRNDM(RNDM,2)
 760       RTHNVE=PI*RNDM(1)
 761       PHINVE=TWPI*RNDM(2)
 762       PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE)
 763       PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE)
 764       PV(3,I)=PP*COS(RTHNVE)
 765       GO TO 8009
 766  8008 CONTINUE
 767       PV(1,I)=PV(1,I)*PP/PP1
 768       PV(2,I)=PV(2,I)*PP/PP1
 769       PV(3,I)=PV(3,I)*PP/PP1
 770  8009 CONTINUE
 771 C
 772       EKIN1=EKIN1+EKIN
 773       CALL ADD(MX7,I,MX7)
 774   602 CONTINUE
 775       CALL ANG(MX7,MX4,COST,TETA)
 776       IF (NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1
 777 C**
 778 C** ROTATE IN DIRECTION OF Z-AXIS, THIS DOES DISTURB IN SOME WAY OUR
 779 C** INCLUSIVE DISTRIBUTIONS, BUT IT IS NESSACARY FOR MOMENTUM CONSER-
 780 C** VATION.
 781 C**
 782   600 PV(1,MX7)=0.
 783       PV(2,MX7)=0.
 784       PV(3,MX7)=0.
 785       PV(4,MX7)=0.
 786       PV(5,MX7)=0.
 787       DO 596 I=1,NT
 788   596 CALL ADD(MX7,I,MX7)
 789 C**
 790 C** SOME SMEARING IN TRANSVERSE DIRECTION FROM FERMI MOTION
 791 C**
 792 *          call rannor(ran1,ran2)
 793       CALL GRNDM(RNDM,2)
 794       RY=RNDM(1)
 795       RZ=RNDM(2)
 796       RX=6.283185*RZ
 797       A1=SQRT(-2.*LOG(RY))
 798       RAN1=A1*SIN(RX)
 799       RAN2=A1*COS(RX)
 800       PV(1,MX7)=PV(1,MX7)+RAN1*0.020*TARG
 801       PV(2,MX7)=PV(2,MX7)+RAN2*0.020*TARG
 802       CALL DEFS(MX4,MX7,MX8)
 803       PV(1,MX7)=0.
 804       PV(2,MX7)=0.
 805       PV(3,MX7)=0.
 806       PV(4,MX7)=0.
 807       PV(5,MX7)=0.
 808       DO 595 I=1,NT
 809       CALL TRAC(I,MX8,I)
 810   595 CALL ADD(MX7,I,MX7)
 811       CALL ANG(MX7,MX4,COST,TETA)
 812       IF (NPRT(4)) WRITE(NEWBCD,2003) TETA
 813 C**
 814 C** ROTATE IN DIRECTION OF PRIMARY PARTICLE, SUBTRACT BINDING ENERGIES
 815 C** AND MAKE SOME FURTHER CORRECTIONS IF REQUIRED (STEEP, STEEQ)
 816 C**
 817       DEKIN=0.
 818       NPIONS=0
 819       EK1=0.
 820       DO 21 I=1,NT
 821       CALL DEFS1(I,MXGKPV-1,I)
 822       IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
 823       IF(ATNO2.LT.1.5) GOTO 21
 824       CALL LENGTX(I,PP)
 825       EKIN=PV(4,I)-ABS(PV(5,I))
 826       CALL NORMAL(RAN)
 827       EKIN=EKIN-CFA*(1.+0.5*RAN)
 828       IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
 829       CALL STEEQ(XXH,I)
 830       DEKIN=DEKIN+EKIN*(1.-XXH)
 831       EKIN=EKIN*XXH
 832       IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) NPIONS=NPIONS+1
 833       IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) EK1=EK1+EKIN
 834       PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
 835       PV(4,I)=EKIN+ABS(PV(5,I))
 836 C
 837       IF (PP .GE. 1.0E-6) GO TO 8010
 838       CALL GRNDM(RNDM,2)
 839       RTHNVE=PI*RNDM(1)
 840       PHINVE=TWPI*RNDM(2)
 841       PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
 842       PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
 843       PV(3,I)=PP1*COS(RTHNVE)
 844       GO TO 8011
 845  8010 CONTINUE
 846       PV(1,I)=PV(1,I)*PP1/PP
 847       PV(2,I)=PV(2,I)*PP1/PP
 848       PV(3,I)=PV(3,I)*PP1/PP
 849  8011 CONTINUE
 850 C
 851    21 CONTINUE
 852       IF(EK1.EQ.0.) GOTO 23
 853       IF(NPIONS.EQ.0) GOTO 23
 854       DEKIN=1.+DEKIN/EK1
 855       DO 22 I=1,NT
 856       IF(ABS(IPA(I)).LT.7.OR.ABS(IPA(I)).GT.9) GOTO 22
 857       CALL LENGTX(I,PP)
 858       EKIN=PV(4,I)-ABS(PV(5,I))
 859       EKIN=EKIN*DEKIN
 860       IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
 861       PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
 862       PV(4,I)=EKIN+ABS(PV(5,I))
 863 C
 864       IF (PP .GE. 1.0E-6) GO TO 8012
 865       CALL GRNDM(RNDM,2)
 866       RTHNVE=PI*RNDM(1)
 867       PHINVE=TWPI*RNDM(2)
 868       PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
 869       PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
 870       PV(3,I)=PP1*COS(RTHNVE)
 871       GO TO 8013
 872  8012 CONTINUE
 873       PV(1,I)=PV(1,I)*PP1/PP
 874       PV(2,I)=PV(2,I)*PP1/PP
 875       PV(3,I)=PV(3,I)*PP1/PP
 876  8013 CONTINUE
 877 C
 878    22 CONTINUE
 879 C**
 880 C** ADD BLACK TRACK PARTICLES, THE TOTAL NUMBER OF PARTICLES PRODUCED
 881 C** IS RESTRICTED TO 198, THIS MAY HAVE INFLUENCE ON VERY HIGH ENERGY
 882 C** FIRST PROTONS AND NEUTRONS
 883 C**
 884    23 IF(ATNO2.LT.1.5) GOTO 40
 885       CALL SELFAB(SPROB)
 886       TEX=ENP(1)
 887       SPALL=TARG
 888       IF(TEX.LT.0.001) GOTO 445
 889       BLACK=(1.5+1.25*TARG)*ENP(1)/(ENP(1)+ENP(3))
 890       CALL POISSO(BLACK,NBL)
 891       IF (NPRT(4)) WRITE(NEWBCD,3003) NBL,TEX
 892       IF(NBL.GT.ATNOL) NBL=ATNOL
 893       IF(NT+NBL.GT.MXGKPV-10) NBL=MXGKPV-10-NT
 894       IF(NBL.LE.0) GOTO 445
 895       EKIN=TEX/NBL
 896       EKIN2=0.
 897       CALL STEEP(XX)
 898       DO 441 I=1,NBL
 899       CALL GRNDM(RNDM,1)
 900       IF(RNDM(1).LT.SPROB) GOTO 441
 901       IF(NT.EQ.MXGKPV-10) GOTO 441
 902       IF(EKIN2.GT.TEX) GOTO 443
 903       CALL GRNDM(RNDM,1)
 904       RAN1=RNDM(1)
 905       CALL NORMAL(RAN2)
 906       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
 907       IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
 908       EKIN1=EKIN1*XX
 909       EKIN2=EKIN2+EKIN1
 910       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
 911       IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
 912       IPA1=16
 913       PNRAT=1.-ZNOL/ATNOL
 914       CALL GRNDM(RNDM,3)
 915       IF(RNDM(1).GT.PNRAT)THEN
 916       IPA1=14
 917       ZNOL = ZNOL -1
 918       ENDIF
 919       ATNOL = ATNOL - 1
 920       NT=NT+1
 921       SPALL=SPALL+1.
 922       COST=-1.+RNDM(2)*2.
 923       SINT=SQRT(ABS(1.-COST*COST))
 924       PHI=TWPI*RNDM(3)
 925       IPA(NT)=-IPA1
 926       SIDE(NT)=-4.
 927       PV(5,NT)=ABS(RMASS(IPA1))
 928       PV(6,NT)=RCHARG(IPA1)
 929       PV(7,NT)=1.
 930       PV(4,NT)=EKIN1+PV(5,NT)
 931       RNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
 932       PP=SQRT(RNVE)
 933       PV(1,NT)=PP*SINT*SIN(PHI)
 934       PV(2,NT)=PP*SINT*COS(PHI)
 935       PV(3,NT)=PP*COST
 936   441 CONTINUE
 937   443 IF(ATNO2.LT.10.) GOTO 445
 938       IF(EK.GT.2.0) GOTO 445
 939       II=NT+1
 940       KK=0
 941       EKA=EK
 942       IF(EKA.GT.1.) EKA=EKA*EKA
 943       IF(EKA.LT.0.1) EKA=0.1
 944       IKA=3.6*EXP((ZNO2**2/ATNO2-35.56)/6.45)/EKA
 945       IF(IKA.LE.0) GO TO 445
 946       DO 444 I=1,NT
 947       II=II-1
 948       IF(IPA(II).NE.-14) GOTO 444
 949       IPA(II)=-16
 950       IPA1  = 16
 951       ZNOL = ZNOL + 1.
 952       PV(5,II)=ABS(RMASS(IPA1))
 953       PV(6,II)=RCHARG(IPA1)
 954       KK=KK+1
 955       IF(KK.GT.IKA) GOTO 445
 956   444 CONTINUE
 957 C**
 958 C** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
 959 C**
 960   445 TEX=ENP(3)
 961       IF(TEX.LT.0.001) GOTO 40
 962       IF(ATNOL .LT. ZNOL + 1.)GOTO 40
 963       IF(ZNOL .LT. 1.)GOTO 40
 964       BLACK=(1.5+1.25*TARG)*ENP(3)/(ENP(1)+ENP(3))
 965       CALL POISSO(BLACK,NBL)
 966       IF(NT+NBL.GT.MXGKPV-10) NBL=MXGKPV-10-NT
 967       IF(NBL.LE.0) GOTO 40
 968       EKIN=TEX/NBL
 969       EKIN2=0.
 970       CALL STEEP(XX)
 971       IF (NPRT(4)) WRITE(NEWBCD,3004) NBL,TEX
 972       DO 442 I=1,NBL
 973       IF(ATNOL .LT. ZNOL + 1.)GOTO 40
 974       IF(ZNOL .LT. 1.)GOTO 40
 975       CALL GRNDM(RNDM,1)
 976       IF(RNDM(1).LT.SPROB) GOTO 442
 977       IF(NT.EQ.MXGKPV-10) GOTO 442
 978       IF(EKIN2.GT.TEX) GOTO 40
 979       CALL GRNDM(RNDM,1)
 980       RAN1=RNDM(1)
 981       CALL NORMAL(RAN2)
 982       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
 983       IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
 984       EKIN1=EKIN1*XX
 985       EKIN2=EKIN2+EKIN1
 986       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
 987       IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
 988       CALL GRNDM(RNDM,3)
 989       COST=-1.+RNDM(1)*2.
 990       SINT=SQRT(ABS(1.-COST*COST))
 991       PHI=TWPI*RNDM(2)
 992       RAN=RNDM(3)
 993       IPA(NT+1)=-30
 994       ATNOL = ATNOL - 2.
 995       ZNOL = ZNOL - 1.
 996       IF(RAN.GT.0.60)THEN
 997       IF(ATNOL .GT. ZNOL + 0.9)THEN
 998       IPA(NT+1)=-31
 999       ATNOL = ATNOL - 1.
1000       IF(RAN.GT.0.90)THEN
1001       IF( (ATNOL .GT. 0.9) .AND. (ZNOL .GT. 0.9))THEN
1002       IPA(NT+1)=-32
1003       ATNOL = ATNOL - 1.
1004       ZNOL = ZNOL - 1.
1005       ENDIF
1006       ENDIF
1007       ENDIF
1008       ENDIF
1009       SIDE(NT+1)=-4.
1010       PV(5,NT+1)=(ABS(IPA(NT+1))-28)*MP
1011       SPALL=SPALL+PV(5,NT+1)*1.066
1012       NT=NT+1
1013       PV(6,NT)=1.
1014       IF(IPA(NT).EQ.-32) PV(6,NT)=2.
1015       PV(7,NT)=1.
1016       PV(4,NT)=PV(5,NT)+EKIN1
1017       RNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
1018       PP=SQRT(RNVE)
1019       PV(1,NT)=PP*SINT*SIN(PHI)
1020       PV(2,NT)=PP*SINT*COS(PHI)
1021       PV(3,NT)=PP*COST
1022   442 CONTINUE
1023 C**
1024 C** STORE ON EVENT COMMON
1025 C**
1026    40 CONTINUE
1027    42 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
1028       EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
1029       EKIN2=0.
1030       CALL TDELAY(TOF1)
1031       CALL GRNDM(RNDM,1)
1032       RAN=RNDM(1)
1033       TOF=TOF-TOF1*LOG(RAN)
1034       DO 44 I=1,NT
1035       EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
1036       IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I)
1037       PV(7,I)=TOF
1038       PV(8,I)=ABS(IPA(I))
1039       PV(9,I)=0.
1040    44 PV(10,I)=0.
1041       IF (NPRT(4)) WRITE(NEWBCD,2006) NT,EKIN,ENP(1),ENP(3),EKIN1,EKIN2
1042       INTCT=INTCT+1.
1043       CALL SETCUR(NT)
1044       NTK=NTK+1
1045       IF(NT.EQ.1) GO TO 9999
1046       DO 50 II=2,NT
1047       I=II-1
1048       IF(NTOT.LT.NSIZE/12) GOTO 43
1049       GO TO 9999
1050    43 CALL SETTRK(I)
1051    50 CONTINUE
1052 C
1053  2002 FORMAT(' *GENXPT* PRODUCTION OF FINAL STATE KINEMATIC AFTER ',I3,
1054      $ ' TRIALS.  KINETIC ENERGIES ',2F6.2,' OUT OF ',2F6.2)
1055  2000 FORMAT(' *GENXPT* CMS PARAMETERS OF FINAL STATE PARTICLES,',
1056      $ ' ENERGIES IN INITIAL AND FINAL STATE ',2F6.2)
1057  2001 FORMAT(' *GENXPT* TRACK',2X,I3,2X,10F8.3,2X,I3,2X,F4.0)
1058  2003 FORMAT(' *GENXPT* TETA,EKIN0,EKIN1,EKIN ',4F10.4)
1059  2006 FORMAT(' *GENXPT* COMP.',1X,I5,1X,5F7.2)
1060  3001 FORMAT(' *GENXPT* NUCLEAR EXCITATION',I5,
1061      $ ' PARTICLES PRODUCED IN ADDITION  TO ',I5,' NORMAL PARTICLES')
1062  3002 FORMAT(' *GENXPT* AVAILABLE ENERGIES ',2F10.4,
1063      $ ' FOR ',2I3,' PARTICLES IN BEAM/TARGET FRAGM. REGION',
1064      $ ' WITH IPA/SIDE ARRAY '/
1065      $ 1H ,5X,10(I3,2X,F3.0,4X))
1066  3003 FORMAT(' *GENXPT* ',I3,' BLACK TRACK PARTICLES PRODUCED',
1067      $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV')
1068  3004 FORMAT(' *GENXPT* ',I5,' HEAVY FRAGMENTS PRODUCED',
1069      $ ' WITH TOTAL ENERGY OF',F8.4,' GEV')
1070 C
1071  9999 CONTINUE
1072 C
1073       END