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 =
- \ g
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