program PMT

C  Author:  Glen Cowan
C  Date:    5-APR-1998

c  The program simulates electron multiplication in a photomultiplier tube.
c  Poisson random numbers are generated with the CERNLIB routine RNPSSN.
c  The program must be linked with the CERN Program Library (including 
c  PACKLIB for HBOOK histogram routines).

      implicit        NONE

c  Constants

      integer         num_dynodes
      parameter       (num_dynodes = 6)

c  Needed for HBOOK routines

      integer         hsize
      parameter       (hsize = 100000)

      integer         hmemor (hsize)
      common  /pawc/  hmemor

c  Local variables

      character*80    outfile

      integer         i
      integer         icycle
      integer         ierror
      integer         istat
      integer         j
      integer         k
      integer         lun
      integer         n
      integer         n_in
      integer         n_out
      integer         num_events

      real            nu(num_dynodes)

c  Initialize HBOOK, open histogram file, book histograms.

      call HLIMIT (hsize)
      lun = 20
      outfile = 'pmt.his'
      call HROPEN (lun, 'histog', outfile, 'N', 1024, istat)
      if (istat .ne. 0) write (*, *) 'HROPEN error, istat = ', istat
      do i = 1, num_dynodes
        call HBOOK1 (i, 'number of electrons', 100, 0., 5000., 0.)
        call HIDOPT (i, 'STAT')
      end do

c  Set the means of the Poisson variables for each dynode.  

      nu(1) = 6.
      do i = 2, num_dynodes
        nu(i) = 3.
      end do

c  Follow each photoelectron through the pmt.  At each dynode,
c  generate a Poisson number representing the number of secondary
c  electrons produced.  Record the total number emerging from each
c  dynode in a histogram.

      write (*, *) 'enter number of events to generate'
      read (*, *) num_events
      do i = 1, num_events
        n_in = 1
        do j = 1, num_dynodes
          n_out = 0
          do k = 1, n_in
            call RNPSSN (nu(j), n, ierror)
            n_out = n_out + n
          end do
          n_in = n_out
          call HF1 (j, FLOAT(n_out), 1.)
        end do
      end do

c  Store histograms

      call HROUT (0, icycle, ' ')
      call HREND ('histog')

      stop
      END