GEANT321/gheisha/twoclu.F

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