;+
; test_stack.pro
;
; PURPOSE: 
; Script to produce plots of stacked counts vs time, and average stacked counts
; vs phase
; 
; EXAMPLE:
;  .RUN test_stack

; HISTORY
;   Sept 01, 2003  Richard Schwartz--object-oriented stacking
;   Sept 28, 2003  Ed Schmahl--test plots 
;-
;OBS_TIME='20-feb-02 '+['11:06','11:07'] & xyoffset=[900.,260.]
;OBS_TIME='20-apr-02 '+['23:26','23:27'] & xyoffset=[935.,-200.]
;OBS_TIME='20-apr-02 '+['23:26','23:27'] & xyoffset=[955.,-200.];shifted 20 asec
;OBS_TIME='01-jun-02 '+['03:56','03:57']  & xyoffset=[-423,-304]
;OBS_TIME='01-jun-02 '+['03:35','03:36']  & xyoffset=[-423,-304] ; only bkgd
OBS_TIME='17-apr-02 '+['00:40','00:41'] & xyoffset=[920,-235]
ob=obs_time(0) & strput,ob,'_',9 & strput,ob,'_',12 & ob=strcompress(ob,/rem)

;if not is_class(o,'hsi_image',/quiet) then $
    o=hsi_image(obs_time=obs_time,energy_band=[12.,25])
o->set,time_range=[0.,43.], xyoffset=xyoffset
o->set, use_auto_time=0
tbd = o->get(/time_bin_def)
o->set,time_bin_min=512, time_bin_def = tbd<64, use_phz_stack=0
o->set,energy_band=[12.,25]
o->plotman, flatfield=0

c0=o->getdata(class='hsi_calib_eventlist',this_det=0)
c1=o->getdata(class='hsi_calib_eventlist',this_det=1)
c2=o->getdata(class='hsi_calib_eventlist',this_det=2)
c3=o->getdata(class='hsi_calib_eventlist',this_det=3)
c4=o->getdata(class='hsi_calib_eventlist',this_det=4)
c5=o->getdata(class='hsi_calib_eventlist',this_det=5)
c6=o->getdata(class='hsi_calib_eventlist',this_det=6)
c7=o->getdata(class='hsi_calib_eventlist',this_det=7)
c8=o->getdata(class='hsi_calib_eventlist',this_det=8)


;
o->set,/use_phz_stacker
;
o->plotman

;
c0s=o->getdata(class='hsi_calib_eventlist',this_det=0)
c1s=o->getdata(class='hsi_calib_eventlist',this_det=1)
c2s=o->getdata(class='hsi_calib_eventlist',this_det=2)
c3s=o->getdata(class='hsi_calib_eventlist',this_det=3)
c4s=o->getdata(class='hsi_calib_eventlist',this_det=4)
c5s=o->getdata(class='hsi_calib_eventlist',this_det=5)
c6s=o->getdata(class='hsi_calib_eventlist',this_det=6)
c7s=o->getdata(class='hsi_calib_eventlist',this_det=7)
c8s=o->getdata(class='hsi_calib_eventlist',this_det=8)

window,0
!p.multi=[0,1,9] & !x.style=1
plot,c0.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 1'
plot,c1.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 2'
plot,c2.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 3'
plot,c3.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 4'
plot,c4.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 5'
plot,c5.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 6'
plot,c6.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 7'
plot,c7.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 8'
plot,c8.count,psym=10,tit='COUNTS VS TIME '+obs_time(0)+' SC 9'
t=tvrd()
write_png,ob+'cts_vs_t_SC1-9.png',255-t

window,1,title='STACKED COUNTS/S '+obs_time(0),xsize=1024,ysize=768
ti0='SC 1 STACKED COUNTS/S' & ti1='SC 2 STACKED COUNTS/S' & ti2='SC 3 STACKED COUNTS/S' & ti3='SC 4 STACKED COUNTS/S' & ti4='SC 5 STACKED COUNTS/S'
ti5='SC 6 STACKED COUNTS/S' & ti6='SC 7 STACKED COUNTS/S' & ti7='SC 8 STACKED COUNTS/S' & ti8='SC 9 STACKED COUNTS/S'

plot,c0s.count/c0s.livetime,psym=10,tit=ti0 
; ; FIND THE ROLL BINS USING UNIQ--SIMPLER THAN REVERSE HISTOGRAM 
u0=uniq(c0s.roll_angle) & n0=n_elements(u0)  ; FIND THE ROLL BINS
for j=0,n0-1 do oplot,[u0[j],u0[j]],[0,1000]

plot,c1s.count/c1s.livetime,psym=10,tit=ti1 
; ; FIND THE ROLL BINS USING UNIQ--SIMPLER THAN REVERSE HISTOGRAM 
u1=uniq(c1s.roll_angle) & n1=n_elements(u1)  ; FIND THE ROLL BINS
for j=0,n1-1 do oplot,[u1[j],u1[j]],[0,1000]

plot,c2s.count/c1s.livetime,psym=10,tit=ti12
; ; FIND THE ROLL BINS USING UNIQ--SIMPLER THAN REVERSE HISTOGRAM 
u2=uniq(c2s.roll_angle) & n2=n_elements(u2)  ; FIND THE ROLL BINS
for j=0,n2-1 do oplot,[u2[j],u2[j]],[0,1000]

plot,c3s.count/c3s.livetime,psym=10,tit=ti3 
; ; FIND THE ROLL BINS USING UNIQ--SIMPLER THAN REVERSE HISTOGRAM 
u3=uniq(c3s.roll_angle) & n3=n_elements(u3)  ; FIND THE ROLL BINS
for j=0,n3-1 do oplot,[u3[j],u3[j]],[0,1000]

plot,c4s.count/c4s.livetime,psym=10,tit=ti4 
; ; FIND THE ROLL BINS USING UNIQ--SIMPLER THAN REVERSE HISTOGRAM 
u4=uniq(c4s.roll_angle) & n4=n_elements(u4)  ; FIND THE ROLL BINS
for j=0,n4-1 do oplot,[u4[j],u4[j]],[0,1000]

plot,c5s.count/c5s.livetime,psym=10,tit=ti5 
; ; FIND THE ROLL BINS USING UNIQ--SIMPLER THAN REVERSE HISTOGRAM 
u5=uniq(c5s.roll_angle) & n5=n_elements(u5)  ; FIND THE ROLL BINS
for j=0,n5-1 do oplot,[u5[j],u5[j]],[0,1000]

plot,c6s.count/c6s.livetime,psym=10,tit=ti6 
u6=uniq(c6s.roll_angle) & n6=n_elements(u6)  ; FIND THE ROLL BINS
for j=0,n6-1 do oplot,[u6[j],u6[j]],[0,1000]


plot,c7s.count/c7s.livetime,psym=10,tit=ti7 
u7=uniq(c7s.roll_angle) & n7=n_elements(u7)  ; FIND THE ROLL BINS
for j=0,n7-1 do oplot,[u7[j],u7[j]],[0,1000]


plot,c8s.count/c8s.livetime,psym=10,tit=ti8 
u8=uniq(c8s.roll_angle) & n8=n_elements(u8)  ; FIND THE ROLL BINS
for j=0,n8-1 do oplot,[u8[j],u8[j]],[0,1000]

t=tvrd()
write_png,ob+'stacks_vs_t_SC1-9.png',255-t
; <><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><>

window,4,title='MEAN STACK PROFILE VS PHASE'
!P.MULTI=[0,3,3] & !x.style=1 & !p.charsize=2
phz=fltarr(8) & dphz=phz

c_phz=fltarr(8)
h=histogram(c0s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do phz(j)=mean((c0s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c0s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
for j=0,7 do  c_phz(j)=mean((c0s.count/c0s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=1',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=1',xtit='PHASE (RADIANS)'

c_phz=fltarr(8)
h=histogram(c1s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do phz(j)=mean((c1s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c1s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
for j=0,7 do  c_phz(j)=mean((c1s.count/c1s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=2',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=2',xtit='PHASE (RADIANS)'

c_phz=fltarr(8)
h=histogram(c2s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do phz(j)=mean((c2s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c2s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
for j=0,7 do  c_phz(j)=mean((c2s.count/c2s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=3',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=3',xtit='PHASE (RADIANS)'

c_phz=fltarr(8)
h=histogram(c3s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do phz(j)=mean((c3s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c3s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
for j=0,7 do  c_phz(j)=mean((c3s.count/c3s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=4',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=4',xtit='PHASE (RADIANS)'

c_phz=fltarr(8)
h=histogram(c4s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do phz(j)=mean((c4s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c4s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
for j=0,7 do  c_phz(j)=mean((c4s.count/c4s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=5',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=5',xtit='PHASE (RADIANS)'

c_phz=fltarr(8)
h=histogram(c5s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do phz(j)=mean((c5s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c5s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
for j=0,7 do  c_phz(j)=mean((c5s.count/c5s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=6',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=6',xtit='PHASE (RADIANS)'

h=histogram(c6s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
for j=0,7 do phz(j)=mean((c6s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c6s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do  c_phz(j)=mean((c6s.count/c6s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=7',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=7',xtit='PHASE (RADIANS)'

h=histogram(c7s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
for j=0,7 do phz(j)=mean((c7s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c7s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do  c_phz(j)=mean((c7s.count/c7s.livetime)[rev(rev(j):rev(j+1)-1)])
;plot,phz,c_phz,psym=6,tit='PHASE PROFILE, SC=8',xtit='PHASE (RADIANS)'
plot,c_phz,psym=10,tit='PHASE PROFILE, SC=8',xtit='PHASE (RADIANS)'

h=histogram(c8s.phase_map_ctr,bin=!pi/4.000001,rev=rev)
for j=0,7 do phz(j)=mean((c8s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
for j=0,7 do dphz(j)=stdev((c8s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])
if max(dphz) GT 0.01 then message,'IMPROPER PHASES'
if n_elements(h) NE 8 then message,'NOT 8 PHASE BINS!',/INFO
for j=0,7 do  c_phz(j)=mean((c8s.count/c8s.livetime)[rev(rev(j):rev(j+1)-1)])
plot,phz,c_phz,psym=10,tit='PHASE PROFILE, SC=9',xtit='PHASE (RADIANS)'

t=tvrd()
write_png,ob+'mean_stack_vs_phz_SC1-9.png',255-t

;plot,c_phz,psym=10,tit='PHASE PROFILE, SC=9',xtit='PHASE (RADIANS)'

; phases = 0.392699 at 0    8   16   24   29   33   41  49
; these indices are rev(9:16) or rev(rev(0):rev(0+1)-1)

; phases = 1.17810 at  1   9  17  25  30  34  42  50
; these indices are rev(9:24) or rev(rev(1):rev(1+1)-1)

;phases = 1.96350 are at rev(rev(2):rev(2+1)-1)

; etc

;c_phz(0)=mean((c8s.count/c8s.livetime)[rev(rev(0):rev(0+1)-1])
;c_phz(1)=mean((c8s.count/c8s.livetime)[rev(rev(1):rev(1+1)-1])
;c_phz(2)=mean((c8s.count/c8s.livetime)[rev(rev(2):rev(2+1)-1])

;phz=fltarr(8)
;for j=0,7 do dphz(j)=mean((c8s.phase_map_ctr)[rev(rev(j):rev(j+1)-1)])


end