c
c filters.f -	computes optimum linear filters for extracting the total
c		charge (Qfilter) and arrival time (Tfilter) from a
c		digitized electronic time series, where the pulse shape
c		is known.  The filter length is chosen by the user.
c author: Richard Jones at uconn.edu
c version: Apr. 1, 2011
c
c notes by R. Jones:
c    1.	The filter is computed by generating a large number of repetitions
c	of the signal pulse, each with a random time offset and a pulse
c	magnitude P[i] drawn from a known pulse-height spectrum.  Let the
c	index i label the repetition and j represent the time series index.
c	The Qfilter is defined by the linear operator,
c                    __
c		Q  = \  g  a    + Q
c                i   /_  j  i,j    0
c		      j
c	One is looking for the choice of filter coefficients g[j] and Q[0]
c	that makes Q a minimum-variance estimator of P.  Form a chisquare
c       <(P-Q)(P-Q)> over the sample and solve for the g[j] and Q[0] that
c	minimize it.  The solution is
c                     __
c               g  =  \   H    C(P,a  )
c                j    /_   j,j'     j'
c                      j'
c	and                __
c		Q  = <P> - \  g <a >
c		 0         /_  j  j
c                           j
c
c       where C(P,a[j]) is the covariance between P and sample j in the
c	time series, and H[j,j'] is the inverse of the covariance matrix
c	between a[j] and a[j'].
c
c    2.	Random noise only enters through terms that are second order in a,
c	so only through H.  This is the most important thing to get from
c	the data themselves, in order to accurately describe the background
c	conditions under which the signal is to be extracted.  I get this
c       by analyzing data where I expect that no signal is found.  In
c	general, I expect it to peak at j=j' and taper off symmetrically
c	on either side as the two samples are further separated in time.
c	I measure this covariance matrix using real data, add it to the
c	one computed using Monte Carlo of the known signal pulse shape,
c	then take the inverse of the sum matrix to obtain H.
c
c    3. Pulse shape fluctuations are included through the simulation
c	methods implemented in the rates study.  These functions are
c	contained in the tresol package, which is included below.

      include 'tresol.f'

      subroutine filters
      implicit none
      include 'filters.inc'
      integer j,jj,jserial,jjserial
      real a0,a1,trise,tfall,R,Rnoise,dtmin
      a0 = 200. ! p.e.
      a1 = 50.  ! p.e.
      trise = 2.5e-9 ! s
      tfall = 5.0e-9 ! s
      R = 1e6  ! /s
      Rnoise = 4e6 ! /s
      dtmin = 10e-9 ! s
      correlated_sampling_noise_rms = 5.0e4 ! makes immune to baseline shifts
      uncorrelated_sampling_noise_rms = 5.0 ! random electronic sampling noise
      jstart = 0 ! samples
      jend = 15 ! samples
      if (jend-jstart+1.gt.max_filter_length) then
        print *, 'max_filter_length too small, please increase.'
        stop
      elseif ((jend-jstart+1)*(jend-jstart+2)/2
     +                       .gt.max_covariance_length) then
        print *, 'max_covariance_length too small, please increase.'
        stop
      endif
      sum1 = 0
      sumP = 0
      sumPt = 0
      jserial = 0
      jjserial = 0
      do j=0,jend-jstart
        jserial = jserial+1
        suma(jserial) = 0
        sumPa(jserial) = 0
        sumPta(jserial) = 0
        do jj=0,j
          jjserial = jjserial+1
          sumaa(jjserial) = 0
        enddo
      enddo
      write(6,*) '**********************************************'
      write(6,*) '* tfilter - computation of linear filters for'
      write(6,*) '*  extracting pulse amplitude and leading-edge'
      write(6,*) '*  time from a time series.'
      write(6,*) '*'
      write(6,*) '*  filter limits: ',jstart,', ',jend,' samples'
      write(6,*) '**********************************************'
      call tresol(a0,a1,trise,tfall,R,Rnoise,dtmin)
      end

      subroutine generate(nrep)
      implicit none
      integer nrep
      include 'filters.inc'
      include 'tresol.inc'
      logical hexist
      external hexist
      real ranlan
      external ranlan
      real tdecay
      real P,t,t0
      integer Npe
      integer ierr
      integer ip
      integer irep
      integer j,j1,j2
      integer jserial
      integer jj,jjserial
      real Pout,Ptout
      real rnd01(max_pixels)
      real unif01(2)
      do irep=1,nrep
        do j=1,max_samples
          sample(j) = 0
        enddo
        call add_noise
        t0 = -tmemory
        do while (t0 .lt. tracelen)
          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
            P = a0*(1+0.17*ranlan(unif01(2)))
            call RNPSSN(P,Npe,ierr)
            Npe = min(Npe,1000) ! truncate the Landau tail
            P = Npe ! renormalize to reduce outlier distortion
            if (verbose) then
              write(6,*) 'pulse amplitude ',Npe
            endif
            call ranlux(rnd01,Npe)
            tdecay = 0
            do ip=1,Npe
              tdecay = tdecay-log(rnd01(ip))*scintdecayt/(Npe-ip+1)
              call add_1photon(t0+tdecay)
            enddo
            j = t0/tstep
            if (j+jstart.gt.0 .and. j+jend.lt.max_samples) then
              t = t0-j*tstep
              sum1 = sum1+1
              sumP = sumP+P
              sumPt = sumPt+P*t
              Pout = g0P
              Ptout = g0Pt
              j1 = j+jstart
              j2 = j+jend
              jserial = 0
              jjserial = 0
              do j=j1,j2
                jserial = jserial+1
                suma(jserial) = suma(jserial)+sample(j)
                sumPa(jserial) = sumPa(jserial)+P*sample(j)
                sumPta(jserial) = sumPta(jserial)+P*t*sample(j)
                if (.not.(sample(j).gt.0 .or. sample(j).le.0)) then
                  print *, 'NaN in sample at position',j
                endif
                do jj=j1,j
                  jjserial = jjserial+1
                  sumaa(jjserial) = sumaa(jjserial)+sample(j)*sample(jj)
                enddo
                Pout = Pout+gP(jserial)*sample(j)
                Ptout = Ptout+gPt(jserial)*sample(j)
              enddo
              if (hexist(idPin)) then
                call hfill(idPin,P,0.,1.)
                call hfill(idPtin,P*t,0.,1.)
                call hfill(idtin,t,0.,1.)
                if (hexist(idPout)) then
                  call hfill(idPout,Pout,0.,1.)
                  call hfill(idPtout,Ptout,0.,1.)
                  call hfill(idtout,Ptout/Pout,0.,1.)
                  call hfill(idPinout,P,Pout,1.)
                  call hfill(idPtinout,P*t,Ptout,1.)
                  call hfill(idtinout,t,Ptout/Pout,1.)
                endif
              endif
            endif
          endif
        enddo
      enddo
      if (.not.hexist(idsum1)) then
        call hbook1(idsum1,'sum1',1,0.,1.,0.)
        call hbook1(idsumP,'sumP',1,0.,1.,0.)
        call hbook1(idsumPt,'sumPt',1,0.,1.,0.)
        call hbook1(idsumPa,'sumPa',jend-jstart+1,0.,jend-jstart+1.,0.)
        call hbook1(idsumPta,'sumPta',jend-jstart+1,0.,
     +                                               jend-jstart+1.,0.)
        call hbook1(idsuma,'suma',jend-jstart+1,0.,jend-jstart+1.,0.)
        call hbook2(idsumaa,'sumaa',jend-jstart+1,0.,jend-jstart+1.,
     +                              jend-jstart+1,0.,jend-jstart+1.,0.)
      endif
      call hfill(idsum1,0.,0.,sngl(sum1))
      call hfill(idsumP,0.,0.,sngl(sumP))
      call hfill(idsumPt,0.,0.,sngl(sumPt))
      jjserial = 0
      do j=0,jend-jstart
        call hfill(idsumPa,float(j),0.,sngl(sumPa(j+1)))
        call hfill(idsumPta,float(j),0.,sngl(sumPta(j+1)))
        call hfill(idsuma,float(j),0.,sngl(suma(j+1)))
        do jj=0,j-1
          jjserial = jjserial+1
          call hfill(idsumaa,float(j),float(jj),sngl(sumaa(jjserial)))
          call hfill(idsumaa,float(jj),float(j),sngl(sumaa(jjserial)))
        enddo
        jjserial = jjserial+1
        call hfill(idsumaa,float(j),float(j),sngl(sumaa(jjserial)))
      enddo
      end

      subroutine makefilters
      implicit none
      include 'filters.inc'
      logical hexist
      external hexist
      integer filter_length
      double precision S(max_filter_length,max_filter_length)
      double precision H(max_filter_length,max_filter_length)
      double precision U(max_filter_length,max_filter_length)
      double precision E(max_filter_length,2)
      double precision F(max_filter_length,2)
      double precision ev(max_filter_length)
      integer j,jj,jjserial
      integer ifail,info
      filter_length=jend-jstart+1
      jjserial = 0
      do j=1,filter_length
        E(j,1) = sumPa(j)-sumP*suma(j)/sum1
        E(j,2) = sumPta(j)-sumPt*suma(j)/sum1
        do jj=1,j
          jjserial = jjserial+1
          S(j,jj) = sumaa(jjserial)-suma(j)*suma(jj)/sum1
     +              + correlated_sampling_noise_rms*sum1
          S(jj,j) = S(j,jj)
          H(j,jj) = S(j,jj)
          H(jj,j) = S(j,jj)
        enddo
        S(j,j) = S(j,j)+uncorrelated_sampling_noise_rms*sum1
        H(j,j) = S(j,j)
      enddo
      call DSYEV('N', !compute eigenvalues only
     +           'L', !passing the lower-triangle of matrix
     +           filter_length, !order of matrix
     +           S,   !matrix to be analyzed (is overwritten)
     +           max_filter_length,
     +           ev,  !returns the eigenvalues in ascending order
     +           U,   !work space for function
     +           max_filter_length**2, !size of work space in U
     +           info)
      do j=1,filter_length
        if (ev(j).lt.0) then
          print *, 'Negative eigenvalue found in the sumaa matrix!'
          print *, 'Eigenvalues listed in ascending order are:'
          print *, (ev(jj),jj=1,filter_length)
          return
        elseif (.not.(ev(j).gt.0)) then
          print *, 'NaN eigenvalue found in the sumaa matrix!'
          print *, 'Eigenvalues listed in ascending order are:'
          print *, (ev(jj),jj=1,filter_length)
          return
        endif
      enddo
      call DSINV(filter_length,H,max_filter_length,ifail)
      if (ifail.ne.0) then
        print *, 'DSINV returns ifail=',ifail,', giving up.'
        return
      endif
      call DMMLT(filter_length,filter_length,filter_length,
     +           S(1,1),   S(1,2),   S(2,1),
     +           H(1,1),   H(1,2),   H(2,1),
     +           U(1,1),   U(1,2),   U(2,1),                  0.)
      call DMMLT(filter_length,filter_length,2,
     +           H(1,1),   H(1,2),   H(2,1),
     +           E(1,1),   E(1,2),   E(2,1),  
     +           F(1,1),   F(1,2),   F(2,1),                  0.)
      if (.not.hexist(idsumaaInv)) then
        call hbook2(idsumaaInv,'sumaaInv',
     +              filter_length,0.,float(filter_length),
     +              filter_length,0.,float(filter_length),0.)
        call hbook2(idsunit,'sunit',
     +              filter_length,0.,float(filter_length),
     +              filter_length,0.,float(filter_length),0.)
        call hbook1(idgP,'gP',filter_length,0.,float(filter_length),0.)
        call hbook1(idgPt,'gPt',filter_length,0.,
     +                                         float(filter_length),0.)
        call hbook1(idg0P,'g0P',1,0.,1.,0.)
        call hbook1(idg0Pt,'g0Pt',1,0.,1.,0.)
      else
        call hreset(idsumaaInv,'')
        call hreset(idsunit,'')
        call hreset(idgP,'')
        call hreset(idgPt,'')
        call hreset(idg0P,'')
        call hreset(idg0Pt,'')
      endif
      g0P = sumP/sum1
      g0Pt = sumPt/sum1
      jjserial = 0
      do j=1,filter_length
        gP(j) = F(j,1)
        gPt(j) = F(j,2)
        g0P = g0P-gP(j)*suma(j)/sum1
        g0Pt = g0Pt-gPt(j)*suma(j)/sum1
        call hfill(idgP,float(j-1),0.,gP(j))
        call hfill(idgPt,float(j-1),0.,gPt(j))
        do jj=1,filter_length
          if (jj.lt.j) then
            jjserial = jjserial + 1
            sumaaInv(jjserial) = H(j,jj)
            sunit(jjserial) = U(j,jj)
          endif
          call hfill(idsumaaInv,float(j-1),float(jj-1),sngl(H(j,jj)))
          call hfill(idsunit,float(j-1),float(jj-1),sngl(U(j,jj)))
        enddo
        jjserial = jjserial+1
        sumaaInv(jjserial) = H(j,j)
        sunit(jjserial) = U(j,j)
      enddo
      call hfill(idg0P,0.,0.,g0P)
      call hfill(idg0Pt,0.,0.,g0Pt)
      end

      subroutine checkfilters(nrep)
      implicit none
      integer nrep
      include 'filters.inc'
      logical hexist
      external hexist
      if (hexist(idPin)) then
        call hdelet(idPin)
        call hdelet(idPout)
        call hdelet(idPinout)
        call hdelet(idPtin)
        call hdelet(idPtout)
        call hdelet(idPtinout)
        call hdelet(idtin)
        call hdelet(idtout)
        call hdelet(idtinout)
      endif
      call hbook1(idPin,'Pin',100,0.,1000.,0)
      call hbook1(idPout,'Pout',100,0.,2000.,0)
      call hbook2(idPinout,'Pout vs Pin',100,0.,1000.,100,0.,1200.,0)
      call hbook1(idPtin,'Ptin',100,0.,4.e-6,0)
      call hbook1(idPtout,'Ptout',100,0.,5.e-6,0)
      call hbook2(idPtinout,'Ptout vs Ptin',100,0.,4.e-6,100,0.,5.e-6,0)
      call hbook1(idtin,'tin',100,0.,4e-9,0)
      call hbook1(idtout,'tout',100,0.,5.e-9,0)
      call hbook2(idtinout,'tout vs tin',100,0.,4.e-9,100,0.,5.e-9,0)
      call generate(nrep)
      end

      subroutine reload
      implicit none
      include 'filters.inc'
      logical hexist
      external hexist
      real hi,hij
      external hi,hij
      integer j,jj,jjserial
      integer filter_length
      filter_length = jend-jstart+1
      if (hexist(idgP)) then
        do j=1,filter_length
          gP(j) = hi(idgP,j)
          gPt(j) = hi(idgPt,j)
        enddo
        g0P = hi(idg0P,1)
        g0Pt = hi(idg0Pt,1)
      endif
      if (hexist(idsumP)) then
        sum1 = hi(idsum1,1)
        sumP = hi(idsumP,1)
        sumPt = hi(idsumPt,1)
        jjserial = 0
        do j=1,filter_length
          sumPa(j) = hi(idsumPa,j)
          sumPta(j) = hi(idsumPta,j)
          suma(j) = hi(idsuma,j)
          do jj=1,j
            jjserial = jjserial+1
            sumaa(jjserial) = hij(idsumaa,j,jj)
            sumaaInv(jjserial) = hij(idsumaaInv,j,jj)
            sunit(jjserial) = hij(idsunit,j,jj)
          enddo
        enddo
      endif
      end