GEANT321/ghrout/higclu.F.ori

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1996/04/18 15:30:15  ravndal
   6 * PCM index overflow protection
   7 *
   8 * Revision 1.1.1.1  1995/10/24 10:21:13  cernlib
   9 * Geant
  10 *
  11 *
  12 #include "geant321/pilot.h"
  13 *CMZ :  3.21/02 29/03/94  15.41.40  by  S.Giani
  14 *-- Author :
  15       SUBROUTINE HIGCLU(IPPP,NFL,AVERN)
  16 C
  17 C *** GENERATION OF X- AND PT- VALUES FOR ALL PRODUCED PARTICLES ***
  18 C *** NVE 01-AUG-1988 CERN GENEVA ***
  19 C
  20 C ORIGIN : H.FESEFELDT (11-OCT-1987)
  21 C
  22 C A SIMPLE TWO CLUSTER MODEL IS USED
  23 C THIS SHOULD BE SUFFICIENT FOR LOW ENERGY INTERACTIONS
  24 C
  25 C
  26 #include "geant321/s_defcom.inc"
  27 #include "geant321/s_genio.inc"
  28 C
  29       REAL NUCSUP
  30       DIMENSION SIDE(MXGKCU),C1PAR(5),G1PAR(5),NUCSUP(6)
  31       DIMENSION RNDM(3)
  32       DATA C1PAR/0.6,0.6,0.35,0.15,0.10/
  33       DATA G1PAR/2.6,2.6,1.8,1.30,1.20/
  34       DATA NUCSUP/1.0,0.7,0.5,0.4,0.35,0.3/
  35       DIMENSION PSUP(6)
  36       DATA PSUP/3.,6.,20.,50.,100.,1000./
  37 C     DATA CB/3.0/
  38       DATA CB/0.01/
  39       BPP(X)=4.000+1.600*LOG(X)
  40 C
  41       MX =MXGKPV-20
  42       MX1=MX+1
  43       MX2=MX+2
  44       MX3=MX+3
  45       MX4=MX+4
  46       MX5=MX+5
  47       MX6=MX+6
  48       MX7=MX+7
  49       MX8=MX+8
  50       EK=ENP(5)
  51       EN=ENP(6)
  52       P=ENP(7)
  53       S=ENP(8)
  54       RS=ENP(9)
  55       CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
  56       IF(P.LT.0.001) GOTO 60
  57       NT=0
  58 C**
  59 C** CHECK MASS-INDICES FOR ALL PARTICLES
  60 C**
  61       DO 1 I=1,100
  62       IF(IPA(I).EQ.0) GOTO 1
  63       NT=NT+1
  64       IPA(NT)=IPA(I)
  65     1 CONTINUE
  66       CALL VZERO(IPA(NT+1),MXGKCU-NT)
  67 C**
  68 C** SET THE EFFECTICE 4-MOMENTUM-VECTOR FOR INTERACTION
  69 C**
  70       PV( 1,MXGKPV-1)=P*PX
  71       PV( 2,MXGKPV-1)=P*PY
  72       PV( 3,MXGKPV-1)=P*PZ
  73       PV( 4,MXGKPV-1)=EN
  74       PV( 5,MXGKPV-1)=AMAS
  75       PV( 6,MXGKPV-1)=NCH
  76       PV( 7,MXGKPV-1)=TOF
  77       PV( 8,MXGKPV-1)=IPART
  78       PV( 9,MXGKPV-1)=0.
  79       PV(10,MXGKPV-1)=USERW
  80       IER(48)=IER(48)+1
  81 C**
  82 C** DISTRIBUTE PARTICLES IN FORWARD AND BACKWARD HEMISPHERE OF CMS
  83 C** OF THE HADRON NUCLEON INTERACTION
  84 C**
  85       SIDE(1)= 1.
  86       SIDE(2)=-1.
  87       TARG=0.
  88       IFOR=1
  89       IBACK=1
  90       DO 3 I=1,NT
  91       IF (I .LE. 2) GO TO 78
  92       SIDE(I)=1.
  93       CALL GRNDM(RNDM,1)
  94       IF (RNDM(1) .LT. 0.5) SIDE(I)=-1.
  95       IF (SIDE(I) .LT. 0.) GO TO 76
  96 C
  97 C --- PARTICLE IN FORWARD HEMISPHERE ---
  98  77   CONTINUE
  99       IFOR=IFOR+1
 100       IF (IFOR .LE. 18) GO TO 78
 101 C
 102 C --- CHANGE IT TO BACKWARD ---
 103       SIDE(I)=-1.
 104       IFOR=IFOR-1
 105       IBACK=IBACK+1
 106       GO TO 78
 107 C
 108 C --- PARTICLE IN BACKWARD HEMISPHERE ---
 109  76   CONTINUE
 110       IBACK=IBACK+1
 111       IF (IBACK .LE. 18) GO TO 78
 112 C
 113 C --- CHANGE IT TO FORWARD ---
 114       SIDE(I)=1.
 115       IBACK=IBACK-1
 116       IFOR=IFOR+1
 117 C**
 118 C** SUPPRESSION OF CHARGED PIONS FOR VARIOUS REASONS
 119 C**
 120    78 IF(IPART.EQ.15.OR.IPART.GE.17) GOTO 3
 121       IF(ABS(IPA(I)).GE.10) GOTO 3
 122       IF(ABS(IPA(I)).EQ. 8) GOTO 3
 123       CALL GRNDM(RNDM,1)
 124       IF(RNDM(1).GT.(10.-P)/6.) GOTO 3
 125       CALL GRNDM(RNDM,1)
 126       IF(RNDM(1).GT.ATNO2/300.) GOTO 3
 127       IPA(I)=14
 128       CALL GRNDM(RNDM,1)
 129       IF(RNDM(1).GT.ZNO2/ATNO2) IPA(I)=16
 130       TARG=TARG+1.
 131     3 CONTINUE
 132       TB=2.*IBACK
 133       CALL GRNDM(RNDM,1)
 134       IF(RS.LT.(2.0+RNDM(1))) TB=(2.*IBACK+NT)/2.
 135 C**
 136 C** NUCLEONS + SOME PIONS FROM INTRANUCLEAR CASCADE
 137 C**
 138       AFC=0.312+0.200*LOG(LOG(S))+S**1.5/6000.
 139       IF(AFC.GT.0.50) AFC= 0.50
 140       XTARG=AFC*(ATNO2**0.33-1.0)*TB
 141       IF(XTARG.LE.0.) XTARG=0.01
 142       DO 881 IPX=1,6
 143          IF(P.LE.PSUP(IPX)) GOTO 882
 144   881 CONTINUE
 145       IPX=6
 146   882 XPNHMF = XTARG*NUCSUP(IPX)
 147       XSHHMF = XTARG - XPNHMF
 148       IF(XSHHMF.LT.0.01) XSHHMF=0.01
 149       IF(XPNHMF.LT.0.01) XPNHMF=0.01
 150       SSHHMF=0.5*XSHHMF
 151       SPNHMF=0.9*XPNHMF
 152       RSHHMF=SSHHMF**2/XSHHMF
 153       RPNHMF=SPNHMF**2/XPNHMF
 154       IF(RSHHMF.LT.1.1) THEN
 155         CALL POISSO(XSHHMF,NSHHMF)
 156         GOTO 541
 157       ELSE
 158         RSHHMF=XSHHMF/(RSHHMF-1.)
 159         IF(RSHHMF.LE.20.) THEN
 160           CALL SVGAM7(RSHHMF,XHMF)
 161         ELSE
 162           KRSHMF=IFIX(RSHHMF+0.5)
 163           CALL SVERL2(KRSHMF,XHMF)
 164         ENDIF
 165         XSHHMF=XHMF*XSHHMF/RSHHMF
 166         CALL POISSO(XSHHMF,NSHHMF)
 167       ENDIF
 168   541 IF(RPNHMF.LE.1.1) THEN
 169         CALL POISSO(XPNHMF,NPNHMF)
 170         GOTO 542
 171       ELSE
 172         RPNHMF=XPNHMF/(RPNHMF-1.)
 173         IF(RPNHMF.LE.20.) THEN
 174            CALL SVGAM7(RPNHMF,XHMF)
 175         ELSE
 176            KRPHMF=IFIX(RPNHMF+0.5)
 177            CALL SVERL2(KRPHMF,XHMF)
 178         ENDIF
 179         XPNHMF=XHMF*XPNHMF/RPNHMF
 180         CALL POISSO(XPNHMF,NPNHMF)
 181       ENDIF
 182   542 NTARG=NSHHMF+NPNHMF
 183       NT2=NT+NTARG
 184       IF(NT2.LE.MXGKPV-30) GOTO 2
 185       NT2=MXGKPV-30
 186       NTARG=NT2-NT
 187     2 CONTINUE
 188       IF(NPRT(4))
 189      *WRITE(NEWBCD,3001) NTARG,NT
 190       NT1=NT+1
 191       IF(NTARG.EQ.0) GOTO 51
 192       DO 4 I=NT1,NT2
 193       IF(NPNHMF.GT.0) GOTO 52
 194       CALL GRNDM(RNDM,1)
 195       IPA(I)=-(7+IFIX(RNDM(1)*3.0))
 196       GOTO 4
 197    52 IPA(I)=-16
 198       PNRAT=1.-ZNO2/ATNO2
 199       CALL GRNDM(RNDM,1)
 200       IF(RNDM(1).GT.PNRAT) IPA(I)=-14
 201       TARG=TARG+1.
 202       SIDE(I)=-2.
 203       NPNHMF=NPNHMF-1
 204     4 CONTINUE
 205       NT=NT2
 206 C**
 207 C** CHOOSE MASSES AND CHARGES FOR ALL PARTICLES
 208 C**
 209    51 DO 5 I=1,NT
 210       IPA1=ABS(IPA(I))
 211       PV(5,I)=RMASS(IPA1)
 212       PV(6,I)=RCHARG(IPA1)
 213       PV(7,I)=1.
 214       IF(PV(5,I).LT.0.) PV(7,I)=-1.
 215       PV(5,I)=ABS(PV(5,I))
 216     5 CONTINUE
 217 C**
 218 C** MARK LEADING STRANGE PARTICLES
 219 C**
 220       LEAD=0
 221       IF(IPART.LT.10.OR.IPART.EQ.14.OR.IPART.EQ.16) GOTO 6
 222       IPA1=ABS(IPA(1))
 223       IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 531
 224       LEAD=IPA1
 225       GOTO 6
 226   531 IPA1=ABS(IPA(2))
 227       IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 6
 228       LEAD=IPA1
 229 C**
 230 C** CHECK AVAILABLE KINETIC ENERGY , CHANGE HEMISPHERE FOR PARTICLES
 231 C** UNTIL IT FITS
 232 C**
 233     6 IF(NT.LE.1) GOTO 60
 234       TAVAI=0.
 235       DO 7 I=1,NT
 236       IF(SIDE(I).LT.-1.5) GOTO 7
 237       TAVAI=TAVAI+ABS(PV(5,I))
 238     7 CONTINUE
 239       IF(TAVAI.LT.RS) GOTO 12
 240       IF(NPRT(4))
 241      *WRITE(NEWBCD,3002) (IPA(I),I=1,20),(SIDE(I),I=1,20),TAVAI,RS
 242  3002 FORMAT(' *HIGCLU* CHECK AVAILABLE ENERGIES'/
 243      *       1H ,20I5/1H ,20F5.0/1H ,'TAVAI,RS ',2F10.3)
 244       DO 10 I=1,NT
 245       II=NT-I+1
 246       IF(SIDE(II).LT.-1.5) GOTO 10
 247       IF(II.EQ.NT) GOTO 11
 248       NT1=II+1
 249       NT2=NT
 250       DO 8 J=NT1,NT2
 251       IPA(J-1)=IPA(J)
 252       SIDE(J-1)=SIDE(J)
 253       DO 8 K=1,10
 254     8 PV(K,J-1)=PV(K,J)
 255       GOTO 11
 256    10 CONTINUE
 257    11 SIDE(NT)=0.
 258       IPA(NT)=0
 259       NT=NT-1
 260       GOTO 6
 261    12 IF(NT.LE.1) GOTO 60
 262       B=BPP(P)
 263       IF(B.LT.CB) B=CB
 264 C**
 265 C** CHOOSE MASSES FOR THE 3 CLUSTER: 1. FORWARD CLUSTER
 266 C**   2. BACKWARD MESON CLUSTER  3. BACKWARD NUCLEON CLUSTER
 267 C**
 268       RMC0=0.
 269       RMD0=0.
 270       RME0=0.
 271       NTC=0
 272       NTD=0
 273       NTE=0
 274       DO 31 I=1,NT
 275       IF(SIDE(I).GT.0.) RMC0=RMC0+ABS(PV(5,I))
 276       IF(SIDE(I).GT.0.) NTC =NTC +1
 277       IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) RMD0=RMD0+ABS(PV(5,I))
 278       IF(                  SIDE(I).LT.-1.5) RME0=RME0+ABS(PV(5,I))
 279       IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) NTD =NTD +1
 280       IF(                  SIDE(I).LT.-1.5) NTE =NTE +1
 281    31 CONTINUE
 282    32 CALL GRNDM(RNDM,1)
 283       RAN=RNDM(1)
 284       RMC=RMC0
 285       IF(NTC.LE.1) GOTO 33
 286       NTC1=NTC
 287       IF(NTC1.GT.5) NTC1=5
 288       RMC=-LOG(1.-RAN)
 289       GPAR=G1PAR(NTC1)
 290       CPAR=C1PAR(NTC1)
 291       DUMNVE=GPAR
 292       IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
 293       RMC=RMC0+RMC**CPAR/DUMNVE
 294    33 RMD=RMD0
 295       IF(NTD.LE.1) GOTO 34
 296       NTD1=NTD
 297       IF(NTD1.GT.5) NTD1=5
 298       CALL GRNDM(RNDM,1)
 299       RAN=RNDM(1)
 300       RMD=-LOG(1.-RAN)
 301       GPAR=G1PAR(NTD1)
 302       CPAR=C1PAR(NTD1)
 303       DUMNVE=GPAR
 304       IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
 305       RMD=RMD0+RMD**CPAR/DUMNVE
 306    34 IF(RMC+RMD.LE.RS) GOTO 35
 307       IF (RMC.LE.RMC0.AND.RMD.LE.RMD0) THEN
 308          HNRMDC = 0.999*RS/(RMC+RMD)
 309          RMD = RMD*HNRMDC
 310          RMC = RMC*HNRMDC
 311       ELSE
 312          RMC=0.1*RMC0+0.9*RMC
 313          RMD=0.1*RMD0+0.9*RMD
 314       END IF
 315       GOTO 34
 316    35 IF(NTE.LE.0) GOTO 38
 317       RME=RME0
 318       IF(NTE.EQ.1) GOTO 38
 319       NTE1=NTE
 320       IF(NTE1.GT.5) NTE1=5
 321       CALL GRNDM(RNDM,1)
 322       RAN=RNDM(1)
 323       RME=-LOG(1.-RAN)
 324       GPAR=G1PAR(NTE1)
 325       CPAR=C1PAR(NTE1)
 326       DUMNVE=GPAR
 327       IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
 328       RME=RME0+RME**CPAR/DUMNVE
 329 C**
 330 C** SET BEAM , TARGET OF FIRST INTERACTION IN CMS
 331 C**
 332    38 PV( 1,MX1) =0.
 333       PV( 2,MX1) =0.
 334       PV( 3,MX1) =P
 335       PV( 5,MX1) =ABS(AMAS)
 336       PV( 4,MX1) =SQRT(P*P+AMAS*AMAS)
 337       PV( 1,MX2) =0.
 338       PV( 2,MX2) =0.
 339       PV( 3,MX2) =0.
 340       PV( 4,MX2) =MP
 341       PV( 5,MX2) =MP
 342
 343 C** TRANSFORM INTO CMS.
 344
 345       CALL ADD(MX1,MX2,MX)
 346       CALL LOR(MX1,MX,MX1)
 347       CALL LOR(MX2,MX,MX2)
 348       PF=(S+RMD*RMD-RMC*RMC)**2 - 4*S*RMD*RMD
 349       IF(PF.LT.0.0001) PF=0.0001
 350       DUMNVE=2.0*RS
 351       IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
 352       PF=SQRT(PF)/DUMNVE
 353       IF(NPRT(4)) WRITE(6,2002) PF,RMC,RMD,RS
 354 C**
 355 C** SET FINAL STATE MASSES AND ENERGIES IN CMS
 356 C**
 357       PV(5,MX3) =RMC
 358       PV(5,MX4) =RMD
 359       PV(4,MX3) =SQRT(PF*PF+RMC*RMC)
 360       PV(4,MX4) =SQRT(PF*PF+RMD*RMD)
 361 C**
 362 C** SET |T| AND |TMIN|
 363 C**
 364       T=-1.0E10
 365       CALL GRNDM(RNDM,1)
 366       IF (B .NE. 0.0) T=LOG(1.-RNDM(1))/B
 367       CALL LENGTX(MX1,PIN)
 368       TACMIN=(PV(4,MX1) -PV(4,MX3))**2 -(PIN-PF)**2
 369 C**
 370 C** CACULATE (SIN(TETA/2.)**2 AND COS(TETA), SET AZIMUTH ANGLE PHI
 371 C**
 372       DUMNVE=4.0*PIN*PF
 373       IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
 374       CTET=-(T-TACMIN)/DUMNVE
 375       CTET=1.0-2.0*CTET
 376       IF (CTET .GT. 1.0) CTET=1.0
 377       IF (CTET .LT. -1.0) CTET=-1.0
 378       DUMNVE=1.0-CTET*CTET
 379       IF (DUMNVE .LT. 0.0) DUMNVE=0.0
 380       STET=SQRT(DUMNVE)
 381       CALL GRNDM(RNDM,1)
 382       PHI=RNDM(1)*TWPI
 383 C**
 384 C** CALCULATE FINAL STATE MOMENTA IN CMS
 385 C**
 386       PV(1,MX3) =PF*STET*SIN(PHI)
 387       PV(2,MX3) =PF*STET*COS(PHI)
 388       PV(3,MX3) =PF*CTET
 389       PV(1,MX4) =-PV(1,MX3)
 390       PV(2,MX4) =-PV(2,MX3)
 391       PV(3,MX4) =-PV(3,MX3)
 392 C**
 393 C** SIMULATE BACKWARD NUCLEON CLUSTER IN LAB. SYSTEM AND TRANSFORM IN
 394 C** CMS.
 395 C**
 396       IF(NTE.EQ.0) GOTO 28
 397       GA=1.2
 398       EKIT1=0.04
 399       EKIT2=0.6
 400       IF(EK.GT.5.) GOTO 666
 401       EKIT1=EKIT1*EK**2/25.
 402       EKIT2=EKIT2*EK**2/25.
 403   666 A=(1.-GA)/(EKIT2**(1.-GA)-EKIT1**(1.-GA))
 404       DO 29 I=1,NT
 405       IF(SIDE(I).GT.-1.5) GOTO 29
 406       CALL GRNDM(RNDM,3)
 407       RAN=RNDM(1)
 408       EKIT=(RAN*(1.-GA)/A+EKIT1**(1.-GA))**(1./(1.-GA))
 409       PV(4,I)=EKIT+PV(5,I)
 410       DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
 411       PP=SQRT(DUMNVE)
 412       RAN=RNDM(2)
 413       COST=LOG(2.23*RAN+0.383)/0.96
 414       IF (COST .LT. -1.0) COST=-1.0
 415       IF (COST .GT. 1.0) COST=1.0
 416       DUMNVE=1.0-COST*COST
 417       IF (DUMNVE .LT. 0.0) DUMNVE=0.0
 418       SINT=SQRT(DUMNVE)
 419       PHI=TWPI*RNDM(3)
 420       PV(1,I)=PP*SINT*SIN(PHI)
 421       PV(2,I)=PP*SINT*COS(PHI)
 422       PV(3,I)=PP*COST
 423       CALL LOR(I,MX,I)
 424    29 CONTINUE
 425 C**
 426 C** FRAGMENTATION OF FORWARD CLUSTER AND BACKWARD MESON CLUSTER
 427 C**
 428    28 PV(1,1)=PV(1,MX3)
 429       PV(2,1)=PV(2,MX3)
 430       PV(3,1)=PV(3,MX3)
 431       PV(4,1)=PV(4,MX3)
 432       PV(1,2)=PV(1,MX4)
 433       PV(2,2)=PV(2,MX4)
 434       PV(3,2)=PV(3,MX4)
 435       PV(4,2)=PV(4,MX4)
 436       DO 17 I=MX5,MX6
 437       DO 16 J=1,3
 438    16 PV(J,I)=-PV(J,I-2)
 439       DO 17 J=4,5
 440    17 PV(J,I)= PV(J,I-2)
 441       KGENEV=1
 442       IF(NTC.LE.1) GOTO 26
 443       TECM=PV(5,MX3)
 444       NPG=0
 445       DO 18 I=1,NT
 446       IF(SIDE(I).LT.0.) GOTO 18
 447       IF(NPG.EQ.18) THEN
 448          SIDE(I)=-SIDE(I)
 449          GOTO 18
 450       ENDIF
 451       NPG=NPG+1
 452       AMASS(NPG)=ABS(PV(5,I))
 453    18 CONTINUE
 454       IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG)
 455       CALL PHASP
 456       NPG=0
 457       DO 19 I=1,NT
 458       IF(SIDE(I).LT.0.OR.NPG.GE.18) GOTO 19
 459       NPG=NPG+1
 460       PV(1,I)=PCM(1,NPG)
 461       PV(2,I)=PCM(2,NPG)
 462       PV(3,I)=PCM(3,NPG)
 463       PV(4,I)=PCM(4,NPG)
 464       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 465       CALL LOR(I,MX5,I)
 466       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
 467    19 CONTINUE
 468    26 IF(NTD.LE.1) GOTO 27
 469       TECM=PV(5,MX4)
 470       NPG=0
 471       DO 20 I=1,NT
 472       IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 20
 473       IF(NPG.EQ.18) THEN
 474          SIDE(I)=-2.
 475          PV(4,I)=ABS(PV(5,I))
 476          DO 48 J=1,3
 477             PV(J,I)=0.
 478    48    CONTINUE
 479          GOTO 20
 480       ENDIF
 481       NPG=NPG+1
 482       AMASS(NPG)=ABS(PV(5,I))
 483    20 CONTINUE
 484       IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG)
 485       CALL PHASP
 486       NPG=0
 487       DO 21 I=1,NT
 488       IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 21
 489       NPG=NPG+1
 490       PV(1,I)=PCM(1,NPG)
 491       PV(2,I)=PCM(2,NPG)
 492       PV(3,I)=PCM(3,NPG)
 493       PV(4,I)=PCM(4,NPG)
 494       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 495       CALL LOR(I,MX6,I)
 496       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
 497    21 CONTINUE
 498 C**
 499 C** LORENTZ TRANSFORMATION IN LAB SYSTEM
 500 C**
 501    27 TARG=0.
 502       DO 36 I=1,NT
 503       IF(PV(5,I).GT.0.5) TARG=TARG+1.
 504       CALL LOR(I,MX2,I)
 505    36 CONTINUE
 506       IF(TARG.LT.0.5) TARG=1.
 507 C**
 508 C** SOMETIMES THE LEADING STRANGE PARTICLES ARE LOST , SET THEM BACK
 509 C**
 510       IF(LEAD.EQ.0) GOTO 6085
 511       DO 6081 I=1,NT
 512       IF(ABS(IPA(I)).EQ.LEAD) GOTO 6085
 513  6081 CONTINUE
 514       I=1
 515       IF(LEAD.GE.14.AND.ABS(IPA(2)).GE.14) I=2
 516       IF(LEAD.LT.14.AND.ABS(IPA(2)).LT.14) I=2
 517       IPA(I)=LEAD
 518       EKIN=PV(4,I)-ABS(PV(5,I))
 519       PV(5,I)=RMASS(LEAD)
 520       PV(7,I)=1.
 521       IF(PV(5,I).LT.0.) PV(7,I)=-1.
 522       PV(5,I)=ABS(PV(5,I))
 523       PV(6,I)=RCHARG(LEAD)
 524       PV(4,I)=PV(5,I)+EKIN
 525       CALL LENGTX(I,PP)
 526       DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
 527       PP1=SQRT(DUMNVE)
 528 C
 529       IF (PP .GE. 1.0E-6) GO TO 8000
 530       CALL GRNDM(RNDM,2)
 531       RTHNVE=PI*RNDM(1)
 532       PHINVE=TWPI*RNDM(2)
 533       PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
 534       PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
 535       PV(3,I)=PP1*COS(RTHNVE)
 536       GO TO 8001
 537  8000 CONTINUE
 538       PV(1,I)=PV(1,I)*PP1/PP
 539       PV(2,I)=PV(2,I)*PP1/PP
 540       PV(3,I)=PV(3,I)*PP1/PP
 541  8001 CONTINUE
 542 C
 543 C** FOR VARIOUS REASONS, THE ENERGY BALANCE IS NOT SUFFICIENT,
 544 C** CHECK THAT,  ENERGY BALANCE, ANGLE OF FINAL SYSTEM E.T.C.
 545  6085 KGENEV=1
 546       PV(1,MX4) =0.
 547       PV(2,MX4) =0.
 548       PV(3,MX4) =P
 549       PV(4,MX4) =SQRT(P*P+AMAS*AMAS)
 550       PV(5,MX4) =ABS(AMAS)
 551       EKIN0=PV(4,MX4) -PV(5,MX4)
 552       PV(1,MX5) =0.
 553       PV(2,MX5) =0.
 554       PV(3,MX5) =0.
 555       PV(4,MX5) =MP*TARG
 556       PV(5,MX5) =PV(4,MX5)
 557       EKIN=PV(4,MX4) +PV(4,MX5)
 558       I=MX4
 559       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 560       I=MX5
 561       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
 562       CALL ADD(MX4,MX5,MX6)
 563       CALL LOR(MX4,MX6,MX4)
 564       CALL LOR(MX5,MX6,MX5)
 565       TECM=PV(4,MX4) +PV(4,MX5)
 566       NPG=NT
 567       PV(1,MX8) =0.
 568       PV(2,MX8) =0.
 569       PV(3,MX8) =0.
 570       PV(4,MX8) =0.
 571       PV(5,MX8) =0.
 572       EKIN1=0.
 573       DO 598 I=1,NPG
 574       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
 575       CALL ADD(MX8,I,MX8)
 576       EKIN1=EKIN1+PV(4,I)-PV(5,I)
 577       EKIN=EKIN-PV(5,I)
 578       IF(I.GT.18) GOTO 598
 579       AMASS(I)=PV(5,I)
 580   598 CONTINUE
 581       IF(NPG.GT.18) GOTO 597
 582       CALL PHASP
 583       EKIN=0.
 584       DO 599 I=1,NPG
 585       PV(1,MX7)=PCM(1,I)
 586       PV(2,MX7)=PCM(2,I)
 587       PV(3,MX7)=PCM(3,I)
 588       PV(4,MX7)=PCM(4,I)
 589       PV(5,MX7)=AMASS(I)
 590       CALL LOR(MX7,MX5,MX7)
 591   599 EKIN=EKIN+PV(4,MX7)-PV(5,MX7)
 592       CALL ANG(MX8,MX4,COST,TETA)
 593       IF(NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1,EKIN
 594 C**
 595 C** MAKE SHURE, THAT  KINETIC ENERGIES ARE CORRECT
 596 C** THE 3. CLUSTER IS NOT PRODUCED WITHIN PROPER KINEMATICS!!!
 597 C** EKIN= KINETIC ENERGY THEORETICALLY
 598 C** EKIN1= KINETIC ENERGY SIMULATED
 599 C**
 600   597 IF(EKIN1.EQ.0.) GOTO 600
 601       PV(1,MX7) =0.
 602       PV(2,MX7) =0.
 603       PV(3,MX7) =0.
 604       PV(4,MX7) =0.
 605       PV(5,MX7) =0.
 606       WGT=EKIN/EKIN1
 607       EKIN1=0.
 608       DO 602 I=1,NT
 609       EKIN=PV(4,I)-PV(5,I)
 610       EKIN=EKIN*WGT
 611       PV(4,I)=EKIN+PV(5,I)
 612       DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
 613       PP=SQRT(DUMNVE)
 614       CALL LENGTX(I,PP1)
 615 C
 616       IF (PP1 .GE. 1.0E-6) GO TO 8002
 617       CALL GRNDM(RNDM,2)
 618       RTHNVE=PI*RNDM(1)
 619       PHINVE=TWPI*RNDM(2)
 620       PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE)
 621       PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE)
 622       PV(3,I)=PP*COS(RTHNVE)
 623       GO TO 8003
 624  8002 CONTINUE
 625       PV(1,I)=PV(1,I)*PP/PP1
 626       PV(2,I)=PV(2,I)*PP/PP1
 627       PV(3,I)=PV(3,I)*PP/PP1
 628  8003 CONTINUE
 629 C
 630       EKIN1=EKIN1+EKIN
 631       CALL ADD(MX7,I,MX7)
 632   602 CONTINUE
 633       CALL ANG(MX7,MX4,COST,TETA)
 634       IF(NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1
 635 C**
 636 C** ROTATE IN DIRECTION OF Z-AXIS, SEE COMMENTS IN 'GENXPT'
 637 C**
 638   600 PV(1,MX7)=0.
 639       PV(2,MX7)=0.
 640       PV(3,MX7)=0.
 641       PV(4,MX7)=0.
 642       PV(5,MX7)=0.
 643       DO 596 I=1,NT
 644   596 CALL ADD(MX7,I,MX7)
 645 *          call rannor(ran1,ran2)
 646       CALL GRNDM(RNDM,2)
 647       RY=RNDM(1)
 648       RZ=RNDM(2)
 649       RX=6.283185*RZ
 650       A1=SQRT(-2.*LOG(RY))
 651       RAN1=A1*SIN(RX)
 652       RAN2=A1*COS(RX)
 653       PV(1,MX7)=PV(1,MX7)+RAN1*0.020*TARG
 654       PV(2,MX7)=PV(2,MX7)+RAN2*0.020*TARG
 655       CALL DEFS(MX4,MX7,MX8)
 656       PV(1,MX7)=0.
 657       PV(2,MX7)=0.
 658       PV(3,MX7)=0.
 659       PV(4,MX7)=0.
 660       PV(5,MX7)=0.
 661       DO 595 I=1,NT
 662       CALL TRAC(I,MX8,I)
 663   595 CALL ADD(MX7,I,MX7)
 664       CALL ANG(MX7,MX4,COST,TETA)
 665       IF(NPRT(4)) WRITE(NEWBCD,2003) TETA
 666 C**
 667 C** ROTATE IN DIRECTION OF PRIMARY PARTICLE
 668 C**
 669       DEKIN=0.
 670       NPIONS=0
 671       EK1=0.
 672       DO 25 I=1,NT
 673       CALL DEFS1(I,MXGKPV-1,I)
 674       IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
 675       IF(ATNO2.LT.1.5) GOTO 25
 676       CALL LENGTX(I,PP)
 677       EKIN=PV(4,I)-ABS(PV(5,I))
 678       CALL NORMAL(RAN)
 679       EKIN=EKIN-CFA*(1.+0.5*RAN)
 680       IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
 681       CALL STEEQ(XXH,I)
 682       DEKIN=DEKIN+EKIN*(1.-XXH)
 683       EKIN=EKIN*XXH
 684       IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) NPIONS=NPIONS+1
 685       IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) EK1=EK1+EKIN
 686       PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
 687       PV(4,I)=EKIN+ABS(PV(5,I))
 688 C
 689       IF (PP .GE. 1.0E-6) GO TO 8004
 690       CALL GRNDM(RNDM,2)
 691       RTHNVE=PI*RNDM(1)
 692       PHINVE=TWPI*RNDM(2)
 693       PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
 694       PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
 695       PV(3,I)=PP1*COS(RTHNVE)
 696       GO TO 8005
 697  8004 CONTINUE
 698       PV(1,I)=PV(1,I)*PP1/PP
 699       PV(2,I)=PV(2,I)*PP1/PP
 700       PV(3,I)=PV(3,I)*PP1/PP
 701  8005 CONTINUE
 702 C
 703    25 CONTINUE
 704       IF(EK1.EQ.0.) GOTO 23
 705       IF(NPIONS.LE.0) GOTO 23
 706       DEKIN=1.+DEKIN/EK1
 707       DO 22 I=1,NT
 708       IF(ABS(IPA(I)).LT.7.OR.ABS(IPA(I)).GT.9) GOTO 22
 709       CALL LENGTX(I,PP)
 710       EKIN=PV(4,I)-ABS(PV(5,I))
 711       EKIN=EKIN*DEKIN
 712       IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
 713       PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
 714       PV(4,I)=EKIN+ABS(PV(5,I))
 715 C
 716       IF (PP .GE. 1.0E-6) GO TO 8006
 717       CALL GRNDM(RNDM,2)
 718       RTHNVE=PI*RNDM(1)
 719       PHINVE=TWPI*RNDM(2)
 720       PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
 721       PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
 722       PV(3,I)=PP1*COS(RTHNVE)
 723       GO TO 8007
 724  8006 CONTINUE
 725       PV(1,I)=PV(1,I)*PP1/PP
 726       PV(2,I)=PV(2,I)*PP1/PP
 727       PV(3,I)=PV(3,I)*PP1/PP
 728  8007 CONTINUE
 729 C
 730    22 CONTINUE
 731    23 IF(ATNO2.LT.1.5) GOTO 40
 732 C**
 733 C** ADD BLACK TRACK PARTICLES
 734 C**
 735       CALL HIGHAB(SPROB)
 736       CALL GRNDM(RNDM,1)
 737       IF(RNDM(1).LT.SPROB) GOTO 40
 738       TEX=ENP(1)
 739       SPALL=TARG
 740       IF(TEX.LT.0.001) GOTO 445
 741       BLACK=(1.5+1.25*TARG)*ENP(1)/(ENP(1)+ENP(3))
 742       CALL POISSO(BLACK,NBL)
 743       IF(NPRT(4))
 744      *WRITE(NEWBCD,3003) NBL,TEX
 745       IF(IFIX(TARG)+NBL.GT.ATNO2) NBL=ATNO2-TARG
 746       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
 747       IF(NBL.LE.0) GOTO 445
 748       EKIN=TEX/NBL
 749       EKIN2=0.
 750       CALL STEEP(XX)
 751       DO 441 I=1,NBL
 752       CALL GRNDM(RNDM,1)
 753       IF(RNDM(1).LT.SPROB) GOTO 441
 754       IF(NT.EQ.MXGKPV-2) GOTO 441
 755       IF(EKIN2.GT.TEX) GOTO 443
 756       CALL GRNDM(RNDM,1)
 757       RAN1=RNDM(1)
 758       CALL NORMAL(RAN2)
 759       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
 760       IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
 761       EKIN1=EKIN1*XX
 762       EKIN2=EKIN2+EKIN1
 763       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
 764       IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
 765       IPA1=16
 766       PNRAT=1.-ZNO2/ATNO2
 767       CALL GRNDM(RNDM,3)
 768       IF(RNDM(1).GT.PNRAT) IPA1=14
 769       NT=NT+1
 770       SPALL=SPALL+1.
 771       COST=-1.0+RNDM(2)*2.0
 772       DUMNVE=1.0-COST*COST
 773       IF (DUMNVE .LT. 0.0) DUMNVE=0.0
 774       SINT=SQRT(DUMNVE)
 775       PHI=TWPI*RNDM(3)
 776       IPA(NT)=-IPA1
 777       SIDE(NT)=-4.
 778       PV(5,NT)=ABS(RMASS(IPA1))
 779       PV(6,NT)=RCHARG(IPA1)
 780       PV(7,NT)=1.
 781       PV(4,NT)=EKIN1+PV(5,NT)
 782       DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
 783       PP=SQRT(DUMNVE)
 784       PV(1,NT)=PP*SINT*SIN(PHI)
 785       PV(2,NT)=PP*SINT*COS(PHI)
 786       PV(3,NT)=PP*COST
 787   441 CONTINUE
 788   443 IF(ATNO2.LT.10.) GOTO 445
 789       IF(EK.GT.2.0) GOTO 445
 790       II=NT+1
 791       KK=0
 792       EKA=EK
 793       IF(EKA.GT.1.) EKA=EKA*EKA
 794       IF(EKA.LT.0.1) EKA=0.1
 795       IKA=3.6*EXP((ZNO2**2/ATNO2-35.56)/6.45)/EKA
 796       IF(IKA.LE.0) GO TO 445
 797       DO 444 I=1,NT
 798       II=II-1
 799       IF(IPA(II).NE.-14) GOTO 444
 800       IPA(II)=-16
 801       IPA1  = 16
 802       PV(5,II)=ABS(RMASS(IPA1))
 803       PV(6,II)=RCHARG(IPA1)
 804       KK=KK+1
 805       IF(KK.GT.IKA) GOTO 445
 806   444 CONTINUE
 807   445 TEX=ENP(3)
 808       IF(TEX.LT.0.001) GOTO 40
 809       BLACK=(1.5+1.25*TARG)*ENP(3)/(ENP(1)+ENP(3))
 810       CALL POISSO(BLACK,NBL)
 811       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
 812       IF(NBL.LE.0) GOTO 40
 813       EKIN=TEX/NBL
 814       EKIN2=0.
 815       CALL STEEP(XX)
 816       IF(NPRT(4))
 817      *WRITE(NEWBCD,3004) NBL,TEX
 818       DO 442 I=1,NBL
 819       CALL GRNDM(RNDM,1)
 820       IF(RNDM(1).LT.SPROB) GOTO 442
 821       IF(NT.EQ.MXGKPV-2) GOTO 442
 822       IF(EKIN2.GT.TEX) GOTO 40
 823       CALL GRNDM(RNDM,1)
 824       RAN1=RNDM(1)
 825       CALL NORMAL(RAN2)
 826       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
 827       IF(EKIN1.LT.0.0) EKIN1=-0.005*LOG(RAN1)
 828       EKIN1=EKIN1*XX
 829       EKIN2=EKIN2+EKIN1
 830       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
 831       IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
 832       CALL GRNDM(RNDM,3)
 833       COST=-1.0+RNDM(1)*2.0
 834       DUMNVE=1.0-COST*COST
 835       IF (DUMNVE .LT. 0.0) DUMNVE=0.0
 836       SINT=SQRT(DUMNVE)
 837       PHI=TWPI*RNDM(2)
 838       RAN=RNDM(3)
 839       IPA(NT+1)=-30
 840       IF(RAN.GT.0.60) IPA(NT+1)=-31
 841       IF(RAN.GT.0.90) IPA(NT+1)=-32
 842       SIDE(NT+1)=-4.
 843       PV(5,NT+1)=(ABS(IPA(NT+1))-28)*MP
 844       SPALL=SPALL+PV(5,NT+1)*1.066
 845       IF(SPALL.GT.ATNO2) GOTO 40
 846       NT=NT+1
 847       PV(6,NT)=1.
 848       IF(IPA(NT).EQ.-32) PV(6,NT)=2.
 849       PV(7,NT)=1.
 850       PV(4,NT)=PV(5,NT)+EKIN1
 851       DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
 852       PP=SQRT(DUMNVE)
 853       PV(1,NT)=PP*SINT*SIN(PHI)
 854       PV(2,NT)=PP*SINT*COS(PHI)
 855       PV(3,NT)=PP*COST
 856   442 CONTINUE
 857 C**
 858 C** STORE ON EVENT COMMON
 859 C**
 860    40 CALL GRNDM(RNDM,1)
 861       IF(RS.GT.(4.+RNDM(1)*1.)) GOTO 42
 862       DO 41 I=1,NT
 863       CALL LENGTX(I,ETB)
 864       IF(ETB.LT.P) GOTO 41
 865       ETF=P
 866       PV(4,I)=SQRT(PV(5,I)**2+ETF**2)
 867       DUMNVE=ETB
 868       IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
 869       ETF=ETF/DUMNVE
 870       PV(1,I)=PV(1,I)*ETF
 871       PV(2,I)=PV(2,I)*ETF
 872       PV(3,I)=PV(3,I)*ETF
 873    41 CONTINUE
 874    42 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
 875       EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
 876       EKIN2=0.
 877       CALL TDELAY(TOF1)
 878       CALL GRNDM(RNDM,1)
 879       RAN=RNDM(1)
 880       TOF=TOF-TOF1*LOG(RAN)
 881       DO 44 I=1,NT
 882       IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I)
 883       PV(7,I)=TOF
 884       PV(8,I)=ABS(IPA(I))
 885       PV(9,I)=0.
 886    44 PV(10,I)=0.
 887       CALL GHETUN(NT)
 888       DO 45 I=1,NT
 889          EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
 890    45 CONTINUE
 891       EKIN2=(EKIN2-EKIN)/EKIN
 892       IF(NPRT(4))
 893      $            WRITE(NEWBCD,2006) NT,EKIN,ENP(1),ENP(3),EKIN1,EKIN2
 894       IF(EKIN2.GT.0.2) GOTO 60
 895       INTCT=INTCT+1.
 896       NMODE=3
 897       IF(SPALL.LT.0.5.AND.ATNO2.GT.1.5) NMODE=14
 898       CALL SETCUR(NT)
 899       NTK=NTK+1
 900       IF(NT.EQ.1) GOTO 300
 901       DO 50 II=2,NT
 902       I=II-1
 903       IF(NTOT.LT.NSIZE/12) GOTO 43
 904       GO TO 9999
 905    43 CALL SETTRK(I)
 906    50 CONTINUE
 907  300  CONTINUE
 908       GO TO 9999
 909 C**
 910 C** IT IS NOT POSSIBLE TO PRODUCE A PROPER TWO CLUSTER FINAL STATE.
 911 C** CONTINUE WITH QUASI ELASTIC SCATTERING
 912 C**
 913    60 IF(NPRT(4)) WRITE(NEWBCD,2005)
 914       DO 61 I=3,MXGKCU
 915    61 IPA(I)=0
 916       IPA(1)=IPART
 917       IPA(2)=14
 918       IF(NFL.EQ.2) IPA(2)=16
 919       CALL TWOB(IPPP,NFL,AVERN)
 920       GO TO 9999
 921 C
 922  2000 FORMAT(' *HIGCLU* CMS PARAMETERS OF FINAL STATE PARTICLES',
 923      $ ' AFTER ',I3,' TRIALS')
 924  2001 FORMAT(' *HIGCLU* TRACK',2X,I3,2X,10F8.2,2X,I3,2X,F3.0)
 925  2002 FORMAT(' *HIGCLU* MOMENTUM ',F8.3,' MASSES ',2F8.4,' RS ',F8.4)
 926  2003 FORMAT(' *HIGCLU* TETA,EKIN0,EKIN1,EKIN ',4F10.4)
 927  2004 FORMAT(' *HIGCLU* TECM,NPB,MASSES: ',F10.4,1X,I3,1X,8F10.4/
 928      $ 1H ,26X,15X,8F10.4)
 929  2005 FORMAT(' *HIGCLU* NUMBER OF FINAL STATE PARTICLES',
 930      $ ' LESS THAN 2 ==> CONTINUE WITH 2-BODY SCATTERING')
 931  2006 FORMAT(' *HIGCLU*  COMP.',1X,I5,1X,5F7.2)
 932  3001 FORMAT(' *HIGCLU* NUCLEAR EXCITATION ',I5,' PARTICLES PRODUCED',
 933      $ ' IN ADDITION TO',I5,' NORMAL PARTICLES')
 934  3003 FORMAT(' *HIGCLU* ',I3,' BLACK TRACK PARTICLES PRODUCED',
 935      $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV')
 936  3004 FORMAT(' *HIGCLU* ',I5,' HEAVY FRAGMENTS WITH TOTAL ENERGY OF ',
 937      $ F8.4,' GEV')
 938 C
 939  9999 CONTINUE
 940       END