Gheisha corrections suggested by Gary Bower (FNAL).

diff --git a/GEANT321/gheisha/caspim.F b/GEANT321/gheisha/caspim.F
index 022367b228aa937ff9cdc5f82b0c5c219d0f2dab..252c0d3b27ab5d8caf42afd39b8675528605c24f 100644 (file)
--- a/GEANT321/gheisha/caspim.F
+++ b/GEANT321/gheisha/caspim.F
@@ -145,10 +145,7 @@ C --- ELASTIC SCATTERING ---
       IF (NFL .EQ. 2) IPA(2)=16
       IF (INT .EQ. 2) GOTO 20
       GOTO 100
-C
-C --- CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION IN REACT.
  20   CONTINUE
-      IF (EAB .LE. RMASS(9)) GO TO 55
 C
 C --- SUPPRESSION OF HIGH MULTIPLICITY EVENTS AT LOW MOMENTUM ---
       IEAB=IFIX(EAB*5.0)+1
@@ -177,6 +174,9 @@ C --- P TARGET ---
       GO TO 100
 C
  23   CONTINUE
+C
+C --- CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION IN REACT.
+      IF (EAB .LE. RMASS(9)) GO TO 55
       N=1.0
 C
       IF (NFL .EQ. 1) GO TO 26

diff --git a/GEANT321/gheisha/caspip.F b/GEANT321/gheisha/caspip.F
index 486f82f74647c5637128bb5823fefafa7f30d47f..63abf0fc8ecc64bd6ed5e125fe71354dc7c938c0 100644 (file)
--- a/GEANT321/gheisha/caspip.F
+++ b/GEANT321/gheisha/caspip.F
@@ -37,8 +37,8 @@ C
       DIMENSION PMUL(2,1200),ANORM(2,60),SUPP(10),CECH(10),B(2)
       DIMENSION RNDM(1)
       SAVE PMUL,ANORM
-      DATA SUPP/0.,0.2,0.45,0.55,0.65,0.75,0.85,0.90,0.94,0.98/
-      DATA CECH/0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.10,0.09,0.07/
+      DATA SUPP/0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98/
+      DATA CECH/1.0,0.95,0.79,0.32,0.19,0.16,0.14,0.12,0.1,0.08/
       DATA B/0.7,0.7/,C/1.25/
 C
 C --- INITIALIZATION INDICATED BY KGINIT(18) ---
@@ -122,26 +122,32 @@ C**  ELASTIC SCATTERING
       IPA(2)=14
       IF(NFL.EQ.2) IPA(2)=16
       IF(INT.EQ.2) GOTO 20
-C**  FOR PI+ N REACTIONS CHANGE SOME OF THE ELASTIC CROSS SECTION
-C**  TO PI+ N --> PI0 P
-      IF(NFL.EQ.1) GOTO 100
-      IPLAB=IFIX(P   *5.)+1
-      IF(IPLAB.GT.10) IPLAB=10
-      CALL GRNDM(RNDM,1)
-      IF(RNDM(1).GT.CECH(IPLAB)/ATNO2**0.42) GOTO 100
-      IPA(1)=8
-      IPA(2)=14
       GOTO 100
-C**  CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION IN REACT.
-  20  IF (EAB .LE. RMASS(7)) GOTO 55
+  20  CONTINUE
 C**  SUPPRESSION OF HIGH MULTIPLICITY EVENTS AT LOW MOMENTUM
       IEAB=IFIX(EAB*5.)+1
       IF(IEAB.GT.10) GOTO 22
       CALL GRNDM(RNDM,1)
       IF(RNDM(1).LT.SUPP(IEAB)) GOTO 22
-      N=1.
-      GOTO (23,24),NFL
+C**  CHARGE EXCHANGE REACTION
+      IPLAB=IFIX(P   *5.)+1
+      IF(IPLAB.GT.10) IPLAB=10
+      CALL GRNDM(RNDM,1)
+      IF(RNDM(1).GT.CECH(IPLAB)) GOTO 23
+      IF(NFL .EQ. 2)GOTO 24
+      INT = 1
+      IPA(1) = 7
+      IPA(2) = 14
+      GOTO 100
+  24  CONTINUE
+      IPA(1)=8
+      IPA(2)=14
+      GOTO 100
  23   CONTINUE
+C**  CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION IN REACT.
+      IF(EAB .LE. RMASS(7))GOTO 55
+      N = 1.
+      IF(NFL .EQ. 2)GOTO 26
       TEST=-(1+B(1))**2/(2.0*C**2)
       IF (TEST .LE. EXPXL) TEST=EXPXL
       IF (TEST .GE. EXPXU) TEST=EXPXU
@@ -157,13 +163,13 @@ C**  SUPPRESSION OF HIGH MULTIPLICITY EVENTS AT LOW MOMENTUM
       NM=0
       NZ=0
       GOTO 50
- 24   CONTINUE
+ 26   CONTINUE
       TEST=-(1+B(2))**2/(2.0*C**2)
       IF (TEST .LE. EXPXL) TEST=EXPXL
       IF (TEST .GE. EXPXU) TEST=EXPXU
       W0=EXP(TEST)
       WP=EXP(TEST)
-      TEST=-(-1+B(2))**2/(2.0*C**2)
+      TEST=-(B(2)-1)**2/(2.0*C**2)
       IF (TEST .LE. EXPXL) TEST=EXPXL
       IF (TEST .GE. EXPXU) TEST=EXPXU
       WM=EXP(TEST)
@@ -178,7 +184,7 @@ C**  SUPPRESSION OF HIGH MULTIPLICITY EVENTS AT LOW MOMENTUM
       NP=1
       NM=0
       NZ=0
-      IF(RAN.LT.WP/WT) GOTO 50
+      IF(RAN.LT.WP/WT)GOTO 50
       NP=0
       NM=1
       NZ=0

diff --git a/GEANT321/gheisha/coranh.F b/GEANT321/gheisha/coranh.F
index 700fc8f487789738e47fdc67cf60692d798502d1..f97c8ac9461294d50d58b4eb696356852a64017f 100644 (file)
--- a/GEANT321/gheisha/coranh.F
+++ b/GEANT321/gheisha/coranh.F
@@ -24,6 +24,7 @@ C
       IF(IPA(1).GE.14) GO TO 9999
       IF(IPA(2).GE.14) GO TO 9999
       NIHIL=1
+      RETURN
 C**
 C**  DO NOT BE CONFUSED, THIS HAS NOTHING TO DO WITH RELATIVISTIC
 C**  KINEMATIC

diff --git a/GEANT321/gheisha/genxpt.F b/GEANT321/gheisha/genxpt.F
index 11c87f667210a4f3f9539e39ca73fe1d2498b701..0b5969b14e3ffa8ce9477e4f0868113f57d1935a 100644 (file)
--- a/GEANT321/gheisha/genxpt.F
+++ b/GEANT321/gheisha/genxpt.F
@@ -23,6 +23,7 @@ C
 #include "geant321/s_defcom.inc"
 #include "geant321/s_genio.inc"
 C
+      REAL ATNOL,ZNOL
       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),
@@ -151,6 +152,9 @@ C**
       TB=2.*NTB
       CALL GRNDM(RNDM,1)
       IF(RS.LT.(2.0+RNDM(1))) TB=(2.*NTB+NT)/2.
+      ZNOL = ZNO2
+      IF(NFL .EQ. 1)ZNOL = ZNOL -1
+      ATNOL = ATNO2 - 1
 C**
 C** ADD PARTICLES FROM INTRANUCLEAR CASCADE
 C**
@@ -175,17 +179,23 @@ C**
   881 CONTINUE
       IPX=6
   882 DO 4 I=NT1,NT2
+      IF(ATNOL .GT. 0.99)THEN
       CALL GRNDM(RNDM,1)
       RAN=RNDM(1)
       IF(RAN.LT.NUCSUP(IPX)) GOTO 52
+      ENDIF
       CALL GRNDM(RNDM,1)
       IPA(I)=-(7+IFIX(RNDM(1)*3.0))
       GOTO 4
    52 IPA(I)=-16
-      PNRAT=1.-ZNO2/ATNO2
+      PNRAT=1.-ZNOL/ATNOL
       CALL GRNDM(RNDM,1)
-      IF(RNDM(1).GT.PNRAT) IPA(I)=-14
+      IF(RNDM(1).GT.PNRAT)THEN
+      IPA(I)=-14
+      ZNOL = ZNOL -1
+      ENDIF
       TARG=TARG+1.
+      ATNOL = ATNOL -1
     4 SIDE(I)=-2.
       NT=NT2
 C**
@@ -651,6 +661,8 @@ C**
       IF(PV(5,I).GT.0.5) TARG=TARG+1.
       CALL LOR(I,MX2,I)
   601 CONTINUE
+      IF(ABS(AMAS) .GT. 0.5)TARG = TARG - 1.
+      IF(NIHIL .GT. 0)TARG = TARG + 2
       IF(TARG.LT.0.5) TARG=1.
       IF(LEAD.EQ.0) GOTO 6085
       DO 6081 I=1,NT
@@ -877,7 +889,7 @@ C**
       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(NBL.GT.ATNOL) NBL=ATNOL
       IF(NT+NBL.GT.MXGKPV-10) NBL=MXGKPV-10-NT
       IF(NBL.LE.0) GOTO 445
       EKIN=TEX/NBL
@@ -898,9 +910,13 @@ C**
       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
       IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
       IPA1=16
-      PNRAT=1.-ZNO2/ATNO2
+      PNRAT=1.-ZNOL/ATNOL
       CALL GRNDM(RNDM,3)
-      IF(RNDM(1).GT.PNRAT) IPA1=14
+      IF(RNDM(1).GT.PNRAT)THEN
+      IPA1=14
+      ZNOL = ZNOL -1
+      ENDIF
+      ATNOL = ATNOL - 1
       NT=NT+1
       SPALL=SPALL+1.
       COST=-1.+RNDM(2)*2.
@@ -932,6 +948,7 @@ C**
       IF(IPA(II).NE.-14) GOTO 444
       IPA(II)=-16
       IPA1  = 16
+      ZNOL = ZNOL + 1.
       PV(5,II)=ABS(RMASS(IPA1))
       PV(6,II)=RCHARG(IPA1)
       KK=KK+1
@@ -942,6 +959,8 @@ C** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
 C**
   445 TEX=ENP(3)
       IF(TEX.LT.0.001) GOTO 40
+      IF(ATNOL .LT. ZNOL + 1.)GOTO 40
+      IF(ZNOL .LT. 1.)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
@@ -951,6 +970,8 @@ C**
       CALL STEEP(XX)
       IF (NPRT(4)) WRITE(NEWBCD,3004) NBL,TEX
       DO 442 I=1,NBL
+      IF(ATNOL .LT. ZNOL + 1.)GOTO 40
+      IF(ZNOL .LT. 1.)GOTO 40
       CALL GRNDM(RNDM,1)
       IF(RNDM(1).LT.SPROB) GOTO 442
       IF(NT.EQ.MXGKPV-10) GOTO 442
@@ -970,12 +991,24 @@ C**
       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
+      ATNOL = ATNOL - 2.
+      ZNOL = ZNOL - 1.
+      IF(RAN.GT.0.60)THEN
+      IF(ATNOL .GT. ZNOL + 0.9)THEN
+      IPA(NT+1)=-31
+      ATNOL = ATNOL - 1.
+      IF(RAN.GT.0.90)THEN
+      IF( (ATNOL .GT. 0.9) .AND. (ZNOL .GT. 0.9))THEN
+      IPA(NT+1)=-32
+      ATNOL = ATNOL - 1.
+      ZNOL = ZNOL - 1.
+      ENDIF
+      ENDIF
+      ENDIF
+      ENDIF
       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.
@@ -990,18 +1023,7 @@ C**
 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
+   40 CONTINUE
    42 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
       EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
       EKIN2=0.

diff --git a/GEANT321/gheisha/twob.F b/GEANT321/gheisha/twob.F
index 6a3a2c55de90e8c3c0537e4205f4e87de0002a65..9e7b90baa9055397247590fdd50c9f056e05b679 100644 (file)
--- a/GEANT321/gheisha/twob.F
+++ b/GEANT321/gheisha/twob.F
@@ -23,6 +23,9 @@ C
 #include "geant321/s_defcom.inc"
       DIMENSION RNDM(3)
 C
+C
+      REAL ZNOL,ATNOL
+C
 C     DATA CB/3./
       DATA CB/0.01/
 C
@@ -33,6 +36,9 @@ C**
 C**  FOR DIFFRACTION SCATTERING ON HEAVY NUCLEI USE BETTER ROUTINE
 C**  "COSCAT"
 C
+      ZNOL = ZNO2
+      IF (NFL .EQ. 1)ZNOL = ZNOL - 1
+      ATNOL = ATNO2 - 1
       TARMAS=RMASS(14)
       IF (NFL .EQ. 2) TARMAS=RMASS(16)
       ENP(8)=RMASS(IPPP)**2+TARMAS**2+2.0*TARMAS*ENP(6)
@@ -195,12 +201,13 @@ C** THE REASON IS, THAT WE HAVE TO SIMULATE ALSO THE NUCLEAR REACTIONS
 C** AT LOW ENERGIES LIKE A(H,P)B, A(H,P P)B, A(H,N)B E.T.C.
 C**
   200 IF(ENP(1).LE.0.0001.AND.ENP(3).LE.0.0001) GOTO 40
+      IF(ATNOL .LT. 0.9)GOTO 445
       SPALL=0.
       TEX=ENP(1)
       IF(TEX.LT.0.0001) GOTO 445
       BLACK=TEX/0.02
       CALL POISSO(BLACK,NBL)
-      IF(NBL.GT.ATNO2) NBL=ATNO2
+      IF(NBL.GT.ATNOL) NBL=ATNOL
       IF(ENP(1).GT.0.0001.AND.NBL.LE.0) NBL=1
       IF (NPRT(4)) WRITE(NEWBCD,3003) NBL,TEX
       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
@@ -221,9 +228,13 @@ C**
       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
       IF(EKIN1.LT.0.) EKIN1=0.0001
       IPA1=16
-      PNRAT=1.-ZNO2/ATNO2
+      PNRAT=1.-ZNOL/ATNOL
       CALL GRNDM(RNDM,3)
-      IF(RNDM(1).GT.PNRAT) IPA1=14
+      IF(RNDM(1).GT.PNRAT)THEN
+      IPA1=14
+      ZNOL = ZNOL -1
+      ENDIF
+      ATNOL = ATNOL -1
       NT=NT+1
       SPALL=SPALL+1.
       COST=-1.+RNDM(2)*2.
@@ -263,11 +274,15 @@ C**
       NBL=IFIX(2.*LOG(ATNO2))
       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
       IF(NBL.LE.0) GOTO 40
+      IF(ATNOL .LT. ZNOL + 1)GOTO 40
+      IF(ZNOL .LT. 1.0)GOTO 40
       EKIN=TEX/NBL
       EKIN2=0.
       CALL STEEP(XX)
       IF (NPRT(4)) WRITE(NEWBCD,3004) NBL,TEX
       DO 442 I=1,NBL
+      IF(ATNOL .LT. ZNOL + 1)GOTO 40
+      IF(ZNOL .LT. 1.0)GOTO 40
       IF(NT.EQ.MXGKPV-2) GOTO 442
       IF(EKIN2.GT.TEX) GOTO 40
       CALL GRNDM(RNDM,1)
@@ -285,12 +300,24 @@ C**
       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
+      ATNOL = ATNOL -2
+      ZNOL = ZNOL -1
+      IF(RAN.GT.0.60)THEN
+      IF(ATNOL .GT. ZNOL + 0.9)THEN
+      ATNOL = ATNOL -1
+      IPA(NT+1)=-31
+      IF(RAN.GT.0.90)THEN
+      IF( (ATNOL .GT. 0.9) .AND. (ZNOL .GT. 0.9) )THEN
+      ATNOL = ATNOL -1
+      ZNOL = ZNOL -1
+      IPA(NT+1)=-32
+      ENDIF
+      ENDIF
+      ENDIF
+      ENDIF
       INVE=ABS(IPA(NT+1))
       PV(5,NT+1)=RMASS(INVE)
       SPALL=SPALL+PV(5,NT+1)*1.066
-      IF(SPALL.GT.ATNO2) GOTO 40
       NT=NT+1
       PV(6,NT)=RCHARG(INVE)
       PV(7,NT)=2.

diff --git a/GEANT321/gheisha/twoclu.F b/GEANT321/gheisha/twoclu.F
index 1ba38dba541eca56ede1f04c884c6bd503bc29e6..86342900c4d1f001ed1ef31986da58ba9e885f14 100644 (file)
--- a/GEANT321/gheisha/twoclu.F
+++ b/GEANT321/gheisha/twoclu.F
@@ -26,6 +26,8 @@ C
 #include "geant321/s_genio.inc"
 C
       REAL NUCSUP
+      INTEGER NIHIL
+      REAL RENORM,ZNOL,ATNOL
       DIMENSION SIDE(MXGKCU),C1PAR(5),G1PAR(5),NUCSUP(5)
       DIMENSION RNDM(3)
       DATA C1PAR/0.6,0.6,0.35,0.15,0.10/
@@ -44,6 +46,7 @@ C
       MX6=MX+6
       MX7=MX+7
       MX8=MX+8
+      CALL CORANH(NIHIL,NFL)
       EK=ENP(5)
       EN=ENP(6)
       P=ENP(7)
@@ -79,6 +82,9 @@ C**
 C** DISTRIBUTE PARTICLES IN FORWARD AND BACKWARD HEMISPHERE OF CMS
 C** OF THE HADRON NUCLEON INTERACTION
 C**
+      ZNOL = ZNO2
+      IF(NFL .EQ. 1)ZNOL = ZNOL -1
+      ATNOL = ATNO2 - 1 
       SIDE(1)= 1.
       SIDE(2)=-1.
       TARG=0.
@@ -121,9 +127,15 @@ C**
       IF(RNDM(1).GT.(10.-P)/6.) GOTO 3
       CALL GRNDM(RNDM,1)
       IF(RNDM(1).GT.ATNO2/300.) GOTO 3
+      IF(ATNOL .LT. 0.9) GOTO 3
       IPA(I)=14
       CALL GRNDM(RNDM,1)
-      IF(RNDM(1).GT.ZNO2/ATNO2) IPA(I)=16
+      IF(RNDM(1).GT.ZNOL/ATNOL)THEN
+      IPA(I)=16
+      ZNOL = ZNOL + 1
+      ENDIF
+      ZNOL = ZNOL -1
+      ATNOL = ATNOL -1
       TARG=TARG+1.
     3 CONTINUE
       TB=2.*IBACK
@@ -148,16 +160,22 @@ C**
       IPX=IFIX(P/3.)+1
       IF(IPX.GT.5) IPX=5
       DO 4 I=NT1,NT2
+      IF(ATNOL .GT. 0.99)THEN
       CALL GRNDM(RNDM,1)
       RAN=RNDM(1)
       IF(RAN.LT.NUCSUP(IPX)) GOTO 52
+      ENDIF
       CALL GRNDM(RNDM,1)
       IPA(I)=-(7+IFIX(RNDM(1)*3.0))
       GOTO 4
    52 IPA(I)=-16
-      PNRAT=1.-ZNO2/ATNO2
+      PNRAT=1.-ZNOL/ATNOL
       CALL GRNDM(RNDM,1)
-      IF(RNDM(1).GT.PNRAT) IPA(I)=-14
+      IF(RNDM(1).GT.PNRAT)THEN
+      IPA(I)=-14
+      ZNOL = ZNOL -1
+      ENDIF
+      ATNOL = ATNOL -1
       TARG=TARG+1.
     4 SIDE(I)=-2.
       NT=NT2
@@ -319,9 +337,6 @@ C**
 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**
@@ -329,6 +344,10 @@ 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
+      RENORM = 1. - EXP(B*(TACMIN - DUMNVE))
+      T=-1.0E10
+      CALL GRNDM(RNDM,1)
+      IF (B .NE. 0.0) T=LOG(1.-RENORM*RNDM(1))/B
       CTET=-(T-TACMIN)/DUMNVE
       CTET=1.0-2.0*CTET
       IF (CTET .GT. 1.0) CTET=1.0
@@ -461,6 +480,8 @@ C**
       IF(PV(5,I).GT.0.5) TARG=TARG+1.
       CALL LOR(I,MX2,I)
    36 CONTINUE
+      IF(ABS(AMAS) .GT. 0.5)TARG = TARG -1
+      IF(NIHIL .GT. 0)TARG = TARG + 2
       IF(TARG.LT.0.5) TARG=1.
 C**
 C** SOMETIMES THE LEADING STRANGE PARTICLES ARE LOST , SET THEM BACK
@@ -698,7 +719,7 @@ C**
       CALL POISSO(BLACK,NBL)
       IF(NPRT(4))
      *WRITE(NEWBCD,3003) NBL,TEX
-      IF(IFIX(TARG)+NBL.GT.ATNO2) NBL=ATNO2-TARG
+      IF(NBL.GT.ATNOL) NBL=ATNOL
       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
       IF(NBL.LE.0) GOTO 445
       EKIN=TEX/NBL
@@ -719,9 +740,13 @@ C**
       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
       IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6
       IPA1=16
-      PNRAT=1.-ZNO2/ATNO2
+      PNRAT=1.-ZNOL/ATNOL
       CALL GRNDM(RNDM,3)
-      IF(RNDM(1).GT.PNRAT) IPA1=14
+      IF(RNDM(1).GT.PNRAT)THEN
+      IPA1=14
+      ZNOL = ZNOL -1
+      ENDIF
+      ATNOL = ATNOL -1
       NT=NT+1
       SPALL=SPALL+1.
       COST=-1.0+RNDM(2)*2.0
@@ -762,6 +787,8 @@ C**
   444 CONTINUE
   445 TEX=ENP(3)
       IF(TEX.LT.0.001) GOTO 40
+      IF(ATNOL .LT. ZNOL + 1)GOTO 40
+      IF(ZNOL .LT. 1)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
@@ -772,6 +799,8 @@ C**
       IF(NPRT(4))
      *WRITE(NEWBCD,3004) NBL,TEX
       DO 442 I=1,NBL
+      IF(ATNOL .LT. ZNOL + 1)GOTO 40
+      IF(ZNOL .LT. 1)GOTO 40
       CALL GRNDM(RNDM,1)
       IF(RNDM(1).LT.SPROB) GOTO 442
       IF(NT.EQ.MXGKPV-2) GOTO 442
@@ -793,12 +822,24 @@ C**
       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
+      ATNOL = ATNOL -2
+      ZNOL = ZNOL -1
+      IF(RAN.GT.0.60)THEN
+      IF(ATNOL .GT. ZNOL + 0.9)THEN
+      IPA(NT+1)=-31
+      ATNOL = ATNOL -1
+      IF(RAN.GT.0.90)THEN
+      IF((ATNOL .GT. 0.9) .AND. (ZNOL .GT. 0.9))THEN
+      IPA(NT+1)=-32
+      ATNOL = ATNOL -1
+      ZNOL = ZNOL -1
+      ENDIF
+      ENDIF
+      ENDIF
+      ENDIF
       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.
@@ -813,20 +854,7 @@ C**
 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
+   40 CONTINUE
    42 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
       EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
       EKIN2=0.