PRO hsi_az_play
;
; This list of IDL commands is used to crudely generate an azimuthal
; histogram of the HETG/LETG alignment data sets (like ring focus).
; Prints to terminal can be catted for a hardcopy.

phase_1_data_dir = '/nfs/atum/d3/dd/dd_d6'

; Read raw HSI events
; Comment in the desired file...
; and its manually determined center:
; and its defocus value:

; old ring focus to test the s/w
;;date_code = '961221'
;;run_id = '106777'
;;cent_Y = -340.0
;;cent_Z = 70.
;;defocus = 65.3
;;hsi_file = '/data/asc13/data/xrcfa/data/'+date_code+'/hsi'+run_id+'i0.fits'

; C-K ring focus at 65.20 mm, Use to characterize HSI response...
; C-K
date_code = '970115'
run_id = '108944'
cent_Y = -360.
cent_Z = 20.0
defocus = 60.5795  ; !?! don't know why this is so small...
hsi_file = phase_1_data_dir+'/'+date_code+'/hsi'+run_id+'i0.fits'

; HETG-reference, new ring focus at 65.20 mm
; Al-K
;;date_code = '961223'
;;run_id = '106856'
;;cent_Y = -290.0
;;cent_Z = 20.0
;;defocus = 60.5795  ; !?! don't know why this is so small...
;;hsi_file = phase_1_data_dir+'/'+date_code+'/hsi'+run_id+'i0.fits'

; HETG alignment at 65.20 mm
;;date_code = '961223'
;;run_id = '106860'
;;cent_Y = -290.0
;;cent_Z = 20.0
;;defocus = 60.5795  ; !?! don't know why this is so small...
;;hsi_file = '/spectra/d6/'+date_code+'/hsi'+run_id+'i0.fits'

; LETG-reference ring focus at 40 mm
;;date_code = '961223'
;;run_id = '106857'
;;cent_Y = -200.0
;;cent_Z = 20.
;;defocus = 60.5795*40./65.2 
;;hsi_file = '/d1/data/'+date_code+'/hsi'+run_id+'i0.fits'

; LETG-inserted Alignment at 40 mm
;;date_code = '961223'
;;run_id = '106858'
;;cent_Y = -200.0
;;cent_Z = 20.
;;defocus = 60.5795*40./65.2 
;;hsi_file = '/d1/data/'+date_code+'/hsi'+run_id+'i0.fits'

; This routine reads in the hsi RAW data and corrects for the
; HSI pixel scale and rotation (/FAC) to get events in microns
; in Facility coordinates.
hsi_dat = xrcf_hsi_read(hsi_file,/FAC,/PBLUR)

print, ' '
print, ' Data from '+date_code+', run id : '+run_id
print, ''
print, ' Center manually determined to be:'
print, '   dY = ',cent_Y, ',    dZ = ',cent_Z, ' (microns)

; Make the centered coords, in mm
Y_wrtc = 1.E-3 * (FLOAT(hsi_dat.Y) - cent_Y)
Z_wrtc = 1.E-3 * (FLOAT(hsi_dat.Z) - cent_Z)

; Plot the events including "cross hairs" to help with
; the manual centering...
plot, Y_wrtc(0:50000), Z_wrtc(0:50000), $
	XRANGE=[-5.,5.],YRANGE=[-5.,5.],PSYM=3, $
	XTITLE='Facility Y (mm)', YTITLE='Facility Z (mm)', $
	TITLE = 'Image from '+date_code+'/hsi'+run_id+'i0.fits'
oplot, [-5.,5.],[0.,0.]
oplot, [0.,0.],[-5.,5.]


; Get the radius for the events
radius = SQRT(Y_wrtc^2 + Z_wrtc^2)

; and the angle of the events
az_angle = ATAN(Z_wrtc,Y_wrtc)/!DTOR

; Make radial distribution
rad_bin = 0.01
lin_hist, radius, rad_bin, rad_val, rad_dist

;;rad_dist = histogram(radius, BINSIZE = rad_bin)
;;rad_val = rad_bin*(0.5+INDGEN(n_elements(rad_dist)))

plot, rad_val, rad_dist, $
	PSYM=10,XRANGE=[0.,5.], TITLE = 'Radial Distribution', $
	XTITLE = 'Radius (mm)', YTITLE='Number of Events'

; Set the rings radii scaled by the defocus values (LEG will be 40 mm)
; and plot these raddi on the above plot
shell_1 = [3.8,4.7] * defocus/65.3
for ir=0,1 do oplot, shell_1(ir)*[1.,1.], [0.,1.E6], LINESTYLE=1
shell_3 = [3.15,3.6] * defocus/65.3
for ir=0,1 do oplot, shell_3(ir)*[1.,1.], [0.,1.E6], LINESTYLE=2
shell_4 = [2.6,3.13] * defocus/65.3
for ir=0,1 do oplot, shell_4(ir)*[1.,1.], [0.,1.E6], LINESTYLE=1
shell_6 = [1.90,2.45] * defocus/65.3
for ir=0,1 do oplot, shell_6(ir)*[1.,1.], [0.,1.E6], LINESTYLE=2

; Make some plots...
pre_print_portrait
!p.multi = [0,1,2]

; Loop over the four shells making plots and rdb files
; This will have to be done when running as a procedure...
for is=0,3 do begin

  CASE is OF
    0: begin
      shell_range = shell_1
      shell_num = 1
    end
    1: begin
      shell_range = shell_3
      shell_num = 3
    end
    2: begin
      shell_range = shell_4
      shell_num = 4
    end
    3: begin
      shell_range = shell_6
      shell_num = 6
    end
  ELSE:
  ENDCASE

  ; Select a shell
  print, ''
  print, ' Select SHELL '+STRCOMPRESS(shell_num)+' events'
  print, '     bin       angle        counts'
  shell_events = WHERE(radius GT shell_range(0) AND $
			radius LT shell_range(1))

  ; See just those events...
  nshow = 20000
  plot, Y_wrtc(shell_events(1:nshow)), Z_wrtc(shell_events(1:nshow)), $
	XRANGE=[-5.,5.],YRANGE=[-5.,5.], XSTY=1, YSTY=1, PSYM=3, $
	XTITLE='Facility Y (mm)', YTITLE='Facility Z (mm)', $
	TITLE = 'Image from '+date_code+'/hsi'+run_id+'i0.fits:' + $
		' Shell '+STRCOMPRESS(shell_num)+' events'
  oplot, [-5.,5.],[0.,0.]
  oplot, [0.,0.],[-5.,5.]
  xyouts, -4.2,4.2, 'Total Events = '+STRCOMPRESS(n_elements(shell_events))+ $
	' (only '+STRCOMPRESS(nshow)+' shown)', $
	CHARSIZE=1.5

  ; Make an azimuthial distribution...
  az_bin = 2.0

  lin_hist, az_angle(shell_events), az_bin, az_bins, az_dist

  plot, az_bins,az_dist, $
	PSYM=10,XRANGE=[-360.,360.], XSTY=1, $
	TITLE = 'Azimuthial Distribution SHELL '+STRCOMPRESS(shell_num), $
	XTITLE = 'Angle from +Y axis (degree)', YTITLE='Number of Events'
  oplot, 360.+az_bins,az_dist, PSYM=10, LINESTYLE=2
  oplot, az_bins - 360.0,az_dist, PSYM=10, LINESTYLE=2

  ; Put the counts data into a structure and write it to an rdb file...
  out_struct = {bin:0.0, count:0.0, err:0.0}
  out_data = REPLICATE(out_struct, n_elements(az_bins))
  out_data.bin = az_bins
  out_data.count = FLOAT(az_dist)
  out_data.err = SQRT(FLOAT(az_dist))

  out_file = 'hsi'+run_id+'_az_hist_shell'+ $
		STRCOMPRESS(shell_num,/REMOVE)+'.rdb'
  hdr = ['# created by xrcf/hsi_az_play.pro, '+SYSTIME(), $
	'#', $
	'#']
  rdb_write, out_file, out_data, HEADER=hdr

end ; of loop over shells

device,/close
SPAWN, 'mv idl.ps hsi_az.'+run_id+'.ps'
set_plot, 'X'

; Put these files in a cmp directory...
SPAWN, 'cp hsi'+run_id+'_az_hist_shell*.rdb '+!DDXRCFCAL+'/cmp/'

RETURN
END