[u/mrichter/AliRoot.git] / GEANT321 / gheisha / genxpt.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:04  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.39  by  S.Giani
*-- Author :
      SUBROUTINE GENXPT(IPPP,NFL,AVERN)
C
C *** GENERATION OF X- AND PT- VALUES FOR ALL PRODUCED PARTICLES ***
C *** NVE 02-MAY-1988 CERN GENEVA ***
C
C ORIGIN : H.FESEFELDT 11-OCT-1987
C
C A SIMPLE SINGLE VARIABLE DESCRIPTION E D3S/DP3= F(Q) WITH
C Q**2 = (M*X)**2 + PT**2 IS USED. FINAL STATE KINEMATIC IS PRODUCED
C BY AN FF-TYPE ITERATIVE CASCADE METHOD
C
#include "geant321/s_defcom.inc"
#include "geant321/s_genio.inc"
C
      REAL MASPAR,LAMB,NUCSUP
      DIMENSION MASPAR(8),BP(8),PTEX(8),C1PAR(5),G1PAR(5),TAVAI(2),
     $          SIDE(MXGKCU),IAVAI(2),BINL(20),DNDL(20),TWSUP(8),
     $          NUCSUP(6),PSUP(6),IPAX(100)
      DIMENSION RNDM(3)
      DATA MASPAR/0.75,0.70,0.65,0.60,0.50,0.40,0.75,0.20/
      DATA     BP/3.50,3.50,3.50,6.00,5.00,4.00,3.50,3.50/
      DATA   PTEX/1.70,1.70,1.50,1.70,1.40,1.20,1.70,1.20/
      DATA  C1PAR/0.6,0.6,0.35,0.15,0.10/
      DATA  G1PAR/2.6,2.6,1.80,1.30,1.20/
      DATA BINL/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.11,1.25
     $         ,1.43,1.67,2.0,2.5,3.33,5.00,10.00/
      DATA TWSUP/1.,1.,0.7,0.5,0.3,0.2,0.1,0.0/
      DATA NUCSUP/1.00,0.7,0.5,0.4,0.35,0.3/
      DATA   PSUP/3.,6.,20.,50.,100.,1000./
C
C**
      CALL HIGSEL(ISEL)
      IF(ISEL.EQ.1) THEN
         CALL HIGXPT(IPPP,NFL,AVERN)
         RETURN
      ENDIF
C**
C**  FOR ANNIHILATION INTERACTIONS INTRODUCE PROPER KINEMATICS
C**
      CALL CORANH(NIHIL,NFL)
C**
C**
C** CHECK FIRST MASS-INDICES
C**
      EK=ENP(5)
      EN=ENP(6)
      P=ENP(7)
      S=ENP(8)
      RS=ENP(9)
      NT=0
      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)
      CALL UCOPY(IPA(1),IPAX(1),100)
C**
C** FOR LOW MULTIPLICITY USE TWO-BODY RESONANCE MODEL OR SINGLE/DOUBLE
C** DIFFRACTION MODEL (--> TWOCLU (--> TWOB (--> COSCAT)))
C**
      CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
      IF(NIHIL.GT.0) GOTO 200
      IF(NT.GE.8) GOTO 200
      IF(EK.LT.1.) GOTO 60
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      IF(IPART.GE.10.AND.IPART.LE.13.AND.RAN.LT.0.5) GOTO 200
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      WSUP=TWSUP(NT)
      IF(RAN.GT.WSUP) GOTO 200
   60 CALL UCOPY(IPAX,IPA,100)
      CALL TWOCLU(IPPP,NFL,AVERN)
      GO TO 9999
C**
C** SET EFFECTIVE 4-MOMENTUM OF PRIMARY PARTICLE
C**
  200 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
      MX9=MX+9
      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(49)=IER(49)+1
C**
C** SOME RANDOMISATION OF ORDER OF FINAL STATE PARTICLES
C**
      DO 201 I=3,NT
      CALL GRNDM(RNDM,1)
      IPX=IFIX(3.+RNDM(1)*(NT-2.))
      IF(IPX.GT.NT) IPX=NT
      IPA1=IPA(IPX)
      IPA(IPX)=IPA(I)
  201 IPA(I)  =IPA1
C**
C** DISTRIBUTE IN FORWARD AND BACKWARD HEMISPHERE IN CMS
C**
      SIDE(1)= 1.
      SIDE(2)=-1.
      NTB=1
      TARG=0.
      IF(IPART.LT.10.OR.IPART.GT.13) GOTO 53
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.0.7) GOTO 53
      IPA1=IPA(1)
      IPA(1)=IPA(2)
      IPA(2)=IPA1
   53 LEAD=0
      IF(IPART.LT.10.OR.IPART.EQ.14.OR.IPART.EQ.16) GOTO 532
      IPA1=ABS(IPA(1))
      IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 531
      LEAD=IPA1
      GOTO 532
  531 IPA1=ABS(IPA(2))
      IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 532
      LEAD=IPA1
  532 DO 3 I=1,NT
      IF(I.LE.2) GOTO 54
      SIDE(I)= 1.
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.0.5) SIDE(I)=-1.
      IF(SIDE(I).LT.-0.5) NTB=NTB+1
   54 CONTINUE
    3 CONTINUE
      TB=2.*NTB
      CALL GRNDM(RNDM,1)
      IF(RS.LT.(2.0+RNDM(1))) TB=(2.*NTB+NT)/2.
C**
C** ADD PARTICLES FROM INTRANUCLEAR CASCADE
C**
      AFC=0.312+0.200*LOG(LOG(S))+S**1.5/6000.
      IF(AFC.GT.0.75) AFC=0.75
      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.MX) GOTO 2
      NT2=MX
      NTARG=NT2-NT
    2 CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,3001) NTARG,NT
      NT1=NT+1
      IF(NTARG.EQ.0) GOTO 51
C**
C** CHECK NUMBER OF EXTRA NUCLEONS AND PIONS
C**
      DO 881 IPX=1,6
      IF(P.LE.PSUP(IPX)) GOTO 882
  881 CONTINUE
      IPX=6
  882 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** CHECK AVAILABLE KINETIC ENERGY, IN THIS MODEL CONSERVATION OF
C** KINETIC ENERGY IN FORWARD AND BACKWARD HEMISPHERE IS ASSUMED
C**
    6 IF(NT.LE.1) GOTO 60
      TAVAI(1)=RS/2.
      TAVAI(2)=(TARG+1.)*RS/2.
      IAVAI(1)=0
      IAVAI(2)=0
      DO 7 I=1,NT
      L=1
      IF(SIDE(I).LT.0.) L=2
      IAVAI(L)=IAVAI(L)+1
      TAVAI(L)=TAVAI(L)-ABS(PV(5,I))
    7 CONTINUE
      NTH=NT
      IF(NTH.GT.10) NTH=10
      IF (NPRT(4))
     $ WRITE(NEWBCD,3002) TAVAI,IAVAI,(IPA(I),SIDE(I),I=1,NTH)
      IF(IAVAI(1).LE.0) GOTO 60
      IF(IAVAI(2).LE.0) GOTO 60
      IF(TAVAI(1).GT.0.) GOTO 11
      CALL GRNDM(RNDM,1)
      ISKIP=IFIX(RNDM(1)*(IAVAI(1)-1))+1
      IS=0
      DO 10  I=1,NT
      II=NT-I+1
      IF(SIDE(II).LT.0.) GOTO 10
      IS=IS+1
      IF(IS.NE.ISKIP) GOTO 10
      IF(II.EQ.NT) GOTO 9
      NT1=II+1
      NT2=NT
      DO 8 J=NT1,NT2
      IPA(J-1)=IPA(J)
      SIDE(J-1)=SIDE(J)
      DO 71 K=1,10
   71 PV(K,J-1)=PV(K,J)
    8 CONTINUE
      GOTO 9
   10 CONTINUE
    9 IPA(NT)=0
      SIDE(NT)=0.
      NT=NT-1
      GOTO 6
   11 IF(TAVAI(2).GT.0.) GOTO 15
      CALL GRNDM(RNDM,1)
      ISKIP=IFIX(RNDM(1)*(IAVAI(2)-1))+1
      IS=0
      DO 14  I=1,NT
      II=NT-I+1
      IF(SIDE(II).GT.0.) GOTO 14
      IS=IS+1
      IF(IS.NE.ISKIP) GOTO 14
      IF(SIDE(II).LT.-1.5) NTARG=NTARG-1
      IF(NTARG.LT.0) NTARG=0
      IF(II.EQ.NT) GOTO 13
      NT1=II+1
      NT2=NT
      DO 12 J=NT1,NT2
      IPA(J-1)=IPA(J)
      SIDE(J-1)=SIDE(J)
      DO 74 K=1,10
   74 PV(K,J-1)=PV(K,J)
   12 CONTINUE
      GOTO 13
   14 CONTINUE
   13 IPA(NT)=0
      SIDE(NT)=0.
      NT=NT-1
      GOTO 6
   15 IF(NT.LE.1) GOTO 60
      IF(NT.EQ.MX) GOTO 29
      NT1=NT+1
      NT2=MX
      DO 28 I=NT1,NT2
   28 IPA(I)=0
   29 CONTINUE
C**
C** NOW THE PREPARATION IS FINISHED.
C** DEFINE INITIAL STATE VECTORS FOR LORENTZ TRANSFORMATIONS.
C**
      PV( 1,MX1)=0.
      PV( 2,MX1)=0.
      PV( 3,MX1)=P
      PV( 4,MX1)=SQRT(P*P+AMAS*AMAS)
      PV( 5,MX1)=ABS(AMAS)
      PV( 1,MX2)=0.
      PV( 2,MX2)=0.
      PV( 3,MX2)=0.
      PV( 4,MX2)=MP
      PV( 5,MX2)=MP
      PV( 1,MX4)=0.
      PV( 2,MX4)=0.
      PV( 3,MX4)=0.
      PV( 4,MX4)=MP*(1.+TARG)
      PV( 5,MX4)=PV(4,MX4)
      PV( 1,MX8)=0.
      PV( 2,MX8)=0.
      PV( 3,MX8)=0.
      PV( 1,MX9)=1.
      PV( 2,MX9)=0.
      PV( 3,MX9)=0.
      CALL ADD(MX1,MX2,MX3)
      CALL ADD(MX4,MX1,MX4)
      CALL LOR(MX1,MX3,MX1)
      CALL LOR(MX2,MX3,MX2)
C**
C** MAIN LOOP FOR 4-MOMENTUM GENERATION , SEE PITHA-REPORT (AACHEN)
C** FOR A DETAILED DESCRIPTION OF THE METHOD.
C**
      CALL GRNDM(RNDM,1)
      PHI=RNDM(1)*TWPI
      EKIN1=0.
      EKIN2=0.
      DO 39 J=1,10
      PV(J,MX5)=0.
   39 PV(J,MX6)=0.
      NPG=0
      TARG1=0.
      DO 16 III=1,NT
      I=NT-III+1
      IPA1=ABS(IPA(I))
C**
C** COUNT NUMBER OF BACKWARD NUCLEONS
C**
      IF(I.EQ.2) GOTO 301
      IF(SIDE(I).LT.-1.5.AND.IPA1.GE.14) GOTO 301
      GOTO 38
  301 NPG=NPG+1
      IF(NPG.GT.18) GOTO 38
      SIDE(I)=-3.
      TARG1=TARG1+1.
      GOTO 16
   38 J=3
      IF(IPA1.LT.14) J=2
      IF(IPA1.LT.10) J=1
      IF(I.LE.2) J=J+3
      IF(SIDE(I).LT.-1.5) J=7
      IF(J.EQ.7.AND.IPA1.GE.14) J=8
C**
C** SET PT - AND PHI VALUES, THEY ARE CHANGED SOMEWHAT IN THE ITERATION
C** LOOP, SET MASS PARAMETER FOR LAMBDA FRAGMENTATION MODEL
C**
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      BPP=BP(J)
      BPE=PTEX(J)
      PT2=-LOG(1.-RAN)/BPP
      ASPAR=MASPAR(J)
      PT2=PT2**BPE
      PT =SQRT(PT2)
      IF(PT.LT.0.001) PT=0.001
      PV(1,I)=PT*COS(PHI)
      PV(2,I)=PT*SIN(PHI)
      PV(10,I)=PT
      BINL(1)=0.
      RLMAX=1./PV(10,I)
      DO 73 J=2,20
   73 BINL(J)=RLMAX*(J-1)/19.
      ET=PV(4,MX1)
      IF(SIDE(I).LT.0.) THEN
         ET=PV(4,MX2)
      ENDIF
      DNDL(1)=0.
      NTRIAL=0
C**
C** START OF BIG ITERATION LOOP
C**
   30 NTRIAL=NTRIAL+1
      IF(NTRIAL.GT. 2) GOTO 169
      DO 17 L=2,20
      DNDL(L)=0.
      X=(BINL(L)+BINL(L-1))/2.
      IF(PV(10,I).LT.0.001) PV(10,I)=0.001
      IF(X.GT.1./PV(10,I)) GOTO 17
      DX=BINL(L)-BINL(L-1)
      DNDL(L)=ASPAR/SQRT((1.+(ASPAR*X)**2)**3)
      DNDL(L)=ET*DNDL(L)/SQRT((X*PV(10,I)*ET)**2+PV(10,I)**2
     *                             +PV(5,I)**2)
      DNDL(L)=DNDL(L)*DX
   17 DNDL(L)=DNDL(L-1)+DNDL(L)
      NTRI=0
   31 CALL GRNDM(RNDM,1)
      RAN=RNDM(1)*DNDL(20)
      DO 18 L=2,20
      IF(RAN.LT.DNDL(L)) GOTO 19
   18 CONTINUE
C**
C** START OF SMALL ITERATION LOOP
C**
   19 NTRI=NTRI+1
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      DX=BINL(L)-BINL(L-1)
      LAMB=BINL(L-1)+RAN*DX/2.
      X=PV(10,I)*LAMB
      IF(X.GT.1.) X=1.
      X=X*SIDE(I)/ABS(SIDE(I))
      PV(3,I)=X*ET
      PV(4,I)=PV(3,I)**2+PV(10,I)**2+PV(5,I)**2
      PV(4,I)=SQRT(PV(4,I))
      IF(SIDE(I).LT.0.) GOTO 165
      IF(I.GT.2) GOTO 20
      EKIN=TAVAI(1)-EKIN1
      CALL NORMAL(RAN)
      IF(EKIN.LT.0.) EKIN=0.04*ABS(RAN)
      PV(4,I)=ABS(PV(5,I))+EKIN
      RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
      PP=SQRT(RNVE)
      CALL LENGTX(I,PP1)
C
      IF (PP1 .GE. 1.0E-6) GO TO 8000
      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 8001
 8000 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
 8001 CONTINUE
C
      CALL ADD(MX5,I,MX5)
      GOTO 16
   20 EKIN=EKIN1+PV(4,I)-ABS(PV(5,I))
      IF(EKIN.LT.0.95*TAVAI(1)) GOTO 161
      IF(NTRI.GT. 5) GOTO 167
      PV(10,I)=PV(10,I)*0.9
      PV( 1,I)=PV( 1,I)*0.9
      PV( 2,I)=PV( 2,I)*0.9
      DNDL(20)=DNDL(20)*0.9
      IF((TAVAI(2)-ABS(PV(5,I))).LT.0.) GOTO 31
      SIDE(I)=-1.
      TAVAI(1)=TAVAI(1)+ABS(PV(5,I))
      TAVAI(2)=TAVAI(2)-ABS(PV(5,I))
      GOTO 31
  161 CALL ADD(MX5,I,MX5)
      EKIN1=EKIN1+PV(4,I)-ABS(PV(5,I))
      GOTO 163
  165 EKIN=EKIN2+PV(4,I)-ABS(PV(5,I))
      XXX=0.95+0.05*TARG/20.
      IF(XXX.GT.0.999) X=0.999
      IF(EKIN.LT.XXX*TAVAI(2)) GOTO 166
      IF(NTRI.GT. 5) GOTO 167
      PV(10,I)=PV(10,I)*0.9
      PV( 1,I)=PV( 1,I)*0.9
      PV( 2,I)=PV( 2,I)*0.9
      DNDL(20)=DNDL(20)*0.9
      IF((TAVAI(1)-ABS(PV(5,I))).LT.0.) GOTO 31
      SIDE(I)=+1.
      TAVAI(1)=TAVAI(1)-ABS(PV(5,I))
      TAVAI(2)=TAVAI(2)+ABS(PV(5,I))
      GOTO 31
  166 CALL ADD(MX6,I,MX6)
      EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
  163 CALL ADD(MX5,MX6,MX7)
      PV(3,MX7)=0.
      CALL ANG(MX7,MX9,COST,PHIS)
      IF(PV(2,MX7).LT.0.) PHIS=TWPI-PHIS
      CALL NORMAL(RAN)
      RAN=RAN*PI/12.
      PHI=PHIS+PI+RAN
      IF(PHI.GT.TWPI) PHI=PHI-TWPI
      IF(PHI.LT.0.) PHI=TWPI-PHI
      GOTO 16
C**
C** PARTICLE MOMENTUM ZERO, REDUCE KINETIC ENERGY OF ALL OTHER
C**
  167 EKIN1=0.
      EKIN2=0.
      DO 162 J=1,10
      PV(J,MX5)=0.
  162 PV(J,MX6)=0.
      II=I+1
      DO 168 L=II,NT
      IF(ABS(IPA(L)).GE.14.AND.SIDE(L).LT.0.) GOTO 168
      PV(4,L)=PV(4,L)*0.95+0.05*ABS(PV(5,L))
      IF(PV(4,L).LT.ABS(PV(5,L))) PV(4,L)=ABS(PV(5,L))
      RNVE=ABS(PV(4,L)**2-PV(5,L)**2)
      PP=SQRT(RNVE)
      CALL LENGTX(L,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,L)=PP*SIN(RTHNVE)*COS(PHINVE)
      PV(2,L)=PP*SIN(RTHNVE)*SIN(PHINVE)
      PV(3,L)=PP*COS(RTHNVE)
      GO TO 8003
 8002 CONTINUE
      PV(1,L)=PV(1,L)*PP/PP1
      PV(2,L)=PV(2,L)*PP/PP1
      PV(3,L)=PV(3,L)*PP/PP1
 8003 CONTINUE
C
      PV(10,L)=SQRT(PV(1,L)**2+PV(2,L)**2)
      IF(SIDE(L).LT.0.) GOTO 164
      EKIN1=EKIN1+PV(4,L)-ABS(PV(5,L))
      CALL ADD(MX5,L,MX5)
      GOTO 168
  164 EKIN2=EKIN2+PV(4,L)-ABS(PV(5,L))
      CALL ADD(MX6,L,MX6)
  168 CONTINUE
C *** NEXT STMT. CHANGED TO PREVENT FROM INFINITE LOOPING ***
C*************      GOTO 38
      GO TO 30
C**
C** SKIP PARTICLE, IF NOT ENOUGH ENERGY
C**
  169 IPA(I)=0
      DO 170 J=1,10
  170 PV(J,I)=0.
      GOTO 163
   16 CONTINUE
      NTRI=0
      II=0
      DO 320 I=1,NT
      IF(IPA(I).EQ.0) GOTO 320
      II=II+1
      IPA(II)=IPA(I)
      SIDE(II)=SIDE(I)
      DO 321 J=1,10
  321 PV(J,II)=PV(J,I)
  320 CONTINUE
      NT=II
C**
C** BACKWARD NUCLEONS PRODUCED WITH A CLUSTER MODEL
C**
      CALL LOR(MX4,MX3,MX7)
      CALL SUB(MX7,MX5,MX7)
      CALL SUB(MX7,MX6,MX7)
      IF(TARG1.GT.1.5) GOTO 310
  322 I=2
      CALL NORMAL(RAN)
      EKIN=TAVAI(2)-EKIN2
      EKINM=RS/2.-MP
      IF(EKIN.GT.EKINM) EKIN=EKINM
      CALL NORMAL(RAN)
      IF(EKIN.LT.0.04) EKIN=0.04*ABS(RAN)
      PV(4,I)=ABS(PV(5,I))+EKIN
      RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
      PP=SQRT(RNVE)
      CALL LENGTX(MX7,PP1)
C
      IF (PP1 .GE. 1.0E-6) GO TO 8004
      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 8005
 8004 CONTINUE
      PV(1,I)=PV(1,MX7)*PP/PP1
      PV(2,I)=PV(2,MX7)*PP/PP1
      PV(3,I)=PV(3,MX7)*PP/PP1
 8005 CONTINUE
C
      CALL ADD(MX6,I,MX6)
      GOTO 330
  310 ITARG1=IFIX(TARG1+0.1)
      IF(ITARG1.GT.5) ITARG1=5
      RMB0=0.
      NPG=0
      DO 311 I=1,NT
      IF(SIDE(I).GT.-2.5) GOTO 311
      NPG=NPG+1
      RMB0=RMB0+ABS(PV(5,I))
  311 CONTINUE
      IF(NPG.LT.2) GOTO 322
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      RMB=-LOG(1.-RAN)
      GPAR=G1PAR(ITARG1)
      CPAR=C1PAR(ITARG1)
      RMB=RMB0+RMB**CPAR/GPAR
      PV(5,MX7)=RMB
      IF(PV(5,MX7).GT.PV(4,MX7)) PV(5,MX7)=PV(4,MX7)
      RNVE=ABS(PV(4,MX7)**2-PV(5,MX7)**2)
      PP=SQRT(RNVE)
      CALL LENGTX(MX7,PP1)
C
      IF (PP1 .GE. 1.0E-6) GO TO 8006
      CALL GRNDM(RNDM,2)
      RTHNVE=PI*RNDM(1)
      PHINVE=TWPI*RNDM(2)
      PV(1,MX7)=PP*SIN(RTHNVE)*COS(PHINVE)
      PV(2,MX7)=PP*SIN(RTHNVE)*SIN(PHINVE)
      PV(3,MX7)=PP*COS(RTHNVE)
      GO TO 8007
 8006 CONTINUE
      PV(1,MX7)=PV(1,MX7)*PP/PP1
      PV(2,MX7)=PV(2,MX7)*PP/PP1
      PV(3,MX7)=PV(3,MX7)*PP/PP1
 8007 CONTINUE
C
      I=MX7
      IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
      PV(1,MX7)=-PV(1,MX7)
      PV(2,MX7)=-PV(2,MX7)
      PV(3,MX7)=-PV(3,MX7)
      KGENEV=1
      TECM=PV(5,MX7)
      NPG=0
      DO 312 I=1,NT
      IF(SIDE(I).GT.-2.5)GOTO 312
      NPG=NPG+1
      AMASS(NPG)=ABS(PV(5,I))
  312 CONTINUE
      CALL PHASP
      NPG=0
      DO 314 I=1,NT
      IF(SIDE(I).GT.-2.5) GOTO 314
      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)
      CALL LOR(I,MX7,I)
      CALL ADD(MX6,I,MX6)
  314 CONTINUE
  330 IF (NPRT(4))
     $ WRITE(NEWBCD,2002) NTRIAL,EKIN1,EKIN2,TAVAI(1),TAVAI(2)
  175 IF (.NOT.NPRT(4)) GOTO 36
      CALL ADD(MX5,MX6,MX7)
      EKIN1=PV(4,MX1)+PV(4,MX2)
      EKIN2=PV(4,MX5)+PV(4,MX6)
      WRITE(NEWBCD,2000) EKIN1,EKIN2
      I=MX1
      WRITE(NEWBCD,2001) I,(PV(J,I),J=1,4)
      I=MX2
      WRITE(NEWBCD,2001) I,(PV(J,I),J=1,4)
      I=MX5
      WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
      I=MX6
      WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5)
      DO 37 I=1,NT
   37 WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I)
C**
C** LORENTZ TRANSFORMATION IN LAB SYSTEM
C**
   36 IF(NT.LE.2) GOTO 60
      TARG=0.
      DO 601 I=1,NT
      IF(PV(5,I).GT.0.5) TARG=TARG+1.
      CALL LOR(I,MX2,I)
  601 CONTINUE
      IF(TARG.LT.0.5) TARG=1.
      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)
      RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
      PP1=SQRT(RNVE)
      PV(1,I)=PP1*PV(1,I)/PP
      PV(2,I)=PP1*PV(2,I)/PP
      PV(3,I)=PP1*PV(3,I)/PP
 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** EKIN= KINETIC ENERGY THEORETICALLY
C** EKIN1= KINETIC ENERGY SIMULATED
C**
  597 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)
      RNVE=ABS(PV(4,I)**2-PV(5,I)**2)
      PP=SQRT(RNVE)
      CALL LENGTX(I,PP1)
C
      IF (PP1 .GE. 1.0E-6) GO TO 8008
      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 8009
 8008 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
 8009 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, THIS DOES DISTURB IN SOME WAY OUR
C** INCLUSIVE DISTRIBUTIONS, BUT IT IS NESSACARY FOR MOMENTUM CONSER-
C** VATION.
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)
C**
C** SOME SMEARING IN TRANSVERSE DIRECTION FROM FERMI MOTION
C**
*          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.
      DO 595 I=1,NT
      CALL TRAC(I,MX8,I)
  595 CALL ADD(MX7,I,MX7)
      CALL ANG(MX7,MX4,COST,TETA)
      IF (NPRT(4)) WRITE(NEWBCD,2003) TETA
C**
C** ROTATE IN DIRECTION OF PRIMARY PARTICLE, SUBTRACT BINDING ENERGIES
C** AND MAKE SOME FURTHER CORRECTIONS IF REQUIRED (STEEP, STEEQ)
C**
      DEKIN=0.
      NPIONS=0
      EK1=0.
      DO 21 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 21
      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 8010
      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 8011
 8010 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
 8011 CONTINUE
C
   21 CONTINUE
      IF(EK1.EQ.0.) GOTO 23
      IF(NPIONS.EQ.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 8012
      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 8013
 8012 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
 8013 CONTINUE
C
   22 CONTINUE
C**
C** ADD BLACK TRACK PARTICLES, THE TOTAL NUMBER OF PARTICLES PRODUCED
C** IS RESTRICTED TO 198, THIS MAY HAVE INFLUENCE ON VERY HIGH ENERGY
C** FIRST PROTONS AND NEUTRONS
C**
   23 IF(ATNO2.LT.1.5) GOTO 40
      CALL SELFAB(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-10) NBL=MXGKPV-10-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-10) 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.+RNDM(2)*2.
      SINT=SQRT(ABS(1.-COST*COST))
      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)
      RNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
      PP=SQRT(RNVE)
      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
C**
C** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
C**
  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-10) NBL=MXGKPV-10-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-10) 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.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
      CALL GRNDM(RNDM,3)
      COST=-1.+RNDM(1)*2.
      SINT=SQRT(ABS(1.-COST*COST))
      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
      RNVE=ABS(PV(4,NT)**2-PV(5,NT)**2)
      PP=SQRT(RNVE)
      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))) 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)
      ETF=ETF/ETB
      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.
      CALL SETCUR(NT)
      NTK=NTK+1
      IF(NT.EQ.1) GO TO 9999
      DO 50 II=2,NT
      I=II-1
      IF(NTOT.LT.NSIZE/12) GOTO 43
      GO TO 9999
   43 CALL SETTRK(I)
   50 CONTINUE
C
 2002 FORMAT(' *GENXPT* PRODUCTION OF FINAL STATE KINEMATIC AFTER ',I3,
     $ ' TRIALS.  KINETIC ENERGIES ',2F6.2,' OUT OF ',2F6.2)
 2000 FORMAT(' *GENXPT* CMS PARAMETERS OF FINAL STATE PARTICLES,',
     $ ' ENERGIES IN INITIAL AND FINAL STATE ',2F6.2)
 2001 FORMAT(' *GENXPT* TRACK',2X,I3,2X,10F8.3,2X,I3,2X,F4.0)
 2003 FORMAT(' *GENXPT* TETA,EKIN0,EKIN1,EKIN ',4F10.4)
 2006 FORMAT(' *GENXPT* COMP.',1X,I5,1X,5F7.2)
 3001 FORMAT(' *GENXPT* NUCLEAR EXCITATION',I5,
     $ ' PARTICLES PRODUCED IN ADDITION  TO ',I5,' NORMAL PARTICLES')
 3002 FORMAT(' *GENXPT* AVAILABLE ENERGIES ',2F10.4,
     $ ' FOR ',2I3,' PARTICLES IN BEAM/TARGET FRAGM. REGION',
     $ ' WITH IPA/SIDE ARRAY '/
     $ 1H ,5X,10(I3,2X,F3.0,4X))
 3003 FORMAT(' *GENXPT* ',I3,' BLACK TRACK PARTICLES PRODUCED',
     $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV')
 3004 FORMAT(' *GENXPT* ',I5,' HEAVY FRAGMENTS PRODUCED',
     $ ' WITH TOTAL ENERGY OF',F8.4,' GEV')
C
 9999 CONTINUE
C
      END