*********************
* file to plot and fit SKAT data, 6/99
* kumac routine for reading and plotting SKAT data from scope at UH,
* original jgl 6/99, 
* much changed by jgl in 1/02

trace on
opt STA
opt fit
opt nfil
opt date

* clean the slate
hi/del *
ve/del *
pic/del *

* now read the data for tubes with water... use correct file
*ve/read tw,ws,wl sk_990608_170615.dat 
*ve/read tw,ws,wl sk_cal.dat '3(1x,F15.12)'
* ve/read tw,ws,wl sk_current.dat '3(1x,F15.12)'
 ve/read tw,ws,wl sk_103150.dat '3(1x,F15.12)'

* get number of points
sigma nwpts = nco(tw)
***>>>
sigma nwpts = 10000

* make vectors for binned data for overlaying
ve/cr wsb(10000) r
ve/cr wlb(10000) r

* convert time in sec to ns
sigma tw = tw*1.0e9
* invert from negative PMT pulses
sigma ws = -ws
sigma wl = -wl


* use last 1000 points for DC offset calculation and subtract
* sumv (vector, not vsum a scalar) makes running integral sum
 sigma wssum = sumv(ws)
 sigma wlsum = sumv(wl)
 sigma wsoff = (wssum(10000)-wssum(9000))/1000.0
 sigma wloff = (wlsum(10000)-wlsum(9000))/1000.0
 sigma ws = ws - wsoff
 sigma wl = wl - wloff
 sigma ws = abs(ws)
 sigma wl = abs(wl)

* set errors at arbitrary 10%
sigma wler = 0.1*wl
sigma wser = 0.1*ws

* get beginning and ending times for data with water
sigma twfirst = tw(1)
sigma dt = tw(2) - tw(1)
* add in extra interval so boundary at upper bin
sigma twlast = tw(nwpts) + dt

1dh 200 'amplitude vs time/ns, 3.54m tube' nwpts twfirst twlast 0.0
1dh 250 'amplitude vs time/ns, 7.41m tube' nwpts twfirst twlast 0.0

opt logy
zone 1 1

hi/put_vect/cont 200 ws
hi/put_vect/err  200 wser
hi/plot 200  

hi/put_vect/cont 250 wl
hi/put_vect/err  250 wler
hi/plot 250 


**********************************************************************
* now do Fourier transform the pulse time series
 
 sigma  twopip = 2.0*3.14159/10000.0
* make vector of frequency values (ie. 1/10000ns - 1000/10000ns )
 ve/cr freq(1000)   r

trace off

****>>temporarily not used, read from file 
 do il=1,1000
  /ve/input freq($eval([il])) $eval($eval([il])*twopip)
 enddo

ve/cre sftl(1000) r
ve/cre cftl(1000) r
ve/cre pftl(1000) r
ve/cre sfts(1000) r
ve/cre cfts(1000) r
ve/cre pfts(1000) r
ve/cre phi(1000) r


*>>> comment out temprarily and read in from file
do idm=1,1000
* get vector of phase angles for this frequency and time
 sigma phi=freq($eval([idm]))*tw
* cos and sin weighted vectors
 sigma sphil=wl*sin(phi)
 sigma cphil=wl*cos(phi)
 sigma sphis=ws*sin(phi)
 sigma cphis=ws*cos(phi)

* and make the sums for the transform, normalized
 sigma ssl=vsum(sphil)
 sigma ccl=vsum(cphil)
 sigma ppl=sqrt(ssl*ssl+ccl*ccl)
 sigma sss=vsum(sphis)
 sigma ccs=vsum(cphis)
 sigma pps=sqrt(sss*sss+ccs*ccs)

 ve/input sftl($eval([idm])) ssl  
 ve/input cftl($eval([idm])) ccl   
 ve/input pftl($eval([idm])) ppl   
 ve/input sfts($eval([idm])) sss  
 ve/input cfts($eval([idm])) ccs   
 ve/input pfts($eval([idm])) pps   

enddo

trace on

* read the precalculated FT instead of calculating here
*VECTOR/READ VLIST=freq,pfts,pftl FNAME=sk_current.pft

1dh 500 'fourier transform, short, water' 1000 0.001 0.1 0
hi/put_vect/cont 500 pfts
hi/pl 500

1dh 510 'fourier transform, long, water' 1000 0.001 0.1 0
hi/put_vect/cont 510 pftl
hi/pl 510

* now find spectral peaks, getting rid of DC part first
 sigma zero = 0.0
 do idm=1,50
  /ve/input pfts($eval([idm])) zero
  /ve/input pftl($eval([idm])) zero
 enddo

sigma ifunds=lvmax(pfts)
sigma ifundl=lvmax(pftl)

* gets the repetition period of pulses
sigma periods = 2.0*3.14159/freq(ifunds)
sigma periodl = 2.0*3.14159/freq(ifundl)
ve/pr periods
ve/pr ifunds
ve/pr periodl
ve/pr ifundl

* gets the index and time of the first (largest) pulse
sigma lpeaks = lvmax(ws)
sigma tpeaks = tw(lpeaks)
sigma lpeakl = lvmax(wl)
sigma tpeakl = tw(lpeakl)

ve/pr lpeaks
ve/pr tpeaks
ve/pr lpeakl
ve/pr tpeakl

* now how many potential peaks in window?
* the last time in the window is
sigma tmax = tw(nwpts)
sigma delts = tmax - tpeaks
sigma npeakss = 1 + int(delts/periods)
sigma deltl = tmax - tpeakl
sigma npeaksl = 1 + int(deltl/periodl)

ve/pr tmax
ve/pr delts
ve/pr npeakss
ve/pr deltl
ve/pr npeaksl

* now make the bins and integrate the histogram
sigma nbinss = npeakss - 1
sigma nbinsl = npeaksl - 1
ve/cr xs(101)   r
ve/cr ys(101)   r
ve/cr zs(101)   r
ve/cr xser(101) r
ve/cr yser(101) r
ve/cr xl(101)   r
ve/cr yl(101)   r
ve/cr zl(101)   r
ve/cr xler(101) r
ve/cr yler(101) r

* the time increment in input data
sigma dt = tw(2)-tw(1)
* the number of input bins between peaks
sigma lbins = int(periods/dt)
sigma lbinl = int(periodl/dt)
* take the start of peak integration as 1/8 period before peak
sigma lbegins = lpeaks - int(0.125*lbins)
sigma lbeginl = lpeakl - int(0.125*lbinl)
* the starting time of integration is then (ns)
sigma tbegins = tw(lbegins)
sigma tbeginl = tw(lbeginl)
* the last input bin will be then
sigma lends = lbegins + nbinss*lbins -1
sigma lendl = lbeginl + nbinsl*lbinl -1
* and the last time used will be (ns)
sigma tends = tw(lends)
sigma tendl = tw(lendl)

ve/pr dt
ve/pr lbins
ve/pr lbegins
ve/pr tbegins
ve/pr lends
ve/pr tends
ve/pr lbinl
ve/pr lbeginl
ve/pr tbeginl
ve/pr lendl
ve/pr tendl

* something seems not to work right about rebin... grrrr
* rebin 250 x y xer yer nbins lbegin lend N
* so do it the long way...

* first short tube
* take the peak bin width as half the peak spacing, rounded down
 sigma lbinws = int(lbins/2.0)
 do ibins=1,nbinss
  sigma lfirsts = lbegins + lbins*($eval([ibins])-1)
  sigma llasts = lfirsts + lbinws -1
  sigma mfirsts = lfirsts + lbinws
  sigma mlasts = llasts + lbinws 
  sigma yyy = 0.0
  sigma zzz = 0.0
  ve/pr lfirsts
  ve/pr llasts
  ve/pr mfirsts
  ve/pr mlasts
trace off
* integrate over the peak
  do jbin=lfirsts,llasts
    sigma yyy = yyy + ws($eval([jbin]))
  enddo
  sigma yyy = yyy/lbinws
  ve/input ys($eval([ibins])) yyy
  do jbin=lfirsts,llasts
    ve/input wsb($eval([jbin])) yyy
  enddo
* integrate the background estimate to be subtracted
  do mbin=mfirsts,mlasts
    sigma zzz = zzz + ws($eval([mbin]))
  enddo
  sigma zzz = zzz/lbinws
  ve/input zs($eval([ibins])) zzz
  do mbin=mfirsts,mlasts
    ve/input wsb($eval([mbin])) zzz
  enddo
 enddo
trace on
 sigma yser = 0.01*ys

* now long
* take the peak bin width as half the peak spacing, rounded down
 sigma lbinwl = int(lbinl/2.0)
trace off
 do ibinl=1,nbinsl
  sigma lfirstl = lbeginl + lbinl*($eval([ibinl])-1)
  sigma llastl = lfirstl + lbinwl -1
  sigma mfirstl = lfirstl + lbinwl
  sigma mlastl = llastl + lbinwl 
  sigma yyy = 0.0
  sigma zzz = 0.0
  ve/pr lfirstl
  ve/pr llastl
  ve/pr mfirstl
  ve/pr mlastl
* integrate over the peak
  do jbin=lfirstl,llastl
    sigma yyy = yyy + wl($eval([jbin]))
  enddo
  sigma yyy = yyy/lbinwl
  ve/input yl($eval([ibinl])) yyy
  do jbin=lfirstl,llastl
    ve/input wlb($eval([jbin])) yyy
  enddo
* integrate the background estimate to be subtracted
  do mbin=mfirstl,mlastl
    sigma zzz = zzz + wl($eval([mbin]))
  enddo
  sigma zzz = zzz/lbinwl
  ve/input zl($eval([ibinl])) zzz
  do mbin=mfirstl,mlastl
    ve/input wlb($eval([mbin])) zzz
  enddo
 enddo
 trace on
 sigma yler = 0.01*yl

* plot the binned data overlaying the raw data
 zone 1 1

 1dh 210 ' ' nwpts twfirst twlast 0.0
 hi/put_vect/cont 210 wsb
 hi/pl 200
 hi/pl 210 s  

 1dh 260 ' ' nwpts twfirst twlast 0.0
 hi/put_vect/cont 260 wlb
 hi/pl 250 
 hi/pl 260 s

 opt liny
 1dh 650 'integrated pulses, water, short' nbinss tbegins tends 0.0
 hi/put_vect/cont 650 ys
 hi/put_vect/err 650 yser

 1dh 652 'integrated pulses, water, long' nbinsl tbeginl tendl 0.0
 hi/put_vect/cont 652 yl
 hi/put_vect/err 652 yler

 sigma zser = 0.1*zs
 sigma zler = 0.1*zl

 1dh 655 'interpulse data, water, short' nbinss tbegins tends 0.0
 hi/put_vect/cont 655 zs
 hi/put_vect/err 655 zser
 hi/fit 655(2:7) e

 1dh 657 'interpulse data, water, long' nbinsl tbeginl tendl 0.0
 hi/put_vect/cont 657 zl
 hi/put_vect/err 657 zler
 hi/fit 657(2:7) e

* bin by bounces, short tube
 1dh 660 'rebin data water short' nbinss 0 nbinss 0
 hi/put_vect/cont 660 ys
 hi/put_vect/err  660 yser
 ve/cr pars(2) r
 ve/cr epars(2) r
 zone 1 1
 opt liny
 hi/fit 660(2:7) e E 2 pars ! ! ! epars
* one bounce transmission is exp{pars(2)}
 ve/pri pars
 ve/pri epars

* bin by bounces, long tube
 1dh 662 'rebin data water long ' nbinsl 0 nbinsl 0
 hi/put_vect/cont 662 yl
 hi/put_vect/err  662 yler
 ve/cr parl(2) r
 ve/cr eparl(2) r
 zone 1 1
 opt liny
 hi/fit 662(2:7) e E 2 parl ! ! ! eparl
 ve/pri parl
 ve/pri eparl

* "calibration" or empty runs give (from data on 6/13/01):
 sigma parcs = -0.239532
 sigma parcl = -0.171707
* make errors arbitrarily smaller for text, by 10x
 sigma eparcs = 0.00244614
 sigma eparcl = 0.00242131

* length of short tube, m
 sigma alengs = 3.54
 sigma alengl = 7.41
 sigma dell = 2.0*(alengl - alengs)

* assuming optical losses/bounce same for each, ditto water, m
 sigma attenws = 2.0*alengs/(parcs-pars(2))
 ve/pri attenws
 sigma atenwse = (attenws**2/(2.0*alengs))*sqrt(eparcs**2 + epars(2)**2)
 ve/pri atenwse

 sigma attenwl = 2.0*alengl/(parcl-parl(2))
 ve/pri attenwl
 sigma atenwle = (attenwl**2/(2.0*alengl))*sqrt(eparcl**2 + eparl(2)**2)
 ve/pri atenwle

 sigma attenw = dell/(pars(2)-parl(2)-(parcs-parcl))
 ve/pri attenw
 sigma attenwe = (attenw**2/dell)*sqrt(eparcl**2+eparcs**2+epars(2)**2+eparl(2))
 ve/pri attenwe

 trace off