PRO toga_full_sim, row_todo, WEB = web, MILLION = million, $
	SPECT_ONLY=spect_only
; Simulation of EIPS-TMA-TOGA-Detector
; using csim
; 8/23/96 dd
; IDL procedure to control csim operation (sorry John!)
; Should be straight forward to translate this to perl, etc.

; 9/29-30 dd Modified to save intermediate files, 
; output IDL save files, and create directories and plots
; for the web.
; 10/1 Modify to do HIREF-W, HIREF-C and DCM as monochromator
; option after an EIPS source
; 10/8 dd Add ACIS 2C OBF as a filter material
; 10/25-26/96 dd Modified to unbundle the source spec generator (xss_sim)
; and works from .trw files too.
;
; Run simulation by  either:
;   editing Manual parameters below (and in xss_sim) and
;   typing  toga_full_sim  to IDL prompt
;  OR
;   set up cmdb_file_path and cmdb_file and
;   type: toga_full_sim, row_number, [/web]
; to simulate the test measurement in row.
; To simulate all measurements in a CMDB type:
;   for ir=0,n_elements(cmdb_lines)-1 do toga_full_sim, ir, /web
; The /web option makes .gif files and puts
; them on the web.


@labx_common
; - - - - - - - Comment out these line if no CMDB driver used - - - - 
@cmdb_common

; Load the CMDB if not already
if n_elements(cmdb_lines) LE 0 then begin
  cmdb_file_path = '/nfs/wiwaxia/h1/dd/idl/cmdb/phase3/'
  cmdb_load, 'phase3_961025.trw'
end

if cmdb_using_cal then begin
@cmdb_cal_integers
end else begin
@cmdb_trw_integers
end

; - - - - - - - 
; Edit parameters below to change them...
;.................................................................

; - - - - - Main parameters to change - - - - - 
if n_elements(row_todo) LE 0 then begin
  ; ** Manual entry - if no ria_todo array given - modify items below
  ;    in this file...
  spect_file = 'use_Dans_labxray_routines' ; 'use_Dans_labxray_routines' OR 'use_Penning'

  ; not using above spect_file string but rather:
  ; Create spectrum file using Dan's X-GEF EIPS model...
  lx_tname = 'Al' ; Select from C,O,Mg,MgO,Al,Si,SiO,Ti,Cr,Fe,Cu,Zr,Mo,W,
  lx_sV = 12.0 ; kV
  lx_sI = 1.0 ; mA

  fmat1 = 'Al'
  fthk1 = 0.0
  fmat2 = 'NONE'
  fthk2 = 0.0

  mono_energy = -1.0  ; Monochromator used if greater than 0.
  mono_res = 300.0    ; resolving power of monochromator
  mono_highorders = 0  ; 0 = no high orders, 1 = HIREFS approx

  det_name = 'HSI'     ; '2C', 'HSI', 'FPC', 'SSD'

  det_defocus = 0.  ; mm from XRCF focus, + is towards TMA/TOGA
  det_dy = 0.0       ; offset of detector center from optical axis in Y

  ; For '2C "0" aim point is at center of third quadrant
  if det_name EQ '2C' then det_dy = det_dy - 0.024*128.

  det_dz = 0.0  ; " in Z
  det_z_size = 18.  ; 2C rows OR detector diameter

  toga_inout = 'TOGA'   ; 'TMA' for no toga, 'TOGA' for toga inserted
  shutters = [0,0,0,0] ; megTop, legNorth, hegBottom, legSouth
                     ; 0 = shutter is out of beam, "open"
                     ; 1 = shutter blocks beam, "closed"

  ; - - - - - - - - - -  Less frequently changed values - - - -
  ; Name for this simulation run... made up from above things...
  run_name = toga_inout + '_' + lx_tname + $
	STRCOMPRESS(STRING(lx_sV,FORMAT='(F4.1)'),/REMOVE) + $
	'kV'+'_'+det_name
  numraystring = '20000' ; rays per quadrant, 10000 per trace iteration
  output_root_dir = '/spectra/d1/toga_sim/'+run_name

end else begin
  ; - - - This whole ELSE is for CMDB usage: it fills the above
  ;        variables based on CMDB line
  ; ** CMDB row index is specified, fill values 
  ; below using CMDB data

  spect_file = 'use_Dans_labxray_routines'
  ; Create spectrum file using Dan's X-GEF EIPS model...
  ; All the sources except Penning start with a EIPS source
  ; setup EIPS source
  CASE STRUPCASE(cmdb_fields(row_todo,c_source)) OF
    'EIPS': begin
      ; Set the EIPS target
      pieces = STR_SEP(cmdb_fields(row_todo,c_target),'-')
      if n_elements(pieces) LE 1 then begin
        print, ' * no valid target value'
        RETURN  ; can't simulate this
      end
      ct = pieces(0)  ; extract the target material
      ; Check if I have this target
      t_list = ['Be','C','O','Mg','MgO','Al', 'Si', $
			'Ti','Cr','Fe','Cu','Mo','Zr','W']
      ; if it is not in the list, return
      if MAX(ct EQ t_list) NE 1 then begin
        print, ' * target not in target list'
        RETURN
      end
      ; if it is 'O' or 'Si' use 'SiO'
      if ct EQ 'O' OR ct EQ 'Si' then ct = 'SiO'
      ; set EIPS current to 1.0 mA
      lx_sI = 1.0
      ; no monochromator
      mono_energy = -1.0
      mono_res = -1.0
      mono_highorders = 0
    end

    'HIREF-W': begin
      ct = 'W'
      lx_sI = 50.0 ; mA
      mono_energy = 0.0
      READS, cmdb_fields(row_todo,c_energy), mono_energy
      mono_res = 300.0
      mono_highorders = 1
    end

    'HIREF-C': begin
      ct = 'C'
      lx_sI = 50.0 ; mA
      mono_energy = 0.0
      READS, cmdb_fields(row_todo,c_energy), mono_energy
      mono_res = 300.0
      mono_highorders = 1
    end

    'DCM': begin
      ct = 'W'
      lx_sI = 50.0 ; mA
      mono_energy = 0.0
      READS, cmdb_fields(row_todo,c_energy), mono_energy
      mono_res = 300.0
      mono_highorders = 1
    end

    'PENNING': begin
      spect_file = 'use_Penning'
    end

    ELSE: begin
      print, ' *** Unrecognized Source: '+cmdb_fields(row_todo,c_source)
      RETURN
    end

  ENDCASE

  if spect_file EQ 'use_Dans_labxray_routines' then begin
    ; EIPS getting source parameters
    lx_tname = ct ;
    ; get voltage
    voltage  = 0.0
    READS, cmdb_fields(row_todo, c_voltage), voltage
    lx_sV =  voltage; kV
  end

  fmat1 = cmdb_fields(row_todo,c_filter_mat_1)
  fthk1 = 0.0
  ; dont need thickness for no filter or OBF
  if NOT(fmat1 EQ 'NONE' OR fmat1 EQ 'N/A' OR fmat1 EQ '' OR $
					fmat1 EQ 'OBF') then begin  
    READS, cmdb_fields(row_todo, c_filter_thick_1), fthk1
  end
  fmat2 = cmdb_fields(row_todo,c_filter_mat_2)
  fthk2 = 0.0
  ; dont need thickness for no filter or OBF
  if NOT(fmat2 EQ 'NONE' OR fmat2 EQ 'N/A' OR fmat2 EQ '' OR $
					fmat2 EQ 'OBF') then begin  
    READS, cmdb_fields(row_todo, c_filter_thick_2), fthk2
  end

  det_name = '2C' ; default
  if STRPOS(cmdb_fields(row_todo, c_detector),'HSI') GE 0 then det_name = 'HSI'
  if STRPOS(cmdb_fields(row_todo, c_detector),'FPC') GE 0 then det_name = 'FPC'
  if STRPOS(cmdb_fields(row_todo, c_detector),'SSD') GE 0 then det_name = 'SSD'

  toga_inout = 'TMA'
  if STRPOS(cmdb_fields(row_todo, c_grating), 'TOGA') GE 0 $
	then toga_inout = 'TOGA'   ; 'TMA' for no toga, 
                                   ; 'TOGA' for toga inserted

  if STRPOS(cmdb_fields(row_todo, c_shutter),'CLOSED') GE 0 then RETURN
  if STRPOS(cmdb_fields(row_todo, c_shutter),'ALL') GE 0 then begin
    ; all shutter open (possibly scanned)
    shutters = [0,0,0,0] ; megTop, legNorth, hegBottom, legSouth
                     ; 0 = shutter is out of beam, "open"
                     ; 1 = shutter blocks beam, "closed"
  end else begin
    if STRPOS(cmdb_fields(row_todo, c_shutter),'DISCRETE') GE 0 then begin
      pieces = str_sep(cmdb_fields(row_todo, c_shutter),',')
      if n_elements(pieces) NE 5 then RETURN  ; mystery DISCRETE
      pieces = pieces(1:4)
      open_shutters=where(STRPOS(pieces,'6') GE 0 OR STRPOS(pieces,'1') GE 0, nopen)
      if nopen LE 0 then RETURN  ; no open shutters?!
      shutters = [1,1,1,1]
      shutters(open_shutters) = 0
    end else RETURN  ; unknown shutter spec.
  end

  ; Set Y and Z locations
  
  ; Start with the CMDB values

  yoffset = 0.0
  yoffstr = cmdb_fields(row_todo, c_offset_y)
  if STRPOS(yoffstr,'FILE') LT 0 then READS, yoffstr, yoffset
  det_dy = yoffset

  zoffset = 0.0
  zoffstr = cmdb_fields(row_todo, c_offset_z)
  if STRPOS(zoffstr,'FILE') LT 0 then READS, zoffstr, zoffset
  det_dz = zoffset  ; " in Z


  ; Modify the Y Z values for HSI aim point offset and
  ; if hsi_pos is specified
  ; the column hsi_pos can contain the phrase "cusp1" to indicate
  ; the cusp as the location...
  if det_name EQ 'HSI' then begin
    if STRPOS(STRUPCASE(cmdb_fields(row_todo,c_hsi_pos)),'CUSP1') GE 0 then begin
      det_dy = det_dy + 0.180*25.4
      det_dz = det_dz + 0.199*25.4
    end else begin
      det_dy = det_dy + 0.0
      det_dz = det_dz + 1.6 ; mm FOA 1.6 mm below
    end
  end

  ; For '2C "0" aim point is at center of third quadrant
  if det_name EQ '2C' then begin
    det_dy = det_dy - 0.024*128.
    det_dz = det_dz
  end

  defoc = 0.0
  defocstr = cmdb_fields(row_todo, c_defocus)
  if STRPOS(defocstr,'FILE') LT 0 then READS, defocstr, defoc
  ; ACIS-S simulation is RC conforming so remove the RC from defocus value
  det_defocus = defoc - det_dy^2/5366.55 ; mm from XRCF focus, + is towards TMA/TOGA


  ; setting det_z_size
  CASE det_name OF
  '2C': begin  
    ; use the acis_frametime to select the z size in rows
    frametime = 0.0
    READS, cmdb_fields(row_todo, c_acis_frametime), frametime
    rows = 1024.
    if frametime LT 3.0 then rows = 200.
    if frametime LT 1.0 then rows = 114.
    if frametime LT 0.4 then rows = 54.
    if frametime LT 0.27 then rows = 38.
    if frametime LT 0.15 then rows = 18.
    if frametime LT 0.08 then rows = 1.
    det_z_size = rows  ; CCD rows
  end
  'HSI': begin
    det_z_size = 18.0 ;mm diameter
  end
  ELSE:  begin   ; assume this is FPC or SSD
    ; look to the hxda_aplist for the information
    aplist = cmdb_fields(row_todo,c_hxda_aplist)
    ; remove {}s
    aplist = STRMID(aplist,1,STRLEN(aplist)-2)
    ; divide it up by commas
    pieces = STR_SEP(aplist, ',')
    ; if it is more than one then do 0.5 mm
    if n_elements(pieces) GT 1 then begin
      det_z_size = 0.5  ; 500 um 
    end else begin
      ; if it is a slit do 200 um
      if STRPOS(pieces(0),'x') GE 0 then begin
        det_z_size = 0.2
      end else begin
        ; get the value
        det_z_size = 500.  ; um
        READS, pieces(0), det_z_size
        det_z_size = det_z_size/1000.
      end
    end
  end
  ENDCASE

  ; - - - - - - - - - -  Less frequently changed values - - - -
  ; Name for this simulation run... made up from above things...
  ; Name it by TRW ID !
  run_name = cmdb_fields(row_todo, c_id)
  ; Set the number of rays traced
  if toga_inout EQ 'TOGA' then begin
    ; With Gratings...
    if KEYWORD_SET(MILLION) then numraystring = '2000000' else $
		numraystring = '200000'
  end else begin
    ; without gratings...
    if KEYWORD_SET(MILLION) then numraystring = '1000000' else $
		numraystring = '40000'
  end
  output_root_dir = '/spectra/d7/phase3_sim/'+run_name

end

; same web directory for all:
; Directory for web output
web_dir = '/nfs/space/web/data/HETG/phase3_sim/'

; same mapping from det_name to det_type and det_ideal
CASE det_name OF
  '2C': begin
    det_type = 'ACIS-S'
    det_ideal = 'no' 
  end
  'HSI': begin
    det_type = 'HRC-S' ; 'ACIS-S' or 'HRC-S'
    det_ideal = 'no' ; 'yes' effic = 1, or 'no' to use nominal efficiency
  end
  'FPC': begin
    det_type = 'HRC-S' ; 'ACIS-S' or 'HRC-S'
    det_ideal = 'yes' ; 'yes' effic = 1, or 'no' to use nominal efficiency
  end
  'SSD': begin
    det_type = 'ACIS-S' ; 'ACIS-S' or 'HRC-S'
    det_ideal = 'no' ; 'yes' effic = 1, or 'no' to use nominal efficiency
  end
  ELSE: begin
    print, 'toga_full_sim: unrecognized detector name: '+det_name
    RETURN
  end

ENDCASE

; report parameters:
  print, ''
  print, ' spect_file    : '+spect_file
if   spect_file EQ 'use_Dans_labxray_routines' then begin
  print, ' . . . . . EIPS and Monos . . . . . . . .'
  print, '  target       : '+lx_tname
  print, '  voltage      : ',lx_sV
  print, '  current      : ',lx_sI
  print, '  mono_energy  : ',mono_energy,' keV'
  print, '  mono_res     : ',mono_res
  print, '  mono_highord : ',mono_highorders+ '0=DCM, 1=HIREFS'
end
  print, ' . . . . . . . . . . . . . . . . . . . . .'
  print, '  filter 1     : '+fmat1, fthk1
  print, '  filter 2     : '+fmat2, fthk2
  print, '  toga_inout   : '+toga_inout
  print, '  det_name     : '+det_name
  print, '  det_type     : '+det_type
  print, '  det_ideal    : '+det_ideal
  print, '  det_dy       : ',det_dy, ' (for 2C includes aim point setting)'
  print, '  det_defocus  : ',det_defocus, ' (from RC)'
  print, '  det_dz       : ',det_dz
  print, "  det_z_size   : ", det_z_size, ' 2C rows or mm diameter'
  print, '  shutters     : ', shutters
  print, '  run_name     : '+run_name
  print, '  output dir   : '+output_root_dir
  print, '  numraystring : '+numraystring
  if KEYWORD_SET(WEB) then print, '  web dir      : '+web_dir
  print, ' '

if NOT(KEYWORD_SET(SPECT_ONLY)) then begin
  ; Remove and re-make the output directories
  SPAWN, 'rm -R '+output_root_dir
  SPAWN, 'mkdir '+output_root_dir
  if KEYWORD_SET(WEB) then begin
    SPAWN, 'rm -R '+web_dir+run_name
    SPAWN, 'mkdir '+web_dir+run_name
  end
end
;.................................................................
; - - - - - - - - - - calculate source spectrum if desired - - - -
if spect_file EQ 'use_Dans_labxray_routines' then begin
  ; Create spectrum
  ;lx_tname already setup
  ;lx_sV already setup
  ;lx_sI = 1.0 ; mA
  lx_sangle = 45.0

  lx_emin = 0.1
  lx_emax = (lx_sV < 9.9)
  lx_estep = 0.0005
  ; reduce the emin for C and Be
  if lx_tname EQ 'Be' OR lx_tname EQ 'C' then begin
    lx_emin = 0.050
    lx_emax = (lx_emax < 5.0)
    lx_estep = (5.0/2.5E4)
  end

  ; fill lx_es and lx_tubespec
  ; lx_tubespec has units of photons/(sec bin steradian)

  ; add C and O contamination lines
  save_target = lx_tname
  save_curr = lx_sI
  ;
  lx_tname = 'C'
  contam_frac_I = 0.1
  lx_sI = contam_frac_I * save_curr
  labx_tubespec
  contam_spec = lx_tubespec
  ;
  lx_tname = 'O'
  contam_frac_I = 0.3
  lx_sI = contam_frac_I * save_curr
  labx_tubespec
  contam_spec = contam_spec + lx_tubespec
  ; put values back
  lx_tname = save_target
  lx_sI = save_curr

  ; real spectrum
  labx_tubespec

  ; add in the contamination
  lx_tubespec = lx_tubespec + contam_spec

  ; Include monochromator
  ; Mono is a filter that transmits 100% at first order, 3% at second,
  ; order 0.5% third order, 0.7% 4th order, 0.3% 5th order
  ; and 1.E-5 at other energies.
  if mono_energy GT 0.0 then begin
    mono_bins = where( approx_equal(lx_es, mono_energy, 0.5/mono_res), $
	nmonobins)
    if nmonobins LE 0 then RETURN
    cont_level = 1.E-6
    lx_tubespec = cont_level * lx_tubespec  ; continuum level
    lx_tubespec(mono_bins) = (1./cont_level) * lx_tubespec(mono_bins)
    if mono_highorders EQ 1 then begin
      ; add higher orders
      order_fracs = [0.03, 0.005, 0.007, 0.003]  ; HIREFS ?
      mono_order = [2.,3.,4.,5.]
      for io=0,n_elements(order_fracs)-1 do begin
        mono_bins = where( approx_equal(lx_es, mono_order(io)*mono_energy, $
  		0.5/mono_res), nmonobins)
        if nmonobins GT 0 then begin
          lx_tubespec(mono_bins) = (1./cont_level) * $
		order_fracs(io)*lx_tubespec(mono_bins)
        end
      end
    end
  end

  ; Include the filters
  fmats = [fmat1,fmat2]
  fthks = [fthk1,fthk2]
  for ifilter = 0,1 do begin  
    CASE fmats(ifilter) OF
      'Al': begin
          ; Catch zero thickness values - put filter in anyway
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=18.38 else fthks(ifilter)=36.76
          end
          lx_tubespec = lx_tubespec * aluminum(lx_es)^fthks(ifilter)
        end
      'Fe': lx_tubespec = lx_tubespec * iron(lx_es)^fthks(ifilter)
      'Be': begin
          ; Catch zero thickness values - put filter in anyway
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=2.*0.903 else fthks(ifilter)=4.*0.903
          end
          lx_tubespec = lx_tubespec * (beryllium(lx_es,/MICRON) > 1.E-10)^fthks(ifilter)
        end
      'C': begin
          ; Catch zero thickness values - put filter in anyway
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=2.*5.16 else fthks(ifilter)=4.*5.16
          end
          lx_tubespec = lx_tubespec * (parylene(lx_es) > 1.E-10)^fthks(ifilter)
        end
      'Mn': lx_tubespec = lx_tubespec * manganese(lx_es)^fthks(ifilter)
      'Mg': lx_tubespec = lx_tubespec * magnesium(lx_es)^fthks(ifilter)
      'Zr': begin
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=2.*1.98 else fthks(ifilter)=4.*1.98
          end
          lx_tubespec = lx_tubespec * zirconium(lx_es)^fthks(ifilter)
        end
      'Ti': lx_tubespec = lx_tubespec * titanium(lx_es)^fthks(ifilter)
      'Cr': lx_tubespec = lx_tubespec * chromium(lx_es)^fthks(ifilter)
      'Ni': lx_tubespec = lx_tubespec * nickel(lx_es)^fthks(ifilter)
      'OBF': lx_tubespec = lx_tubespec * acis2c_obf(lx_es)
      'NONE':
      'N/A':
      ELSE: print, '*** Unrecognized filter : '+fmats(ifilter)+'.'
    ENDCASE
  end
  ; convert to photons/cm^2/s/keV
  ; Real flux at XRCF at TMA input
  lx_tubespec = ( lx_tubespec/(lx_es(1)-lx_es(0)) ) * (1./53175.4)^2

  ; write out the data in photons/cm^2/s/ke
  spect_file = output_root_dir + '/spec.txt'
  OPENW, spect_unit, spect_file, /GET
  for ie=0,n_elements(lx_es)-1 do printf, spect_unit, $
	lx_es(ie), lx_tubespec(ie)
  CLOSE, spect_unit
  FREE_LUN, spect_unit
end else begin
  ; Penning source?
  if spect_file EQ 'use_Penning' then begin
    ; use a Kramer continuum 'cause there is some!
    lx_tZ = 0.1 ; Set to get continuum about right
    lx_tname = 'Penning'
    lx_sV = 0.9  ; kV
    lx_sI = 800.0 ; mA
    ; setup lx_ variables ...
    lx_emin = 0.050
    lx_emax = lx_sV
    lx_estep = 0.0001
    ; generate a grid
    lx_es = lx_emin+lx_estep*INDGEN(1+(lx_emax-lx_emin)/lx_estep)    
    ; and call the continuum model
    lx_tubespec = labx_tubeinuum(lx_es)
    if lx_sV GT 0.53 then begin
      ; add some O
      aboveO = where(lx_es GE 0.524)
      lx_tubespec(aboveO(0)) = lx_tubespec(aboveO(0))*1000.0
    end

    ; convert to photons/cm^2/s/keV
    ; Real flux at XRCF at TMA input
    lx_tubespec = ( lx_tubespec/(lx_es(1)-lx_es(0)) ) * (1./53175.4)^2

    ; Add in the Penning lines
    ; put in three place holders
    ; 160A , 10,000 ph/cm2sec
    ; 130A , 50 ph/cm2sec
    ; 56A , 4 ph/cm2sec
    ; p_lams = [160.,130.,56.]
    ; p_flux = [10000., 50., 4.]
; Use the energies in the PIGS-9 plot:
;    p_lams = lx_akeV/[0.0775, 0.095, 0.1, 0.105, $
;	0.111, 0.115, 0.119, 0.125, 0.132, 0.1375, 0.140, 0.145, $
;	0.210, 0.221, 0.258]
    ; wavelengths from JeffK from Finley SPIE vol.18 no 5 (1979)p.649
    p_lams = [170., 160., 130., 124.2, 117., $
	111.5, 108., 104., 99.5, 94., 90.8, 88.6, 85.6, $
	59., 56.]    
    p_flux = [2000.,8000., 50., 15., 12., $
	15., 15., 7., 1.5, 1.5, 1.5, 0.8, 0.8, $
	1.3, 4.]
    for il=0,n_elements(p_lams)-1 do begin
      e_lam = lx_akeV/p_lams(il)
      above_line = where(lx_es GE e_lam)
      lx_tubespec(above_line(0)) = lx_tubespec(above_line(0)) + $
		p_flux(il)/(lx_es(1)-lx_es(0))
    end
  ; ALWAYS have a 0.5 um parylene window
  lx_tubespec = lx_tubespec * (parylene(lx_es) > 1.E-10)^0.5

  ; Include the filters
  fmats = [fmat1,fmat2]
  fthks = [fthk1,fthk2]
  for ifilter = 0,1 do begin  
    CASE fmats(ifilter) OF
      'Al': begin
          ; Catch zero thickness values - put filter in anyway
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=18.38 else fthks(ifilter)=36.76
          end
          lx_tubespec = lx_tubespec * aluminum(lx_es)^fthks(ifilter)
        end
      'Fe': lx_tubespec = lx_tubespec * iron(lx_es)^fthks(ifilter)
      'Be': begin
          ; Catch zero thickness values - put filter in anyway
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=2.*0.903 else fthks(ifilter)=4.*0.903
          end
          lx_tubespec = lx_tubespec * (beryllium(lx_es,/MICRON) > 1.E-10)^fthks(ifilter)
        end
      'C': begin
          ; Catch zero thickness values - put filter in anyway
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=2.*5.16 else fthks(ifilter)=4.*5.16
          end
          lx_tubespec = lx_tubespec * (parylene(lx_es) > 1.E-10)^fthks(ifilter)
        end
      'Mn': lx_tubespec = lx_tubespec * manganese(lx_es)^fthks(ifilter)
      'Mg': lx_tubespec = lx_tubespec * magnesium(lx_es)^fthks(ifilter)
      'Zr': begin
          if fthks(ifilter) EQ 0 then begin
            if ifilter EQ 0 then fthks(ifilter)=2.*1.98 else fthks(ifilter)=4.*1.98
          end
          lx_tubespec = lx_tubespec * zirconium(lx_es)^fthks(ifilter)
        end
      'Ti': lx_tubespec = lx_tubespec * titanium(lx_es)^fthks(ifilter)
      'Cr': lx_tubespec = lx_tubespec * chromium(lx_es)^fthks(ifilter)
      'Ni': lx_tubespec = lx_tubespec * nickel(lx_es)^fthks(ifilter)
      'OBF': lx_tubespec = lx_tubespec * acis2c_obf(lx_es)
      'NONE':
      'N/A':
      ELSE: print, '*** Unrecognized filter : '+fmats(ifilter)+'.'
    ENDCASE
  end

    ; write out the data in photons/cm^2/s/keV
    spect_file = output_root_dir + '/spec.txt'
    OPENW, spect_unit, spect_file, /GET
    for ie=0,n_elements(lx_es)-1 do printf, spect_unit, $
	lx_es(ie), lx_tubespec(ie)
    CLOSE, spect_unit
    FREE_LUN, spect_unit
  end else begin
  ; something else?!?
    print, ' spect file is unrecognized: '+spect_file
    RETURN
  end
end


; - - - - - - - - - - 
if NOT(KEYWORD_SET(SPECT_ONLY)) then begin

; Now if the simulation is TMA then only need one pass of csim
; otherwise need to run csim up to 8 times to do the 4 grating quadrants
; and to do the direct beam (clear aperture) that leaks between
; the gratings into zero-order.

; check shutters are OK - need at least one 0 entry
open_shutters = where(shutters EQ 0, nopen)
if nopen LE 0 then begin
  print, ' * No open shutters : ', shutters
  return
  end
if toga_inout EQ 'TMA' then begin
  par_files = ['csim.par.tma']
  shutter_string = STRING(shutters, FORMAT='(4I1)')
  print, 'Shutter string =', shutter_string,'.'
end else begin
  all_par_files = ['csim.par.toga.megT', 'csim.par.toga.legN', $
		'csim.par.toga.hegB', 'csim.par.toga.legS']
  all_clr_files = ['csim.par.toga.clrT', 'csim.par.toga.clrN', $
		'csim.par.toga.clrB', 'csim.par.toga.clrS']
  ; just use the ones with open shutters:
  par_files = [ all_par_files(open_shutters) , all_clr_files(open_shutters) ]
end

; Loop through the par files and concatenate the simulations
for ip = 0, n_elements(par_files)-1 do begin

  ; Simulate
  par_file_name = par_files(ip)
  ; Because csim can modify some parameters in the par file
  ;  get a fresh copy from a write-protected file each time
  SPAWN, 'cp '+par_file_name+ ' csim.par'

  ; Make a sub-directory for each sub-configuration ray-traced
  pieces = str_sep(par_file_name,'.')
  npieces = n_elements(pieces)
  output_sub_dir = output_root_dir + '/sim.' + pieces(npieces-1)

  if toga_inout EQ 'TMA' then begin
    ; only need one csim call with shutters set as desired
    csim_line =  'csim @@csim.par '+ $
	'NumRays=' + numraystring + ' ' + $
	"SpectrumFile="+ spect_file + " " + $
	"Shutters6="+ shutter_string + " " + $
	"DetectorType="+ det_type + " " + $
	"ACISIdeal="+ det_ideal + " " + $
        ; -69.2 is the on-axis focus offset for operation of the TMA
        ; at finite source dist (see par files).
	"ACISOffsetX="+ STRCOMPRESS(STRING((det_defocus-69.2), $
		FORMAT='(F7.3)'),/REMOVE) + " " + $
	"HRCIdeal="+ det_ideal + " " + $
	"HRCOffsetX="+ STRCOMPRESS(STRING((det_defocus-69.2), $
		FORMAT='(F7.3)'),/REMOVE) + " " + $
	"OutputDir="+ output_sub_dir
  end else begin
    ; here for toga in beam - no need to set shutters
    ; because each par file has correct shutter open
    csim_line =  'csim @@csim.par '+ $
	'NumRays=' + numraystring + ' ' + $
	"SpectrumFile="+ spect_file + " " + $
	"DetectorType="+ det_type + " " + $
	"ACISIdeal="+ det_ideal + " " + $
        ; -69.2 is the on-axis focus offset for operation of the TMA
        ; at finite source dist (see par files).
	"ACISOffsetX="+ STRCOMPRESS(STRING((det_defocus-69.2), $
		FORMAT='(F7.3)'),/REMOVE) + " " + $
	"HRCIdeal="+ det_ideal + " " + $
	"HRCOffsetX="+ STRCOMPRESS(STRING((det_defocus-69.2), $
		FORMAT='(F7.3)'),/REMOVE) + " " + $
	"OutputDir="+ output_sub_dir
  end  ; selecting csim line
  print, '* * * * *'
  print, ' csim.par is from: '+par_file_name
  print, csim_line
  print, '* * * * *'
  SPAWN, csim_line
  if det_ideal EQ 'yes' then begin  ; ideal option does not create
                                    ; pulseheight.dat so make one...
    SPAWN, 'cp '+output_sub_dir+'/energy.dat '+output_sub_dir+'/pulseheight.dat'
  end

  ; reset variables to be filled
  ypos = FLTARR(1)
  zpos = FLTARR(1)
  chipx = FLTARR(1)
  chipy = FLTARR(1)
  pha = FLTARR(1)

  ; Simulate  ;Read from simulation directory
  read_csim_data, output_sub_dir, orders, $
		ypos, zpos, chipx, chipy, pha
  ; also get the times
  time = read_csim_file(output_sub_dir+'/time.dat')

  ; offset the values to put them into chip coord and
  ; not facility coords
  ypos = ypos - det_dy
  zpos = zpos - det_dz

  ; For FPC detector we need to add efficiency and energy blur
  if det_name EQ 'FPC' then begin
    ; include the FPC efficiency, hxda_fpc(E).
    randvals = RANDOMU(SEED,n_elements(pha))
    effics = hxda_fpc(pha)
    keep = where(randvals LE effics, nkeep)
    if nkeep GE 1 then begin
      ypos = ypos(keep)
      zpos = zpos(keep)
      pha = pha(keep)
      time = time(keep)
    end else begin
      ypos = [0.]
      zpos = [0.]
      pha = [0.]
      time = [0.]
    end

    ; Add escape effect
    e_esc = 2.957 ; Ar
    f_esc = 0.08  ; escape fraction
    randvals = RANDOMU(SEED,n_elements(pha))
    ie=LONG(0)
    while ie LE n_elements(pha)-1 do begin
      if pha(ie) GT e_esc then begin
        ; escape candidate
        if randvals(ie) LT f_esc then begin
          ; it is escaped!
          pha(ie) = pha(ie) - e_esc
        end
      end
      ie = ie+1
    end


    ; Blur the FPC energies using the FWHM(keV)=labx_smd_res(E) function
    esigma = labx_smd_res(pha)/2.35
    ie=LONG(0)
    while ie LE n_elements(pha)-1 do begin
      pha(ie) = pha(ie) + esigma(ie)*g_rand()
      ie = ie+1
    end
  end
 
  ; save/concatenate the simulations
  if ip EQ 0 then begin
    save_y = ypos
    save_z = zpos
    save_pha = pha
    save_time = time
  end else begin
    save_y = [save_y, ypos]
    save_z = [save_z, zpos]
    save_pha = [save_pha, pha]
    save_time = [save_time, time]
  end

; write to IDL save file too
SAVE, ypos, zpos, pha, time, FILENAME=output_sub_dir+'/sim_idl_save.dat'

end  ; loop over par files
print, ' - - - - - - - Simulation complete - - - - - - -'
print, ''

;................................................................
; write IDL save file here
ypos = save_y
zpos = save_z
pha = save_pha
time = save_time
SAVE, ypos, zpos, pha, time, FILENAME=output_root_dir+'/sim_idl_save.dat'
;................................................................

; Make some plots...

integ_time = MAX(save_time) - MIN(save_time)

end ; SPECT_ONLY

!p.multi = 0

; output a plot of the spectrum
; in photons/cm^2/s/keV
sourcexrange = [(MIN(lx_es) > 0.02),(MAX(lx_es) < 10.0)]
max_flux = MAX(lx_tubespec)
pre_print_sqr
plot_oo, lx_es, lx_tubespec, $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XTITLE = 'Energy (keV)', XRANGE = sourcexrange, $
	YTITLE='Flux (photons/cm^2/s/keV)', $
	YRANGE = [1.E-8*max_flux,max_flux]
device, /close
set_plot, 'X'
SPAWN, 'cp idl.ps '+output_root_dir+'/source_spectrum.ps'
SPAWN, 'cd '+output_root_dir+' ; ps2gif source_spectrum.ps'
if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+'/source_spectrum.gif '+$
	web_dir+run_name

; output a plot of the spectrum as seen by BND in 30 sec
; lx_tubespec is in photons/cm^2/s/keV
bnd_spec = lx_tubespec * 48.0 * 1000.0 * (lx_es(1)-lx_es(0))  ; photons per bin
; 
; For BND detector we need to add efficiency and energy blur
; include the FPC efficiency, hxda_fpc(E).
bnd_temp = hxda_fpc(lx_es) * bnd_spec
bnd_es = lx_es

; Add the blur
;  ... first ...
; Increase the bin size by factor of 20
bnd_estep = 20.0 * lx_estep
ipoint = LONG(0)
iout = 0
while ipoint+19 LE n_elements(bnd_temp)-1 do begin
  bnd_es(iout) = TOTAL(lx_es(ipoint:ipoint+19))/20.0
  bnd_spec(iout) = TOTAL(bnd_temp(ipoint:ipoint+19))
  ipoint = ipoint + 20
  iout = iout + 1
end
bnd_es = bnd_es(0:iout-1)
bnd_spec = bnd_spec(0:iout-1)

; add escape too!
save_tube = lx_tubespec
save_es = lx_es
lx_tubespec = bnd_spec
lx_es = bnd_es
lx_dname = 'SMD'
labx_escape
bnd_spec = lx_tubespec
lx_tubespec = save_tube
lx_es = save_es

; detector resolution vs energy, FWHM in keV
res = labx_smd_res(bnd_es)
res = res/(2.35*bnd_estep) ; one-sigma in bins

npts = n_elements(bnd_spec)
nptsm = npts-1
out_spec = 0. * bnd_spec
dout = 0. * bnd_spec

; go through the input spectrum
for in = 0, nptsm do begin
  res_e = res(in)  ; one_sigma in bins at this energy
  amp = bnd_spec(in)
  dout = 0. * dout
  for io = 0, nptsm do begin
    sigmas = (in-io)/res_e
    if(sigmas GT -6.0 AND sigmas LT 6.0) then $
	      dout(io) = EXP(-(sigmas^2)/2.)
  end
  out_spec = out_spec + (amp/res_e)*dout
end

; Add Poisson counting stats
bnd_spec = g_poiss(out_spec/SQRT(2.*!PI))

max_flux = MAX(bnd_spec)
pre_print_sqr
plot_oo, bnd_es, bnd_spec, PSYM=10, $
	TITLE = run_name + $
	' ( BND Counts,  1000 sec, '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XTITLE = 'Energy (keV)', XRANGE = sourcexrange, $
	YTITLE='Flux (BND events/bin, bin size ='+$
		STRING(bnd_estep*1000.0,FORMAT='(F6.1)')+' eV)', $
	YRANGE = [0.1,max_flux]
device, /close
set_plot, 'X'
SPAWN, 'cp idl.ps '+output_root_dir+'/bnd_spectrum.ps'
SPAWN, 'cd '+output_root_dir+' ; ps2gif bnd_spectrum.ps'
if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+'/bnd_spectrum.gif '+$
	web_dir+run_name

if NOT(KEYWORD_SET(SPECT_ONLY)) then begin

; event plot - in Facility coords and no spatial limits, but
; include the detector where ever it is
; or sub-frame aliasing
pre_print_sqr

; Find out the range of things to overplot...
; 
full_fac_y = save_y + det_dy
full_fac_z = save_z + det_dz
; Add in the detector outlines
if det_name EQ '2C' then begin
  ; over plot CCD and aliased-CCD sizes
  full_fac_y = [det_dy+1024.*0.024*[-0.5,0.5,0.5,-0.5,-0.5],full_fac_y]
  full_fac_z = [det_dz+1024.*0.024*[0.5, 0.5, -0.5, -0.5, 0.5],full_fac_z]
end else begin
  ; other detectors are circular, so draw the circle
  theta_start = 0.0
  theta_end = 360.  
  n_pts=17 * ABS(theta_end - theta_start)/90.
  X_circ = FLTARR(n_pts)
  Y_circ = FLTARR(n_pts)
  for ip=0, n_pts-1 do begin
    theta = theta_start + FLOAT(ip)*(theta_end - theta_start)/FLOAT(n_pts-1)
    X_circ(ip) = det_z_size/2. * COS(theta*!DTOR)
    Y_circ(ip) = det_z_size/2. * SIN(theta*!DTOR)
  end
  full_fac_y = [det_dy + X_circ, full_fac_y]
  full_fac_z = [det_dz + Y_circ, full_fac_z]
end

; Make it have same scale on each axis
all_fac = [full_fac_y, full_fac_z]
full_range = [MIN(all_fac),MAX(all_fac)]
plot, full_range, full_range, PSYM=3, $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(n_elements(save_y),FORMAT='(I8)'),/REMOVE_ALL)+ $
	' events  /  '+ $
	STRCOMPRESS(STRING(integ_time,FORMAT='(F8.2)'),/REMOVE_ALL) +' seconds '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XTITLE = 'Facility Y (mm)', $
	YTITLE='Facility Z (mm)', $
	CHARSIZE = 1.0, /NODATA
; Here are the events
oplot, save_y + det_dy, save_z + det_dz, PSYM=3
; and over plots
if det_name EQ '2C' then begin
  oplot, det_dy+1024.*0.024*[-0.5,0.5,0.5,-0.5,-0.5], $
	det_dz+1024.*0.024*[0.5, 0.5, -0.5, -0.5, 0.5]
  oplot, det_dy+1024.*0.024*[-0.5,0.5,0.5,-0.5,-0.5], $
	det_dz+det_z_size*0.024*[0.5, 0.5, -0.5, -0.5, 0.5]
end else begin
  ; other detectors are circular, so draw the circle
  oplot, det_dy + X_circ, det_dz + Y_circ
  ; if it is HSI get the mask outline too
  if det_name EQ 'HSI' then begin
    hsi_mask, x_hsi, y_hsi
    oplot, det_dy + x_hsi, det_dz + y_hsi
    hsi_bad_oplot, det_dy, det_dz
  end
end
  
device, /close
set_plot, 'X'
SPAWN, 'cp idl.ps '+output_root_dir+'/facility_events.ps'
SPAWN, 'cd '+output_root_dir+' ; ps2gif facility_events.ps'
if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+'/facility_events.gif '+$
	web_dir+run_name

; energy vs y plot
pre_print_sqr
plot, save_y + det_dy, save_pha, PSYM=3, $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(n_elements(save_y),FORMAT='(I8)'),/REMOVE_ALL)+ $
	' events  /  '+ $
	STRCOMPRESS(STRING(integ_time,FORMAT='(F8.2)'),/REMOVE_ALL) +' seconds '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XTITLE = 'Facility Y (mm)', YTITLE='Energy (keV)', $
	CHARSIZE = 1.0
device, /close
set_plot, 'X'
SPAWN, 'cp idl.ps '+output_root_dir+'/e_vs_fac_y.ps'
SPAWN, 'cd '+output_root_dir+' ; ps2gif e_vs_fac_y.ps'
if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+'/e_vs_fac_y.gif '+$
	web_dir+run_name

; event plot - use events that are on/in the detector
if det_name EQ '2C' then begin
  ;chip coords limited to one chip spatial limit in Y
  onchip = where(save_y GT -12.29 AND save_y LT 12.29 AND save_z GT -12.29 $
	AND save_z LT 12.29, nonchip)
end else begin
  onchip = where(SQRT(save_y^2 + save_z^2) LT det_z_size/2.0, nonchip)
end

if nonchip GT 1 then begin

  ; use only the events in/on the detector
  chip_y = save_y(onchip)
  chip_z = save_z(onchip)
  chip_pha = save_pha(onchip)

  ; nominal range 
  height = det_z_size
  width = det_z_size

  if det_name EQ '2C' then begin
    ; now alias the Z values
    ; det_z_size is full height of region in rows (pixels)
    ; convert it to millimeters
    height = det_z_size*0.024
    width = 1024.*0.024 ; mm
    above = where(chip_z GT -1.*height/2.,nabove)
    if nabove GE 1 then begin
      ratio = (chip_z(above) - height/2.0)/height
      frac_ratio = ratio - FLOAT(FIX(ratio))
      chip_z(above) = height * (frac_ratio - 0.5)
    end
    below = where(chip_z LT height/2.,nbelow)
    if nbelow GE 1 then begin
      ratio = (height/2. - chip_z(below) )/height
      frac_ratio = ratio - FLOAT(FIX(ratio))
      chip_z(below) = height * ((1.-frac_ratio) - 0.5)
    end
  end

  pre_print_sqr
  plot, chip_y, chip_z, PSYM=3, $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(n_elements(chip_y),FORMAT='(I8)'),/REMOVE_ALL)+ $
	' events  /  '+ $
	STRCOMPRESS(STRING(integ_time,FORMAT='(F8.2)'),/REMOVE_ALL) +' seconds '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XTITLE = 'Y '+det_name+' (mm)', XRANGE=width*[-0.5,0.5], $
	YTITLE='Z '+det_name+' (mm)', YRANGE = height*[-0.5,0.5], $
	CHARSIZE = 1.0
  ; overplot the circle for non 2C detectors
  if det_name NE '2C' then oplot, X_circ, Y_circ
  if det_name EQ 'HSI' then begin
    oplot, x_hsi, y_hsi
    ; add bad regions
    hsi_bad_oplot
  end
  device, /close
  set_plot, 'X'
  SPAWN, 'cp idl.ps '+output_root_dir+'/'+det_name+'_events.ps'
  SPAWN, 'cd '+output_root_dir+' ; ps2gif '+det_name+'_events.ps'
  if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+'/'+det_name+'_events.gif '+$
	web_dir+run_name

if det_name EQ '2C' OR det_name EQ 'HSI' then begin
  pre_print_sqr
  plot, chip_y, chip_pha, PSYM=3, $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(n_elements(chip_y),FORMAT='(I8)'),/REMOVE_ALL)+ $
	' events  /  '+ $
	STRCOMPRESS(STRING(integ_time,FORMAT='(F8.2)'),/REMOVE_ALL) +' seconds '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )' , $
	XRANGE=width*[-0.5,0.5], XTITLE = 'Y '+det_name+' (mm)', $
	YTITLE='Energy (keV)', $
	CHARSIZE = 1.0
  device, /close
  set_plot, 'X'
  SPAWN, 'cp idl.ps '+output_root_dir+'/e_vs_'+det_name+'_y.ps'
  SPAWN, 'cd '+output_root_dir+' ; ps2gif e_vs_'+det_name+'_y.ps'
  if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+ $
	'/e_vs_'+det_name+'_y.gif '+web_dir+run_name

  ; position histogram
  bin_size = 0.05
  histmax = 0.0
  histmin = 0.0
  hist_array = histogram(chip_y, BINSIZE=bin_size,MIN=-15., MAX=15., $
	OMAX = histmax, OMIN = histmin)
  nhist = n_elements(hist_array)
  xarray = (histmax-histmin)*(indgen(nhist)/FLOAT(nhist-1))+histmin
  pre_print_sqr
  plot_io, xarray, hist_array,YRANGE=[1.0,MAX(hist_array)],PSYM=10, $
	XTITLE=det_name+' Y (mm)', YTITLE='Counts/bin (bin size = '+ $
	'50 um)', $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(n_elements(chip_y),FORMAT='(I8)'),/REMOVE_ALL)+ $
	' events  /  '+ $
	STRCOMPRESS(STRING(integ_time,FORMAT='(F8.2)'),/REMOVE_ALL) +' seconds '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )'
  device, /close
  set_plot, 'X'
  SPAWN, 'cp idl.ps '+output_root_dir+'/'+det_name+'_y_hist.ps'
  SPAWN, 'cd '+output_root_dir+' ; ps2gif '+det_name+'_y_hist.ps'
  if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+ $
	'/'+det_name+'_y_hist.gif '+web_dir+run_name
end  ; 2C and HSI only

  ; energy histogram
  bin_size = 0.01
  if det_name EQ 'FPC' then bin_size = 0.05
  if lx_emin LT 0.09 then begin  ; Penning and low E anodes
    if det_name EQ 'HSI' then bin_size = 0.001
    if det_name EQ 'FPC' then bin_size = 0.010
  end
  histmax = 0.0
  histmin = 0.0
  hist_array = histogram(chip_pha, BINSIZE=bin_size,MIN=0., MAX=10., $
	OMAX = histmax, OMIN = histmin)
  nhist = n_elements(hist_array)
  xarray = (histmax-histmin)*(indgen(nhist)/FLOAT(nhist-1))+histmin
  ; set the energy range based on the lx_emin value
  erange = [0.1,10.]
  if lx_emin LT 0.09 then erange = [0.01, 1.0]
  pre_print_sqr
  plot_oo, xarray, hist_array, $
	XRANGE= erange, YRANGE=[1.,MAX(hist_array)],PSYM=10, $ 
	XTITLE='Energy (keV)', YTITLE='Counts detected/bin (bin size = '+ $
	STRCOMPRESS(STRING(bin_size*1.E3,FORMAT='(F8.1)'),/REMOVE_ALL)+' eV)', $
	TITLE = run_name + $
	' ( '+ $
	STRCOMPRESS(STRING(n_elements(chip_pha),FORMAT='(I8)'),/REMOVE_ALL)+ $
	' events  /  '+ $
	STRCOMPRESS(STRING(integ_time,FORMAT='(F8.2)'),/REMOVE_ALL) +' seconds '+ $
	STRCOMPRESS(STRING(lx_sI,FORMAT='(F7.2)'),/REMOVE_ALL)+' mA )'
  device, /close
  set_plot, 'X'
  SPAWN, 'cp idl.ps '+output_root_dir+'/'+det_name+'_e_hist.ps'
  SPAWN, 'cd '+output_root_dir+' ; ps2gif '+det_name+'_e_hist.ps'
  if KEYWORD_SET(WEB) then SPAWN, 'cp '+output_root_dir+'/'+det_name+'_e_hist.gif '+$
	web_dir+run_name

  ; If the maximum value of the spectrum and histogram is a single line then
  ; this will calculate the effective area in that line...
  line_counts = MAX(hist_array)
  print, 'Max counts in histogram (line)= ', line_counts
  print, '  line is at ',xarray(where(hist_array EQ line_counts)), ' keV'
  print, 'Integration time = ', integ_time, ' s'

  ; lx_tubespec is in photons/cm^2/s/keV
  ; Real flux in line at XRCF at TMA input in phots/cm^2/s:
  line_flux = MAX(lx_tubespec)*(lx_es(1)-lx_es(0))
  print, 'Line flux = ',line_flux,' phots/cm^2 s'
  line_area = line_counts/(line_flux*integ_time)
  print, 'Effective area = ',  line_area, ' cm^2'

end  ; checking for photons on chip

end ; SPECT_ONLY

!p.multi = 0


RETURN
END