c
c tresol.f - simulation of the pulse train in a single channel of the
c            gluex tagger microscope.  The pulses are generated with an
c            exponential t0 and a Landau amplitude distribution, with
c            fixed rise and fall times.  A leading-edge discriminator
c            is used to read off crossing times.  Use of this package
c            is for extracting time resolution and efficiency as a
c            function of count rate in a single channel.
c author: Richard Jones at uconn.edu
c version: Sept. 28, 2009
c
c modifications by R. Jones:
c    1.	Noise was added to simulate dark rate in a silicon photomultiplier.
c	The noise is modeled as a random rate of pixel clusters firing
c	with the same pulse shape as signal pulses.  The pixel cluster
c	size is distributed as Poisson with a mean set by 1+fcrosstalk
c	where fcrosstalk is set based on multi-poisson fits to measured
c	noise pulse-height spectra.
c    2.	Whereas before the most-probable pulse height of the signal pulses
c	was specified in arbitrary units, now the units must be normalized
c	by some convention so that the relative size of noise and signal
c	pulses can be properly adjusted.  Pulse height shall henceforth be
c	specified in units of photoelectrons (SiPM pixels).
c version: Sept. 20, 2010
c
c modifications by R. Jones:
c    3. Before I was generating the entire pulse corresponding to a tag
c       electron by rescaling the single-pixel pulse by an amplitude
c       factor.  What this misses is the contribution to the spread in
c       the leading edge timing coming from photoelectron statistics.
c       I had tried to incorporate the broadening effects of the scintillation
c       decay time constant on the pulse shape by introducing a broadened
c       pulse called scint_pulseN (N is the extra attenuation introduced
c       to simulate the effect of N meters of RG-58 cable between the preamp
c       and the discriminator) but that has very little effect, only slightly
c       increasing the rise time of the pulse.
c                 * sipm_pulse: rise time 2.45 ns
c                 * wiki_pulse: rise time 2.43 ns
c                 * scint_pulse0: rise time 2.70 ns
c       The change made here is to switch from generating the signal as
c       single pulses of varying magnitude, I now generate them as a cluster
c       of many individual pulses, each corresponding to the pulse height of
c       a single pixel.  This increases the run time a great deal, but it is
c       the only way to incorporate statistical time resolution effects on
c       the same basis as electronic noise due to dark rate contributes.
c    4. Read the ranlux random number generator seed from a file randomseed
c       if it exists in the current working directory and is readable.
c version: Dec. 14, 2010
c
c modifications by R. Jones:
c    5. Added after-pulsing to the simulation.  Multiple exponentials are
c	supported in the afterpulsing model, through parameters
c		* nafterpulse - number of exponentials
c		* tafterpulse - after-pulsing time constant
c		* fafterpulse - mean after-pulses per pixel
c    6. Serialized the pixels generated in a single pulse.  Before I had
c	been generating the individual pixels in the simulation independently
c	according to the exponential decay distribution for the scintillator.
c       This has the disadvantage that the dead-pixel correction is applied to
c	the overall yield in the pulse, and not preferentially to the later
c	part of the pulse as required by the physical model.  Now they are
c	generated in time order, and the dead pixel correction updated as
c	each occurs.
c    7. Added cross-talk hits to the signal pulse as well as to the noise
c	pulses.  Before I had included these only for the noise pulses.
c version: Dec. 29, 2010

      subroutine tresol(a0_,a1_,trise_,tfall_,R_,Rnoise_,dtmin_)
      implicit none
      real a0_,a1_,trise_,tfall_,R_,Rnoise_,dtmin_
      include 'tresol.inc'
      real gamma
      external gamma
      integer n
      data idtrace,idtres,idgen,iddet/10,11,12,13/
      data tracelen,tmemory/1e-6,2e-6/
      data scintdecayt/2.7e-9/ ! BCF-20 fast-green plastic scintillator
c========================================================================
c This section describes the CPTA 2x2 mm^2 SSPM from Photonique 
c operated at 2V above the reverse-bias breakdown potential
      data Rnoise/9e6/ !Hz
      data trecovery/1e-6/ !s
      data fcrosstalk/0.13/
      data Npixels/1700/
      data nafterpulse/0/
c========================================================================
c This section describes the Hamamatus S10931 3x3 mm^2 MPPC (25 micron pixels)
c operated at 1V above the reverse-bias breakdown potential
c     data Rnoise/4e6/ !Hz
c     data trecovery/15e-9/ !s
c     data fcrosstalk/0.02/
c     data Npixels/14400/
c     data nafterpulse/2/
c     data fafterpulse/0.06,0.05,7*0/
c     data tafterpulse/18e-9,70e-9,7*0/
c========================================================================
      a0 = a0_
      a1 = a1_
      trise = trise_
      tfall = tfall_
      R = R_
      Rnoise = Rnoise_
      dtmin = dtmin_
      tstep = tracelen/max_samples
      ppar(1) = tfall/5
      ppar(2) = 10*trise/tfall
      ppar(3) = exp(ppar(2))/(1e9*ppar(1)*ppar(2))**ppar(2)
      ppar(4) = 1.5e-6
      ppar(5) = 1.0e-7
      ppar(6) = gamma(ppar(2)+1)*(ppar(1)*1e9)**(ppar(2)+1)/
     +          (ppar(4)-ppar(5))/1e9
      if (Npixels.gt.max_pixels) then
        write(6,*) 'tresol error - Npixels',Npixels,
     +             ' must be less than max_pixels',max_pixels
        return
      endif
      call get_pulse_characteristics
      write(6,*) 'tresol - simulation of pulse train in a single',
     +           ' microscope channel'
      write(6,*) '   MP amplitude = ',a0,' photoelectrons'
      write(6,*) '   scintillator time constant = ', scintdecayt,' s'
      write(6,*) '   discriminator threshold = ',a1, ' photoelectrons'
      write(6,*) '   minimum pulse separation = ',dtmin,' s'
      write(6,*) '   pulse rise time = ',trise,' s'
      write(6,*) '   pulse fall time = ',tfall,' s'
      write(6,*) '   pulse peak time = ',dtmax,' s'
      write(6,*) '   pulse 50% delay = ',dtcross,' s'
      write(6,*) '   mean pulse rate = ',R,' Hz'
      write(6,*) '   full trace length = ',tracelen,' s'
      write(6,*) '   time step length = ',tstep,' s'
      write(6,*) '   dark rate = ', Rnoise,' Hz'
      write(6,*) '   pixel recovery time = ',trecovery, ' s'
      write(6,*) '   cross-talk fraction = ',fcrosstalk
      do n=1,nafterpulse
        write(6,1000) fafterpulse(n),tafterpulse(n)
 1000   format(f8.2,' after-pulses with time constant = ',e12.5,' s')
      enddo
      call seedranlux
      end

      subroutine simulate(nrep)
      implicit none
      integer nrep
      include 'tresol.inc'
      logical hexist
      external hexist
      real ranlan
      external ranlan
      real tdecay
      real t0
      integer Npe
      integer ip
      integer ierr
      integer irep,i
      integer next_pulse
      real rnd01(max_pixels)
      real unif01(2)
      call deftrace
      call deftres
      call defeffic
      do irep=1,nrep
        do i=1,max_samples
          sample(i) = 0
        enddo
        t0 = -tmemory
        last_event = 0
        next_pulse = 0
        do while (t0 .lt. tracelen)
          next_pulse = next_pulse+1
          call ranlux(unif01,2)
          t0 = t0-log(unif01(1))/R
          if (t0 .lt. tracelen) then
c
c The factor 0.17 has been adjusted to obtain a good description of the
c fluctations in energy deposition by electrons in 2cm of scintillator.
c
            Npe = a0*(1+0.17*ranlan(unif01(2)))
            Npe = min(Npe,max_pixels) ! truncate the Landau tail
            if (verbose) then
              write(6,*) 'pulse amplitude ',Npe
            endif
            call ranlux(rnd01,Npe)
            pcount = 0
            tdecay = 0
            do ip=1,Npe
              tdecay = tdecay-log(rnd01(ip))*scintdecayt/(Npe-ip+1)
              call add_1photon(t0+tdecay)
            enddo
            if (t0+dtcross.gt.0 .and. t0+dtmax.lt.tracelen) then
              last_event = last_event+1
              if (last_event .ge. max_events) then
                write(6,*) 'tresol::simulate error - ',
     +                     'overflow of event arrays, ',
     +                     'please increase the value of max_events'
                return
              endif
              event_t0(last_event) = t0
              event_amp(last_event) = Npe
              call hfill(idgen,float(Npe),0.,1.)
              if (verbose) then
                print 1000, last_event,t0,Npe,pcount
 1000           format(' event',i6,' with t0,Npe,pcount =',
     +                 e12.5,' s,',i5,',',f5.0)
              endif
            endif
          endif
        enddo
        event_t0(last_event+1) = -9999
        event_amp(last_event+1) = -9999
        call add_noise
        call hpak(idtrace,sample)
        call analyze_trace
      enddo
      end

      subroutine add_1photon(t0)
      implicit none
      real t0
      include 'tresol.inc'
      real pulse
      external pulse
      real deadpixels
      real lastt0
      common /deadtime/deadpixels,lastt0
      data deadpixels/0/
      data lastt0/0/
      integer ierr
      integer i,i1,i2
      integer Npe,nafter,n
      real rnd01(max_pixels)
      real lambda
      real t,amp
      real t1
      if (t0.lt.lastt0) then
        deadpixels = Rnoise*trecovery
      else
        deadpixels = Rnoise*trecovery + (deadpixels-Rnoise*trecovery)
     +               *exp(-(t0-lastt0)/trecovery)
      endif
      deadpixels = min(nint(deadpixels),Npixels-1)
      lambda = (1+fcrosstalk)*(1-deadpixels/Npixels)
      call RNPSSN(lambda,Npe,ierr)
      if (ierr .ne. 0) then
        write(6,*) 'tresol::add_1photon - ',
     +             'error ',ierr,' returned from RNPSSN, ',
     +             'called with mean',lambda,', cannot continue!'
        return
      elseif (Npe .gt. Npixels) then
        Npe = Npixels
      elseif (Npe .le. 0) then
        return
      endif
      deadpixels = deadpixels+Npe
      i1 = t0/tstep+0.9999
      i1 = max(i1,1)
      i2 = i1+tmemory/tstep
      i2 = min(i2,int(tracelen/tstep))
      do i=i1,i2
        t = tstep*(i-0.5)
        sample(i) = sample(i)+Npe*pulse(t-t0)
      enddo
      lastt0 = t0
      pcount = pcount+Npe
      do nafter=1,nafterpulse
        if (fafterpulse(nafter) .gt. 0) then
          call RNPSSN(fafterpulse(nafter),Npe,ierr)
          if (ierr .ne. 0) then
            write(6,*) 'tresol::add_1photon - ',
     +                 'error ',ierr,' returned from RNPSSN, ',
     +                 'called with mean',fafterpulse(nafter),
     +                 'cannot continue!'
          endif
          Npe = min(Npe,max_pixels)
          call ranlux(rnd01,Npe)
          do n=1,Npe
            t1 = t0-log(rnd01(n))*tafterpulse(nafter)/(Npe-n+1)
            i1 = t1/tstep+0.9999
            i1 = max(i1,1)
            i2 = i1+tmemory/tstep
            i2 = min(i2,int(tracelen/tstep))
            amp = 1-exp(-trecovery*(t1-t0))
            do i=i1,i2
              t = tstep*(i-0.5)
              sample(i) = sample(i)+amp*pulse(t-t1)
            enddo
            t0 = t1
          enddo
          pcount = pcount+Npe
        endif
      enddo
      end

      subroutine add_noise
      implicit none
      include 'tresol.inc'
      real pulse
      external pulse
      real urandom(max_events)
      integer prandom
      integer ierr
      real sample_noise(max_samples)
      save sample_noise
      real twidth
      integer Nnoise
      data Nnoise/0/
      save Nnoise
      integer n
      integer i,i1,i2
      real t,t0
      real amp
      twidth = tracelen+tmemory
      if (Nnoise .ne. int(Rnoise*twidth+0.5)) then
        print *, 'generating noise trace...'
        Nnoise = Rnoise*twidth+0.5
        if (Nnoise .ge. max_events) then
          write(6,*) 'tresol::add_noise error - ',
     +               'overflow of event arrays, ',
     +               'please increase the value of max_events'
          return
        endif
        call ranlux(urandom,Nnoise)
        do i=1,max_samples
          sample_noise(i) = 0
        enddo
        do n=1,Nnoise
          t0 = urandom(n)*twidth-tmemory
          call RNPSSN(fcrosstalk,prandom,ierr)
          if (ierr .ne. 0) then
             write(6,*) 'tresol::add_noise error - ',
     +                  'error ',ierr,' returned from RNPSSN, ',
     +                  'called with mean',fcrosstalk,
     +                  'cannot continue!'
             return
          endif
          amp = 1+prandom
          i1 = t0/tstep+0.9999
          i1 = max(i1,1)
          i2 = i1+tmemory/tstep
          i2 = min(i2,int(tracelen/tstep))
          do i=i1,i2
            t = tstep*(i-0.5)
            sample_noise(i) = sample_noise(i)+amp*pulse(t-t0)
          enddo
        enddo
        print *, '...done generating noise trace with',Nnoise,' pulses.'
      endif
      call ranlux(urandom,1)
      t0 = urandom(1)*tracelen
      i1 = t0/tstep+0.9999
      i2 = tracelen/tstep+0.9999
      do i=1,i2
        sample(i) = sample(i)+sample_noise(mod(i1+i,i2)+1)
      enddo
      end

      subroutine analyze_trace
      implicit none
      include 'tresol.inc'
      integer this_event
      integer i,imax
      real tcross
      real f
      this_event = 1
      i = 1
      do while (i .lt. max_samples-1)
        if (sample(i).lt.a1) then
          if (sample(i+1).gt.a1) then
            f = (a1-sample(i))/(sample(i+1)-sample(i))
            tcross = (i+f-0.5)*tstep
            do while (abs(tcross-event_t0(this_event)) .gt.
     +                abs(tcross-event_t0(this_event+1)))
              this_event = this_event + 1
            enddo
            call hfill(iddet,event_amp(this_event),0.,1.)
            imax = i+nint((dtmax-dtcross)/tstep)
            call hfill(idtres,sample(min(imax,max_samples)),
     +                 tcross-event_t0(this_event)-dtcross,1.)
            if (verbose) then
              print *, 'pulse with amplitude',
     +                 sample(min(imax,max_samples)),
     +                 ' has tcross gen,rec =',
     +                 event_t0(this_event)+dtcross,tcross
            endif
            i = i+dtmin/tstep
          endif
        endif
        i = i+1
      enddo
      end

      subroutine deftrace
      implicit none
      include 'tresol.inc'
      logical hexist
      external hexist
      character*80 title
      if (hexist(idtrace)) then
        call hdelet(idtrace)
      endif
      write(title,*) 'trace of pulse train'
      call hbook1(idtrace,title,max_samples,0,tracelen,0)
      end

      subroutine deftres
      implicit none
      include 'tresol.inc'
      logical hexist
      external hexist
      character*80 title
      if (hexist(idtres)) then
        call hdelet(idtres)
      endif
      write(title,*) 'time resolution vs pulse height'
      call hbook2(idtres,title,100,0.,a0*10,200,-5e-9,5e-9,0)
      end

      subroutine defeffic
      implicit none
      include 'tresol.inc'
      logical hexist
      external hexist
      character*80 title
      if (hexist(idgen)) then
        call hdelet(idgen)
        call hdelet(iddet)
      endif
      write(title,*) 'pulses generated'
      call hbook1(idgen,title,200,0,a0*10,0)
      write(title,*) 'pulses detected'
      call hbook1(iddet,title,200,0,a0*10,0)
      end

      real function dflandau(x)
      implicit none
      real x
      include 'tresol.inc'
      real denlan
      external denlan
      if (x .gt.-5*a0) then
        dflandau = denlan((x-a0)/(0.17*a0))
      else
        dflandau = 0
      endif
      end

      subroutine get_pulse_characteristics
      implicit none
      include 'tresol.inc'
      real pulse
      external pulse
      real t10up,t50up,t90up,tpeak
      real t90down,t10down
      real t,f
      integer i
      real s0,s1,s2
      s0 = 0
      do i=1,max_samples
        t = (i-0.5)*tstep
        s2 = s1
        s1 = s0
        s0 = pulse(t)
        if (s1 .lt. 0.10 .and. s0 .gt. 0.10) then
           f = (s0-0.10)/(s0-s1)
           t10up = t-f*tstep
           goto 1
        endif
      enddo
    1 continue
      do i=i+1,max_samples
        t = (i-0.5)*tstep
        s2 = s1
        s1 = s0
        s0 = pulse(t)
        if (s1 .lt. 0.50 .and. s0 .gt. 0.50) then
           f = (s0-0.50)/(s0-s1)
           t50up = t-f*tstep
           goto 2
        endif
      enddo
    2 continue
      do i=i+1,max_samples
        t = (i-0.5)*tstep
        s2 = s1
        s1 = s0
        s0 = pulse(t)
        if (s1 .lt. 0.90 .and. s0 .gt. 0.90) then
           f = (s0-0.90)/(s0-s1)
           t90up = t-f*tstep
           goto 3
        endif
      enddo
    3 continue
      do i=i+1,max_samples
        t = (i-0.5)*tstep
        s2 = s1
        s1 = s0
        s0 = pulse(t)
        if (s0 .lt. s1) then
           f = 0.5*(s0-s2)/(2*s1-s0-s2)
           tpeak = t+(f-1)*tstep
           goto 4
        endif
      enddo
    4 continue
      do i=i+1,max_samples
        t = (i-0.5)*tstep
        s2 = s1
        s1 = s0
        s0 = pulse(t)
        if (s0 .lt. 0.90 .and. s1 .gt. 0.90) then
           f = (s0-0.90)/(s0-s1)
           t90down = t-f*tstep
           goto 5
        endif
      enddo
    5 continue
      do i=i+1,max_samples
        t = (i-0.5)*tstep
        s2 = s1
        s1 = s0
        s0 = pulse(t)
        if (s0 .lt. 0.10 .and. s1 .gt. 0.10) then
           f = (s0-0.10)/(s0-s1)
           t10down = t-f*tstep
           goto 6
        endif
      enddo
    6 continue
      trise = t90up-t10up
      tfall = t10down-t90down
      dtmax = tpeak
      dtcross = t50up
      end

      real function pulse(t)
      implicit none
      real t
      real sipm_pulse
      real wiki_pulse
      real scint_pulse0
      real scint_pulse30
      real scint_pulse55
*     pulse = sipm_pulse(t)
      pulse = wiki_pulse(t)
*     pulse = scint_pulse0(t)
*     pulse = scint_pulse30(t)
*     pulse = scint_pulse55(t)
      end

      real function sipm_pulse(t)
      implicit none
      real t
      include 'tresol.inc'
      if (t .gt. 0) then
        sipm_pulse = ppar(3)*(exp(-t/ppar(1))*(t*1e9)**ppar(2)
     +             - ppar(6)*(exp(-t/ppar(4))-exp(-t/ppar(5)))) 
      else
        sipm_pulse = 0
      endif
      end

      real function wiki_pulse(t)
      implicit none
      real t
      include 'tresol.inc'
      real para(3)
      data para/1.00,9.50e-9,1.44e-9/
      real parb(4)
      data parb/-56.380,1.6168e10,-1.4996e18,4.5625e25/
      real parc(4)
      data parc/4.66,-4.48e8,-4.15,-5.13e7/
      real pard(6)
      data pard/-3.9944,0.43269E+07,49196.,
     +          -0.42559E-02,0.60271E+06,-0.74869E+06/
      if (t.lt.para(2)) then
        wiki_pulse = para(1)*exp(-0.5*((t-para(2))/para(3))**2)
      else if (t.lt.11e-9) then
        wiki_pulse = parb(1)+parb(2)*t+parb(3)*t**2+parb(4)*t**3
      else if (t.lt.32e-9) then
        wiki_pulse = exp(parc(1)+parc(2)*t)+exp(parc(3)+parc(4)*t)
      else
        wiki_pulse = pard(1)*exp(-t*pard(2))*sin((t-5e-8)*pard(3)) +
     +               pard(4)*exp(-t*pard(5))*sin((t-5e-8)*pard(6))
      endif
c     wiki_pulse = wiki_pulse - (1+tanh(t*1e8-5))*
c    +             (exp(-t/ppar(4))-exp(-t/ppar(5)))
c    +             /(ppar(4)-ppar(5))*2.58e-9
      end

      real function scint_pulse0(t)
      implicit none
      real t
      include 'tresol.inc'
      real para(3)
      data para/2.50e-3,22.1e-9,1.53e-9/
      real parb(4)
      data parb/-0.64881,0.81319e8,-0.33661e16,0.46161e23/
      real parc(4)
      data parc/3.16,-3.99e8,-1.15,-2.54e8/
      real pard(6)
      data pard/-1.146e-4,0.23818e+7,-0.17975e+8,
     +           0.57134e-1,0.10459e+9,-0.55833e+6/
      if (t.lt.para(2)) then
        scint_pulse0 = para(1)*exp(-0.5*((t-para(2))/para(3))**2)
      else if (t.lt.25e-9) then
        scint_pulse0 = parb(1)+parb(2)*t+parb(3)*t**2+parb(4)*t**3
      else if (t.lt.42e-9) then
        scint_pulse0 = exp(parc(1)+parc(2)*t)+exp(parc(3)+parc(4)*t)
      else
        scint_pulse0 = pard(1)*exp(-t*pard(2))
     +                        *sin(abs((t-5e-8)*pard(3))**0.5)
     +               + pard(4)*exp(-t*pard(5))
     +                        *sin(abs((t-5e-8)*pard(6))**0.5)
      endif
      scint_pulse0 = scint_pulse0/0.0026
      end

      real function scint_pulse30(t)
      implicit none
      real t
      include 'tresol.inc'
      real para(3)
      data para/2.21e-3,22.4e-9,1.60e-9/
      real parb(4)
      data parb/-0.41921,0.51662e8,-0.20968e16,0.28141e23/
      real parc(4)
      data parc/3.37,-4.06e8,-2.84,-1.94e8/
      real pard(6)
      data pard/-1.0555e-4,0.23649e+7,-0.17836e+8,
     +           0.75701e-1,0.77241e+8,-4.9260e+4/
      if (t.lt.para(2)) then
        scint_pulse30 = para(1)*exp(-0.5*((t-para(2))/para(3))**2)
      else if (t.lt.25e-9) then
        scint_pulse30 = parb(1)+parb(2)*t+parb(3)*t**2+parb(4)*t**3
      else if (t.lt.42e-9) then
        scint_pulse30 = exp(parc(1)+parc(2)*t)+exp(parc(3)+parc(4)*t)
      else
        scint_pulse30 = pard(1)*exp(-t*pard(2))
     +                         *sin(abs((t-4.8e-8)*pard(3))**0.5)
     +                + pard(4)*exp(-t*pard(5))
     +                         *sin(abs((t-4.8e-8)*pard(6))**0.5)
      endif
      scint_pulse30 = scint_pulse30/0.0022
      end

      real function scint_pulse55(t)
      implicit none
      real t
      include 'tresol.inc'
      real para(3)
      data para/1.96e-3,22.6e-9,1.62e-9/
      real parb(4)
      data parb/-0.41487,0.50456e8,-0.20235e16,0.26860e23/
      real parc(4)
      data parc/2.23,-3.49e8,-5.15,-1.30e8/
      real pard(6)
      data pard/-1.0523e-4,0.24084e+7,-0.17922e+8,
     +           0.46463e-1,0.53160e+8,-2.0966e+4/
      if (t.lt.para(2)) then
        scint_pulse55 = para(1)*exp(-0.5*((t-para(2))/para(3))**2)
      else if (t.lt.27e-9) then
        scint_pulse55 = parb(1)+parb(2)*t+parb(3)*t**2+parb(4)*t**3
      else if (t.lt.42e-9) then
        scint_pulse55 = exp(parc(1)+parc(2)*t)+exp(parc(3)+parc(4)*t)
      else
        scint_pulse55 = pard(1)*exp(-t*pard(2))
     +                         *sin(abs((t-5e-8)*pard(3))**0.5)
     +                + pard(4)*exp(-t*pard(5))
     +                         *sin(abs((t-5e-8)*pard(6))**0.5)
      endif
      scint_pulse55 = scint_pulse55/0.0020
      end

      subroutine changeid(id)
      implicit none
      integer id
      include 'tresol.inc'
      if (id .gt. 0) then
        idtrace = id
        idtres = id+1
        idgen = id+2
        iddet = id+3
      else
        write(6,*) 'idtrace = ',idtrace
        write(6,*) 'idtres = ',idtres
        write(6,*) 'idgen = ',idgen
        write(6,*) 'iddet = ',iddet
      endif
      end

      subroutine seedranlux
      implicit none
      include 'tresol.inc'
      integer iseq
      open(unit=1,file='randomseed',status='old',err=9)
      read(1,*,err=8) iseq
      call rluxgo(2,iseq,0,0)
    8 close(unit=1)
    9 continue
      end
*
* $Id: rnpssn.F,v 1.2 2005/04/19 08:34:14 mclareni Exp $
*
* $Log: rnpssn.F,v $
* Revision 1.2  2005/04/19 08:34:14  mclareni
* Protection against possible invalid negative number submitted by Piotr Niezurawski
*
* Revision 1.1.1.1  1996/04/01 15:02:56  mclareni
* Mathlib gen
*
*
*#include "gen/pilot.h"
      SUBROUTINE RNPSSN(AMU,N,IERR)
      EXTERNAL RANLUX
 
      SAVE EMU,AMU0,AMAX
 
      DATA AMU0 /-12345.67/, AMAX /88/
      PARAMETER (AMXA = 88)
 
      IERR=0
      IF(AMU .LE. 0) THEN
       IERR=1
       J=0
      ELSEIF(AMU .GT. AMAX) THEN
       CALL RNORMX(R,1,RANLUX)
       J=R*SQRT(AMU)+AMU+0.5
      ELSE
       IF(AMU .NE. AMU0) THEN
        AMU0=AMU
        EMU=EXP(-AMU)
       ENDIF
       P=1
       J=-1
    1  J=J+1
       CALL RANLUX(R,1)
       P=P*R
       IF(P .GT. EMU) GO TO 1
      ENDIF
* PN
      IF (J.LT.0) THEN
       PRINT *,' RNPSSN: Warning: J<0; J=',J
       PRINT *,'         Correction: J=0'
       PRINT *,'         Increase AMAX value!'
       J=0
      ENDIF		       		
* PN
      N=J
      RETURN
 
*     ENTRY RNPSET(AMX)
*     AMAX=MIN(AMX,AMXA)
*     WRITE(6,'(/7X,''+++++ CERN V136 RNPSSN :  SWITCH TO '',
*    1 ''NORMAL APPROXIMATION FOR      AMU > '',F7.2/)') AMAX
*     RETURN
      END