GEANT321/ghrout/higclu.F

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