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