GEANT321/ghrout/higxpt.F

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