PRO cip_weights_and_errors
;+
; PRO cip_weights_and_errors
; -----------------------------------------
;  CALDB critical
; -----------------------------------------
; Procedure to read in the multi-grating created shell efficiencies
; and create :
;  - vignetting corrected shell effic.s
;  - the (Flight, 2pi) mirror-area-weighted efficiencies
;  - the error files for the shell areas and weighted effic.s.
; ----------- Summary of actions ---------
;  - read in needed data
;  - create/write plus-minus order averaged shell files that include
;       the shell vignetting factor
;  - create/write HEG MEG and HETG effective efficiency files
;  - create/write shell error files
;  - create/write HEG MEG and HETG effective efficiency error files
; ----------------------------------------
; CALDB files used as input:
;     hetg/fcp/multiout_981227_rdb/HETG_shell[N]_multiout.rdb
;     hetg/cip/hetgD1996-11-01basicN0001.rdb
;     hrma/cip/hrmaD1996-11-01effareaN0004.rdb
;-----------------------------------------
;-
; 1997-09-17 dd Modified to create HETG weighted output as well.
; (970917 version saved as cip_weights_and_errors.970917.pro)
; 1998-12-28 dd Overhauled for CALDB use; use new HRMA areas;
;  use new grating multiout files (improved polyimide edges);
;  use 5% for HRMA shell area error (instead of 10%).
; 1998-12-29 dd Use Herman's modified HRMA areas to avoid
;  2 keV discontinuity.
; 1999-01-09 dd Arg: remove unneeded "grating" prefix on files...
;   no other changes.
; 1999-07-14 dd Modified to create N0003 set of efficiencies:
;               changes:
;                       - modify the shell efficiencies
;                         zero and first orders
;                         with the function defined in 
;                         n0003_effic_corr.990714.pro
;                       - use HRMA 'N0004 flight values
;                       - change the error to a flat 10% for first order,
;                         30% for 3rd order
; ----------------------------------------


; --- preliminary: Get the data ------
; The efficiency data ('multiout.rdb files) are in
fcp_dir = !DDHETGCAL+'/fcp/' 
; The CIP directory is:
cip_dir = !DDHETGCAL+'/cip/'
; hrma cip dir
hrma_dir = !DDHRMACAL+'/cip/'

; --- version of the products ---
version_str = 'N0003'

; Read all the shell efficiencies in
multi_dir = 'multiout_981227_rdb/'
sh1 = rdb_read(fcp_dir+multi_dir+'HETG_shell1_multiout.rdb')
sh3 = rdb_read(fcp_dir+multi_dir+'HETG_shell3_multiout.rdb')
sh4 = rdb_read(fcp_dir+multi_dir+'HETG_shell4_multiout.rdb')
sh6 = rdb_read(fcp_dir+multi_dir+'HETG_shell6_multiout.rdb')

; Plot the plus and minus first orders...
plot_oi, sh1.energy, sh1.op1, XRANGE=[0.1,10.], $
	TITLE="HETG multiout Efficiencies"
oplot, sh1.energy, sh1.om1
oplot, sh3.energy, sh3.op1
oplot, sh3.energy, sh3.om1
oplot, sh4.energy, sh4.op1
oplot, sh4.energy, sh4.om1
oplot, sh6.energy, sh6.op1
oplot, sh6.energy, sh6.om1

; Read the HRMA areas
; rdb columns are: energy,ea_1,err_1, etc., ea_hrma, err_hrma
hrma_ea = rdb_read(hrma_dir+'hrmaD1996-11-01effareaN0004.rdb')
; put these into the arrays: es = energy, hrmaN = mirror pair N area:
es = hrma_ea.energy
hrma1 = hrma_ea.ea_1
hrma3 = hrma_ea.ea_3
hrma4 = hrma_ea.ea_4
hrma6 = hrma_ea.ea_6

; get the basic header string array that will describe these
; efficiency data
cip_effic_header, effic_hdr

; read the vignetting terms for each shell, HETGvign1 etc.
basic = rdb_read(cip_dir+'hetgD1996-11-01basicN0001.rdb')
shells = [1,3,4,6]
vigns = FLTARR(4)
for iv=0,n_elements(vigns)-1 do begin
  vign_name = 'HETGvign'+STRCOMPRESS(shells(iv),/REMOVE_ALL)
  this_vign = where(STRPOS(basic.name, vign_name) GE 0, nfound)
  if nfound EQ 1 then begin
    vigns(iv) = basic(this_vign(0)).value
  end else begin
    print, "* Couldn't find "+vign_name
    vigns(iv) = 0.935
  end
end
print, shells
print, vigns

;------ Write out the plus-minus averaged shell efficiencies ---
;       including vignetting

;------------ Write out the errors in these as well ------------

; Loop over the shells
shells = [1,3,4,6]
for is=0,3 do begin
  CASE is OF
  0: sh = sh1
  1: sh = sh3
  2: sh = sh4
  3: sh = sh6
  ENDCASE

  ; Go through all the tags (except the first one, energy)
  ; and put the average values into the output
  sh_out = sh
  ; Zero order tag index
  zind = n_tags(sh)/2
  ; for zero-order just include vignetting and the N0003 correction
  sh_out.(zind) = vigns(is) * sh.(zind) * $
	n0003_effic_corr(sh.energy, is, /ZERO)
  ; average the non-zero orders and include shell vignetting factor
  ; treat first order separately
  for io=1, 1 do begin
    sh_out.(zind+io) = vigns(is) * n0003_effic_corr(sh.energy, is) * $
	( sh.(zind+io) + sh.(zind-io) )/2.0
    sh_out.(zind-io) = sh_out.(zind+io)
  end
  for io=2, zind-1 do begin
    sh_out.(zind+io) = vigns(is) * ( sh.(zind+io) + sh.(zind-io) )/2.0
    sh_out.(zind-io) = sh_out.(zind+io)
  end

  ; Write the rdb file
  filename = 'hetgmp' + $
	STRCOMPRESS( STRING(shells(is)),/REMOVE_ALL ) $
	+'D1996-11-01greff'+ version_str +'.rdb'
  hdr_strings = STRARR(5)
  hdr_strings(0) = '#filename:	'+filename
  hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
  hdr_strings(2) = '#'
  hdr_strings(3) = '# HETG Shell '+ $
	STRCOMPRESS(STRING(shells(is)),/REMOVE_ALL) + $
	' Efficiency including vignetting'
  hdr_strings(4) = '#'
  ; append generic effic header comments
  hdr_strings =[ hdr_strings, effic_hdr ]
  ; and comments for marx use
  marx_hdr = [ $
  '# These lines are for John Davis and MARX use', $
  '# when they are commented/moved.', $
  '#', $
  '#Max_Order:	11', $
  '#Num_Energies:	5001', $
  '#................................' ]
  hdr_strings =[ hdr_strings, marx_hdr ]
  rdb_write, cip_dir+filename, sh_out, HEADER=hdr_strings

  ; Create the error values
  err_out = sh_out

  ; Zero order tag index
  zind = n_tags(err_out)/2
  ; 10% zero-order error everywhere
  err_array = 0.1 + 0.0*err_out.energy
  err_out.(zind) = err_array*sh_out.(zind)

  ; First order
  ; 5% for first-order above 2.5 keV,
  ; growing to 20% at and below 1 keV
;;  err_array = 0.05 + 0.0*err_out.energy
;;  below = where(err_out.energy LE 1.0)
;;  err_array(below) = 0.20
;;  between = where(err_out.energy GT 1.0 and $
;;		err_out.energy LT 2.5)
;;  err_array(between) = 0.05 + 0.15*(2.5 - err_out(between).energy)/(2.5-1.0)
;;  plot_oi, err_out.energy, 1.+err_array, LINESTYLE = 2, $
;;	XRANGE=[0.1,10.], XTITLE = 'Energy (keV)', $
;;	TITLE = 'First-order Error Band vs. Energy'
;;  plot, err_out.energy, 1.-err_array, LINESTYLE = 2
  ; --> fix it at 10%
  err_array = 0.10 + 0.0*err_out.energy
  err_out.(zind-1) = err_array*sh_out.(zind-1)
  err_out.(zind+1) = err_array*sh_out.(zind+1)
  ; 20% for 2 order
  for io= 2, 2 do begin
    err_array = 0.2 + 0.0*err_out.energy
    err_out.(zind+io) = err_array*sh_out.(zind+io)
    err_out.(zind-io) = err_array*sh_out.(zind-io)
  end
  ; 3 order error was 50% for N0002,
  ; reduce it to 30% given XRCF Phase H results on MEG.
  for io= 3, 3 do begin
    err_array = 0.3 + 0.0*err_out.energy
    err_out.(zind+io) = err_array*sh_out.(zind+io)
    err_out.(zind-io) = err_array*sh_out.(zind-io)
  end
  ; 90% for higher orders
  for io=4, zind-1 do begin
    err_array = 0.9 + 0.0*err_out.energy
    err_out.(zind+io) = err_array*sh_out.(zind+io)
    err_out.(zind-io) = err_array*sh_out.(zind-io)
  end

  ; Write the rdb file
  filename = 'hetgmp' + $
	STRCOMPRESS( STRING(shells(is)),/REMOVE_ALL ) $
	+'D1996-11-01grerr'+ version_str +'.rdb'
  hdr_strings = STRARR(5)
  hdr_strings(0) = '#filename:	'+filename
  hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
  hdr_strings(2) = '#'
  hdr_strings(3) = '# HETG Shell '+ $
	STRCOMPRESS(STRING(shells(is)),/REMOVE_ALL) + $
	' Efficiency Errors'
  hdr_strings(4) = '#'
  ; append generic effic header comments
  hdr_strings =[ hdr_strings, effic_hdr ]
  rdb_write, cip_dir+filename, err_out, HEADER=hdr_strings

; For each shell save the vignetting-included, averaged
; shell efficiency and the error for use below
  CASE is OF
  0: vsh1 = sh_out
  1: vsh3 = sh_out
  2: vsh4 = sh_out
  3: vsh6 = sh_out
  ENDCASE
  CASE is OF
  0: err1 = err_out
  1: err3 = err_out
  2: err4 = err_out
  3: err6 = err_out
  ENDCASE

end

;------ Make the mirror weighted values -------
; Algorithm:
;  - weight by mirror areas
;----------------------------------------------
; MEG
; Setup the output structure, same as shell efficiency
meg = vsh1
; Assume that shell 1 and 3 have same energy spacing (they do)
; and interpolate the mirror data to it
h1 = INTERPOL(hrma1, es, vsh1.energy )
h3 = INTERPOL(hrma3, es, vsh3.energy )
; Create the weighting coefficients
w1 = h1/(h1+h3)
w3 = h3/(h1+h3)
; Go through all the tags (except the first one, energy)
; and put the weighted values into the output
for it=1, n_tags(meg)-1 do begin
  meg.(it) = ( w1*vsh1.(it) + w3*vsh3.(it) )
end
; Write the rdb file
filename = 'hetgmegD1996-11-01greff'+ version_str +'.rdb'
hdr_strings = STRARR(5)
hdr_strings(0) = '#filename:	'+filename
hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
hdr_strings(2) = '#'
hdr_strings(3) = '#   MEG Mirror Weighted Efficiency'
hdr_strings(4) = '#'
; append generic effic header comments
hdr_strings =[ hdr_strings, effic_hdr ]
rdb_write, cip_dir+filename, meg, HEADER=hdr_strings

; HEG
; Setup the output structure, same as shell efficiency
heg = vsh4
; Assume that shell 4 and 6 have same energy spacing (they do)
; and interpolate the mirror data to it
h4 = INTERPOL(hrma4, es, vsh4.energy )
h6 = INTERPOL(hrma6, es, vsh6.energy )
; Create the weighting coefficients
w4 = h4/(h4+h6)
w6 = h6/(h4+h6)
; Go through all the tags (except the first one, energy)
; and put the weighted values into the output
for it=1, n_tags(heg)-1 do begin
  heg.(it) = ( w4*vsh4.(it) + w6*vsh6.(it) )
end
; Write the rdb file
filename = 'hetghegD1996-11-01greff'+ version_str +'.rdb'
hdr_strings = STRARR(5)
hdr_strings(0) = '#filename:	'+filename
hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
hdr_strings(2) = '#'
hdr_strings(3) = '#   HEG Mirror Weighted Efficiency'
hdr_strings(4) = '#'
; append generic effic header comments
hdr_strings =[ hdr_strings, effic_hdr ]
rdb_write, cip_dir+filename, heg, HEADER=hdr_strings

; HETG
; Setup the output structure, same as heg and meg
; Do this full HETG by combining the heg and meg
hetg = meg
; We already have the h1 - h6
; Create the weighting coefficients
;  wh, wm = weighting for heg and meg
wh = (h4 + h6)/(h1+h3+h4+h6)
wm = (h1 + h3)/(h1+h3+h4+h6)
; Go through all the tags (except the first one, energy)
; and put the weighted values into the output
for it=1, n_tags(hetg)-1 do begin
  hetg.(it) = ( wm*meg.(it) + wh*heg.(it) )
end
; Write the rdb file
filename = 'hetgD1996-11-01greff'+ version_str +'.rdb'
hdr_strings = STRARR(5)
hdr_strings(0) = '#filename:	'+filename
hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
hdr_strings(2) = '#'
hdr_strings(3) = '#   HETG Mirror Weighted Efficiency'
hdr_strings(4) = '#'
; append generic effic header comments
hdr_strings =[ hdr_strings, effic_hdr ]
rdb_write, cip_dir+filename, hetg, HEADER=hdr_strings

;------ Make the mirror weighted ERROR values -------
; Algorithm:
;   o shell errors are combined just as shell effic.s
;     'cause the errors are systematic mostly
;   o add 5% of the shell difference to account for
;     possible weighting errors (HRMA eff area errors)
;----------------------------------------------
hrma_shell_error = 0.05

; MEG
; Setup the output structure, same as shell efficiency
meg_err = meg
; Go through all the tags (except the first one, energy)
; calculate the error...
for it=1, n_tags(meg_err)-1 do begin
  meg_err.(it) = w1*err1.(it) + w3*err3.(it)
  wgt_err = ABS(hrma_shell_error*(vsh1.(it) - vsh3.(it)))
  meg_err.(it) = meg_err.(it) + wgt_err
end
; Write the rdb file
filename = 'hetgmegD1996-11-01grerr'+ version_str +'.rdb'
hdr_strings = STRARR(5)
hdr_strings(0) = '#filename:	'+filename
hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
hdr_strings(2) = '#'
hdr_strings(3) = '#   MEG Mirror Weighted Efficiency Error'
hdr_strings(4) = '#'
; append generic effic header comments
hdr_strings =[ hdr_strings, effic_hdr ]
rdb_write, cip_dir+filename, meg_err, HEADER=hdr_strings

; HEG
; Setup the output structure, same as shell efficiency
heg_err = heg
; Go through all the tags (except the first one, energy)
; calculate the error...
for it=1, n_tags(heg_err)-1 do begin
  heg_err.(it) = w4*err4.(it) + w6*err6.(it)
  wgt_err = ABS(hrma_shell_error*(vsh4.(it) - vsh6.(it)))
  heg_err.(it) = heg_err.(it) + wgt_err
end
; Write the rdb file
filename = 'hetghegD1996-11-01grerr'+ version_str +'.rdb'
hdr_strings = STRARR(5)
hdr_strings(0) = '#filename:	'+filename
hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
hdr_strings(2) = '#'
hdr_strings(3) = '#   HEG Mirror Weighted Efficiency Error'
hdr_strings(4) = '#'
; append generic effic header comments
hdr_strings =[ hdr_strings, effic_hdr ]
rdb_write, cip_dir+filename, heg_err, HEADER=hdr_strings

; HETG
; Setup the output structure, same as hetg
hetg_err = hetg
; Go through all the tags (except the first one, energy)
; calculate the error...
for it=1, n_tags(hetg_err)-1 do begin
  hetg_err.(it) = wm*meg_err.(it) + wh*heg_err.(it)
  wgt_err = ABS(hrma_shell_error*(meg.(it) - heg.(it)))
  hetg_err.(it) = hetg_err.(it) + wgt_err
end
; Write the rdb file
filename = 'hetgD1996-11-01grerr'+ version_str +'.rdb'
hdr_strings = STRARR(5)
hdr_strings(0) = '#filename:	'+filename
hdr_strings(1) = '#data produced by cip_weights_and_errors.pro '+SYSTIME()
hdr_strings(2) = '#'
hdr_strings(3) = '#   HETG Mirror Weighted Efficiency Error'
hdr_strings(4) = '#'
; append generic effic header comments
hdr_strings =[ hdr_strings, effic_hdr ]
rdb_write, cip_dir+filename, hetg_err, HEADER=hdr_strings


RETURN
END