subroutine gheishp
     x(ipart,pvec,mprod,nprod,iprod,tprod,pprod)
c
c |gheisha| includes
#include "gelhad/ghcdes/mxgkgh.inc"
C
C |geant| commons
#include "geant321/gcbank.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gccuts.inc"
#include "geant321/gcflag.inc"
#include "geant321/gcking.inc"
#include "geant321/gcmate.inc"
#include "geant321/gcphys.inc"
#include "geant321/gctmed.inc"
#include "geant321/gctrak.inc"
#include "geant321/gcunit.inc"
#include "geant321/gsecti.inc"
C
C /gheisha/ commons
#include "gelhad/ghcdes/blankp.inc"
#include "gelhad/ghcdes/consts.inc"
#include "gelhad/ghcdes/event.inc"
      integer *4 nprod,mprod,code,stop
      integer *4 iprod(mprod)
      real *4 pprod(3,mprod)
      real *4 tprod(mprod)
      integer *4 ipart
      real *4 pvec(3)
      character *(*) todo
      real *4 arg
      real *4 value
      integer *4 ivalue
      equivalence (value,ivalue)
c
      integer *4 intforce/0/
c
c *** main steering for hadron shower development ***
c *** nve 15-jun-1988 cern geneva ***
c
c called by : guhadr (user routine)
c origin : f.carminati, h.fesefeldt
c                       routines : calim  16-sep-1987
c                                  setres 19-aug-1985
c                                  intact 06-oct-1987
      integer *4 ind
c
c
c
c
c
c --- "nevent" changed to "kevent" in common /curpar/ due to clash ---
c --- with variable "nevent" in geant common ---
c
      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
c ipart->kpart because of name conflict with /geant321/
c
      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
c --- "absl(21)" changed to "abslth(21)" in common /mat/ due to clash ---
c --- with variable "absl" in geant common ---
c
      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)
c
      dimension ipelos(35)
      save ideol
c
c --- random number array                          --
      dimension rndm(1)
c
c --- dimension stmts. for geant321/gheisha particle code conversions ---
c --- kipart(i)=gheisha code corresponding to geant   code i ---
c --- ikpart(i)=geant   code corresponding to gheisha code i ---
c
      integer kipart(48),ikpart(35)
c
c --- data stmts. for geant321/gheisha particle code conversions ---
c --- kipart(i)=gheisha code corresponding to geant   code i ---
c --- ikpart(i)=geant   code corresponding to gheisha code i ---
c
      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/
c
      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
c
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
c --- lowerbound of kinetic energy bin in n cross-section tables ---
      data teklow /0.0001/
c
c --- kinetic energy to switch from "casn" to "gnslwd" for n cascade ---
      data swtekn /0.05/
c
      data ideol/0/
c
c --- set the interaction mechanism to "hadr" ---
c
c --- init output ---
      stop=0
      code=0
      nprod=0
      destep=0.0
c
c
 9004 continue
      kpart=kipart(ipart)
      kkpart=kpart
c
c --- transport the track number to gheisha and initialise some numbers
      ntk=0
      intct=0.0
      next=1
      ntot=0
      int=0
      tof=0.0
c
c --- fill result common for this track with geant values ---
c --- calim code ---
      xend=0.0
      yend=0.0
      zend=0.0
      amas=rmass(kpart)
      nch=rcharg(kpart)
      charge=rcharg(kpart)
      p=sqrt(pvec(1)**2+pvec(2)**2+pvec(3)**2)
      if(p.lt.1.0e-10) then
       px=0.0
       py=0.0
       pz=0.0
       stop=1
      else
      px=pvec(1)/p
      py=pvec(2)/p
      pz=pvec(3)/p
      endif !p.lt.1.0e-10)
c --- setres code ---
      amasq=amas*amas
      en=sqrt(amasq+p*p)
      ek=abs(en-abs(amas))
      enold=en
c
c
      sinl=0.0
      cosl=1.0
      sinp=0.0
      cosp=1.0
c
      if (abs(p) .le. 1.0e-10) go to 1
      sinl=pz
      cosl=sqrt(abs(1.0-sinl**2))
c
 1    continue
      call grndm(rndm,1)
      phi=rndm(1)*twpi
      if ((px .eq. 0.0) .and. (py .eq. 0.0)) goto 3
      if (abs(px) .lt. 1.e-10) goto 2
      phi=atan2(py,px)
      goto 3
c
 2    continue
      if (py .gt. 0.0) phi=pi/2.0
      if (py .le. 0.0) phi=3.0*pi/2.0
c
 3    continue
      sinp=sin(phi)
      cosp=cos(phi)
c
c --- set gheisha index for the current medium always to 1 ---
      ind=1
c
c --- transfer global material constants for current medium ---
c --- detailed data for compounds is obtained via routine compo ---
      atno(ind+1)=a
      zno(ind+1)=z
      den(ind+1)=dens
      radlth(ind+1)=radl
      abslth(ind+1)=absl
c
c --- setup parmat for physics steering ---
      parmat(ind+1,5)=0.0
      parmat(ind+1,8)=ipfis
      parmat(ind+1,9)=0.0
      parmat(ind+1,10)=0.0
      if(jtm.gt.0) then !in some media ?
         jtmn=lq(jtm)
         if(jtmn.ge.1) parmat(ind+1,5)=q(jtmn+26)
      endif !jtm.gt.0
c
c --- check for stopping track ---
      if(stop.ne.0) then
       call ghstopp(ipart,code,stop)
       if(code.eq.5) go  to 9999
       if(ihadr.ne.2) go to 40
       if(ipelos(kpart).eq.0) then
        destep=destep+en !unstable deposit all energy
       else
        destep=destep+ek !stable deposit kin energy only
       endif !ipelos(kpart).eq.0
       go to 9999
      endif !stop.ne.0
      stop=0
c
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
c *** division into various interaction channels denoted by "int" ***
c the convention for "int" is the following
c
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
c --- intact code ---
      kk=abs(q(jma+11)) !number of material components
      alam1=0.0
      call grndm(rndm,1)
      rat=rndm(1)*alam
      atno2=a
      zno2 =z
c
      do 6 k=1,kk
      if (kk .le. 0) go to 6
c
      if (kk .eq. 1) go to 7
      atno2=q(jmixt+k)
      zno2=q(jmixt+kk+k)
c
 7    continue
c
c force selected interaction type
      if(intforce.ne.0) then
       int=intforce
       go to 8
      endif !intforce
c
c --- try for elastic scattering ---
      int=1
      code=13
      alam1=alam1+aiel(k)
      if (rat .lt. alam1) go to 8
c
c --- try for inelastic scattering ---
      int=2
      code=20
      alam1=alam1+aiin(k)
      if (rat .lt. alam1) go to 8
c
c --- try for nuclear fission with inelastic scattering ---
      int=3
      code=15
      alam1=alam1+aifi(k)
      if (rat .lt. alam1) go to 8
c
c --- try for neutron capture ---
      int=4
      code=18
      alam1=alam1+aica(k)
      if (rat .lt. alam1) go to 8
c
 6    continue
c --- no reaction selected ==> elastic scattering ---
      int=1
      code=13
c
c *** take action according to selected reaction channel ***
c --- following code is a translation of "calim" into geant jargon ---
c
 8    continue
c
c --- in case of no interaction or unknown cross sections ==> done ---
      if (int .le. 0) go to 40
c
c --- in case of non-elastic scattering and no generation of sec. ---
c --- particles deposit total particle energy and return ---
      if ((int .eq. 1) .or. (ihadr .ne. 2)) go to 9
      stop=2
      destep=destep+en
      go to 9999
c
 9    continue
      if (int .ne. 4) go to 10
c
c --- neutron capture ---
      stop=1
      call captur(nopt)
      go to 40
c
 10   continue
      if (int .ne. 3) go to 11
c --- nuclear fission ---
      stop=1
      tkin=fissio(ek)
      int=0
      go to 40
c
11    continue
c
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
c --- additional parameters to simulate fermi motion and evaporation ---
      do 111 jenp=1,10
       enp(jenp)=0.0
111   continue
      enp(5)=ek
      enp(6)=en
      enp(7)=p
c
      if (int .ne. 1) go to 12
c
c *** elastic scattering processes ***
c
c --- only nuclear interactions for heavy fragments ---
      if ((kpart .ge. 30) .and. (kpart .le. 32)) go to 35
c
c --- normal elastic scattering for light media ---
      if (atno2 .lt. 1.5) go to 35
c
c --- coherent elastic scattering for heavy media ---
      call coscat
      go to 40
c
c *** non-elastic scattering processes ***
 12   continue
c
c --- only nuclear interactions for heavy fragments ---
      if ((kpart .ge. 30) .and. (kpart .le. 32)) go to 35
c
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
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
c *** in case of energy above cut-off let the particle cascade ***
      test=abs(charge)
      if ((test .gt. 1.0e-10) .and. (enp(5) .gt. cuthad)) go to 35
      if ((test .le. 1.0e-10) .and. (enp(5) .gt. cutneu)) go to 35
c
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
c *** particle with energy below cut-off ***
c --- ==> only nuclear evaporation and quasi-elastic scattering ---
 13   continue
c
      stop=3
c
c
      if ((kpart .ne. 14) .and. (kpart .ne. 16)) go to 14
      if (kpart .eq. 16) go to 86
c
c --- slow proton ---
 85   continue
      call nucrec(nopt,2)
c
      if (nopt .ne. 0) go to 50
c
      call coscat
      go to 40
c
c --- slow neutron ---
 86   continue
      nucflg=0
      call gnslwd(nucflg,int,nfl,teklow)
      if (nucflg .ne. 0) go to 50
      go to 40
c
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
      call twob(ippp,nfl,avern)
      goto 40
c
c --- initialisation of cascade quantities ---
 35   continue
c
c *** cascade generation ***
c --- calculate final state multiplicity and longitudinal and ---
c --- transverse momentum distributions ---
c
c --- fixed particle type to steer the cascade ---
      kkpart=kpart
c
c --- no cascade for leptons ---
      if (kkpart .le. 6) go to 9999
c
c *** what to do with "new particles" for gheisha ?????? ***
c --- return for the time being ---
      if (kkpart .ge. 35) go to 9999
c
c --- cascade of heavy fragments
      if ((kkpart .ge. 30) .and. (kkpart .le. 32)) go to 390
c
c --- initialize the ipa array ---
      call vzero(ipa(1),100)
c
c --- cascade of omega - and omega - bar ---
      if (kkpart .eq. 33) go to 330
      if (kkpart .eq. 34) go to 331
c
      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
c
 15   continue
      nvepar=kkpart-6
      go to (307,308,309,310,311,312,313,314,315,316,317,318),nvepar
c
c --- pi+ cascade ---
 307  continue
      call caspip(j,int,nfl)
      go to 40
c
c --- pi0 ==> no cascade ---
 308  continue
      go to 40
c
c --- pi- cascade ---
 309  continue
      call caspim(j,int,nfl)
      go to 40
c
c --- k+ cascade ---
 310  continue
      call caskp(j,int,nfl)
      go to 40
c
c --- k0 cascade ---
 311  continue
      call cask0(j,int,nfl)
      go to 40
c
c --- k0 bar cascade ---
 312  continue
      call cask0b(j,int,nfl)
      go to 40
c
c --- k- cascade ---
 313  continue
      call caskm(j,int,nfl)
      go to 40
c
c --- proton cascade ---
 314  continue
      call casp(j,int,nfl)
      go to 40
c
c --- proton bar cascade ---
 315  continue
c      if (nprt(9)) print 2013
 2013 format(' *gheish* routine caspb will be called')
      call caspb(j,int,nfl)
      go to 40
c
c --- neutron cascade ---
 316  continue
      nucflg=0
      if (ek .gt. swtekn) call casn(j,int,nfl)
      if (ek .le. swtekn) call gnslwd(nucflg,int,nfl,teklow)
      if (nucflg .ne. 0) go to 50
      go to 40
c
c --- neutron bar cascade ---
 317  continue
      call casnb(j,int,nfl)
      go to 40
c
c --- lambda cascade ---
 318  continue
      call casl0(j,int,nfl)
      go to 40
c
c --- lambda bar cascade ---
 319  continue
c      if (nprt(9)) print 2018
 2018 format(' *gheish* routine casal0 will be called')
      call casal0(j,int,nfl)
      go to 40
c
c --- sigma + cascade ---
 320  continue
      call cassp(j,int,nfl)
      go to 40
c
c --- sigma 0 ==> no cascade ---
 321  continue
      go to 40
c
c --- sigma - cascade ---
 322  continue
      call cassm(j,int,nfl)
      go to 40
c
c --- sigma + bar cascade ---
 323  continue
      call casasp(j,int,nfl)
      go to 40
c
c --- sigma 0 bar ==> no cascade ---
 324  continue
      go to 40
c
c --- sigma - bar cascade ---
 325  continue
      call casasm(j,int,nfl)
      go to 40
c
c --- xi 0 cascade ---
 326  continue
      call casx0(j,int,nfl)
      go to 40
c
c --- xi - cascade ---
 327  continue
      call casxm(j,int,nfl)
      go to 40
c
c --- xi 0 bar cascade ---
 328  continue
      call casax0(j,int,nfl)
      go to 40
c
c --- xi - bar cascade ---
 329  continue
c      if (nprt(9)) print 2026
 2026 format(' *gheish* routine casaxm will be called')
      call casaxm(j,int,nfl)
      go to 40
c
c --- omega - cascade ---
 330  continue
c      if (nprt(9)) print 2027
 2027 format(' *gheish* routine casom will be called')
      call casom(j,int,nfl)
      go to 40
c
c --- omega - bar cascade ---
 331  continue
c      if (nprt(9)) print 2028
 2028 format(' *gheish* routine casaom will be called')
      call casaom(j,int,nfl)
      go to 40
c
c --- heavy fragment cascade ---
 390  continue
      nucflg=0
      call casfrg(nucflg,int,nfl)
      if (nucflg .ne. 0) go to 50
c
c *** check whether there are new particles generated ***
 40   continue
      if ((ntot .ne. 0) .or. (kkpart .ne. kpart)) go to 50
       nprod=1
       iprod(1)=ikpart(kpart)
       pprod(1,1)=px*p
       pprod(2,1)=py*p
       pprod(3,1)=pz*p
       tprod(1)=tprod(1)+tof*0.5e-10
       edep=abs(enold-en)
       if(stop.eq.0) destep=destep+edep
      go to 9999
c
c *** current particle is not the same as in the beginning or/and ***
c *** one or more secondaries have been generated ***
 50   continue
c
      nvedum=kipart(ipart)
c
c --- initial particle type has been changed ==> put new type on ---
c --- the geant temporary stack ---
c
c --- chose beteen ks/kl
      if(kpart.eq.11.or.kpart.eq.12) then
       kpart=1
       call grndm(rndm,1)
       if(rndm(1).gt.0.5) kpart=12
      endif !kpart.eq.11.or.kpart.eq.12
c
c
c --- put particle on the stack ---
       nprod=1
       iprod(1)=ikpart(kpart)
       pprod(1,1)=px*p
       pprod(2,1)=py*p
       pprod(3,1)=pz*p
       tprod(1)=tprod(1)+tof*0.5e-10
c
c
c *** check whether secondaries have been generated and copy them ***
c *** also on the geant stack ***
 60   continue
c
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
c
      if (ntot .le. 0) go to 9999
c
c --- one or more secondaries have been generated ---
      do 61 l=1,ntot
      index=(l-1)*12
      jnd=eve(index+12)
c
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
c --- forget about neutrinos ---
 63   continue
      if (jnd .eq. 2) go to 61
c
c --- swith to geant quantities ---
      ity=ikpart(jnd)
      plx=eve(index+9)
      ply=eve(index+10)
      plz=eve(index+11)
c
c
c --- add particle to the stack if stack not yet full ---
      fail=1301
      if (nprod.ge.mprod) go to 1313
      nprod=nprod+1
      iprod(nprod)=ity
      pprod(1,nprod)=plx
      pprod(2,nprod)=ply
      pprod(3,nprod)=plz
      tprod(nprod)=tprod(nprod)+eve(index+8)*0.5e-10
c
c
 61   continue
c
 9999 continue
      return
c
1313  continue
      write(6,*) 'fail=',fail,' in |gheish|'
      if(fail.eq.1301) then
       write(6,*) 'produced particle array overflowed'
       write(6,*) 'results will be truncated'
      endif !fail.eq.1301
      return
c
      entry gheiset(todo,arg)
      value=arg
      if(todo.eq.'ihadr') then
       ihadr=ivalue
      else if (todo.eq.'ipfis') then
       ipfis=ivalue
      else if(todo.eq.'cuthad') then
       cuthad=value
      else if(todo.eq.'cutneu') then
       cutneu=value
      else if(todo.eq.'intforce') then
       intforce=ivalue
      endif !todo.eq.'ihadr'
      end