*
* $Id: caspb.F,v 1.1.1.1 1995/10/24 10:21:01 cernlib Exp $
*
* $Log: caspb.F,v $
* Revision 1.1.1.1  1995/10/24 10:21:01  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.39  by  S.Giani
*-- Author :
      SUBROUTINE CASPB(K,INT,NFL)
C
C *** CASCADE OF ANTI PROTON ***
C *** NVE 04-MAY-1988 CERN GENEVA ***
C
C ORIGIN : H.FESEFELDT (13-SEP-1987)
C
C PB  UNDERGOES INTERACTION WITH NUCLEON WITHIN NUCLEUS.
C CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE PIONS/KAONS.
C IF NOT ASSUME NUCLEAR EXCITATION OCCURS AND INPUT PARTICLE
C IS DEGRADED IN ENERGY.    NO OTHER PARTICLES PRODUCED.
C IF REACTION IS POSSIBLE FIND CORRECT NUMBER OF PIONS/PROTONS/
C NEUTRONS PRODUCED USING AN INTERPOLATION TO MULTIPLICITY DATA.
C REPLACE SOME PIONS OR PROTONS/NEUTRONS BY KAONS OR STRANGE BARYONS
C ACCORDING TO AVERAGE MULTIPLICITY PER INELASTIC REACTIONS.
C
#include "geant321/mxgkgh.inc"
#include "geant321/s_consts.inc"
#include "geant321/s_curpar.inc"
#include "geant321/s_result.inc"
#include "geant321/s_prntfl.inc"
#include "geant321/s_kginit.inc"
#include "geant321/limits.inc"
C
      REAL N
      DIMENSION PMUL1(2,1200),PMUL2(2,400),ANORM1(2,60),ANORM2(2,60),
     $          SUPP(10),CECH(20),ANHL(29),B(2)
      DIMENSION RNDM(1)
      SAVE PMUL1,ANORM1,PMUL2,ANORM2
      DATA SUPP/0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98/
      DATA CECH/0.14,0.17,0.18,0.18,0.18,0.17,0.17,0.16,0.155,0.145,
     *          0.11,0.082,0.065,0.050,0.041,0.035,0.028,0.024,0.010
     *         ,0.0/
      DATA ANHL/1.00,1.00,1.00,1.00,1.0,1.00,1.0,1.00,1.00,0.90
     *         ,0.6,0.52,0.47,0.44,0.41,0.39,0.37,0.35,0.34,0.24
     *         ,0.19,0.15,0.12,0.10,0.09,0.07,0.06,0.05,0./
      DATA B/0.7,0.7/,C/1.25/
C
C --- INITIALIZATION INDICATED BY KGINIT(11) ---
      IF (KGINIT(11) .NE. 0) GO TO 10
      KGINIT(11)=1
C
C --- INITIALIZE PMUL AND ANORM ARRAYS ---
      DO 9000 J=1,1200
      DO 9001 I=1,2
      PMUL1(I,J)=0.0
      IF (J .LE. 400) PMUL2(I,J)=0.0
      IF (J .LE. 60) ANORM1(I,J)=0.0
      IF (J .LE. 60) ANORM2(I,J)=0.0
 9001 CONTINUE
 9000 CONTINUE
C
C** COMPUTE NORMALIZATION CONSTANTS
C** FOR P AS TARGET
C
      L=0
      DO 1 NP1=1,20
      NP=NP1-1
      NMM1=NP1-1
      IF(NMM1.LE.1) NMM1=1
      NPP1=NP1+1
      DO 1 NM1=NMM1,NPP1
      NM=NM1-1
      DO 1 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.1200) GOTO 1
      NT=NP+NM+NZ
      IF(NT.LE.0.OR.NT.GT.60) GOTO 1
      PMUL1(1,L)=PMLTPC(NP,NM,NZ,NT,B(1),C)
      ANORM1(1,NT)=ANORM1(1,NT)+PMUL1(1,L)
    1 CONTINUE
C** FOR N AS TARGET
      L=0
      DO 2 NP1=1,20
      NP=NP1-1
      NPP1=NP1+2
      DO 2 NM1=NP1,NPP1
      NM=NM1-1
      DO 2 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.1200) GOTO 2
      NT=NP+NM+NZ
      IF(NT.LE.0.OR.NT.GT.60) GOTO 2
      PMUL1(2,L)=PMLTPC(NP,NM,NZ,NT,B(2),C)
      ANORM1(2,NT)=ANORM1(2,NT)+PMUL1(2,L)
    2 CONTINUE
      DO 3 I=1,60
      IF(ANORM1(1,I).GT.0.) ANORM1(1,I)=1./ANORM1(1,I)
      IF(ANORM1(2,I).GT.0.) ANORM1(2,I)=1./ANORM1(2,I)
    3 CONTINUE
      IF(.NOT.NPRT(10)) GOTO 9
      WRITE(NEWBCD,2001)
      DO 4 NFL=1,2
      WRITE(NEWBCD,2002) NFL
      WRITE(NEWBCD,2003) (ANORM1(NFL,I),I=1,60)
      WRITE(NEWBCD,2003) (PMUL1(NFL,I),I=1,1200)
    4 CONTINUE
C** DO THE SAME FOR ANNIHILATION CHANNELS
C** FOR P AS TARGET
C
    9 L=0
      DO 5 NP1=1,20
      NP=NP1-1
      NM=NP
      DO 5 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.400) GOTO 5
      NT=NP+NM+NZ
      IF(NT.LE.1.OR.NT.GT.60) GOTO 5
      PMUL2(1,L)=PMLTPC(NP,NM,NZ,NT,B(1),C)
      ANORM2(1,NT)=ANORM2(1,NT)+PMUL2(1,L)
    5 CONTINUE
C** FOR N AS TARGET
      L=0
      DO 6 NP1=1,20
      NP=NP1-1
      NM=NP+1
      DO 6 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.400) GOTO 6
      NT=NP+NM+NZ
      IF(NT.LE.1.OR.NT.GT.60) GOTO 6
      PMUL2(2,L)=PMLTPC(NP,NM,NZ,NT,B(2),C)
      ANORM2(2,NT)=ANORM2(2,NT)+PMUL2(2,L)
    6 CONTINUE
      DO 7 I=1,60
      IF(ANORM2(1,I).GT.0.) ANORM2(1,I)=1./ANORM2(1,I)
      IF(ANORM2(2,I).GT.0.) ANORM2(2,I)=1./ANORM2(2,I)
    7 CONTINUE
      IF(.NOT.NPRT(10)) GOTO 10
      WRITE(NEWBCD,3001)
      DO 8 NFL=1,2
      WRITE(NEWBCD,3002) NFL
      WRITE(NEWBCD,3003) (ANORM2(NFL,I),I=1,60)
      WRITE(NEWBCD,3003) (PMUL2(NFL,I),I=1,400)
    8 CONTINUE
C** CHOOSE PROTON OR NEUTRON AS TARGET
   10 NFL=2
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.ZNO2/ATNO2) NFL=1
      TARMAS=RMASS(14)
      IF (NFL .EQ. 2) TARMAS=RMASS(16)
      S=AMASQ+TARMAS**2+2.0*TARMAS*EN
      RS=SQRT(S)
      ENP(8)=AMASQ+TARMAS**2+2.0*TARMAS*ENP(6)
      ENP(9)=SQRT(ENP(8))
      EAB=RS-TARMAS-ABS(RMASS(15))
C**  ELASTIC SCATTERING
      NCECH=0
      NP=0
      NM=0
      NZ=0
      N=0.
      IF(INT.EQ.2) GOTO 20
C** INTRODUCE CHARGE EXCHANGE REACTION PB P --> NB N
      IF(NFL.EQ.2) GOTO 100
      IPLAB=IFIX(P*10.)+1
      IF(IPLAB.GT.10) IPLAB=IFIX(P)+10
      IF(IPLAB.GT.20) IPLAB=20
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).GT.CECH(IPLAB)/ATNO2**0.75) GOTO 100
      NCECH=1
      GOTO 100
C** ANNIHILATION CHANNELS
   20 IPLAB=IFIX(P*10.)+1
      IF(IPLAB.GT.10) IPLAB=IFIX(P)+10
      IF(IPLAB.GT.19) IPLAB=IFIX(P/10.)+19
      IF(IPLAB.GT.28) IPLAB=29
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).GT.ANHL(IPLAB)) GOTO 19
      EAB=RS
      IF (EAB .LE. 2.0*RMASS(7)) GOTO 55
      GOTO 222
C**  CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION IN REACT.
   19 IF (EAB .LE. RMASS(7)) GOTO 55
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 (24,23),NFL
 23   CONTINUE
      TEST=-(1+B(1))**2/(2.0*C**2)
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      W0=EXP(TEST)
      TEST=-(-1+B(1))**2/(2.0*C**2)
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      WM=EXP(TEST)
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      NP=0
      NM=0
      NZ=1
      IF(RAN.LT.W0/(W0+WM)) GOTO 100
      NP=0
      NM=1
      NZ=0
      GOTO 100
 24   CONTINUE
      TEST=-(1+B(2))**2/(2.0*C**2)
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      W0=EXP(TEST)
      WP=EXP(TEST)
      TEST=-(-1+B(2))**2/(2.0*C**2)
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      WM=EXP(TEST)
      WT=W0+WP+WM
      WP=W0+WP
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      NP=0
      NM=0
      NZ=1
      IF(RAN.LT.W0/WT) GOTO 100
      NP=1
      NM=0
      NZ=0
      IF(RAN.LT.WP/WT) GOTO 100
      NP=0
      NM=1
      NZ=0
      GOTO 100
C
   22 ALEAB=LOG(EAB)
C** NO. OF TOTAL PARTICLES VS SQRT(S)-2*MP
      N=3.62567+0.665843*ALEAB+0.336514*ALEAB*ALEAB
     * +0.117712*ALEAB*ALEAB*ALEAB+0.0136912*ALEAB*ALEAB*ALEAB*ALEAB
      N=N-2.
C** NORMALIZATION CONSTANT FOR  KNO-DISTRIBUTION
      ANPN=0.
      DO 21 NT=1,60
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=PI*NT/(2.0*N*N)
      DUM2=ABS(DUM1)
      DUM3=EXP(TEST)
      ADDNVE=0.0
      IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
      IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GT. 1.0E-10)) ADDNVE=DUM1*DUM3
      ANPN=ANPN+ADDNVE
   21 CONTINUE
      ANPN=1./ANPN
C** P OR N AS TARGET
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      EXCS=0.
      GOTO (30,40),NFL
C** FOR P AS TARGET
   30 L=0
      DO 31 NP1=1,20
      NP=NP1-1
      NMM1=NP1-1
      IF(NMM1.LE.1) NMM1=1
      NPP1=NP1+1
      DO 31 NM1=NMM1,NPP1
      NM=NM1-1
      DO 31 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.1200) GOTO 31
      NT=NP+NM+NZ
      IF(NT.LE.0.OR.NT.GT.60) GOTO 31
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=ANPN*PI*NT*PMUL1(1,L)*ANORM1(1,NT)/(2.0*N*N)
      DUM2=ABS(DUM1)
      DUM3=EXP(TEST)
      ADDNVE=0.0
      IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
      IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GT. 1.0E-10)) ADDNVE=DUM1*DUM3
      EXCS=EXCS+ADDNVE
      IF(RAN.LT.EXCS) GOTO 100
   31 CONTINUE
      GOTO 80
C** FOR N AS TARGET
   40 L=0
      DO 41 NP1=1,20
      NP=NP1-1
      NPP1=NP1+2
      DO 41 NM1=NP1,NPP1
      NM=NM1-1
      DO 41 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.1200) GOTO 41
      NT=NP+NM+NZ
      IF(NT.LE.0.OR.NT.GT.60) GOTO 41
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=ANPN*PI*NT*PMUL1(2,L)*ANORM1(2,NT)/(2.0*N*N)
      DUM2=ABS(DUM1)
      DUM3=EXP(TEST)
      ADDNVE=0.0
      IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
      IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GT. 1.0E-10)) ADDNVE=DUM1*DUM3
      EXCS=EXCS+ADDNVE
      IF(RAN.LT.EXCS) GOTO 100
   41 CONTINUE
      GOTO 80
C** ANNIHILATION CHANNELS
  222 IPA(1)=0
      IPA(2)=0
      ALEAB=LOG(EAB)
C** NO. OF TOTAL PARTICLES VS SQRT(S)
      N=3.62567+0.665843*ALEAB+0.336514*ALEAB*ALEAB
     * +0.117712*ALEAB*ALEAB*ALEAB+0.0136912*ALEAB*ALEAB*ALEAB*ALEAB
      N=N-2.
C** NORMALIZATION CONSTANT FOR  KNO-DISTRIBUTION
      ANPN=0.
      DO 221 NT=2,60
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=PI*NT/(2.0*N*N)
      DUM2=ABS(DUM1)
      DUM3=EXP(TEST)
      ADDNVE=0.0
      IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
      IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GT. 1.0E-10)) ADDNVE=DUM1*DUM3
      ANPN=ANPN+ADDNVE
  221 CONTINUE
      ANPN=1./ANPN
C** P OR N AS TARGET
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      EXCS=0.
      GOTO (230,240),NFL
C** FOR P AS TARGET
  230 L=0
      DO 231 NP1=1,20
      NP=NP1-1
      NM=NP
      DO 231 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.400) GOTO 231
      NT=NP+NM+NZ
      IF(NT.LE.1.OR.NT.GT.60) GOTO 231
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=ANPN*PI*NT*PMUL2(1,L)*ANORM2(1,NT)/(2.0*N*N)
      DUM2=ABS(DUM1)
      DUM3=EXP(TEST)
      ADDNVE=0.0
      IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
      IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GT. 1.0E-10)) ADDNVE=DUM1*DUM3
      EXCS=EXCS+ADDNVE
      IF(RAN.LT.EXCS) GOTO 120
  231 CONTINUE
      GOTO 80
C** FOR N AS TARGET
  240 L=0
      DO 241 NP1=1,20
      NP=NP1-1
      NM=NP+1
      DO 241 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF(L.GT.400) GOTO 241
      NT=NP+NM+NZ
      IF(NT.LE.1.OR.NT.GT.60) GOTO 241
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=ANPN*PI*NT*PMUL2(2,L)*ANORM2(2,NT)/(2.0*N*N)
      DUM2=ABS(DUM1)
      DUM3=EXP(TEST)
      ADDNVE=0.0
      IF (DUM2 .GE. 1.0) ADDNVE=DUM1*DUM3
      IF ((DUM2 .LT. 1.0) .AND. (DUM3 .GT. 1.0E-10)) ADDNVE=DUM1*DUM3
      EXCS=EXCS+ADDNVE
      IF(RAN.LT.EXCS) GOTO 120
  241 CONTINUE
      GOTO 80
   50 IF(NPRT(4))
     *WRITE(NEWBCD,1003) EAB,N,NFL,NP,NM,NZ
      CALL STPAIR
      IF(INT.EQ.1) CALL TWOB(15,NFL,N)
      IF(INT.EQ.2) CALL GENXPT(15,NFL,N)
      GO TO 9999
   55 IF(NPRT(4))
     *WRITE(NEWBCD,1001)
      GOTO 53
C** EXCLUSIVE REACTION NOT FOUND,ASSUME ELASTIC SCATTERING
   80 IF(NPRT(4))
     *WRITE(NEWBCD,1004)EAB,N
   53 INT=1
      NP=0
      NM=0
      NZ=0
  100 DO 101 I=1,60
  101 IPA(I)=0
      IF(INT.LE.0) GOTO 131
      GOTO (112,102),NFL
  102 GOTO (103,104),INT
  103 IPA(1)=15
      IPA(2)=16
      NT=2
      GOTO 130
  104 IF(NP.EQ.-1+NM) GOTO 105
      IF(NP.EQ.   NM) GOTO 106
      IPA(1)=17
      IPA(2)=14
      GOTO 120
  105 IPA(1)=15
      IPA(2)=14
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.0.5) GOTO 120
      IPA(1)=17
      IPA(2)=16
      GOTO 120
  106 IPA(1)=15
      IPA(2)=16
      GOTO 120
  112 GOTO (113,114),INT
  113 IPA(1)=15
      IPA(2)=14
      NT=2
      IF(NCECH.EQ.0) GOTO 130
      IPA(1)=17
      IPA(2)=16
      GOTO 130
  114 IF(NP.EQ.  NM) GOTO 115
      IF(NP.EQ.1+NM) GOTO 116
      IPA(1)=17
      IPA(2)=14
      GOTO 120
  115 IPA(1)=17
      IPA(2)=16
      CALL GRNDM(RNDM,1)
      IF(RNDM(1).LT.0.33) GOTO 120
      IPA(1)=15
      IPA(2)=14
      GOTO 120
  116 IPA(1)=15
      IPA(2)=16
  120 NT=2
      IF(NP.EQ.0) GOTO 122
      DO 121 I=1,NP
      NT=NT+1
  121 IPA(NT)=7
  122 IF(NM.EQ.0) GOTO 124
      DO 123 I=1,NM
      NT=NT+1
  123 IPA(NT)=9
  124 IF(NZ.EQ.0) GOTO 130
      DO 125 I=1,NZ
      NT=NT+1
  125 IPA(NT)=8
  130 IF(NPRT(4))
     *WRITE(NEWBCD,2004) NT,(IPA(I),I=1,20)
      GOTO 50
  131 IF(NPRT(4))
     *WRITE(NEWBCD,2005)
C
1001  FORMAT('0*CASPB* CASCADE ENERGETICALLY NOT POSSIBLE',
     $ ' CONTINUE WITH QUASI-ELASTIC SCATTERING')
1003  FORMAT(' *CASPB* ANTIPROTON-INDUCED CASCADE,',
     $ ' AVAIL. ENERGY',2X,F8.4,
     $ 2X,'<NTOT>',2X,F8.4,2X,'FROM',4(2X,I3),2X,'PARTICLES')
1004  FORMAT(' *CASPB* ANTIPROTON-INDUCED CASCADE,',
     $ ' EXCLUSIVE REACTION',
     $ ' NOT FOUND  TRY ELASTIC SCATTERING  AVAIL. ENERGY',2X,F8.4,2X,
     $ ' <NTOT>',2X,F8.4)
2001  FORMAT('0*CASPB* TABLES FOR MULT. DATA ANTIPROTON INDUCED ',
     $ 'REACTION  FOR DEFINITION OF NUMBERS SEE FORTRAN CODING')
2002  FORMAT(' *CASPB* TARGET PARTICLE FLAG',2X,I5)
2003  FORMAT(1H ,10E12.4)
2004  FORMAT(' *CASPB* ',I3,2X,'PARTICLES , MASS INDEX ARRAY',2X,20I4)
2005  FORMAT(' *CASPB* NO PARTICLES PRODUCED')
3001  FORMAT('0*CASPB* TABLES FOR MULT. DATA ANTIPROTON INDUCED ',
     $ ' ANNIHILATION REACTION  FOR DEFINITION OF NUMBERS SEE FORTRAN',
     $ ' CODING')
3002  FORMAT(' *CASPB* TARGET PARTICLE FLAG',2X,I5)
3003  FORMAT(1H ,10E12.4)
C
 9999 CONTINUE
      END