      SUBROUTINE CGHEI
 
C-----------------------------------------------------------------------
C  C(ORSIKA) GHE(ISHA) I(NTERFACE)
C
C  MAIN STEERING ROUTINE FOR HADRON PACKAGE GHEISHA ***
C  THIS SUBROUTINE IS CALLED FROM NUCINT
C
C  ORIGIN  : F.CARMINATI, H.FESEFELDT (ROUTINE GHESIG)
C  REDESIGN: P. GABRIEL IK1  FZK KARLSRUHE
C-----------------------------------------------------------------------
 
*KEEP,CGCOMP.
      PARAMETER (KK=3)
      COMMON/CGCOMP/ ACOMP,ZCOMP,WCOMP
      REAL           ACOMP(KK),ZCOMP(KK),WCOMP(KK)
*KEEP,ELABCT.
      COMMON /ELABCT/  ELCUT
      DOUBLE PRECISION ELCUT(4)
*KEEP,ELADPM.
      COMMON /ELADPM/  ELMEAN,ELMEAA,IELDPM,IELDPA
      DOUBLE PRECISION ELMEAN(37),ELMEAA(37)
      INTEGER          IELDPM(37,13),IELDPA(37,13)
*KEEP,ELASTY.
      COMMON /ELASTY/  ELAST,IELIS,IELHM,IELNU,IELPI
      DOUBLE PRECISION ELAST
      INTEGER          IELIS(20),IELHM(20),IELNU(20),IELPI(20)
*KEEP,GENER.
      COMMON /GENER/   GEN,ALEVEL
      DOUBLE PRECISION GEN,ALEVEL
*KEEP,MULT.
      COMMON /MULT/    EKINL,MSMM,MULTMA,MULTOT
      DOUBLE PRECISION EKINL
      INTEGER          MSMM,MULTMA(37,13),MULTOT(37,13)
*KEEP,PAM.
      COMMON /PAM/     PAMA,SIGNUM
      DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
*KEEP,PARPAR.
      COMMON /PARPAR/  CURPAR,SECPAR,PRMPAR,OUTPAR,C,
     *                 E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
      DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
     *                 C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
      INTEGER          ITYPE,LEVL
*KEEP,RUNPAR.
      COMMON /RUNPAR/  FIXHEI,THICK0,HILOECM,HILOELB,
     *                 STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
     *                 MONIOU,MDEBUG,NUCNUC,
     *                 CETAPE,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE,
     *                 DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
      COMMON /RUNPAC/  DSN,HOST,USER
      DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
      REAL             STEPFC
      INTEGER          NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE
      INTEGER          CETAPE
      CHARACTER*79     DSN
      CHARACTER*20     HOST,USER
 
      LOGICAL          DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
*KEND.
 
      DOUBLE PRECISION ELASTI,ELABOR,PLX,PLY,PLZ,PLSQ,PLTOT,RMASSK
 
      COMMON/GSECTI/ AIEL(20),AIIN(20),AIFI(20),AICA(20),ALAM,K0FLAG
      INTEGER K0FLAG
      REAL AIEL,AIIN,AIFI,AICA,ALAM
 
C --- GHEISHA COMMONS ---
      PARAMETER (MXGKGH=100)
      PARAMETER (MXGKPV=MXGKGH)
      COMMON /VECUTY/ PV(10,MXGKPV)
 
      COMMON/CONSTS/ PI,TWPI,PIBTW,MP,MPI,MMU,MEL,MKCH,MK0,SMP,SMPI,
     $               SMU,CT,CTKCH,CTK0,
     $               ML0,MSP,MS0,MSM,MX0,MXM,CTL0,CTSP,CTSM,CTX0,CTXM,
     $               RMASS(35),RCHARG(35)
 
                     REAL MP,MPI,MMU,MEL,MKCH,MK0,
     *                    ML0,MSP,MS0,MSM,MX0,MXM
 
      PARAMETER (MXEVEN=12*MXGKGH)
      COMMON/EVENT / NSIZE,NCUR,NEXT,NTOT,EVE(MXEVEN)
 
      COMMON/PRNTFL/INBCD,NEWBCD,INBIN,NEWBIN,NPEVT,NEVTP,LPRT,NPRT(10)
                    LOGICAL LPRT,NPRT
 
 
C --- "NEVENT" CHANGED TO "KEVENT" IN COMMON /CURPAR/ DUE TO CLASH ---
C --- WITH VARIABLE "NEVENT" IN GEANT COMMON ---
 
      PARAMETER (MXGKCU=MXGKGH)
      COMMON /CURPAR/ WEIGHT(10),DDELTN,IFILE,IRUN,NEVT,KEVENT,SHFLAG,
     $                ITHST,ITTOT,ITLST,IFRND,TOFCUT,CMOM(5),CENG(5),
     $                RS,S,ENP(10),NP,NM,NN,NR,NO,NZ,IPA(MXGKCU),
     $                ATNO2,ZNO2
 
C --- "IPART" CHANGED TO "KPART" IN COMMON /RESULT/ DUE TO CLASH ---
C --- WITH VARIABLE "IPART" IN GEANT COMMON ---
 
      COMMON /RESULT/ XEND,YEND,ZEND,RCA,RCE,AMAS,NCH,TOF,PX,PY,PZ,
     $                USERW,INTCT,P,EN,EK,AMASQ,DELTN,ITK,NTK,KPART,IND,
     $                LCALO,ICEL,SINL,COSL,SINP,COSP,
     $                XOLD,YOLD,ZOLD,POLD,PXOLD,PYOLD,PZOLD,
     $                XSCAT,YSCAT,ZSCAT,PSCAT,PXSCAT,PYSCAT,PZSCAT
                      REAL NCH,INTCT
 
C --- "ABSL(21)" CHANGED TO "ABSLTH(21)" IN COMMON /MAT/ DUE TO CLASH ---
C --- WITH VARIABLE "ABSL" IN GEANT COMMON ---
 
      COMMON /MAT/ LMAT,
     $             DEN(21),RADLTH(21),ATNO(21),ZNO(21),ABSLTH(21),
     $             CDEN(21),MDEN(21),X0DEN(21),X1DEN(21),RION(21),
     $             MATID(21),MATID1(21,24),PARMAT(21,10),
     $             IFRAT,IFRAC(21),FRAC1(21,10),DEN1(21,10),
     $             ATNO1(21,10),ZNO1(21,10)
 
      DIMENSION IPELOS(35)
      REAL EMAX,EEESQ
      SAVE IDEOL
 
      DIMENSION RNDM(1)
 
C --- DIMENSION STMTS. FOR GEANT/GHEISHA PARTICLE CODE CONVERSIONS ---
C --- KIPART(I)=GHEISHA CODE CORRESPONDING TO GEANT   CODE I ---
C --- IKPART(I)=GEANT   CODE CORRESPONDING TO GHEISHA CODE I ---
 
      DIMENSION KIPART(48),IKPART(35)
 
C --- DATA STMTS. FOR GEANT/GHEISHA PARTICLE CODE CONVERSIONS ---
C --- KIPART(I)=GHEISHA CODE CORRESPONDING TO GEANT   CODE I ---
C --- IKPART(I)=GEANT   CODE CORRESPONDING TO GHEISHA CODE I ---
 
      DATA KIPART/
     $               1,   3,   4,   2,   5,   6,   8,   7,
     $               9,  12,  10,  13,  16,  14,  15,  11,
     $              35,  18,  20,  21,  22,  26,  27,  33,
     $              17,  19,  23,  24,  25,  28,  29,  34,
     $              35,  35,  35,  35,  35,  35,  35,  35,
     $              35,  35,  35,  35,  30,  31,  32,  35/
 
      DATA IKPART/
     $               1,   4,   2,   3,   5,   6,   8,   7,
     $               9,  11,  16,  10,  12,  14,  15,  13,
     $              25,  18,  26,  19,  20,  21,  27,  28,
     $              29,  22,  23,  30,  31,  45,  46,  47,
     $              24,  32,  48/
 
 
C --- DENOTE STABLE PARTICLES ACCORDING TO GHEISHA CODE ---
C --- STABLE : GAMMA, NEUTRINO, ELECTRON, PROTON AND HEAVY FRAGMENTS ---
C --- WHEN STOPPING THESE PARTICLES ONLY LOOSE THEIR KINETIC ENERGY ---
      DATA IPELOS/
     $             1,   1,   0,   1,   0,   0,   0,   0,
     $             0,   0,   0,   0,   0,   1,   0,   0,
     $             0,   0,   0,   0,   0,   0,   0,   0,
     $             0,   0,   0,   0,   0,   1,   1,   1,
     $             0,   0,   1/
 
C --- LOWERBOUND OF KINETIC ENERGY BIN IN N CROSS-SECTION TABLES ---
      DATA TEKLOW /0.0001/
 
C --- KINETIC ENERGY TO SWITCH FROM "CASN" TO "GNSLWD" FOR N CASCADE ---
      DATA SWTEKN /0.05/
 
      DATA IDEOL/0/
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,*) 'CGHEI :'
 
C --- DEFINE PARTICLE TYPE
      IF     ( ITYPE .LE.  48 ) THEN
         IPART = ITYPE
      ELSEIF ( ITYPE .EQ. 201 ) THEN
         IPART = 45
      ELSEIF ( ITYPE .EQ. 301 ) THEN
         IPART = 46
      ELSEIF ( ITYPE .EQ. 402 ) THEN
         IPART = 47
      ELSE
         WRITE (MONIOU,7795) ITYPE
 7795    FORMAT (//,' *CGHEI*  ILLEGAL PARTICLE TYPE OCCURS =',I5)
         IPART = 48
      ENDIF
 
      NETEST=IKPART(KPART)
      IF ( NETEST .EQ. IPART ) GO TO 9004
 
      WRITE(MONIOU,8881) IPART,KPART
 8881 FORMAT(' *CGHEI* IPART,KPART = ',2(I3,1X)/
     $ ' *CGHEI* ======> PARTICLE TYPES DO NOT MATCH <=======')
      STOP
 
 9004 CONTINUE
      KPART=KIPART(IPART)
      KKPART=KPART
 
C --- TRANSPORT THE TRACK NUMBER TO GHEISHA AND INITIALISE SOME NUMBERS
C --- NTK=ITRA   ITRA = CURRENT TRACK NUMBER IN GEANT (GCKINE)
      NTK=0
      INTCT=0.0
      NEXT=1
      NTOT=0
      INT=0
      TOF=0.0
 
C --- STORE COORDINATES FOR SECONDARIES AND RESET ITYPE
      SECPAR(1) = 0.
      DO 7001 LK = 5, 8
         SECPAR(LK) = CURPAR(LK)
 7001 CONTINUE
 
 
C --- FILL RESULT COMMON FOR THIS TRACK WITH CORSIKA VALUES ---
 
      AMAS=RMASS(KPART)
      NCH=RCHARG(KPART)
  107 XEND = CURPAR(7)
      YEND = CURPAR(8)
      ZEND = CURPAR(5)
      SINL = -CURPAR(3)
      PHI  = CURPAR(4)
      USERW=0.0
 
      AMASQ=AMAS*AMAS
      EN = CURPAR(2) * ABS(AMAS)
      EK = ABS ( EN - ABS(AMAS) )
      ENOLD=EN
      EMAX = 0.
      P  = SQRT ( EN*EN - AMASQ )
      ELABOR = EN
 
      SINP = SIN(PHI)
      COSP = COS(PHI)
      COSL = SQRT ( ABS(1.-SINL**2) )
 
      PX = COSL * COSP
      PY = COSL * SINP
      PZ = SINL
 
C --- SET GHEISHA INDEX FOR THE CURRENT MEDIUM ALWAYS TO 1 ---
      IND=1
 
C --- TRANSFER GLOBAL MATERIAL CONSTANTS FOR CURRENT MEDIUM ---
C --- DETAILED DATA FOR COMPOUNDS IS OBTAINED VIA ROUTINE COMPO ---
      ATNO(IND+1) = 14.56
      ZNO(IND+1)  =  7.265
      DEN(IND+1)  =  0.0
      RADLTH(IND+1)= 0.0
      ABSLTH(IND+1)= 0.0
 
C --- SETUP PARMAT FOR PHYSICS STEERING ---
      PARMAT(IND+1,10)=0.0
 
  5   CONTINUE
 
C --- INDICATE LIGHT (<= PI) AND HEAVY PARTICLES (HISTORICALLY) ---
C --- CALIM CODE ---
      J=2
      TEST=RMASS(7)-0.001
      IF (ABS(AMAS) .LT. TEST) J=1
 
C *** DIVISION INTO VARIOUS INTERACTION CHANNELS DENOTED BY "INT" ***
C THE CONVENTION FOR "INT" IS THE FOLLOWING
 
C INT  = -1 REACTION CROSS SECTIONS NOT YET TABULATED/PROGRAMMED
C      =  0 NO INTERACTION
C      =  1 ELEASTIC SCATTERING
C      =  2 INELASTIC SCATTERING
C      =  3 NUCLEAR FISSION WITH INELEASTIC SCATTERING
C      =  4 NEUTRON CAPTURE
C INT  =  3, 4  SHOULD BE DELETED FOR AIR TARGET
 
C --- INTACT CODE ---
      ALAM1=0.0
      CALL GRNDM(RNDM,1)
      RAT=RNDM(1)*ALAM
      ATNO2 = 14.56
      ZNO2  = 7.265
 
      DO 6 K=1,KK
      ATNO2 = ACOMP(K)
      ZNO2 =  ZCOMP(K)
 
C --- TRY FOR ELASTIC SCATTERING ---
      INT=1
      ALAM1=ALAM1+AIEL(K)
      IF (RAT .LT. ALAM1) GO TO 8
 
C --- TRY FOR INELASTIC SCATTERING ---
      INT=2
      ALAM1=ALAM1+AIIN(K)
      IF (RAT .LT. ALAM1) GO TO 8
 
C --- TRY FOR NEUTRON CAPTURE ---
      INT=4
      ALAM1=ALAM1+AICA(K)
      IF (RAT .LT. ALAM1) GO TO 8
 
 6    CONTINUE
 
C --- NO REACTION SELECTED ==> ELASTIC SCATTERING ---
      INT=1
 
C *** TAKE ACTION ACCORDING TO SELECTED REACTION CHANNEL ***
C --- FOLLOWING CODE IS A TRANSLATION OF "CALIM" INTO GEANT JARGON ---
 
 8    CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,1001) INT
 1001 FORMAT(' *CGHEI* INTERACTION TYPE CHOSEN INT = ',I3)
 
      IF (INT .NE. 4) GO TO 10
 
C --- NEUTRON CAPTURE ---
      IF (NPRT(9)) WRITE(MDEBUG,2000)
 2000 FORMAT(' *CGHEI* ROUTINE CAPTUR WILL BE CALLED')
      CALL CAPTUR(NOPT)
      GO TO 40
 
 10   CONTINUE
 
C --- ELASTIC AND INELASTIC SCATTERING ---
      PV(1,MXGKPV)=P*PX
      PV(2,MXGKPV)=P*PY
      PV(3,MXGKPV)=P*PZ
      PV(4,MXGKPV)=EN
      PV(5,MXGKPV)=AMAS
      PV(6,MXGKPV)=NCH
      PV(7,MXGKPV)=TOF
      PV(8,MXGKPV)=KPART
      PV(9,MXGKPV)=0.
      PV(10,MXGKPV)=USERW
 
C --- ADDITIONAL PARAMETERS TO SIMULATE FERMI MOTION AND EVAPORATION ---
      DO 111 JENP=1,10
         ENP(JENP)=0.
 111  CONTINUE
      ENP(5)=EK
      ENP(6)=EN
      ENP(7)=P
 
      IF (INT .NE. 1) GO TO 12
 
C *** ELASTIC SCATTERING PROCESSES ***
 
C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS ---
      IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35
 
C --- NORMAL ELASTIC SCATTERING FOR LIGHT MEDIA ---
      IF (ATNO2 .LT. 1.5) GO TO 35
 
C --- COHERENT ELASTIC SCATTERING FOR HEAVY MEDIA ---
      IF (NPRT(9)) WRITE(MDEBUG,2002)
 2002 FORMAT(' *CGHEI* ROUTINE COSCAT WILL BE CALLED')
      CALL COSCAT
      GO TO 40
 
C *** NON-ELASTIC SCATTERING PROCESSES ***
 12   CONTINUE
 
C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS ---
      IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35
 
C *** USE SOMETIMES NUCLEAR REACTION ROUTINE "NUCREC" FOR LOW ENERGY ***
C *** PROTON AND NEUTRON SCATTERING ***
      CALL GRNDM(RNDM,1)
      TEST1=RNDM(1)
      TEST2=4.5*(EK-0.01)
      IF ((KPART .EQ. 14) .AND. (TEST1 .GT. TEST2)) GO TO 85
      IF ((KPART .EQ. 16) .AND. (TEST1 .GT. TEST2)) GO TO 86
 
C *** FERMI MOTION AND EVAPORATION ***
      TKIN=CINEMA(EK)
      PV(9,MXGKPV) = TKIN
      ENP(5)=EK+TKIN
C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
      IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
      ENP(6)=ENP(5)+ABS(AMAS)
      ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
      ENP(7)=SQRT(ABS(ENP(7)))
      TKIN=FERMI(ENP(5))
      ENP(5)=ENP(5)+TKIN
C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
      IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
      ENP(6)=ENP(5)+ABS(AMAS)
      ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
      ENP(7)=SQRT(ABS(ENP(7)))
      TKIN=EXNU(ENP(5))
      ENP(5)=ENP(5)-TKIN
C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES ---
      IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW
      ENP(6)=ENP(5)+ABS(AMAS)
      ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS)
      ENP(7)=SQRT(ABS(ENP(7)))
 
C *** IN CASE OF ENERGY ABOVE CUT-OFF LET THE PARTICLE CASCADE ***
      IF ( ENP(5) .GT. ELCUT(1)) GO TO 35
 
C --- SECOND CHANCE FOR ANTI-BARYONS DUE TO POSSIBLE ANNIHILATION ---
      IF ((AMAS .GE. 0.0) .OR. (KPART .LE. 14)) GO TO 13
      ANNI=1.3*P
      IF (ANNI .GT. 0.4) ANNI=0.4
      CALL GRNDM(RNDM,1)
      TEST=RNDM(1)
      IF (TEST .GT. ANNI) GO TO 35
 
C *** PARTICLE WITH ENERGY BELOW CUT-OFF ***
C --- ==> ONLY NUCLEAR EVAPORATION AND QUASI-ELASTIC SCATTERING ---
 13   CONTINUE
 
      IF (NPRT(9)) WRITE(MDEBUG,1002) KPART,EK,EN,P,ENP(5),ENP(6),ENP(7)
 1002 FORMAT(' *CGHEI* ENERGY BELOW CUT-OFF FOR GHEISHA PARTICLE ',I3/
     $ ' EK,EN,P,ENP(5),ENP(6),ENP(7) = ',6(G12.5,1X))
 
      IF ((KPART .NE. 14) .AND. (KPART .NE. 16)) GO TO 14
      IF (KPART .EQ. 16) GO TO 86
 
C --- SLOW PROTON ---
 85   CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2003) EK,KPART
 2003 FORMAT(' *CGHEI* ROUTINE NUCREC WILL BE CALLED',
     $ ' EK = ',G12.5,' GEV  KPART = ',I3)
      CALL NUCREC(NOPT,2)
 
      IF (NOPT .NE. 0) GO TO 50
 
      IF (NPRT(9)) WRITE(MDEBUG,2004)EK,KPART
 2004 FORMAT(' *CGHEI* ROUTINE COSCAT WILL BE CALLED',
     $ ' EK = ',G12.5,' GEV  KPART = ',I3)
      CALL COSCAT
      GO TO 40
 
C --- SLOW NEUTRON ---
 86   CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2015)
      NUCFLG=0
      CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW)
      IF (NUCFLG .NE. 0) GO TO 50
      GO TO 40
 
C --- OTHER SLOW PARTICLES ---
 14   CONTINUE
      IPA(1)=KPART
C --- DECIDE FOR PROTON OR NEUTRON TARGET ---
      IPA(2)=16
      CALL GRNDM(RNDM,1)
      TEST1=RNDM(1)
      TEST2=ZNO2/ATNO2
      IF (TEST1 .LT. TEST2) IPA(2)=14
      AVERN=0.0
      NFL=1
      IF (IPA(2) .EQ. 16) NFL=2
      IPPP=KPART
      IF (NPRT(9)) WRITE(MDEBUG,2005)
 2005 FORMAT(' *CGHEI* ROUTINE TWOB WILL BE CALLED')
      CALL TWOB(IPPP,NFL,AVERN)
      GOTO 40
 
C --- INITIALISATION OF CASCADE QUANTITIES ---
 35   CONTINUE
 
C *** CASCADE GENERATION ***
C --- CALCULATE FINAL STATE MULTIPLICITY AND LONGITUDINAL AND ---
C --- TRANSVERSE MOMENTUM DISTRIBUTIONS ---
 
C --- FIXED PARTICLE TYPE TO STEER THE CASCADE ---
      KKPART=KPART
 
C --- NO CASCADE FOR LEPTONS ---
      IF (KKPART .LE. 6) GO TO 9999
 
C *** WHAT TO DO WITH "NEW PARTICLES" FOR GHEISHA ?????? ***
C --- RETURN FOR THE TIME BEING ---
      IF (KKPART .GE. 35) GO TO 9999
 
C --- CASCADE OF HEAVY FRAGMENTS
      IF ((KKPART .GE. 30) .AND. (KKPART .LE. 32)) GO TO 390
 
C --- INITIALIZE THE IPA ARRAY ---
      CALL VZERO(IPA(1),MXGKCU)
 
C --- CASCADE OF OMEGA - AND OMEGA - BAR ---
      IF (KKPART .EQ. 33) GO TO 330
      IF (KKPART .EQ. 34) GO TO 331
 
      NVEPAR=KKPART-17
      IF (NVEPAR .LE. 0) GO TO 15
      GO TO (318,319,320,321,322,323,324,325,326,327,328,329),NVEPAR
 
 15   CONTINUE
      NVEPAR=KKPART-6
      GO TO (307,308,309,310,311,312,313,314,315,316,317,318),NVEPAR
 
C --- PI+ CASCADE ---
 307  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2006)
 2006 FORMAT(' *CGHEI* ROUTINE CASPIP WILL BE CALLED')
      CALL CASPIP(J,INT,NFL)
      GO TO 40
 
C --- PI0 ==> NO CASCADE ---
 308  CONTINUE
      GO TO 40
 
C --- PI- CASCADE ---
 309  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2007)
 2007 FORMAT(' *CGHEI* ROUTINE CASPIM WILL BE CALLED')
      CALL CASPIM(J,INT,NFL)
      GO TO 40
 
C --- K+ CASCADE ---
 310  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2008)
 2008 FORMAT(' *CGHEI* ROUTINE CASKP WILL BE CALLED')
      CALL CASKP(J,INT,NFL)
      GO TO 40
 
C --- K0 CASCADE ---
 311  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2009)
 2009 FORMAT(' *CGHEI* ROUTINE CASK0 WILL BE CALLED')
      CALL CASK0(J,INT,NFL)
      GO TO 40
 
C --- K0 BAR CASCADE ---
 312  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2010)
 2010 FORMAT(' *CGHEI* ROUTINE CASK0B WILL BE CALLED')
      CALL CASK0B(J,INT,NFL)
      GO TO 40
 
C --- K- CASCADE ---
 313  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2011)
 2011 FORMAT(' *CGHEI* ROUTINE CASKM WILL BE CALLED')
      CALL CASKM(J,INT,NFL)
      GO TO 40
 
C --- PROTON CASCADE ---
 314  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2012)
 2012 FORMAT(' *CGHEI* ROUTINE CASP WILL BE CALLED')
      CALL CASP(J,INT,NFL)
      GO TO 40
 
C --- PROTON BAR CASCADE ---
 315  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2013)
 2013 FORMAT(' *CGHEI* ROUTINE CASPB WILL BE CALLED')
      CALL CASPB(J,INT,NFL)
      GO TO 40
 
C --- NEUTRON CASCADE ---
 316  CONTINUE
      NUCFLG=0
      IF (EK .GT. SWTEKN) THEN
         CALL CASN(J,INT,NFL)
         IF (NPRT(9)) WRITE(MDEBUG,2014)
 2014 FORMAT(' *CGHEI* ROUTINE CASN WILL BE CALLED')
      ELSE
         CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW)
         IF (NPRT(9)) WRITE(MDEBUG,2015)
 2015 FORMAT(' *CGHEI* ROUTINE GNSLWD WILL BE CALLED')
      ENDIF
      IF (NUCFLG .NE. 0) GO TO 50
      GO TO 40
 
C --- NEUTRON BAR CASCADE ---
 317  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2016)
 2016 FORMAT(' *CGHEI* ROUTINE CASNB WILL BE CALLED')
      CALL CASNB(J,INT,NFL)
      GO TO 40
 
C --- LAMBDA CASCADE ---
 318  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2017)
 2017 FORMAT(' *CGHEI* ROUTINE CASL0 WILL BE CALLED')
      CALL CASL0(J,INT,NFL)
      GO TO 40
 
C --- LAMBDA BAR CASCADE ---
 319  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2018)
 2018 FORMAT(' *CGHEI* ROUTINE CASAL0 WILL BE CALLED')
      CALL CASAL0(J,INT,NFL)
      GO TO 40
 
C --- SIGMA + CASCADE ---
 320  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2019)
 2019 FORMAT(' *CGHEI* ROUTINE CASSP WILL BE CALLED')
      CALL CASSP(J,INT,NFL)
      GO TO 40
 
C --- SIGMA 0 ==> NO CASCADE ---
 321  CONTINUE
      GO TO 40
 
C --- SIGMA - CASCADE ---
 322  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2020)
 2020 FORMAT(' *CGHEI* ROUTINE CASSM WILL BE CALLED')
      CALL CASSM(J,INT,NFL)
      GO TO 40
 
C --- SIGMA + BAR CASCADE ---
 323  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2021)
 2021 FORMAT(' *CGHEI* ROUTINE CASASP WILL BE CALLED')
      CALL CASASP(J,INT,NFL)
      GO TO 40
 
C --- SIGMA 0 BAR ==> NO CASCADE ---
 324  CONTINUE
      GO TO 40
 
C --- SIGMA - BAR CASCADE ---
 325  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2022)
 2022 FORMAT(' *CGHEI* ROUTINE CASASM WILL BE CALLED')
      CALL CASASM(J,INT,NFL)
      GO TO 40
 
C --- XI 0 CASCADE ---
 326  CONTINUE
      IF (NPRT(9)) PRINT 2023
 2023 FORMAT(' *CGHEI* ROUTINE CASX0 WILL BE CALLED')
      CALL CASX0(J,INT,NFL)
      GO TO 40
 
C --- XI - CASCADE ---
 327  CONTINUE
      IF (NPRT(9)) PRINT 2024
 2024 FORMAT(' *CGHEI* ROUTINE CASXM WILL BE CALLED')
      CALL CASXM(J,INT,NFL)
      GO TO 40
 
C --- XI 0 BAR CASCADE ---
 328  CONTINUE
      IF (NPRT(9)) PRINT 2025
 2025 FORMAT(' *CGHEI* ROUTINE CASAX0 WILL BE CALLED')
      CALL CASAX0(J,INT,NFL)
      GO TO 40
 
C --- XI - BAR CASCADE ---
 329  CONTINUE
      IF (NPRT(9)) PRINT 2026
 2026 FORMAT(' *CGHEI* ROUTINE CASAXM WILL BE CALLED')
      CALL CASAXM(J,INT,NFL)
      GO TO 40
 
C --- OMEGA - CASCADE ---
 330  CONTINUE
      IF (NPRT(9)) PRINT 2027
 2027 FORMAT(' *CGHEI* ROUTINE CASOM WILL BE CALLED')
      CALL CASOM(J,INT,NFL)
      GO TO 40
 
C --- OMEGA - BAR CASCADE ---
 331  CONTINUE
      IF (NPRT(9)) PRINT 2028
 2028 FORMAT(' *CGHEI* ROUTINE CASAOM WILL BE CALLED')
      CALL CASAOM(J,INT,NFL)
      GO TO 40
 
C --- HEAVY FRAGMENT CASCADE ---
 390  CONTINUE
      IF (NPRT(9)) WRITE(MDEBUG,2090)
 2090 FORMAT(' *CGHEI* ROUTINE CASFRG WILL BE CALLED')
      NUCFLG=0
      CALL CASFRG(NUCFLG,INT,NFL)
      IF (NUCFLG .NE. 0) GO TO 50
 
C *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED ***
 40   CONTINUE
      IF ((NTOT .NE. 0) .OR. (KKPART .NE. KPART)) GO TO 50
 
 50   CONTINUE
 
      NVEDUM=KIPART(IPART)
      IF (NPRT(9)) WRITE(MDEBUG,1004)NTOT,IPART,KPART,KKPART,NVEDUM
 1004 FORMAT(' *CGHEI* SEC. GEN. NTOT,IPART,KPART,KKPART,KIPART = ',
     $ 5(I3,1X))
 
C --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON ---
C --- THE TEMPORARY STACK ---
 
C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT ---
      IF ((KPART .NE. 11) .AND. (KPART .NE. 12)) GO TO 52
      CALL GRNDM(RNDM,1)
      KPART=11.5+RNDM(1)
 
 52   CONTINUE
 
C --- IN CASE THE NEW PARTICLE IS A NEUTRINO ==> FORGET IT ---
      IF (KPART .EQ. 2) GO TO 60
 
C --- PUT CURRENT GHEISHA PARTICLE ON THE CORSIKA STACK
C --- ( IF SURVIVING ANGLE CUT ! )
      NGKINE = 1
      SECPAR(3) = -PZ
 
C --- CALCULATE ELASTICITY
      IF ( EN .GT. EMAX ) THEN
        EMAX = EN
      ENDIF
 
      IF ( SECPAR(3) .GT. C(29) ) THEN
 
         ITY=IKPART(KPART)
         IF     ( ITY .LT. 45 ) THEN
            SECPAR(1) = DBLE(ITY)
         ELSEIF ( ITY .EQ. 45 ) THEN
            SECPAR(1) = 201.D0
         ELSEIF ( ITY .EQ. 46 ) THEN
            SECPAR(1) = 301.D0
         ELSEIF ( ITY .EQ. 47 ) THEN
            SECPAR(1) = 402.D0
         ENDIF
         IF ( ABS(AMAS) .LT. 1.E-9 ) THEN
            SECPAR(2) = EN
         ELSE
            SECPAR(2) = DBLE(EN) / DBLE(ABS(AMAS))
         ENDIF
         IF ( PX .NE. 0. .OR. PY .NE. 0. )  THEN
            SECPAR(4) = ATAN2 ( DBLE(PY) , DBLE(PX) )
         ELSE
            SECPAR(4) = 0.D0
         ENDIF
 
        CALL TSTACK
 
      ENDIF
 
C *** CHECK WHETHER SECONDARIES HAVE BEEN GENERATED AND COPY THEM ***
C *** ALSO ON THE GEANT STACK ***
 60   CONTINUE
 
C --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE ---
C --- CONVENTION IS THE FOLLOWING ---
C
C EVE(INDEX+ 1)= X
C EVE(INDEX+ 2)= Y
C EVE(INDEX+ 3)= Z
C EVE(INDEX+ 4)= NCAL
C EVE(INDEX+ 5)= NCELL
C EVE(INDEX+ 6)= MASS
C EVE(INDEX+ 7)= CHARGE
C EVE(INDEX+ 8)= TOF
C EVE(INDEX+ 9)= PX
C EVE(INDEX+10)= PY
C EVE(INDEX+11)= PZ
C EVE(INDEX+12)= TYPE
 
      IF (NTOT .LE. 0) GO TO 9999
 
C --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED ---
      DO 61 L=1,NTOT
      INDEX=(L-1)*12
      JND=EVE(INDEX+12)
 
C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT ---
      IF ((JND .NE. 11) .AND. (JND .NE. 12)) GO TO 63
      CALL GRNDM(RNDM,1)
      JND=11.5+RNDM(1)
 
C --- FORGET ABOUT NEUTRINOS ---
 63   CONTINUE
      IF (JND .EQ. 2) GO TO 61
 
C --- SWITCH TO CORSIKA QUANTITIES ---
      ITY=IKPART(JND)
      IF (NPRT(9)) WRITE(MDEBUG,1006) ITY,NGKINE,L,(EVE(INDEX+J),J=1,12)
 1006 FORMAT(' *CGHEI* GEANT PART. ',I3,' ALSO PUT ONTO STACK AT',
     $ ' POS. ',I3/
     $ ' EVE(',I2,') = ',(' ',10G12.5))
 
      PLX = EVE(INDEX+9)
      PLY = EVE(INDEX+10)
      PLZ = EVE(INDEX+11)
      PLSQ = PLX**2 + PLY**2 + PLZ**2
      PLTOT = SQRT (PLSQ)
      RMASSK = ABS(RMASS(JND))
 
C  FIND HIGHEST ENERGY PARTICLE FOR ELASTICITY
      EEESQ = PLSQ + RMASSK**2
      IF ( EEESQ .GT. EMAX**2 ) THEN
        EMAX = SQRT(EEESQ)
      ENDIF
 
C --- APPLY ANGLE CUT AND
C --- ADD PARTICLE TO THE CORSIKA STACK (RESTRICTED TO 100) ---
      IF ( PLTOT .LE. 1.D-10 ) GOTO 61
      SECPAR(3) = -PLZ / PLTOT
      IF ( SECPAR(3) .LE. C(29) )  GOTO 61
 
      IF (NGKINE .GE. MXGKGH) GO TO 9999
      NGKINE=NGKINE+1
      IF     ( ITY .LT. 45 ) THEN
         SECPAR(1) = DBLE(ITY)
      ELSEIF ( ITY .EQ. 45 ) THEN
         SECPAR(1) = 201.D0
      ELSEIF ( ITY .EQ. 46 ) THEN
         SECPAR(1) = 301.D0
      ELSEIF ( ITY .EQ. 47 ) THEN
         SECPAR(1) = 402.D0
      ELSE
         SECPAR(1) = 0.D0
         WRITE(MONIOU,*) '*CGHEI*  ILLEGAL PARTICLE TYPE'
      ENDIF
      IF ( RMASSK .LT. 1.D-9 ) THEN
         SECPAR(2) = PLTOT
      ELSE
         SECPAR(2) = SQRT (PLSQ+RMASSK**2) / RMASSK
      ENDIF
      IF ( PLX .NE. 0.D0 .OR. PLY .NE. 0.D0 )  THEN
         SECPAR(4) = ATAN2 ( PLY,PLX )
      ELSE
         SECPAR(4) = 0.D0
      ENDIF
 
      CALL TSTACK
 
 61   CONTINUE
 
C --- COUNTER FOR ENERGY-MULTIPLICITY MATRIX
      MSMM = MSMM + NTOT
 
C --- FILL ELASTICITY IN MATRICES
      ELASTI = EMAX/ENOLD
      MELL = MIN ( 1.D0+10.D0*      MAX( 0.D0, ELASTI ) , 11.D0 )
      MEN  = MIN ( 4.D0+ 3.D0*LOG10(MAX( .1D0, EKINL )), 37.D0 )
      IELDPM(MEN,MELL) = IELDPM(MEN,MELL) + 1
      IELDPA(MEN,MELL) = IELDPA(MEN,MELL) + 1
      IF ( ELASTI .LT. 1. ) THEN
        ELMEAN(MEN) = ELMEAN(MEN) + ELASTI
        ELMEAA(MEN) = ELMEAA(MEN) + ELASTI
      ENDIF
 
 9999 CONTINUE
C --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW ---
      NGKINE=MIN(NGKINE,MXGKGH)
 
      RETURN
      END
