*
* $Id$
*
* $Log$
* Revision 1.1  1995/10/24 10:21:04  cernlib
* Initial revision
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.40  by  S.Giani
*-- Author :
      SUBROUTINE CASX0(K,INT,NFL)
C
C *** CASCADE OF XI0 ***
C *** NVE 20-JAN-1989 CERN GENEVA ***
C
C XI0  UNDERGOES INTERACTION WITH NUCLEON WITHIN NUCLEUS.
C CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE PIONS/KAONS.
C IF NOT, ASSUME NUCLEAR EXCITATION OCCURS, DEGRADE INPUT PARTICLE
C IN ENERGY AND NO OTHER PARTICLES ARE 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 PMUL(2,1200),ANORM(2,60),CECH(10),IIPA(12,2),B(2)
      DIMENSION RNDM(1)
      SAVE PMUL,ANORM
      DATA CECH/0.50,0.45,0.40,0.35,0.30,0.25,0.06,0.04,0.005,0./
C --- ARRAY IIPA DENOTES THE STRANGENESS AND CHARGE EXCHAGE REACTIONS ---
C XI0 P --> S+ S0,  XI0 P --> S0 S+
C XI0 P --> S+ L0,  XI0 P --> L0 S+
C XI0 P --> P XI0
C XI0 N --> S0 S0
C XI0 N --> L0 L0
C XI0 N --> XI- P,  XI0 N --> P XI-
C XI0 N --> S+ S-,  XI0 N --> S- S+
C XI0 N --> N XI0
      DATA IIPA/20,21,20,18,14, 21,18,27,14,20,22,16,
     *          21,20,18,20,26, 21,18,14,27,22,20,26/
      DATA B/0.7,0.7/,C/1.25/
C
C --- INITIALIZATION INDICATED BY KGINIT(20) ---
      IF (KGINIT(20) .NE. 0) GO TO 10
      KGINIT(20)=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 *** COMPUTE NORMALIZATION CONSTANTS ***
C
C --- FOR P TARGET ---
      L=0
      DO 1 NP1=1,20
      NP=NP1-1
      NMM1=NP1-2
      IF (NMM1 .LE. 0) 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) GO TO 1
      NT=NP+NM+NZ
      IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 1
      PMUL(1,L)=PMLTPC(NP,NM,NZ,NT,B(1),C)
      ANORM(1,NT)=ANORM(1,NT)+PMUL(1,L)
 1    CONTINUE
C --- FOR N TARGET ---
      L=0
      DO 2 NP1=1,20
      NP=NP1-1
      NMM1=NP1-1
      IF (NMM1 .LE. 0) NMM1=1
      NPP1=NP1+2
      DO 2 NM1=NMM1,NPP1
      NM=NM1-1
      DO 2 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF (L .GT. 1200) GO TO 2
      NT=NP+NM+NZ
      IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 2
      PMUL(2,L)=PMLTPC(NP,NM,NZ,NT,B(2),C)
      ANORM(2,NT)=ANORM(2,NT)+PMUL(2,L)
 2    CONTINUE
C
      DO 3 I=1,60
      IF (ANORM(1,I) .GT. 0.) ANORM(1,I)=1./ANORM(1,I)
      IF (ANORM(2,I) .GT. 0.) ANORM(2,I)=1./ANORM(2,I)
 3    CONTINUE
C
      IF (.NOT. NPRT(10)) GO TO 10
C
      WRITE(NEWBCD,2001)
 2001 FORMAT('0*CASX0* TABLES FOR MULT. DATA XI0 INDUCED REACTION',
     $ ' FOR DEFINITION OF NUMBERS SEE FORTRAN CODING')
      DO 4 NFL=1,2
      WRITE(NEWBCD,2002) NFL
 2002 FORMAT(' *CASX0* TARGET PARTICLE FLAG',2X,I5)
      WRITE(NEWBCD,2003) (ANORM(NFL,I),I=1,60)
      WRITE(NEWBCD,2003) (PMUL(NFL,I),I=1,1200)
 2003 FORMAT(1H ,10E12.4)
 4    CONTINUE
C
C --- SELECT TARGET NUCLEON ---
 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(26)
C
C --- RESET STRANGENESS FIXING FLAG ---
      NVEFIX=0
C
C *** ELASTIC SCATTERING ***
      NP=0
      NM=0
      NZ=0
      N=0.
      IPA(1)=26
      IPA(2)=14
      IF (NFL .EQ. 2) IPA(2)=16
C
      IF (INT .EQ. 2) GO TO 20
C
C *** INTRODUCE CHARGE AND STRANGENESS EXCHANGE REACTIONS ***
      IPLAB=IFIX(P*2.5)+1
      IF (IPLAB .GT. 10) IPLAB=10
      CALL GRNDM(RNDM,1)
      IF (RNDM(1) .GT. (CECH(IPLAB)/ATNO2**0.42)) GO TO 120
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      IRN=IFIX(RAN*5.)+1
      IF (NFL .EQ. 2) IRN=5+IFIX(RAN*7.)+1
      IF (NFL .EQ. 1) IRN=MAX(IRN,5)
      IF (NFL .EQ. 2) IRN=MAX(IRN,12)
      IPA(1)=IIPA(IRN,1)
      IPA(2)=IIPA(IRN,2)
      GO TO 120
C
C --- CHECK IF ENERGETICALLY POSSIBLE TO PRODUCE ONE EXTRA PION ---
 20   CONTINUE
      IF (EAB .LE. RMASS(7)) GO TO 55
C
C --- NO. OF TOTAL PARTICLES VS SQRT(S)-MP-MSM ---
      ALEAB=LOG(EAB)
      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
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 .GE. 1.0E-10)) ADDNVE=DUM1*DUM3
      ANPN=ANPN+ADDNVE
 21   CONTINUE
      ANPN=1./ANPN
C
C --- CHECK FOR TARGET NUCLEON TYPE ---
      CALL GRNDM(RNDM,1)
      RAN=RNDM(1)
      EXCS=0.
      GO TO (30,40),NFL
C
C --- PROTON TARGET ---
 30   CONTINUE
      L=0
      DO 31 NP1=1,20
      NP=NP1-1
      NMM1=NP1-2
      IF (NMM1 .LE. 0) 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) GO TO 31
      NT=NP+NM+NZ
      IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 31
      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) GO TO 100
   31 CONTINUE
      GO TO 80
C
C --- NEUTRON TARGET ---
 40   CONTINUE
      L=0
      DO 41 NP1=1,20
      NP=NP1-1
      NMM1=NP1-1
      IF (NMM1 .LE. 0) NMM1=1
      NPP1=NP1+2
      DO 41 NM1=NMM1,NPP1
      NM=NM1-1
      DO 41 NZ1=1,20
      NZ=NZ1-1
      L=L+1
      IF (L .GT. 1200) GO TO 41
      NT=NP+NM+NZ
      IF ((NT .LE. 0) .OR. (NT .GT. 60)) GO TO 41
      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) GO TO 100
   41 CONTINUE
      GO TO 80
C
 50   CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,1003) EAB,N,NFL,NP,NM,NZ
 1003 FORMAT(' *CASX0* XI0 -INDUCED CASCADE,',
     $ ' AVAIL. ENERGY',2X,F8.4,
     $ 2X,'<NTOT>',2X,F8.4,2X,'FROM',4(2X,I3),2X,'PARTICLES')
      IF (INT .EQ. 1) CALL TWOB(27,NFL,N)
      IF (INT .EQ. 2) CALL GENXPT(27,NFL,N)
      GO TO 9999
C
C *** ENERGETICALLY NOT POSSIBLE TO PRODUCE ONE EXTRA PION ***
 55   CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,1001)
 1001 FORMAT('0*CASX0* CASCADE ENERGETICALLY NOT POSSIBLE',
     $ ' CONTINUE WITH QUASI-ELASTIC SCATTERING')
      GO TO 53
C
C *** EXCLUSIVE REACTION NOT FOUND ***
 80   CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,1004) RS,N
 1004 FORMAT(' *CASX0* XI0 -INDUCED CASCADE,',
     $ ' EXCLUSIVE REACTION NOT FOUND',
     $ ' TRY ELASTIC SCATTERING  AVAIL. ENERGY',2X,F8.4,2X,
     $ '<NTOT>',2X,F8.4)
C
 53   CONTINUE
      INT=1
      NP=0
      NM=0
      NZ=0
      IPA(1)=26
      IPA(2)=14
      IF (NFL .EQ. 2) IPA(2)=16
      GO TO 120
C
C *** INELASTIC INTERACTION HAS OCCURRED ***
C *** NUMBER OF SECONDARY MESONS DETERMINED BY KNO DISTRIBUTION ***
 100  CONTINUE
      DO 101 I=1,60
      IPA(I)=0
 101  CONTINUE
C
      IF (INT .LE. 0) GO TO 131
C
C --- TAKE TARGET NUCLEON TYPE INTO ACCOUNT ---
      GO TO (102,112),NFL
C
C --- PROTON TARGET ---
 102  CONTINUE
C --- CHECK FOR TOTAL CHARGE OF FINAL STATE MESONS TO DETERMINE ---
C --- THE KIND OF BARYONS TO BE PRODUCED TAKING INTO ACCOUNT    ---
C --- CHARGE AND STRANGENESS CONSERVATION                       ---
      NCHT=NP-NM
      IF (NCHT .LT. 1) GO TO 103
      IF (NCHT .EQ. 1) GO TO 104
      IF (NCHT .GT. 1) GO TO 105
C
 103  CONTINUE
C --- XI0 P ---
      IPA(1)=26
      IPA(2)=14
      IF (NCHT .EQ. 0) GO TO 120
C --- CHARGE MISMATCH ==> TAKE A S+ AND CORRECT THE STRANGENESS ---
C --- BY REPLACING A PI- BY K- ---
C --- S+ P ---
      IPA(1)=20
      IPA(2)=14
      NVEFIX=1
      GO TO 120
C
 104  CONTINUE
C --- XI0 N ---
      IPA(1)=26
      IPA(2)=16
      CALL GRNDM(RNDM,1)
      IF (RNDM(1) .LT. 0.5) GO TO 120
C --- XI- P ---
      IPA(1)=27
      IPA(2)=14
      GO TO 120
C
 105  CONTINUE
C --- XI- N ---
      IPA(1)=27
      IPA(2)=16
      GO TO 120
C
C --- NEUTRON TARGET ---
 112  CONTINUE
C --- CHECK FOR TOTAL CHARGE OF FINAL STATE MESONS TO DETERMINE ---
C --- THE KIND OF BARYONS TO BE PRODUCED TAKING INTO ACCOUNT    ---
C --- CHARGE AND STRANGENESS CONSERVATION                       ---
      NCHT=NP-NM
      IF (NCHT .LT. 0) GO TO 113
      IF (NCHT .EQ. 0) GO TO 114
      IF (NCHT .GT. 0) GO TO 115
C
 113  CONTINUE
C --- XI0 P ---
      IPA(1)=26
      IPA(2)=14
      IF (NCHT .EQ. -1) GO TO 120
C --- CHARGE MISMATCH ==> TAKE A S+ AND CORRECT THE STRANGENESS ---
C --- BY REPLACING A PI- BY K- ---
C --- S+ P ---
      IPA(1)=20
      IPA(2)=14
      NVEFIX=1
      GO TO 120
C
 114  CONTINUE
C --- XI0 N ---
      IPA(1)=26
      IPA(2)=16
      CALL GRNDM(RNDM,1)
      IF (RNDM(1) .LT. 0.5) GO TO 120
C --- XI- P ---
      IPA(1)=27
      IPA(2)=14
      GO TO 120
C
 115  CONTINUE
C --- XI- N ---
      IPA(1)=27
      IPA(2)=16
C
C --- TAKE PIONS FOR ALL SECONDARY MESONS ---
C --- REPLACE PI BY K IN CASE OF STRANGENESS TO BE FIXED ---
 120  CONTINUE
      NT=2
C
      IF (NP .EQ. 0) GO TO 122
C
C --- PI+ ---
      DO 121 I=1,NP
      NT=NT+1
      IPA(NT)=7
 121  CONTINUE
C
 122  CONTINUE
      IF (NM .EQ. 0) GO TO 124
C
C --- PI- ---
      DO 123 I=1,NM
      NT=NT+1
      IPA(NT)=9
      IF (NVEFIX .GE. 1) IPA(NT)=13
      IF (NPRT(4) .AND. (NVEFIX .GE. 1)) PRINT 3000
 3000 FORMAT(' *CASX0* K- INTRODUCED')
      NVEFIX=NVEFIX-1
 123  CONTINUE
C
 124  CONTINUE
      IF (NZ .EQ. 0) GO TO 130
C
C --- PI0 ---
      DO 125 I=1,NZ
      NT=NT+1
      IPA(NT)=8
 125  CONTINUE
C
C --- ALL SECONDARY PARTICLES HAVE BEEN DEFINED ---
C --- NOW GO FOR MOMENTA AND X VALUES ---
 130  CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,2004) NT,(IPA(I),I=1,60)
 2004 FORMAT(' *CASX0* ',I3,' PARTICLES PRODUCED. MASS INDEX ARRAY : '/
     $ 3(1H ,20(I3,1X)/))
      GO TO 50
C
 131  CONTINUE
      IF (NPRT(4)) WRITE(NEWBCD,2005)
 2005 FORMAT(' *CASX0* NO PARTICLES PRODUCED')
C
 9999 CONTINUE
      END
