; egs analysis again ; mwt/ajb 19/11/97 pro egsagain files = ['cs3015b.chn', 'cs3016b.chn', 'cs3017.chn', 'cs3018.chn' ,'cs3019.chn' ] gains= [ $ ; centroids of the pulser from manual fittings 1671.4, $ 1640.8, $ 1641.9, $ 1644.2, $ 1650.4 ] labels = [ $ 'Air only', $ '27.9 mm water', $ '52.4 mm water', $ '76.3 mm water', $ '87.9 mm water' ] depths = [0.0, 27.9, 52.4, 76.3, 87.9] ;depths = depths + 8.7 hard_roi = [1005.0, 1324.00] ;gains=REPLICATE(1.0, 5) n = n_elements(files) n_mca = 2048 bigdata = fltarr(n,n_mca) ;bigx = fltarr(n,2048 == n_mca) ;0, 1, 2, ...., 2047 ;0, 1, 2, ...., 2047 ;(n times) ;0, 1, 2, ...., 2047 little_x = FINDGEN(n_mca) big_x = FLTARR(n, n_mca) for i=0, n-1 do begin big_x(i,*) = gains(0)/gains(i) * little_x endfor ;plot, y ; plots aginst index for ifile = 0, n-1 do begin read_chn, files(ifile), data, info print, ifile, files(ifile), gains(0)/gains(ifile), ' ', info.sample, info.tlive if (ifile EQ 0) then biginfo=info else biginfo=[biginfo,info] bigdata(ifile,*)=data endfor ; now count the counts in a certain window ; the last nonzero count is 1344 n_windows=8 in_window=fltarr(n, n_windows) lastbin = 1360 ; by inspection for i_window = 0,n_windows-1 do begin firstbin = (1.0*i_window/n_windows)*lastbin ; print, firstbin in_window_0 = TOTAL(bigdata(0, firstbin:lastbin)) for ifile = 0,n-1 do begin lt_norm = float(biginfo(0).tlive)/float(biginfo(ifile).tlive) in_win = 1.0*TOTAL(bigdata(ifile, firstbin:lastbin)) in_win_norm = in_win * (float(biginfo(0).tlive)/float(biginfo(ifile).tlive)) / in_window_0 in_window(ifile, i_window) = in_win_norm endfor endfor log_in_win = alog10(in_window) degree = 1 poly_fits = fltarr(degree+1, n_windows) for i_window=0,n_windows-1 do begin poly_fits(*, i_window) = POLY_FIT(depths, log_in_win(*, i_window), degree) print, i_window, poly_fits(*, i_window) endfor ; now count: a) the entire spectrum; b) the hard region ; as done before for Jottings.xls ; i wonder if entire gets the pulser, causing a problem ; normalise to the reference data set == bigdata(0) count_entire = fltarr(n) count_almost_entire = fltarr(n) count_hard = fltarr(n) entire_0 = float(total(bigdata(0, *)) ) almost_0 = float(total(bigdata(0, 0:1360))) hard_0 = float(total(bigdata(0, hard_roi(0):hard_roi(1)))) for ifile = 0,n-1 do begin t_norm = float(biginfo(0).tlive)/float(biginfo(ifile).tlive) n_entire = float(total(bigdata(ifile, *))) n_almost = float(total(bigdata(ifile, 0:1360))) n_hard = float(total(bigdata(ifile, hard_roi(0):hard_roi(1)))) count_entire(ifile) = t_norm * n_entire $ / entire_0 count_almost_entire(ifile) = t_norm * n_almost $ / almost_0 count_hard(ifile) = t_norm * n_hard $ / hard_0 print, ifile, n_almost, n_hard, t_norm, count_almost_entire(ifile) , count_hard(ifile) endfor ; plot of everything IF !d.name NE 'PS' then begin window, /free, title='EGS again, spectra' w_spectra = !d.window window, /FREE, title='EGS again, in-Window counts' w_in_window = !d.window ; window, /FREE, title='EGS again == Jottings.xls' ; w_jottings = !d.window wset, w_spectra endif ; spectrum plot if !d.name EQ 'PS' then begin device, file='egs_spect.ps', /color help, /device endif n_bins = 256 ; 2048 for raw data smoothing = 5 plot_io, big_x(0,*), smooth( rebin(bigdata(0,*)>1, 1, n_bins) , smoothing) , $ xrange=[0,170], xstyle=1, $ ; 1.0*(n_bins/n_mca)*[0.,1600.], yrange=[ 40, max(rebin (bigdata(0,*), 1, n_bins ))], ystyle=1, $ psym=10, $ xtitle = 'channel', ytitle='count rate' xyouts, 0.6, 0.8 , /norm, labels(0); biginfo(0).sample for i=1,n-1 do begin coleur = !D.n_colors*0.5*(2.-float(i)/n) norml = (1.0*biginfo(0).tlive/biginfo(i).tlive) oplot, big_x(i,*), SMOOTH( rebin( bigdata(i,*), 1, n_bins )*norml, smoothing), $ psym=10, color=coleur xyouts, 0.6, 0.8 - (i* 0.05 ), labels(i), $ ; biginfo(i).sample /norm, color=coleur endfor for i_window = 1, n_windows-1 do begin i_color=!D.n_colors*0.5*(1.0-float(i_window)/n_windows) firstbin = (1.0*i_window/n_windows)*lastbin ;nextbin = (1.0*(i_window+1)/n_windows)*lastbin oplot, 1.0* n_bins / n_mca * [firstbin, firstbin] , [600, 400 ] ,$ color=i_color endfor oplot, [1.0, 1.0 ]* n_bins / n_mca * hard_roi(0), [500, 300] if !d.name ne 'PS' then begin bitmap = tvrd() write_gif, 'egs_spect.gif', bitmap endif if !d.name eq 'PS' then device, /close ; depths plot IF !d.name NE 'PS' then begin wset, w_in_window endif if !d.name eq 'PS' then device, file='egs_inwin.ps', /color xdepths = (findgen(55) / 54.) * depths(4) ; array of depths plot_io, depths, in_window(0:4, 0), yrange = [0.4,1.2], ystyle=1, $ xtitle='Depth of Water', ytitle='Log(in-window rates), relative to Air Only' for i_window = 1, n_windows-1 do begin i_color=!D.n_colors*0.5*(1.-float(i_window)/n_windows) oplot, depths, in_window(0:4, i_window), color=i_color fitted_absorptions = $ 10^(poly_fits(0, i_window) + xdepths * poly_fits(1, i_window) );+ xdepths^2 * poly_fits(2, i_window)) oplot, xdepths, fitted_absorptions, color=i_color, psym=3 endfor for i_file = 0, n-1 do begin coleur = !D.n_colors*0.5*(2.-float(i_file)/n) oplot, replicate(depths(i_file), n_windows), in_window(i_file, *), $ color=coleur, psym=4 endfor if !d.name ne 'PS' then begin bitmap = tvrd() write_gif, 'egs_inwin.gif', bitmap endif if !d.name eq 'PS' then device, /close ;IF !d.name NE 'PS' then begin ; wset, w_jottings ;endif ;plot, depths, count_entire(0:4), yrange = [0,1.2] ;oplot, depths, count_entire(0:4), psym=1 ;oplot, depths, count_almost_entire(0:4) ;oplot, depths, count_almost_entire(0:4), psym=1 ;oplot, depths, count_hard(0:4) ;oplot, depths, count_hard(0:4), psym=1 end 0.739233 0.808418 39 1060.80 1.00000 0.827003 0.790306 0.726962 0.667051 40 1088.00 1.00000 0.793736 0.773013 0.659374 0.661559 41 1115.20 1.00000 0.887691 0.734423 0.607891 0.642624 4