PRO letg_cip_etc
; -----------------------------------------
;  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 LETG effective efficiency files
;  - create/write shell error files
;  - create/write LETG effective efficiency error files
; ----------------------------------------
; 1998-12-29 Created from letg_cip_etc.pro and the
;   new version of cip_weights_and_errors.981229.pro
;   Creating N0002 product, changes:
;    - use newly calculated grating multiout files (should
;        not differ for LETG from previous because gold consts
;        are the same.  These new multiout files do NOT
;        include ANY vignetting factors.)
;    - the LETGbasic file is updated with new values/interpretation
;        of the vignetting terms.
;    - use 5% for HRMA shell area error (instead of 10%).
;        and Herman's modified HRMA areas to avoid
;        2 keV discontinuity in weighting.
; 1999-01-09 dd Arg: remove unneeded "grating" prefix on files...
;   Change CIP directories for basic and output products to
;   be letg/cip .
; 1999-07-17 dd Modified for 'N0003 efficiency calculation from
;               hetg/fcp/multiout_990716_rdb/ files, from 990510 optical
;               constants.
; ----------------------------------------
; CALDB files used as input:
;     hetg/fcp/multiout_990716_rdb/LETG_shell[N]_multiout.rdb
;     letg/cip/letgD1996-11-01basicN0003.rdb
;     hrma/cip/hrmaD1996-11-01effareaN0004.rdb
;-----------------------------------------

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

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

; Read all the shell efficiencies in
sh1 = rdb_read(fcp_dir+'LETG_shell1_multiout.rdb')
sh3 = rdb_read(fcp_dir+'LETG_shell3_multiout.rdb')
sh4 = rdb_read(fcp_dir+'LETG_shell4_multiout.rdb')
sh6 = rdb_read(fcp_dir+'LETG_shell6_multiout.rdb')

; Plot the plus and minus first orders...
plot_oi, sh1.energy, sh1.op1, XRANGE=[0.1,10.], $
	TITLE="LETG 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 describes these
; efficiency data
cip_effic_header, effic_hdr

; read the vignetting terms for each shell, LETGvign1 etc.
basic = rdb_read(cip_dir+'letgD1996-11-01basicN0003.rdb')
shells = [1,3,4,6]
vigns = FLTARR(4)
for iv=0,n_elements(vigns)-1 do begin
  vign_name = 'LETGvign'+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
    print, " Using the average value of 0.8346"
    print, ""
    vigns(iv) = 0.8346
  end
end
print, shells

; Get the values of LETGvigncoarse and LETGvignfine
; These are the values of fraction of X-rays that are not
; absorbed by the structure.
vign_name = 'LETGvigncoarse'
this_vign = where(STRPOS(basic.name, vign_name) GE 0, nfound)
if nfound EQ 1 then begin
    LETGvigncoarse = basic(this_vign(0)).value
end else begin
    print, "* Couldn't find "+vign_name+'; using 0.90'
    print, ""
    LETGvigncoarse = 0.90
end
vign_name = 'LETGvignfine'
this_vign = where(STRPOS(basic.name, vign_name) GE 0, nfound)
if nfound EQ 1 then begin
    LETGvignfine = basic(this_vign(0)).value
end else begin
    print, "* Couldn't find "+vign_name+'; using 0.90'
    print, ""
    LETGvignfine = 0.90
end

vigns = vigns * LETGvigncoarse * LETGvignfine

print, 'Total vignetting factors for the shells are:'
print, vigns
print, ''

;------ 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
  sh_out.(zind) = vigns(is) * sh.(zind)
  ; average the non-zero orders and include shell vignetting factor
  for io=1, 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 = 'letgmp' + $
	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 letg_cip_etc.pro '+SYSTIME()
  hdr_strings(2) = '#'
  hdr_strings(3) = '# LETG 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:	25', $
  '#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)
  ; 10% for first-order everywhere
  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
  ; 50% for 3 order
  for io= 3, 3 do begin
    err_array = 0.5 + 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 = 'letgmp' + $
	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 letg_cip_etc.pro '+SYSTIME()
  hdr_strings(2) = '#'
  hdr_strings(3) = '# LETG 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
;----------------------------------------------
; 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
meg = vsh1
for it=1, n_tags(meg)-1 do begin
  meg.(it) = ( w1*vsh1.(it) + w3*vsh3.(it) )
end

; 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

; LETG
; Setup the output structure, same as heg and meg
; Do this full LETG by combining the heg and meg
LETG = 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(LETG)-1 do begin
  LETG.(it) = ( wm*meg.(it) + wh*heg.(it) )
end
; Write the rdb file
filename = 'letgD1996-11-01greff'+ version_str +'.rdb'
hdr_strings = STRARR(5)
hdr_strings(0) = '#filename:	'+filename
hdr_strings(1) = '#data produced by letg_cip_etc.pro '+SYSTIME()
hdr_strings(2) = '#'
hdr_strings(3) = '#   LETG Mirror Weighted Efficiency'
hdr_strings(4) = '#'
; append generic effic header comments
hdr_strings =[ hdr_strings, effic_hdr ]
rdb_write, cip_dir+filename, LETG, 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 10% 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
; 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

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


RETURN
END