pro pileup,table
;
common arrays,flx,enmean,ftot,e0a
common profile,p150,pofn
common probdist,probspectrum
;
;on_matherror,crash
on_ioerror,crash
;
;
;read in spectrum and energy values
dread,'bonomo2idl',flx,enmean ;read flux and energy bins calculated by bonomo
esfn=findgen(1000)+10 ;create 1 keV energy bins
sfn=interpol(flx,enmean,esfn) ;interpolate bonomo flux to linear bins in esfn
;
;normalize input spectrum to 1.
totsfn=total(sfn)
;put 1.e-30 where division by totsfn might be too small (< 1.e-30)
logtot=alog10(totsfn) & sw=where(alog10(sfn)-logtot lt -30)
first=sw(0)
if first eq -1 then sfn=sfn/totsfn else $
sfn=[sfn(0:first-1)/totsfn,fltarr(1000-first)+1.e-30]
;
;compute counting rate rt
rt=68.62*totsfn ;68.62 cm is the detector area, bins are 1 keV wide
rtfirst=rt ;save this original counting rate
print,'Counting rate (rt) = ',rt
;
; for table, compute a1 such that rt=400000 for starters.
if table eq 1 then begin a1new=400000/rt & rt=400000 & endif else a1new=1
print,'a1new=',a1new,' rt=',rt
;
;get r in units of tau=.6e-6 s, tc is coincidence time in tau's
tau=.6e-6 & tc=3 & r=rt*tau
;
;
e0a=findgen(1000)+.01 ; energy bins for spectrum starting at zero keV
f0a=[fltarr(10),sfn] ; put 10's in front of sfn array to make it match e0a.
f0b=f0a ; f0b will contain baseline corrections
;
; do baseline corrections if counting rate > than 100,000 counts/sec (rt>1.e5)
; results are in f0b.
;print,'before dobaseline'
if rt gt 1.e5 then dobaseline,r,sfn,esfn,f0a,f0b
;print,'safely finished dobaseline'
;
;
; Do convolution for 8 orders of pileup
norder=8
;
;(PERM$DATA:WMATRIX2) was originally obtained from the pulse shape using
;pro convolvemat in "convolve.pro". WMATRIX2 contains matrix elements from
;the voltage pulse shape distribution which connect the input pulse height
;distribution with the energy to be added due to pileup.
;
; compute Datlowe's w function
ww=fltarr(1000)
;print,'before opening w:matrix2'
openr,1,'perm$data:wmatrix2'
awm=assoc(1,fltarr(1000))
for i=0,999 do ww(i)=total(awm(i)*f0a(i:999)) ;f0a is always used, not f0b
;print,'finished ww calculation.'
close,1
;print,'closed wmatrix2'
;
; ww is the w function for the normalized input spectrum sfn
wwinvert=ww(999-indgen(1000))
;help,f0b
efc=fltarr(1000,norder)
efc(0,0)=f0b(0:999) ;f0b includes baseline correction if it was done, otherwise f0b=f0a
;print,'after defined efc'
;help,efc
for i=1,norder-1 do begin; &$
for j=0,999 do begin; &$
scale=f0b(j) ; scale factor so don't get underflow in multiplication
if j lt 20 then scale=1. ; guard against 0's at beginning of f0b
efc(j,i)=total(wwinvert(999-j:999)*(efc(0:j,i-1)/SCALE)) * SCALE; &$
endfor
endfor
;
; compute coefficients (pp) for combining the different orders of pileup
;print,'before pcoef'
pcoef,r,tc,pp,norder
;
; ftot is resulting spectrum
;print,'before flot=(efc#pp)'
ftot=(efc#pp) * totsfn * a1new ; combine pileup orders weighted by coeffs.,
; mult. by normalizing factor and scale to a1.
;
if table ne 1 then begin
dwrite,'idl2bonomo',ftot,e0a ; write resulting spectrum and energy bins
; to transfer them into BONOMO.
print,'pileup done.' & return & endif
;
;--------------------------------------------------------------------------
; the following is done only when writing the conversion factor table.
;
openu,15,'shelby$data:convfac.cvf',256
outrec=assoc(15,fltarr(64))
record = get_symbol('record')
energies=10.*exp((findgen(64)+1)*alog(100)/64.)
print,'record=',record
;
; save select elements of ftot for counting rate of 400000 c/s in fsave array
fsave=fltarr(64,10)
fsave(0,0)=ftot(energies)
;
a1save=fltarr(10) ; array for saving a1's to send back to bonomo
a1save(0)=a1new
;
; set up array of rates at which to compute convolved flux for table
rates=fltarr(9)
rates(0)=325000 & rates(1)=250000 & rates(2)=175000 & rates(3)=100000
rates(4)=50000 & rates(5)=25000 & rates(6)=12500 & rates(7)=6000
rates(8)=1000
;
; for table output, loop setting a1 for each counting rate in rates array,
; compute ftot, write ftot(energies)
;
efcbase=1 ; flag = 0/1 efc's have/don't have baseline correction
;
;print,'before large efc loop'
for i=0,8 do begin
a1new=rates(i)/rtfirst & r=rates(i)*tau &$
; print,'a1new=',a1new,' rt new = ',rates(i)
pcoef,r,tc,pp,norder ; compute coefficients pp
if r/tau gt 100000 then begin
dobaseline,r,sfn,esfn,f0a,f0b
efc=fltarr(1000,norder)
efc(0,0)=f0b ;f0b includes baseline correction
for ii=1,norder-1 do begin &$
for j=0,999 do begin &$
scale=f0b(j) ; scale factor so don't get underflow in multiplication
if j lt 20 then scale=1. ; guard against 0's at beginning of f0b
efc(j,ii)=total(wwinvert(999-j:999)*(efc(0:j,ii-1)/SCALE)) * SCALE
endfor
endfor
endif
;
;print,'before 2nd efc if in large loop'
if (efcbase eq 1) and (r/tau le 100000) then begin
efc=fltarr(1000,norder) ; recompute efc's without baseline corr.
efc(0,0)=f0a ;f0a doesn't include baseline correction
for ii=1,norder-1 do begin
for j=0,999 do begin
scale=f0b(j) ; scale factor so don't get underflow in multiplication
if j lt 20 then scale=1. ; guard against 0's at beginning of f0b
efc(j,ii)=total(wwinvert(999-j:999)*(efc(0:j,ii-1)/SCALE)) * SCALE &$
endfor
endfor
efcbase=0 ; set flag so won't recalculate efc's without baseline
endif
;
;print,'at second flot statement'
ftot=(efc#pp) * totsfn * a1new
; save select elements of ftot for writing in table
; print,ftot(energies)
fsave(0,i+1)=ftot(energies)
a1save(i+1)=a1new
endfor
;
;invert order of a1's and fsave's so that a1's in table will in be ascending
;order, and write inverted fsave's in table
a1save=[a1save(9-indgen(10)),fltarr(54)]
for i=0,9 do outrec(record+i+1)=fsave(0:63,9-i)
;
close,15 ; close table disk file
;
dwrite,'idl2bonomo',a1save,fix(energies); return a1's used to
; BONOMO to write in header of table
return
;
crash:
print,'Math or I/O error in pileup.'
print,'Error was:',strmessage(!err_string)
exit
;
end