subroutine smoother(NBIN1,NBIN2,NBIN3,NBIN4,NBIN5,
     +      inarray,outarray)
      implicit none
      integer NPASS_MAX,DMAX,NPASS_MIN
      parameter (NPASS_MAX = 6)
      parameter (NPASS_MIN = 1)
      parameter (DMAX = (NPASS_MAX + 1)**2 -1)
      integer NBIN1,NBIN2,NBIN3,NBIN4,NBIN5
      integer nreq
      parameter (nreq=20)
      logical o1,o2,o3,o4
      integer ij1,ij2,ij3,ij4,d1,d2,d3,d4,c
      integer i1,i2,i3,i4,i5,n,ntot,maxbin(5),max_n,mxb(5),mnb(5)
      integer i(5),j1,j2,j3,j4,j5,d,nlocal,like_lo,like_hi,j
      integer ipass,nbin_tot,id
      equivalence (i(1),i1),(i(2),i2),(i(3),i3),(i(4),i4),(i(5),i5)
      integer*4 outarray(NBIN1,NBIN2,NBIN3,NBIN4,NBIN5)
      integer*4 inarray(NBIN1,NBIN2,NBIN3,NBIN4,NBIN5)
      integer ncum(0:DMAX),nbcum(0:DMAX)
      real t,v
      real*8 v_tot

c       loop through the bins to count the total number of events
c       Also check for non zero values
c
        ntot = 0
        do i1 = 1,NBIN1
         do i2 = 1,NBIN2
          do i3 = 1,NBIN3
           do i4 = 1,NBIN4
            do i5 = 1,NBIN5
             n = inarray(i1,i2,i3,i4,i5)
              if (n.lt.0) then
               print *,'Inarray value less than zero',n,i1,i2,i3,i4,i5
               endif

               ntot=ntot+n

               enddo
               enddo
               enddo
               enddo
               enddo
               

c    begin smoothing, loop over bins
c    
               v_tot = 0.
               do i1 = 1,NBIN1
                 do i2 = 1,NBIN2
                   do i3 = 1,NBIN3
                    do i4 = 1,NBIN4 
                     do i5 = 1,NBIN5
                      call vzero(ncum(0),DMAX+1)
                      call vzero(nbcum(0),DMAX+1)

c
c      set up for a maximum of NPASS_MAX passes                      
c
              do ipass = 0,NPASS_MAX
               do j1 = max(-ipass,1-i1),min(ipass,NBIN1-i1)
                ij1 = i1 + j1
                d1 = j1**2
                o1 = iabs(j1).eq.ipass
                do j2 = max(-ipass,1-i2),min(ipass,NBIN2-i2)
                 d2 = d1 + j2**2
                 if (d2.gt.DMAX) cycle
                 ij2 = i2 + j2
                 o2 = o1.or.iabs(j2).eq.ipass
                 do j3 = max(-ipass,1-i3),min(ipass,NBIN3-i3)
                  d3 = d2 + j3**2
                  if (d3.gt.DMAX) cycle
                  ij3 = i3 + j3
                  o3 = o2.or.iabs(j3).eq.ipass
                  do j4 = max(-ipass,1-i4),min(ipass,NBIN4-i4)
                   d4 = d3 + j4**2
                   if (d4.gt.DMAX) cycle
                   ij4 = i4 + j4
                   o4 = o3.or.iabs(j4).eq.ipass
c              Check for valid combination, must have at
c          least one value of -ipass or +ipass
c         (otherwise,it was done in a previous pass)

           if (o4) then
         do j5= max(-ipass,1-i5),min(ipass,NBIN5-i5)
           d = d4 + j5**2
           if (d.gt.DMAX) cycle
           ncum(d) = ncum(d) +inarray(ij1,ij2,ij3,ij4,i5+j5)
           nbcum(d) = nbcum(d) + 1
           enddo
        else
                d = d4 + ipass**2
                if (d.le.DMAX) then
                   if (i5 - ipass.ge.1) then
              ncum(d) = ncum(d) +inarray(ij1,ij2,ij3,ij4,i5-ipass)
              nbcum(d) = nbcum(d) + 1
              endif
              if (i5 + ipass.le.NBIN5) then
              ncum(d) = ncum(d) +inarray(ij1,ij2,ij3,ij4,i5+ipass)
              nbcum(d) = nbcum(d) + 1
              endif
              endif
              endif
              enddo
              enddo
              enddo
              enddo

c
c       End of pass, check if enough events
c
          nlocal = 0
          nbin_tot = 0
          id = 0
          do while ((nlocal.lt.nreq).and.(id.lt.(ipass + 1)**2))
           nlocal = nlocal + ncum(id)
           nbin_tot = nbin_tot + nbcum(id)
           id = id + 1
           enddo

           if (nlocal.ge.nreq.and.ipass.ge.NPASS_MIN) exit
           enddo

c    Calsulate the soothed ln likelihood; events are normalized to one 
c    per bin, the result is stored as 2000*ln
           
           if (nlocal.eq.0) nlocal = 1
           v = float(nlocal)*NBIN1*NBIN2*NBIN3*NBIN4*NBIN5
     +      /float(nbin_tot)/float(ntot)

           
           v_tot = v_tot +v
           outarray(i1,i2,i3,i4,i5) =2000.*alog(v)
          
           
           enddo
           enddo
           enddo
           enddo
           enddo


  

c  Unfortunately, the output array is not perfectly normalized.
c
           c=2000.*dlog(NBIN1*NBIN2*NBIN3*NBIN4*NBIN5/v_tot)
           
             do i1=1,NBIN1
              do i2=1,NBIN2
               do i3=1,NBIN3
                do i4=1,NBIN4
                 do i5=1,NBIN5
          outarray(i1,i2,i3,i4,i5)=outarray(i1,i2,i3,i4,i5)+c
              enddo
              enddo
              enddo
              enddo
              enddo

           return
           end