SUBROUTINE HIST_2 ( NPTS, DATA, PEAK, DELTA, TITLE, PLOT )
C     Make a 'smoothed' histogram of the data in (DATA(i),i=1,NPTS).
C     The histogram ranges from XMIN to XMAX, in NB bins.
C     A bin is incremented if the bin is within DELTA of the value.
C
C     This subroutine returns the peak in the smoothed histogram.
C
      IMPLICIT UNDEFINED (A-Z)
      INTEGER *4 NB
      PARAMETER ( NB = 200 )
      INTEGER *4 NPTS, I, J, NB, IPEAK
      REAL *4 DATA(*), XMIN, XMAX, PEAK, DELTA
      REAL *4 YMIN, YMAX, XHI, XLO, YHI, YLO, XBAR(2), YBAR(2)
      REAL *4 HIST(NB), XDATA(NB), YPEAK
      CHARACTER*(*) TITLE
      LOGICAL PLOT
      real*8 matrixa(3,3),matrixb(3,1),coeffa,coeffb,coeffc
      real*8 sn,sy,sx,sxx,sxxx,sxxxx,sxy,sxxy,lhist,shxdata
      integer*4 gjn,gjnp,gjm,gjmp
      real*4 xplot(nb),yplot(nb),chdelta,shpeak
C-----------------------------------------------------------------------
      XMIN = DATA(1)
      XMAX = DATA(1)
      DO I = 1, NPTS
          XMIN = MIN ( XMIN, DATA(I) )
          XMAX = MAX ( XMAX, DATA(I) )
      END DO
C-----------------------------------------------------------------------
      DO I = 1, NB
          HIST(I) = 0.
          XDATA(I) = XMIN + (XMAX-XMIN)*FLOAT(I-1)/FLOAT(NB-1)
      END DO
C-----------------------------------------------------------------------
C      write(6,*) ' hist_2 peak    = ', peak
C      write(6,*) ' hist_2 spacing = ', delta
C      write(6,*) ' histog range   = ', xmin, xmax
C-----------------------------------------------------------------------
C     Accumulate the histogram.
C
      DO I = 1, NPTS
          DO J = 1, NB
              IF ( ABS(XDATA(J)-DATA(I)) .LT. DELTA ) THEN
                  HIST(J) = HIST(J) + 1
              END IF
          END DO
      END DO
C-----------------------------------------------------------------------
C     Find the peak
C
      IPEAK = 1
      DO I = 1, NB
          IF ( HIST(I) .GT. HIST(IPEAK) ) THEN
              IPEAK = I
          END IF
      END DO
C     WRITE(6,*) ' hist_2 peak found at ', IPEAK, XDATA(IPEAK)
      YPEAK = XDATA(IPEAK)
c
      if(hist(ipeak).gt.5)then
chummel Now determine exact peak position of highest peak
c
      if(index(title,'dD2C').eq.0)then ! otherwise we plot the dD2C histogram
      sy=0
      sxy=0
      sxxy=0
      sn=0
      sx=0
      sxx=0
      sxxx=0
      sxxxx=0
      chdelta=4.0*delta
      shpeak=ypeak
      if(index(title,'CHAN 3').ne.0)chdelta=chdelta*2.0
      do i=ipeak,nb
       if(abs(xdata(i)-ypeak).lt.chdelta
     $    .and.hist(i).gt.hist(ipeak)/3.)then
        lhist=log(hist(i))
        shxdata=xdata(i)-shpeak
        sy=sy+lhist
        sxy=sxy+shxdata*lhist
        sxxy=sxxy+shxdata*shxdata*lhist
        sn=sn+1
        sx=sx+shxdata
        sxx=sxx+shxdata*shxdata
        sxxx=sxxx+shxdata*shxdata*shxdata
        sxxxx=sxxxx+shxdata*shxdata*shxdata*shxdata
       else
        goto 10
       endif
      enddo
   10 do i=ipeak-1,1,-1
       if(abs(xdata(i)-ypeak).lt.chdelta
     $    .and.hist(i).gt.hist(ipeak)/3.)then
        lhist=log(hist(i))
        shxdata=xdata(i)-shpeak
        sy=sy+lhist
        sxy=sxy+shxdata*lhist
        sxxy=sxxy+shxdata*shxdata*lhist
        sn=sn+1
        sx=sx+shxdata
        sxx=sxx+shxdata*shxdata
        sxxx=sxxx+shxdata*shxdata*shxdata
        sxxxx=sxxxx+shxdata*shxdata*shxdata*shxdata
       else
        goto 20
       endif
      enddo
   20 matrixa(1,1)=dble(sn)
      matrixa(1,2)=dble(sx)
      matrixa(1,3)=dble(sxx)
      matrixa(2,1)=dble(sx)
      matrixa(2,2)=dble(sxx)
      matrixa(2,3)=dble(sxxx)
      matrixa(3,1)=dble(sxx)
      matrixa(3,2)=dble(sxxx)
      matrixa(3,3)=dble(sxxxx)
      matrixb(1,1)=dble(sy)
      matrixb(2,1)=dble(sxy)
      matrixb(3,1)=dble(sxxy)
      gjn=3
      gjnp=3
      gjm=1
      gjmp=1
      call gaussj(matrixa,gjn,gjnp,matrixb,gjm,gjmp)
      coeffb=-matrixb(3,1)
      coeffc=matrixb(2,1)/2.D0/coeffb
      coeffa=exp(matrixb(1,1)+coeffc*coeffc*coeffb)
      if(coeffb.lt.0)then
       write(12,'(a,a)')' *** Warning: Gaussfit failed!',title
      else
       ypeak=real(coeffc)+shpeak
      endif
c     write(12,*)' GAUSSJ: peak found at ',ypeak
      j=0
      do i=1,nb
       if(abs(xdata(i)-ypeak).lt.chdelta)then
        j=j+1
        xplot(j)=xdata(i)
        shxdata=xdata(i)-shpeak-real(coeffc)
        yplot(j)=real(coeffa)*exp(-real(coeffb)*shxdata*shxdata)
       endif
      enddo
      endif  ! endif for if(index(title,'dD2C').eq.0
      endif  ! endif for if(hist(ipeak).gt.5)
C-----------------------------------------------------------------------
C    Plot the histogram
C
      IF ( PLOT ) THEN
          YMIN = 0
          YMAX = HIST(IPEAK)
          CALL PGSCI ( 5 )
          CALL PGSCH  ( 2.0 )
          CALL PGRNGE ( XMIN, XMAX, XLO, XHI )
          CALL PGRNGE ( YMIN, YMAX, YLO, YHI )
          CALL PGENV  ( XLO, XHI, YLO, YHI, 0, 1 )
          call pgsch(1.5)
          CALL PGMTEXT ( 'T', 1.5, 0.5, 0.5, TITLE )
          CALL PGMTEXT ( 'B', 3.5, 0.5, 0.5, 'TWO COLOR DELAY' )
          CALL PGMTEXT ( 'L', 3.0, 0.5, 0.5, 'HISTOGRAM' )
          CALL PGSCH ( 1.0 )
          CALL PGSCI ( 3 )
          CALL PGBIN ( NB, XDATA, HIST, .TRUE. )
          CALL PGSCI ( 5 )
C
C     Plot the position of the fitted peak with a solid line,
C     the position of the one color peak with a dashed line.
C
          call pgsci(5)
          XBAR(1) = YPEAK
          XBAR(2) = YPEAK
          YBAR(1) = YLO
          YBAR(2) = YHI
      if(index(title,'dD2C').eq.0)then ! otherwise we plot the dD2c histogram
          CALL PGSLS ( 1 )
          CALL PGLINE( 2, XBAR, YBAR )
          call pgsci(11)
          call pgline(j,xplot,yplot)
          call pgsci(5)
      endif
          XBAR(1) = PEAK
          XBAR(2) = PEAK
          CALL PGSLS ( 2 )
          CALL PGLINE( 2, XBAR, YBAR )
          CALL PGSLS ( 1 )
      END IF
C-----------------------------------------------------------------------
      PEAK = YPEAK
      RETURN
      END