[u/mrichter/AliRoot.git] / GEANT321 / ghrout / cohert.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:13  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.40  by  S.Giani
*-- Author :
      SUBROUTINE COHERT(IPPP,NFL,AVERN)
C
C *** GENERATION OF X- AND PT- VALUES FOR ALL PRODUCED PARTICLES ***
C
C
C GENERATION OF DIFFRACTION DISSOCIATION AT HIGH ENERGIES
C (NOT USED IN STANDARD VERSION)
C
#include "geant321/s_defcom.inc"
#include "geant321/s_genio.inc"
C
      REAL NUCSUP
      DIMENSION SIDE(200),C1PAR(5),G1PAR(5),NUCSUP(5)
      DIMENSION RNDM(3)
      DATA C1PAR/0.6,0.6,0.35,0.15,0.10/
      DATA G1PAR/2.6,2.6,1.8,1.30,1.20/
      DATA NUCSUP/1.0,0.8,0.6,0.5,0.4/
C     DATA CB/3.0/
      DATA CB/0.01/
      BPP(X)=5.000+0.300*LOG(X)
 
C
      MX =MXGKPV-20
      MX1=MX+1
      MX2=MX+2
      MX3=MX+3
      MX4=MX+4
      MX5=MX+5
      MX6=MX+6
      MX7=MX+7
      MX8=MX+8
      EK=ENP(5)
      EN=ENP(6)
      P=ENP(7)
      S=ENP(8)
      RS=ENP(9)
      CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
      IF(P.LT.0.001) GOTO 60
      NT=0
C**
      IREHMF=4
      IF(IABS(IPA(1)).NE.IPART) IREHMF=5
 
C** CHECK MASS-INDICES FOR ALL PARTICLES
C**
      DO 1 I=1,100
      IF(IPA(I).EQ.0) GOTO 1
      NT=NT+1
      IPA(NT)=IPA(I)
    1 CONTINUE
      CALL VZERO(IPA(NT+1),MXGKCU-NT)
C**
C** SET THE EFFECTICE 4-MOMENTUM-VECTOR FOR INTERACTION
C**
      PV( 1,MXGKPV-1)=P*PX
      PV( 2,MXGKPV-1)=P*PY
      PV( 3,MXGKPV-1)=P*PZ
      PV( 4,MXGKPV-1)=EN
      PV( 5,MXGKPV-1)=AMAS
      PV( 6,MXGKPV-1)=NCH
      PV( 7,MXGKPV-1)=TOF
      PV( 8,MXGKPV-1)=IPART
      PV( 9,MXGKPV-1)=0.
      PV(10,MXGKPV-1)=USERW
      IER(48)=IER(48)+1
C**
C** DISTRIBUTE PARTICLES IN FORWARD AND BACKWARD HEMISPHERE OF CMS
C** OF THE HADRON NUCLEON INTERACTION
C**
      SIDE(1)= 1.
      SIDE(2)=-1.
      TARG=0.
      IFOR=1
      IBACK=1
      DO 3 I=1,NT
      IF (I .LE. 2) GO TO 78
      SIDE(I)= -1.
      IF (SIDE(I) .LT. 0.) GO TO 76
C
C --- PARTICLE IN FORWARD HEMISPHERE ---
 77   CONTINUE
      IFOR=IFOR+1
      IF (IFOR .LE. 18) GO TO 78
C
C --- CHANGE IT TO BACKWARD ---
      SIDE(I)=-1.
      IFOR=IFOR-1
      IBACK=IBACK+1
      GO TO 78
C
C --- PARTICLE IN BACKWARD HEMISPHERE ---
 76   CONTINUE
      IBACK=IBACK+1
      IF (IBACK .LE. 18) GO TO 78
C
C --- CHANGE IT TO FORWARD ---
      SIDE(I)=1.
      IBACK=IBACK-1
      IFOR=IFOR+1
C**
C** SUPPRESSION OF CHARGED PIONS FOR VARIOUS REASONS
C**
   78 IF(IPART.EQ.15.OR.IPART.GE.17) GOTO 3
      IF(ABS(IPA(I)).GE.10) GOTO 3
      IF(ABS(IPA(I)).EQ. 8) GOTO 3
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).GT.(10.-P)/6.) GOTO 3
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).GT.ATNO2/300.) GOTO 3
      IPA(I)=14
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).GT.ZNO2/ATNO2) IPA(I)=16
      TARG=TARG+1.
    3 CONTINUE
      TB=2.*IBACK
      CALL GRNDM(RNDM,1)
      IF(RS.LT.(2.0+RNDM(1))) TB=(2.*IBACK+NT)/2.
C**
C** NUCLEONS + SOME PIONS FROM INTRANUCLEAR CASCADE
C**
      AFC=0.312+0.200*LOG(LOG(S))
      XTARG=AFC*(ATNO2**0.33-1.0)*TB
      IF(XTARG.LE.0.) XTARG=0.01
      CALL POISSO(XTARG,NTARG)
      NT2=NT+NTARG
      IF(NT2.LE.MXGKPV-30) GOTO 2
      NT2=MXGKPV-30
      NTARG=NT2-NT
    2 CONTINUE
      IF(NPRT(4))
     *WRITE(NEWBCD,3001) NTARG,NT
      NT1=NT+1
      IF(NTARG.EQ.0) GOTO 51
      IPX=IFIX(P/3.)+1
      IF(IPX.GT.5) IPX=5
      DO 4 I=NT1,NT2
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      IF(RAN.LT.NUCSUP(IPX)) GOTO 52
      CALL GRNDM(RNDM,1)
      IPA(I)=-(7+IFIX(RNDM(1)*3.0))
      GOTO 4
   52 IPA(I)=-16
      PNRAT=1.-ZNO2/ATNO2
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).GT.PNRAT) IPA(I)=-14
      TARG=TARG+1.
    4 SIDE(I)=-2.
      NT=NT2
C**
C** CHOOSE MASSES AND CHARGES FOR ALL PARTICLES
C**
   51 DO 5 I=1,NT
      IPA1=ABS(IPA(I))
      PV(5,I)=RMASS(IPA1)
      PV(6,I)=RCHARG(IPA1)
      PV(7,I)=1.
      IF(PV(5,I).LT.0.) PV(7,I)=-1.
      PV(5,I)=ABS(PV(5,I))
    5 CONTINUE
C**
C** MARK LEADING STRANGE PARTICLES
C**
      LEAD=0
      IF(IPART.LT.10.OR.IPART.EQ.14.OR.IPART.EQ.16) GOTO 6
      IPA1=ABS(IPA(1))
      IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 531
      LEAD=IPA1
      GOTO 6
  531 IPA1=ABS(IPA(2))
      IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 6
      LEAD=IPA1
C**
C** CHECK AVAILABLE KINETIC ENERGY , CHANGE HEMISPHERE FOR PARTICLES
C** UNTIL IT FITS
C**
    6 IF(NT.LE.1) GOTO 60
      TAVAI=0.
      DO 7 I=1,NT
      IF(SIDE(I).LT.-1.5) GOTO 7
      TAVAI=TAVAI+ABS(PV(5,I))
    7 CONTINUE
      IF(TAVAI.LT.RS) GOTO 12
      IF(NPRT(4))
     *WRITE(NEWBCD,3002) (IPA(I),I=1,20),(SIDE(I),I=1,20),TAVAI,RS
 3002 FORMAT(' *COHERT* CHECK AVAILABLE ENERGIES'/
     *       1H ,20I5/1H ,20F5.0/1H ,'TAVAI,RS ',2F10.3)
      DO 10 I=1,NT
      II=NT-I+1
      IF(SIDE(II).LT.-1.5) GOTO 10
      IF(II.EQ.NT) GOTO 11
      NT1=II+1
      NT2=NT
      DO 8 J=NT1,NT2
      IPA(J-1)=IPA(J)
      SIDE(J-1)=SIDE(J)
      DO 8 K=1,10
    8 PV(K,J-1)=PV(K,J)
      GOTO 11
   10 CONTINUE
   11 SIDE(NT)=0.
      IPA(NT)=0
      NT=NT-1
      GOTO 6
   12 IF(NT.LE.1) GOTO 60
      B=BPP(ATNO2)
      IF(B.LT.CB) B=CB
C**
C** CHOOSE MASSES FOR THE 3 CLUSTER: 1. FORWARD CLUSTER
C**   2. BACKWARD MESON CLUSTER  3. BACKWARD NUCLEON CLUSTER
C**
      RMC0=0.
      RMD0=0.
      RME0=0.
      NTC=0
      NTD=0
      NTE=0
      DO 31 I=1,NT
      IF(SIDE(I).GT.0.) RMC0=RMC0+ABS(PV(5,I))
      IF(SIDE(I).GT.0.) NTC =NTC +1
      IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) RMD0=RMD0+ABS(PV(5,I))
      IF(                  SIDE(I).LT.-1.5) RME0=RME0+ABS(PV(5,I))
      IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) NTD =NTD +1
      IF(                  SIDE(I).LT.-1.5) NTE =NTE +1
   31 CONTINUE
   32 CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      RMC=RMC0
      IF(NTC.LE.1) GOTO 33
      NTC1=NTC
      IF(NTC1.GT.5) NTC1=5
      RMC=-LOG(1.-RAN)
      GPAR=G1PAR(NTC1)
      CPAR=C1PAR(NTC1)
      DUMNVE=GPAR
      IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
      RMC=RMC0+RMC**CPAR/DUMNVE
   33 RMD=RMD0
      IF(NTD.LE.1) GOTO 34
      NTD1=NTD
      IF(NTD1.GT.5) NTD1=5
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      RMD=-LOG(1.-RAN)
      GPAR=G1PAR(NTD1)
      CPAR=C1PAR(NTD1)
      DUMNVE=GPAR
      IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
      RMD=RMD0+RMD**CPAR/DUMNVE
   34 IF(RMC+RMD.LT.RS) GOTO 35
      IF (RMC.LE.RMC0.AND.RMD.LE.RMD0) THEN
         HNRMDC = 0.999*RS/(RMC+RMD)
         RMD = RMD*HNRMDC
         RMC = RMC*HNRMDC
      ELSE
         RMC=0.1*RMC0+0.9*RMC
         RMD=0.1*RMD0+0.9*RMD
      END IF
      GOTO 34
   35 IF(NTE.LE.0) GOTO 38
      RME=RME0
      IF(NTE.EQ.1) GOTO 38
      NTE1=NTE
      IF(NTE1.GT.5) NTE1=5
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      RME=-LOG(1.-RAN)
      GPAR=G1PAR(NTE1)
      CPAR=C1PAR(NTE1)
      DUMNVE=GPAR
      IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
      RME=RME0+RME**CPAR/DUMNVE
C**
C** SET BEAM , TARGET OF FIRST INTERACTION IN CMS
C**
   38 PV(1,MX1)=0.
      PV(2,MX1)=0.
      PV(3,MX1)=P
      PV(5,MX1)=ABS(AMAS)
      PV(4,MX1)=SQRT(P*P+AMAS*AMAS)
      PV(1,MX2)=0.
      PV(2,MX2)=0.
      PV(3,MX2)=0.
      PV(4,MX2)=MP
      PV(5,MX2)=MP
 
C** TRANSFORM INTO CMS.
 
      CALL ADD(MX1,MX2,MX )
      CALL LOR(MX1,MX ,MX1)
      CALL LOR(MX2,MX ,MX2)
      PF=(S+RMD*RMD-RMC*RMC)**2 - 4*S*RMD*RMD
      IF(PF.LT.0.0001) PF=0.0001
      DUMNVE=2.0*RS
      IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
      PF=SQRT(PF)/DUMNVE
      IF(NPRT(4)) WRITE(6,2002) PF,RMC,RMD,RS
C**
C** SET FINAL STATE MASSES AND ENERGIES IN CMS
C**
      PV(5,MX3)=RMC
      PV(5,MX4)=RMD
      PV(4,MX3)=SQRT(PF*PF+PV(5,MX3)*PV(5,MX3))
      PV(4,MX4)=SQRT(PF*PF+PV(5,MX4)*PV(5,MX4))
C**
C** SET |T| AND |TMIN|
C**
      T=-1.0E10
      CALL GRNDM(RNDM,1)
      IF (B .NE. 0.0) T=LOG(1.-RNDM(1))/B
      CALL LENGTX(MX1,PIN)
      TACMIN=(PV(4,MX1)-PV(4,MX3))**2-(PIN-PF)**2
C**
C** CACULATE (SIN(TETA/2.)**2 AND COS(TETA), SET AZIMUTH ANGLE PHI
C**
      DUMNVE=4.0*PIN*PF
      IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
      CTET=-(T-TACMIN)/DUMNVE
      CTET=1.0-2.0*CTET
      IF (CTET .GT. 1.0) CTET=1.0
      IF (CTET .LT. -1.0) CTET=-1.0
      DUMNVE=1.0-CTET*CTET
      IF (DUMNVE .LT. 0.0) DUMNVE=0.0
      STET=SQRT(DUMNVE)
      CALL GRNDM(RNDM,1)
      PHI=RNDM(1)*TWPI
C**
C** CALCULATE FINAL STATE MOMENTA IN CMS
C**
      PV(1,MX3)=PF*STET*SIN(PHI)
      PV(2,MX3)=PF*STET*COS(PHI)
      PV(3,MX3)=PF*CTET
      PV(1,MX4)=-PV(1,MX3)
      PV(2,MX4)=-PV(2,MX3)
      PV(3,MX4)=-PV(3,MX3)
C**
C** SIMULATE BACKWARD NUCLEON CLUSTER IN LAB. SYSTEM AND TRANSFORM IN
C** CMS.
C**
      IF(NTE.EQ.0) GOTO 28
      GA=1.2
      EKIT1=0.04
      EKIT2=0.6
      IF(EK.GT.5.) GOTO 666
      EKIT1=EKIT1*EK**2/25.
      EKIT2=EKIT2*EK**2/25.
  666 A=(1.-GA)/(EKIT2**(1.-GA)-EKIT1**(1.-GA))
      DO 29 I=1,NT
      IF(SIDE(I).GT.-1.5) GOTO 29
      CALL GRNDM(RNDM,3)
      RAN=RNDM(1)
      EKIT=(RAN*(1.-GA)/A+EKIT1**(1.-GA))**(1./(1.-GA))
      PV(4,I)=EKIT+PV(5,I)
      DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
      PP=SQRT(DUMNVE)
      RAN=RNDM(2)
      COST=LOG(2.23*RAN+0.383)/0.96
      IF (COST .LT. -1.0) COST=-1.0
      IF (COST .GT. 1.0) COST=1.0
      DUMNVE=1.0-COST*COST
      IF (DUMNVE .LT. 0.0) DUMNVE=0.0
      SINT=SQRT(DUMNVE)
      PHI=TWPI*RNDM(3)
      PV(1,I)=PP*SINT*SIN(PHI)
      PV(2,I)=PP*SINT*COS(PHI)
      PV(3,I)=PP*COST
      CALL LOR(I,MX ,I)
   29 CONTINUE
C**
C** FRAGMENTATION OF FORWARD CLUSTER AND BACKWARD MESON CLUSTER
C**
   28 PV(1,1)=PV(1,MX3)
      PV(2,1)=PV(2,MX3)
      PV(3,1)=PV(3,MX3)
      PV(4,1)=PV(4,MX3)
      PV(1,2)=PV(1,MX4)
      PV(2,2)=PV(2,MX4)
      PV(3,2)=PV(3,MX4)
      PV(4,2)=PV(4,MX4)
      DO 17 I=MX5,MX6
      DO 16 J=1,3
   16 PV(J,I)=-PV(J,I-2)
      DO 17 J=4,5
   17 PV(J,I)= PV(J,I-2)
      KGENEV=1
      IF(NTC.LE.1) GOTO 26
      TECM= PV(5,MX3)
      NPG=0
      DO 18 I=1,NT
      IF(SIDE(I).LT.0.) GOTO 18
      NPG=NPG+1
      AMASS(NPG)=ABS(PV(5,I))
   18 CONTINUE
      IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG)
      CALL PHASP
      NPG=0
      DO 19 I=1,NT
      IF(SIDE(I).LT.0.) GOTO 19
      NPG=NPG+1
      PV(1,I)=PCM(1,NPG)
      PV(2,I)=PCM(2,NPG)
      PV(3,I)=PCM(3,NPG)
      PV(4,I)=PCM(4,NPG)
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
      CALL LOR(I,MX5,I)
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
   19 CONTINUE
   26 IF(NTD.LE.1) GOTO 27
      TECM= PV(5,MX4)
      NPG=0
      DO 20 I=1,NT
      IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 20
      NPG=NPG+1
      AMASS(NPG)=ABS(PV(5,I))
   20 CONTINUE
      IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG)
      CALL PHASP
      NPG=0
      DO 21 I=1,NT
      IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 21
      NPG=NPG+1
      PV(1,I)=PCM(1,NPG)
      PV(2,I)=PCM(2,NPG)
      PV(3,I)=PCM(3,NPG)
      PV(4,I)=PCM(4,NPG)
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
      CALL LOR(I,MX6,I)
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
   21 CONTINUE
C**
C** LORENTZ TRANSFORMATION IN LAB SYSTEM
C**
   27 TARG=0.
      DO 36 I=1,NT
      IF(PV(5,I).GT.0.5) TARG=TARG+1.
      CALL LOR(I,MX2,I)
   36 CONTINUE
      IF(TARG.LT.0.5) TARG=1.
C**
C** SOMETIMES THE LEADING STRANGE PARTICLES ARE LOST , SET THEM BACK
C**
      IF(LEAD.EQ.0) GOTO 6085
      DO 6081 I=1,NT
      IF(ABS(IPA(I)).EQ.LEAD) GOTO 6085
 6081 CONTINUE
      I=1
      IF(LEAD.GE.14.AND.ABS(IPA(2)).GE.14) I=2
      IF(LEAD.LT.14.AND.ABS(IPA(2)).LT.14) I=2
      IPA(I)=LEAD
      EKIN=PV(4,I)-ABS(PV(5,I))
      PV(5,I)=RMASS(LEAD)
      PV(7,I)=1.
      IF(PV(5,I).LT.0.) PV(7,I)=-1.
      PV(5,I)=ABS(PV(5,I))
      PV(6,I)=RCHARG(LEAD)
      PV(4,I)=PV(5,I)+EKIN
      CALL LENGTX(I,PP)
      DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
      PP1=SQRT(DUMNVE)
C
      IF (PP .GE. 1.0E-6) GO TO 8000
      CALL GRNDM(RNDM,2)
      RTHNVE=PI*RNDM(1)
      PHINVE=TWPI*RNDM(2)
      PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
      PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
      PV(3,I)=PP1*COS(RTHNVE)
      GO TO 8001
 8000 CONTINUE
      PV(1,I)=PV(1,I)*PP1/PP
      PV(2,I)=PV(2,I)*PP1/PP
      PV(3,I)=PV(3,I)*PP1/PP
 8001 CONTINUE
C
C** FOR VARIOUS REASONS, THE ENERGY BALANCE IS NOT SUFFICIENT,
C** CHECK THAT,  ENERGY BALANCE, ANGLE OF FINAL SYSTEM E.T.C.
 6085 KGENEV=1
      PV(1,MX4)=0.
      PV(2,MX4)=0.
      PV(3,MX4)=P
      PV(4,MX4)=SQRT(P*P+AMAS*AMAS)
      PV(5,MX4)=ABS(AMAS)
      EKIN0=PV(4,MX4)-PV(5,MX4)
      PV(1,MX5)=0.
      PV(2,MX5)=0.
      PV(3,MX5)=0.
      PV(4,MX5)=MP*TARG
      PV(5,MX5)=PV(4,MX5)
      EKIN=PV(4,MX4)+PV(4,MX5)
      I=MX4
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
      I=MX5
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
      CALL ADD(MX4,MX5,MX6)
      CALL LOR(MX4,MX6,MX4)
      CALL LOR(MX5,MX6,MX5)
      TECM=PV(4,MX4)+PV(4,MX5)
      NPG=NT
      PV(1,MX8)=0.
      PV(2,MX8)=0.
      PV(3,MX8)=0.
      PV(4,MX8)=0.
      PV(5,MX8)=0.
      EKIN1=0.
      DO 598 I=1,NPG
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
      CALL ADD(MX8,I,MX8)
      EKIN1=EKIN1+PV(4,I)-PV(5,I)
      EKIN=EKIN-PV(5,I)
      IF(I.GT.18) GOTO 598
      AMASS(I)=PV(5,I)
  598 CONTINUE
      IF(NPG.GT.18) GOTO 597
      CALL PHASP
      EKIN=0.
      DO 599 I=1,NPG
      PV(1,MX7)=PCM(1,I)
      PV(2,MX7)=PCM(2,I)
      PV(3,MX7)=PCM(3,I)
      PV(4,MX7)=PCM(4,I)
      PV(5,MX7)=AMASS(I)
      CALL LOR(MX7,MX5,MX7)
  599 EKIN=EKIN+PV(4,MX7)-PV(5,MX7)
      CALL ANG(MX8,MX4,COST,TETA)
      IF(NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1,EKIN
C**
C** MAKE SHURE, THAT  KINETIC ENERGIES ARE CORRECT
C** THE 3. CLUSTER IS NOT PRODUCED WITHIN PROPER KINEMATICS!!!
C** EKIN= KINETIC ENERGY THEORETICALLY
C** EKIN1= KINETIC ENERGY SIMULATED
C**
  597 EKIN1=0.
      IF(EKIN1.EQ.0.) GOTO 600
      PV(1,MX7)=0.
      PV(2,MX7)=0.
      PV(3,MX7)=0.
      PV(4,MX7)=0.
      PV(5,MX7)=0.
      WGT=EKIN/EKIN1
      EKIN1=0.
      DO 602 I=1,NT
      EKIN=PV(4,I)-PV(5,I)
      EKIN=EKIN*WGT
      PV(4,I)=EKIN+PV(5,I)
      DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2)
      PP=SQRT(DUMNVE)
      CALL LENGTX(I,PP1)
C
      IF (PP1 .GE. 1.0E-6) GO TO 8002
      CALL GRNDM(RNDM,2)
      RTHNVE=PI*RNDM(1)
      PHINVE=TWPI*RNDM(2)
      PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE)
      PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE)
      PV(3,I)=PP*COS(RTHNVE)
      GO TO 8003
 8002 CONTINUE
      PV(1,I)=PV(1,I)*PP/PP1
      PV(2,I)=PV(2,I)*PP/PP1
      PV(3,I)=PV(3,I)*PP/PP1
 8003 CONTINUE
C
      EKIN1=EKIN1+EKIN
      CALL ADD(MX7,I,MX7)
  602 CONTINUE
      CALL ANG(MX7,MX4,COST,TETA)
      IF(NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1
C**
C** ROTATE IN DIRECTION OF Z-AXIS, SEE COMMENTS IN 'GENXPT'
C**
  600 PV(1,MX7)=0.
      PV(2,MX7)=0.
      PV(3,MX7)=0.
      PV(4,MX7)=0.
      PV(5,MX7)=0.
      DO 596 I=1,NT
  596 CALL ADD(MX7,I,MX7)
*          call rannor(ran1,ran2)
      CALL GRNDM(RNDM,2)
      RY=RNDM(1)
      RZ=RNDM(2)
      RX=6.283185*RZ
      A1=SQRT(-2.*LOG(RY))
      RAN1=A1*SIN(RX)
      RAN2=A1*COS(RX)
      PV(1,MX7)=PV(1,MX7)+RAN1*0.020*TARG
      PV(2,MX7)=PV(2,MX7)+RAN2*0.020*TARG
      CALL DEFS(MX4,MX7,MX8)
      PV(1,MX7)=0.
      PV(2,MX7)=0.
      PV(3,MX7)=0.
      PV(4,MX7)=0.
      PV(5,MX7)=0.
C     DO 595 I=1,NT
C     CALL TRAC(I,MX8,I)
C 595 CALL ADD(MX7,I,MX7)
C     CALL ANG(MX7,MX4,COST,TETA)
      IF(NPRT(4)) WRITE(NEWBCD,2003) TETA
C**
C** ROTATE IN DIRECTION OF PRIMARY PARTICLE
C**
      DEKIN=0.
      NPIONS=0
      EK1=0.
      DO 25 I=1,NT
      CALL DEFS1(I,MXGKPV-1,I)
      IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
      IF(ATNO2.LT.1.5) GOTO 25
      CALL LENGTX(I,PP)
      EKIN=PV(4,I)-ABS(PV(5,I))
      CALL NORMAL(RAN)
      EKIN=EKIN-CFA*(1.+0.5*RAN)
      IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
      CALL STEEQ(XXH,I)
      DEKIN=DEKIN+EKIN*(1.-XXH)
      EKIN=EKIN*XXH
      IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) NPIONS=NPIONS+1
      IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) EK1=EK1+EKIN
      PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
      PV(4,I)=EKIN+ABS(PV(5,I))
C
      IF (PP .GE. 1.0E-6) GO TO 8004
      CALL GRNDM(RNDM,2)
      RTHNVE=PI*RNDM(1)
      PHINVE=TWPI*RNDM(2)
      PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
      PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
      PV(3,I)=PP1*COS(RTHNVE)
      GO TO 8005
 8004 CONTINUE
      PV(1,I)=PV(1,I)*PP1/PP
      PV(2,I)=PV(2,I)*PP1/PP
      PV(3,I)=PV(3,I)*PP1/PP
 8005 CONTINUE
C
   25 CONTINUE
      EK1=0.
      IF(EK1.EQ.0.) GOTO 23
      IF(NPIONS.LE.0) GOTO 23
      DEKIN=1.+DEKIN/EK1
      DO 22 I=1,NT
      IF(ABS(IPA(I)).LT.7.OR.ABS(IPA(I)).GT.9) GOTO 22
      CALL LENGTX(I,PP)
      EKIN=PV(4,I)-ABS(PV(5,I))
      EKIN=EKIN*DEKIN
      IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6
      PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
      PV(4,I)=EKIN+ABS(PV(5,I))
C
      IF (PP .GE. 1.0E-6) GO TO 8006
      CALL GRNDM(RNDM,2)
      RTHNVE=PI*RNDM(1)
      PHINVE=TWPI*RNDM(2)
      PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE)
      PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE)
      PV(3,I)=PP1*COS(RTHNVE)
      GO TO 8007
 8006 CONTINUE
      PV(1,I)=PV(1,I)*PP1/PP
      PV(2,I)=PV(2,I)*PP1/PP
      PV(3,I)=PV(3,I)*PP1/PP
 8007 CONTINUE
C
   22 CONTINUE
   23 IGEN=0
      IF(ATNO2.LT.1.5) GOTO 40
C**
C** ADD BLACK TRACK PARTICLES
C**
      CALL HIGHAB(SPROB)
      TEX=ENP(1)
      SPALL=TARG
      IF(TEX.LT.0.001) GOTO 445
      BLACK=(1.5+1.25*TARG)*ENP(1)/(ENP(1)+ENP(3))
      CALL POISSO(BLACK,NBL)
      IF(NPRT(4))
     *WRITE(NEWBCD,3003) NBL,TEX
      IF(IFIX(TARG)+NBL.GT.ATNO2) NBL=ATNO2-TARG
      IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
      IF(NBL.LE.0) GOTO 445
      EKIN=TEX/NBL
      EKIN2=0.
      CALL STEEP(XX)
      DO 441 I=1,NBL
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.SPROB) GOTO 441
      IF(NT.EQ.MXGKPV-2) GOTO 441
      IF(EKIN2.GT.TEX) GOTO 443
      CALL GRNDM(RNDM,1)
      RAN1=RNDM(1)
      CALL NORMAL(RAN2)
      EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
      IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
      EKIN1=EKIN1*XX
      EKIN2=EKIN2+EKIN1
      IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
      IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
      IPA1=16
      PNRAT=1.-ZNO2/ATNO2
      CALL GRNDM(RNDM,3)
      IF(RNDM(1).GT.PNRAT) IPA1=14
      NT=NT+1
      SPALL=SPALL+1.
      COST=-1.0+RNDM(2)*2.0
      DUMNVE=1.0-COST*COST
      IF (DUMNVE .LT. 0.0) DUMNVE=0.0
      SINT=SQRT(DUMNVE)
      PHI=TWPI*RNDM(3)
      IPA(NT)=-IPA1
      SIDE(NT)=-4.
      PV(5,NT)=ABS(RMASS(IPA1))
      PV(6,NT)=RCHARG(IPA1)
      PV(7,NT)=1.
      PV(4,NT)=EKIN1+PV(5,NT)
      DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
      PP=SQRT(DUMNVE)
      PV(1,NT)=PP*SINT*SIN(PHI)
      PV(2,NT)=PP*SINT*COS(PHI)
      PV(3,NT)=PP*COST
  441 CONTINUE
  443 IF(ATNO2.LT.10.) GOTO 445
      IF(EK.GT.2.0) GOTO 445
      II=NT+1
      KK=0
      EKA=EK
      IF(EKA.GT.1.) EKA=EKA*EKA
      IF(EKA.LT.0.1) EKA=0.1
      IKA=3.6*EXP((ZNO2**2/ATNO2-35.56)/6.45)/EKA
      IF(IKA.LE.0) GO TO 445
      DO 444 I=1,NT
      II=II-1
      IF(IPA(II).NE.-14) GOTO 444
      IPA(II)=-16
      IPA1  = 16
      PV(5,II)=ABS(RMASS(IPA1))
      PV(6,II)=RCHARG(IPA1)
      KK=KK+1
      IF(KK.GT.IKA) GOTO 445
  444 CONTINUE
  445 TEX=ENP(3)
      IF(TEX.LT.0.001) GOTO 40
      BLACK=(1.5+1.25*TARG)*ENP(3)/(ENP(1)+ENP(3))
      CALL POISSO(BLACK,NBL)
      IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
      IF(NBL.LE.0) GOTO 40
      EKIN=TEX/NBL
      EKIN2=0.
      CALL STEEP(XX)
      IF(NPRT(4))
     *WRITE(NEWBCD,3004) NBL,TEX
      DO 442 I=1,NBL
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.SPROB) GOTO 442
      IF(NT.EQ.MXGKPV-2) GOTO 442
      IF(EKIN2.GT.TEX) GOTO 40
      CALL GRNDM(RNDM,1)
      RAN1=RNDM(1)
      CALL NORMAL(RAN2)
      EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
      IF(EKIN1.LT.0.0) EKIN1=-0.005*LOG(RAN1)
      EKIN1=EKIN1*XX
      EKIN2=EKIN2+EKIN1
      IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
      IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
      CALL GRNDM(RNDM,3)
      COST=-1.0+RNDM(1)*2.0
      DUMNVE=1.0-COST*COST
      IF (DUMNVE .LT. 0.0) DUMNVE=0.0
      SINT=SQRT(DUMNVE)
      PHI=TWPI*RNDM(2)
      RAN=RNDM(3)
      IPA(NT+1)=-30
      IF(RAN.GT.0.60) IPA(NT+1)=-31
      IF(RAN.GT.0.90) IPA(NT+1)=-32
      SIDE(NT+1)=-4.
      PV(5,NT+1)=(ABS(IPA(NT+1))-28)*MP
      SPALL=SPALL+PV(5,NT+1)*1.066
      IF(SPALL.GT.ATNO2) GOTO 40
      NT=NT+1
      PV(6,NT)=1.
      IF(IPA(NT).EQ.-32) PV(6,NT)=2.
      PV(7,NT)=1.
      PV(4,NT)=PV(5,NT)+EKIN1
      DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
      PP=SQRT(DUMNVE)
      PV(1,NT)=PP*SINT*SIN(PHI)
      PV(2,NT)=PP*SINT*COS(PHI)
      PV(3,NT)=PP*COST
  442 CONTINUE
C**
C** STORE ON EVENT COMMON
C**
   40 CALL GRNDM(RNDM,1)
      IF(RS.GT.(4.+RNDM(1)*1.)) GOTO 42
      DO 41 I=1,NT
      CALL LENGTX(I,ETB)
      IF(ETB.LT.P) GOTO 41
      ETF=P
      PV(4,I)=SQRT(PV(5,I)**2+ETF**2)
      DUMNVE=ETB
      IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10
      ETF=ETF/DUMNVE
      PV(1,I)=PV(1,I)*ETF
      PV(2,I)=PV(2,I)*ETF
      PV(3,I)=PV(3,I)*ETF
   41 CONTINUE
   42 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
      EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
      EKIN2=0.
      CALL TDELAY(TOF1)
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      TOF=TOF-TOF1*LOG(RAN)
      DO 44 I=1,NT
      EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
      IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I)
      PV(7,I)=TOF
      PV(8,I)=ABS(IPA(I))
      PV(9,I)=0.
   44 PV(10,I)=0.
      IF(NPRT(4)) WRITE(NEWBCD,2006) NT,EKIN,ENP(1),ENP(3),EKIN1,EKIN2
      INTCT=INTCT+1.
      NMODE=3
      IF(SPALL.LT.0.5.AND.ATNO2.GT.1.5) NMODE=14
      CALL SETCUR(NT)
      NTK=NTK+1
      IF(NT.EQ.1) GOTO 300
      DO 50 II=2,NT
      I=II-1
      IF(NTOT.LT.NSIZE/12) GOTO 43
      GO TO 9999
   43 CALL SETTRK(I)
   50 CONTINUE
 300  CONTINUE
      GO TO 9999
C**
C** IT IS NOT POSSIBLE TO PRODUCE A PROPER TWO CLUSTER FINAL STATE.
C** CONTINUE WITH QUASI ELASTIC SCATTERING
C**
   60 IF(NPRT(4)) WRITE(NEWBCD,2005)
      DO 61 I=3,MXGKCU
   61 IPA(I)=0
      IPA(1)=IPART
      IPA(2)=14
      IF(NFL.EQ.2) IPA(2)=16
      CALL TWOB(IPPP,NFL,AVERN)
      GO TO 9999
C
 2000 FORMAT(' *COHERT* CMS PARAMETERS OF FINAL STATE PARTICLES',
     $ ' AFTER ',I3,' TRIALS')
 2001 FORMAT(' *COHERT* TRACK',2X,I3,2X,10F8.2,2X,I3,2X,F3.0)
 2002 FORMAT(' *COHERT* MOMENTUM ',F8.3,' MASSES ',2F8.4,' RS ',F8.4)
 2003 FORMAT(' *COHERT* TETA,EKIN0,EKIN1,EKIN ',4F10.4)
 2004 FORMAT(' *COHERT* TECM,NPB,MASSES: ',F10.4,1X,I3,1X,8F10.4/
     $ 1H ,26X,15X,8F10.4)
 2005 FORMAT(' *COHERT* NUMBER OF FINAL STATE PARTICLES',
     $ ' LESS THAN 2 ==> CONTINUE WITH 2-BODY SCATTERING')
 2006 FORMAT(' *COHERT*  COMP.',1X,I5,1X,5F7.2)
 3001 FORMAT(' *COHERT* NUCLEAR EXCITATION ',I5,' PARTICLES PRODUCED',
     $ ' IN ADDITION TO',I5,' NORMAL PARTICLES')
 3003 FORMAT(' *COHERT* ',I3,' BLACK TRACK PARTICLES PRODUCED',
     $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV')
 3004 FORMAT(' *COHERT* ',I5,' HEAVY FRAGMENTS WITH TOTAL ENERGY OF ',
     $ F8.4,' GEV')
C
 9999 CONTINUE
      RETURN
      END