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