subroutine tres(iseed)
      integer iseed
      call rluxgo(2,iseed,0,0)
      end

      real function tresgen(id,nphot,ncast)
      integer id,ncast,nphot
      common /ntup/np,fdisc,tdisc,time(10000)
      integer np
      real fdisc,tdisc,time
      character*80 ntform
      data ntform/'nphot[0,10000]:I,fdisc:R,tdisc:R,time(nphot):R'/
      integer n,i
      real t(10000)
      integer st(10000)
      real rnd(10000)
      logical hexist
      external hexist
      if (nphot.gt.10000) then
        print *, 'tres only compiled for nphot up to 10000'
        tres=0
        return
      endif
      if (hexist(id)) then
        call hdelet(id)
      endif
      call hbnt(id,'discriminator simulation',' ')
      call hbname(id,'ntup',np,ntform)
      fdisc=0.30
      np=nphot
      do n=1,ncast
        call ranlux(rnd,nphot)
        do i=1,nphot
          t(i)=-log(rnd(i))
        enddo
        tdisc=discriminate(fdisc,t,nphot)
        call sortzv(t,st,nphot,1,0,0)
        do i=1,nphot
          time(i)=t(st(i))
        enddo
        call hfnt(id)
      enddo
      call hrout(id,icycle,' ')
      end

      real function discriminate(fdisc,tlist,nt)
      real fdisc,tlist(*)
      integer nt
      real tstep
      integer tbins
      parameter (tstep=5e-3,tbins=1000)
      real pulse(tbins)
      real f
      integer i,j
      do j=1,tbins
        pulse(j)=0
      enddo
      do i=1,nt
        do j=1,tbins
          pulse(j)=pulse(j)+fpulse(j*tstep-tlist(i))
        enddo
      enddo
      do j=1,tbins-1
        if (pulse(j+1).gt.nt*fdisc) goto 10
      enddo
   10 continue
      f=(nt*fdisc-pulse(j))/(pulse(j+1)-pulse(j))
      discriminate=(j+f)*tstep
      end

      real function fpulse(t)
      real t
      real trise,tfall
      data trise,tfall/1.0,5.0/
      real tmax
      if (t.lt.0) then
        fpulse=0
        return
      endif
      fpulse=exp(-t/tfall)-exp(-t/trise)
      tmax=(log(tfall)-log(trise))/(1/trise-1/tfall)
      fpulse=fpulse/(exp(-tmax/tfall)-exp(-tmax/trise))
      end