*
* $Id: caspim.F,v 1.1.1.1 1995/10/24 10:21:00 cernlib Exp $
*
* $Log: caspim.F,v $
* Revision 1.1.1.1  1995/10/24 10:21:00  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
*-- Author :
      SUBROUTINE CASPIM(K,INT,NFL)
C
C *** CASCADE OF PI- ***
C *** NVE 04-MAY-1988 CERN GENEVA ***
C
C ORIGIN : H.FESEFELDT 13-SEP-1987
C
C PI-  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/limits.inc"
#include "geant321/s_kginit.inc"
C
      REAL N
      DIMENSION PMUL(2,1200),ANORM(2,60),SUPP(10),CECH(10),B(2)
      DIMENSION RNDM(1)
      SAVE PMUL,ANORM
      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.95,0.79,0.32,0.19,0.16,0.14,0.12,0.10,0.08/
      DATA B/0.7,0.7/,C/1.25/
C
C --- INITIALIZATION INDICATED BY KGINIT(16) ---
      IF (KGINIT(16) .NE. 0) GO TO 10
      KGINIT(16)=1
C
C --- INITIALIZE PMUL AND ANORM ARRAYS ---
      DO 9000 J=1,1200
      DO 9001 I=1,2
      PMUL(I,J)=0.0
      IF (J .LE. 60) ANORM(I,J)=0.0
 9001 CONTINUE
 9000 CONTINUE
C
C *** COMPUTATION OF NORMALIZATION CONSTANTS ***
C
C --- P TARGET ---
      L=0
      DO 1100 NP1=1,20
      NP=NP1-1
      NMM1=NP1-1
      IF (NMM1 .LE. 1) NMM1=1
      NPP1=NP1+1
C
      DO 1101 NM1=NMM1,NPP1
      NM=NM1-1
C
      DO 1102 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF (L .GT. 1200) GOTO 1199
      NT=NP+NM+NZ
      IF (NT .LE. 0) GO TO 1102
      IF (NT .GT. 60) GO TO 1102
      PMUL(1,L)=PMLTPC(NP,NM,NZ,NT,B(1),C)
      ANORM(1,NT)=ANORM(1,NT)+PMUL(1,L)
 1102 CONTINUE
C
 1101 CONTINUE
C
 1100 CONTINUE
C
 1199 CONTINUE
C
C --- N TARGET ---
      L=0
      DO 1200 NP1=1,20
      NP=NP1-1
      NPP1=NP1+2
C
      DO 1201 NM1=NP1,NPP1
      NM=NM1-1
C
      DO 1202 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF (L .GT. 1200) GO TO 1299
      NT=NP+NM+NZ
      IF (NT .LE. 0) GO TO 1202
      IF (NT .GT. 60) GO TO 1202
      PMUL(2,L)=PMLTPC(NP,NM,NZ,NT,B(2),C)
      ANORM(2,NT)=ANORM(2,NT)+PMUL(2,L)
 1202 CONTINUE
C
 1201 CONTINUE
C
 1200 CONTINUE
C
 1299 CONTINUE
C
      DO 3 I=1,60
      IF (ANORM(1,I) .GT. 0.0) ANORM(1,I)=1.0/ANORM(1,I)
      IF (ANORM(2,I) .GT. 0.0) ANORM(2,I)=1.0/ANORM(2,I)
    3 CONTINUE
C
      IF (.NOT. NPRT(10)) GO TO 10
      WRITE(NEWBCD,2001)
      DO 4 NFL=1,2
      WRITE(NEWBCD,2002) NFL
      WRITE(NEWBCD,2003) (ANORM(NFL,I),I=1,60)
      WRITE(NEWBCD,2003) (PMUL(NFL,I),I=1,1200)
    4 CONTINUE
C
C --- CHOOSE PROTON OR NEUTRON AS TARGET ---
 10   CONTINUE
      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-RMASS(9)
C
C --- ELASTIC SCATTERING ---
      NP=0
      NM=0
      NZ=0
      N=0.0
      IPA(1)=9
      IPA(2)=14
      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
      IF (IEAB .GT. 10) GO TO 22
      CALL GRNDM(RNDM,1)
      IF (RNDM(1) .LT. SUPP(IEAB)) GO TO 22
C
C --- CHARGE EXCHANGE REACTION (IS INCLUDED IN INELASTIC CROSS SECTION)
      IPLAB=IFIX(P*5.0)+1
      IF (IPLAB .GT. 10) IPLAB=10
      CALL GRNDM(RNDM,1)
      IF (RNDM(1) .GT. CECH(IPLAB)) GO TO 23
C
      IF (NFL .EQ. 1) GOTO 24
C
C --- N TARGET ---
      INT=1
      IPA(1)=9
      IPA(2)=16
      GO TO 100
C
C --- P TARGET ---
 24   CONTINUE
      IPA(1)=8
      IPA(2)=16
      GO TO 100
C
 23   CONTINUE
      N=1.0
C
      IF (NFL .EQ. 1) GO TO 26
C
C --- N TARGET ---
      DUM=-(1+B(2))**2/(2.0*C**2)
      IF (DUM .LT. EXPXL) DUM=EXPXL
      IF (DUM .GT. EXPXU) DUM=EXPXU
      W0=EXP(DUM)
      DUM=-(-1+B(2))**2/(2.0*C**2)
      IF (DUM .LT. EXPXL) DUM=EXPXL
      IF (DUM .GT. EXPXU) DUM=EXPXU
      WM=EXP(DUM)
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      NP=0
      NM=0
      NZ=1
      IF (RAN .LT. W0/(W0+WM)) GO TO 50
      NP=0
      NM=1
      NZ=0
      GO TO 50
C
C --- P TARGET ---
 26   CONTINUE
      DUM=-(1+B(1))**2/(2.0*C**2)
      IF (DUM .LT. EXPXL) DUM=EXPXL
      IF (DUM .GT. EXPXU) DUM=EXPXU
      W0=EXP(DUM)
      WP=EXP(DUM)
      DUM=-(-1+B(1))**2/(2.0*C**2)
      IF (DUM .LT. EXPXL) DUM=EXPXL
      IF (DUM .GT. EXPXU) DUM=EXPXU
      WM=EXP(DUM)
      WP=WP*10.
      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) GO TO 50
      NP=1
      NM=0
      NZ=0
      IF (RAN .LT. WP/WT) GO TO 50
      NP=0
      NM=1
      NZ=0
      GOTO 50
C
 22   CONTINUE
      ALEAB=LOG(EAB)
C
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.0
C
C --- NORMALIZATION CONSTANT FOR KNO-DISTRIBUTION ---
      ANPN=0.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 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
      ANPN=ANPN+ADDNVE
   21 CONTINUE
      ANPN=1.0/ANPN
C
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      EXCS=0.0
      IF (NFL .EQ. 2) GO TO 40
C
C --- P TARGET ---
      L=0
      DO 310 NP1=1,20
      NP=NP1-1
      NMM1=NP1-1
      IF (NMM1 .LE. 1) NMM1=1
      NPP1=NP1+1
C
      DO 311 NM1=NMM1,NPP1
      NM=NM1-1
C
      DO 312 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF (L .GT. 1200) GO TO 80
      NT=NP+NM+NZ
      IF (NT .LE. 0) GO TO 312
      IF (NT .GT. 60) GO TO 312
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=ANPN*PI*NT*PMUL(1,L)*ANORM(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 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
      EXCS=EXCS+ADDNVE
      IF (RAN .LT. EXCS) GOTO 50
 312  CONTINUE
C
 311  CONTINUE
C
 310  CONTINUE
      GOTO 80
C
C --- N TARGET ---
 40   CONTINUE
      L=0
      DO 410 NP1=1,20
      NP=NP1-1
      NPP1=NP1+2
C
      DO 411 NM1=NP1,NPP1
      NM=NM1-1
C
      DO 412 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF (L .GT. 1200) GO TO 80
      NT=NP+NM+NZ
      IF (NT .LE. 0) GO TO 412
      IF (NT .GT. 60) GO TO 412
      TEST=-(PI/4.0)*(NT/N)**2
      IF (TEST .LT. EXPXL) TEST=EXPXL
      IF (TEST .GT. EXPXU) TEST=EXPXU
      DUM1=ANPN*PI*NT*PMUL(2,L)*ANORM(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 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
      EXCS=EXCS+ADDNVE
      IF (RAN .LT. EXCS) GOTO 50
 412  CONTINUE
C
 411  CONTINUE
C
 410  CONTINUE
      GO TO 80
C
 50   CONTINUE
      IF (NFL .EQ. 2) GO TO 65
C
C --- P TARGET ---
      IF (NP .EQ. NM) GO TO 61
      IF (NP .EQ. 1+NM) GO TO 63
      IPA(1)=8
      IPA(2)=14
      GO TO 100
C
 61   CONTINUE
      CALL GRNDM(RNDM,1)
      IF (RNDM(1) .LT. 0.75) GO TO 62
      IPA(1)=8
      IPA(2)=16
      GO TO 100
C
 62   CONTINUE
      IPA(1)=9
      IPA(2)=14
      GO TO 100
C
 63   CONTINUE
      IPA(1)=9
      IPA(2)=16
      GO TO 100
C
C --- N TARGET ---
 65   CONTINUE
      IF (NP .EQ. -1+NM) GO TO 66
      IF (NP .EQ. NM) GO TO 68
      IPA(1)=8
      IPA(2)=16
      GO TO 100
C
 66   CONTINUE
      CALL GRNDM(RNDM,1)
      IF (RNDM(1) .LT. 0.50) GO TO 67
      IPA(1)=8
      IPA(2)=16
      GO TO 100
C
 67   CONTINUE
      IPA(1)=9
      IPA(2)=14
      GO TO 100
C
 68   CONTINUE
      IPA(1)=9
      IPA(2)=16
      GO TO 100
C
 70   CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,1003) EAB,N,NFL,NP,NM,NZ
      CALL STPAIR
      IF (INT .EQ. 1) CALL TWOB(9,NFL,N)
      IF (INT .EQ. 2) CALL GENXPT(9,NFL,N)
      GO TO 9999
C
C --- ENERGETICALLY NOT POSSIBLE TO PRODUCE CASCADE-PARTICLES ---
C --- CONTINUE WITH QUASI-ELASTIC SCATTERING ---
 55   CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,1001)
      GO TO 53
C
C --- EXCLUSIVE REACTION NOT FOUND ---
 80   CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,1004) RS,N
C
 53   CONTINUE
      INT=1
      NP=0
      NM=0
      NZ=0
      N=0.0
      IPA(1)=9
      IPA(2)=14
      IF (NFL .EQ. 2) IPA(2)=16
C
 100  CONTINUE
      DO 101 I=3,60
      IPA(I)=0
 101  CONTINUE
      IF (INT .LE. 0) GO TO 131
C
 120  CONTINUE
      NT=2
      IF (NP .EQ. 0) GO TO 122
      DO 121 I=1,NP
      NT=NT+1
      IPA(NT)=7
 121  CONTINUE
C
 122  CONTINUE
      IF (NM .EQ. 0) GO TO 124
      DO 123 I=1,NM
      NT=NT+1
      IPA(NT)=9
 123  CONTINUE
C
 124  CONTINUE
      IF (NZ .EQ. 0) GO TO 130
      DO 125 I=1,NZ
      NT=NT+1
      IPA(NT)=8
 125  CONTINUE
C
 130  CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,2004) NT,(IPA(I),I=1,20)
      IF (IPA(1) .EQ. 7) NP=NP+1
      IF (IPA(1) .EQ. 8) NZ=NZ+1
      IF (IPA(1) .EQ. 9) NM=NM+1
      GO TO 70
C
 131  CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,2005)
C
1001  FORMAT('0*CASPIM* CASCADE ENERGETICALLY NOT POSSIBLE',
     $ ' CONTINUE WITH QUASI-ELASTIC SCATTERING')
1003  FORMAT(' *CASPIM* PION- -INDUCED CASCADE, AVAIL. ENERGY',2X,F8.4,
     $ 2X,'<NTOT>',2X,F8.4,2X,'FROM',4(2X,I3),2X,'PARTICLES')
1004  FORMAT(' *CASPIM* PION- -INDUCED CASCADE, EXCLUSIVE REACTION',
     $ ' NOT FOUND TRY ELASTIC SCATTERING  AVAIL. ENERGY',2X,F8.4,2X,
     * '<NTOT>',2X,F8.4)
2001  FORMAT('0*CASPIM* TABLES FOR MULTIPLICITY DATA PION- INDUCED',
     $ 'REACTION FOR DEFINITION OF NUMBERS SEE FORTRAN CODING')
2002  FORMAT(' *CASPIM* TARGET PARTICLE FLAG',2X,I5)
2003  FORMAT(1H ,10E12.4)
2004  FORMAT(' *CASPIM* ',I3,2X,'PARTICLES, MASS INDEX ARRAY',2X,20I4)
2005  FORMAT(' *CASPIM* NO PARTICLES PRODUCED')
C
 9999 CONTINUE
      END