forked from foxsi/calsoft
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmakehist.pro
executable file
·73 lines (56 loc) · 3.16 KB
/
makehist.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FUNCTION MAKEHIST, FILE, $
SUBTRACT_COMMON = SUBTRACT_COMMON, CMN_AVERAGE = CMN_AVERAGE, $
CMN_MEDIAN = CMN_MEDIAN, STOP = STOP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;if not keyword_set(savefile) then savefile = 'pedestal.sav'
; call preceding function to read in the data file(s).
;data = read_data_struct(file, subtract = subtract_common, stop = stop)
restore,file
n_evts = n_elements(data)
print, n_evts, ' total events'
; Split data up into ASICs. subtract common mode noise if desired.
i0 = where(rebin(data.packet_error,1,n_evts) lt 2.5)
cmn = intarr(n_evts,4)
hist = fltarr(1124, 64, 4, 2)
for as=0, 3 do begin
; if keyword_set(subtract_common) then cmn[*,as] = data[*].common_mode[as] + randomu(seed,n_evts*4+as) - 0.5
if keyword_set(subtract_common) then cmn[*,as] = data[*].common_mode[as]
if keyword_set(cmn_average) then cmn[*,as] = data[*].cmn_average[as]
if keyword_set(cmn_median) then cmn[*,as] = data[*].cmn_median[as] + randomu(seed,n_evts*4+as) - 0.5
endfor
; debugging
print, 'events for the spectra', n_elements(i0)
; Make histograms of data for each channel.
for as=0, 3 do begin
for k = 0, 63 do begin
if i0[0] gt -1 then gooddata = data[i0].data[as, k]-cmn[i0, as] else gooddata = 0
hist[*, k, as, 0] = findgen(1124)-100
hist[*, k, as, 1] = histogram([gooddata], min = -100, max = 1023)
endfor
endfor
; save data in current directory with chosen name, for easy recall.
;save, a0, a1, a2, a3, file = savefile
window, xsize = 800, ysize = 800
;x_axis = indgen(1124)-100
!p.multi = [0, 2, 2]
mm = max((rebin(hist, 1124, 1, 4, 2))(*,0,*,1))
plot, (rebin(hist, 1124, 1, 4, 2))(*, 0, 0, 0), (rebin(hist, 1124, 1, 4, 2))(*, 0, 0, 1), xrange = [-100, 800], yrange = [0.01, mm], $
xtitle = 'ASIC 0 raw counts', psym = 10,/ylog
plot, (rebin(hist, 1124, 1, 4, 2))(*, 0, 1, 0), (rebin(hist, 1124, 1, 4, 2))(*, 0, 1, 1), xrange = [-100, 800], yrange = [0.01, mm], $
xtitle = 'ASIC 1 raw counts', psym = 10,/ylog
plot, (rebin(hist, 1124, 1, 4, 2))(*, 0, 2, 0), (rebin(hist, 1124, 1, 4, 2))(*, 0, 2, 1), xrange = [-100, 800], yrange = [0.01, mm], $
xtitle = 'ASIC 2 raw counts', psym = 10,/ylog
plot, (rebin(hist, 1124, 1, 4, 2))(*, 0, 3, 0), (rebin(hist, 1124, 1, 4, 2))(*, 0, 3, 1), xrange = [-100, 800], yrange = [0.01, mm], $
xtitle = 'ASIC 3 raw counts', psym = 10,/ylog
; also store a copy of the save file in data storage directory
; if n_elements(file) eq 1 then save_name = file else save_name = file[0]
; if keyword_set(subtract_common) then save_name = save_name + '_sub'
; save_name = save_name + '.sav'
; save, a0, a1, a2, a3, file = save_name
return, hist
END