subroutine isim(chdet,chmeth)
c     implicit none
      character*(*) chdet,      ! 'hybrid', 'lg'
     &              chmeth      ! 'simfast', 'libs', 'libr' - seq. or random lib. reading
      real          zmcplane, x, y, nx, ny
      logical sim_inited
      data    sim_inited /.false./
      save    sim_inited
#include "geant321/gckine.inc"
#include "geant321/gctrak.inc"
#include "cphoto.inc"
*
*
      zmcplane  = 625.0506
*
      if(.not.sim_inited) then
        sim_inited = .true.
        call init_sim_geom(chdet,zmcplane)
      endif
*
      if (VECT(6).ne.0.) then
        nx  = VECT(4)/VECT(6)
        ny  = VECT(5)/VECT(6)
      else
        nx  = 0.
        ny  = 0.
      endif
      x = VECT(1) + nx*(VECT(3)-zmcplane)
      y = VECT(2) + ny*(VECT(3)-zmcplane)
*
      if(ipart.ge.1.and.ipart.le.3) then
        if(chmeth.eq.'simfast') then
          icase = 0
          call addhits_sim(x,y,nx,ny,getot)
          if(igeometry.eq.1) then
            icase = 1
            call addhits_sim(x,y,nx,ny,getot)
          endif
        else
          icase = 0
          call addhits_sim(x,y,nx,ny,getot)
          if(igeometry.eq.1) then
            icase = 1
            if(chmeth.eq.'libr') then
              call addhits_lib(x,y,nx,ny,getot,'r')
            elseif(chmeth.eq.'libs') then
              call addhits_lib(x,y,nx,ny,getot,'s')
            else
              print*, 'unknown option ', chmeth
              call exit(1)
            endif
          endif
        endif
      else if(ipart.eq.4) then
      else if(ipart.eq.7.or.ipart.eq.17) then
      else if(ipart.ge.5.and.ipart.le.6) then
          icase = 0
          call addhits_mip(x,y,nx,ny,getot)
          if(igeometry.eq.1) then
            icase = 1
            call addhits_mip(x,y,nx,ny,getot)
          endif
      else if(ipart.le.32.and.charge.eq.0) then
          icase = 0
          call addhits_had(x,y,nx,ny,getot)
          if(igeometry.eq.1) then
            icase = 1
            call addhits_had(x,y,nx,ny,getot)
         endif
      else if(ipart.le.32.and.charge.ne.0) then
          icase = 0
          call addhits_mip(x,y,nx,ny,getot)
          call addhits_had(x,y,nx,ny,getot)
          if(igeometry.eq.1) then
            icase = 1
            call addhits_mip(x,y,nx,ny,getot)
            call addhits_had(x,y,nx,ny,getot)
         endif
      endif
*
      return
      end
*
*
* define geometry and cuts here
*
      subroutine init_sim_geom(chdet,z)
      implicit none
*
      character*(*) chdet       ! 'hybrid', 'lg'
      real z        ! [cm] distance from target to calorimeter
*
#include "cphoto.inc"
#include "phot_cell_cbk_mc.inc"
#include "adcgam_bk.inc"
*
        if(chdet.eq.'hybrid') then
          igeometry = 1
        elseif(chdet.eq.'lg') then
          igeometry = 0
        else
          print*, 'init_sim_geom err: Unknowon option chdet=  ', chdet
          call exit(1)
        endif
*
        dtarg(1,1) = 2.077/2.   ! HyCal PWO semisizes
        dtarg(1,2) = 2.075/2.
        dtarg(1,3) = 18./2.
        NROW(1)    = 36         ! 73x73 modules (151cmx151cm)
        NCOL(1)    = 36
c       icase      = 1          ! PWO
*
        dtarg(0,1) = 3.815/2.   ! HyCal LG semisizes
        dtarg(0,2) = 3.815/2.
        dtarg(0,3) = 45./2.
        NROW(0)    = 30         ! R = 120cm
        NCOL(0)    = 30
c       icase      = 0          ! LG
*
        a(1) = 0.942003         ! profile parameters and sigma reduction factor
        b(1) = 0.107323
        c(1) = -1.0674
        d(1) = 1.42019
        e(1) = 1.1254
        f(1) = 1.21292
        sfactor(1) = 0.007      ! 0.022
*
        a(0) = .796
        b(0) = 0.21
        c(0) = -.096
        d(0) = 2.
        e(0) = .3
        f(0) = 0.0516
        sfactor(0) = 0.1
*
*     cut and parameter values:
*
      zhycal     = z          ! [cm]
      min_energy = 0.1        ! [GeV] discard clusters with energy below this threshold
      global_gain= 1.000      ! global energy gain factor ("calibration")
      delch      = 4.         ! chi2 improvement cut for tgamma
      xm2cut     = 1.2        ! inv. mass over e cut to split hits
      minpk      = 9          ! [10MeV] min cell energy to be good for the peak
      chisq1     = 3          ! min chi2 to try to split at prelimnirary stage
      chisq2     = 1.4 ! 0.4  ! min chi2 to try to split at final stage
      event_flag = 1          ! break event loop if 0
      event_display = 0       ! whether to run event display
*
*     Notes:
*      1) hole size is hardcoded to 5x5 cells for PWO and 3x3 cells for LG
*      2) id_adcgam now is 10*col+row of central id of a cluster
*      3) type_adcgam now: 0 signle; +10 from multimaxima island; +1 or +2 from split
*
      return
      end
*