PRO hsi_source_spectra, trw_id, LPR = lpr, NO_ANAL=no_anal, $ COMPARE=compare, PS_COMPARE = ps_compare, HEAE_ANAL=heae_anal ;+ ; Analyze HSI images to extract source spectra ; /LPR creates output plot files ; /NO_ANAL does no processing: useful to see the various parameters ; of the measurement printed out. ; /COMPARE converts counts/keV to photons/cm^2 keV and compares ; result with labxray source model. ;- ; 12/15/97 dd Added /COMPARE: conversion of counts/keV spectra to ; photons/cm^2 keV using HRMA-Grating-HSI CIPs ; 4/8/99 dd Update the HRMA, Grating, and HSI QE values from CALDBcal. ; 4/13/99 dd Do a QE depression correction when reading in the HSI data, ; burried below are the lines: ; qed_value = hsi_qe_depress(run_id_numerical) ; hsi_data = xrcf_hsi_read(file_name,/FAC,EXPOSURE=hsi_time, $ ; RAW_HSI= raw_hsi, $ ; QE_DEPRESS=qed_value, QE_INFO_STR=qe_info_str) ; And also make "tile plots" for central HSI tile, using: ; hsi_tile_plots, raw_hsi, plot_title, QE_INFO_STR = qe_info_str ; 4/23/99 dd Add the HEAE_ANAL option to add counts info to the ; heae_whcounts.rdb file. ; - The /HEAE_ANAL keyword to hsi_source_spectra will ; have it read heae_whcounts.rdb, modify the appropriate line, ; by adding focal_counts, focal_err_plus, focal_err_minus, ; focal_time, focal_analysis=ROI,ALL , etc. ; and then write it back out. ; 5/18,31/99 dd Use new HRMA 'N0004 QEs and HSI QE ('990530). ; - - - - - - - - - - - - - @cmdb_common @labx_common ; Make a common for big, long-to-read-in things COMMON internal, trw_run_ids, cip_hrma, cip_heg, cip_meg, cip_leg, cip_hsi ; - - - - - - - - - - - - - - - - - - ; Preparation work... ; Use TRW ID to derrive all other info... ; lookup table to map TRW-ID to runid and date: if n_elements(trw_run_ids) LT 2 then begin trw_run_ids = rdb_read(!DDIDL+'/xrcf/eae_trw_mod.rdb') end ; Get cal products from CALDB cal areas ; Load them if not already if (KEYWORD_SET(COMPARE) AND n_elements(cip_hrma) LT 10) then begin print, 'Loading CIPs...' print, ' HSI .' cip_hsi = rdb_read(!DDHXDSDIR+'/hsi_qe_N19990530.rdb') ; Use the HRMA N0003 XRCF value to remove artificial 2 keV jump... ; Use the HRMA N0004 XRCF values (artificial 0.94 keV jump) print, ' HRMA at XRCF ...' cip_hrma = rdb_read(!DDHRMACAL+'/cip/hrmaD1996-11-01XeffareaN0004.rdb') print, ' MEG .....' cip_meg = rdb_read(!DDHETGCAL+'/cip/hetgmegD1996-11-01greffpr005N0002.rdb') print, ' HEG .....' cip_heg = rdb_read(!DDHETGCAL+'/cip/hetghegD1996-11-01greffpr005N0002.rdb') print, ' LETG ..........' cip_leg = rdb_read(!DDLETGCAL+'/cip/letgD1996-11-01greffpr005N0002.rdb') endif ; Setup where the data come from data_prefix = !DDHXDSDATA+'/phase1/data/' ; Assume CMDB is available if cmdb_using_cal GT 0 then begin @cmdb_cal_integers end else begin @cmdb_trw_integers end ; Get XRCF parameters xrcf_params, ANGLES=angles, PERIODS=periods, HC=hc, $ ROWLANDS=rowlands, NAMES=names ; Spectral extraction width spec_width = 2.0 ; mm ; Where to keep the resulting files ; Not really something to fill up CALDB area with... save_em_dir = !DDHSISPECTRA ; - - - - - - - - - - - - - - - - - - ; Get and Report Parameters needed for analysis print, '..............................' print, ' TRW ID : '+trw_id this_meas = where(trw_run_ids.trwid EQ trw_id AND $ STRPOS(trw_run_ids.irig_start,'NULL') LT 0, nfound) if nfound GE 1 then begin run_id = STRCOMPRESS(STRING( $ LONG(trw_run_ids(this_meas(0)).runid) ),/REMOVE_ALL) run_id_numerical = trw_run_ids(this_meas(0)).runid irig_str = trw_run_ids(this_meas(0)).irig_start pieces = STR_SEP(irig_str,' ') date_pieces = STR_SEP(pieces(0),'/') if n_elements(date_pieces) EQ 3 then begin d_m = date_pieces(0) if STRLEN(d_m) EQ 1 then d_m = '0'+date_pieces(0) d_d = date_pieces(1) if STRLEN(d_d) EQ 1 then d_d = '0'+date_pieces(1) d_y = STRMID(date_pieces(2),2,2) date_str = d_y + d_m + d_d end else begin print, 'For '+trwid+ $ " didn't get 3 pieces from irig date: "+irig_str date_str = 'NULL' end date_code = date_str end else begin print, 'TRW ID '+trw_id+' not found in xrcf/eae_trw_mod.rdb' RETURN end print, 'date code : '+date_code print, ' run id : '+run_id ; Use CMDB to find out the grating type, source, and energy this_meas = where(cmdb_fields(*,c_id) EQ trw_id, nfound) grating = -1 if nfound EQ 1 then begin if cmdb_fields(this_meas(0),c_grating) EQ 'LETG' then grating = 0 if cmdb_fields(this_meas(0),c_grating) EQ 'HETG' then begin shutter_str = cmdb_fields(this_meas(0),c_shutters) if STRPOS(shutter_str,'MEG') GE 0 then grating = 1 if STRPOS(shutter_str,'HEG') GE 0 then grating = 2 if STRPOS(shutter_str,'ALL,ALL') GE 0 then grating = -1 end ; get a nice source string source_str = cmdb_source_str(this_meas(0)) print, ' source : '+source_str ; and the CMDB line energy src_energy = cmdb_data(this_meas(0)).energy print, ' energy : '+STRING(src_energy,FORMAT='(F8.4)')+' keV' end else begin print, trw_id+' found no or multiple times in cmdb :'+cmdb_file RETURN end ; Find the HSI location in the focal plane if STRPOS(trw_id, '-3D-') GT 0 OR STRPOS(trw_id, '-dF-') GT 0 then begin ; Read the loc file to know where the HSI images are xrcf_loc_file_read, trw_id, Xloc, Yloc, Zloc end else begin ; Use the CMDB location ir = this_meas(0) if STRPOS(cmdb_fields(ir,c_measurement_config),'HXDA') $ GE 0 then begin ; Using HXDA Y_str = cmdb_fields(ir,c_HXDA_Y) Z_str = cmdb_fields(ir,c_HXDA_Z) X_str = cmdb_fields(ir,c_HXDA_defocus) end else begin ; Using FAM etc. Y_str = cmdb_fields(ir,c_five_axis_y) Z_str = cmdb_fields(ir,c_five_axis_Z) X_str = cmdb_fields(ir,c_five_axis_X) end temp = 0.0 READS, X_str, temp Xloc = [temp] READS, Y_str, temp Yloc = [temp] READS, Z_str, temp Zloc = [temp] end ; Some kludges ; For D-HXH-3D-11.043(44) the grating is MEG(HEG) and it is offset correctly ; .loc file was calculated at zero angle because of "ALL,ALL" shutters. ; So the Yloc is really the dispersion hypooneus... if trw_id EQ 'D-HXH-3D-11.043' then begin grating = 1 Zloc = Yloc * SIN(angles(grating)*!DTOR) Yloc = Yloc * COS(angles(grating)*!DTOR) end if trw_id EQ 'D-HXH-3D-11.044' then begin grating = 2 Zloc = Yloc * SIN(angles(grating)*!DTOR) Yloc = Yloc * COS(angles(grating)*!DTOR) end if trw_id EQ 'D-HXH-3D-11.037' then begin grating = 1 Zloc = Yloc * SIN(angles(grating)*!DTOR) Yloc = Yloc * COS(angles(grating)*!DTOR) end if trw_id EQ 'D-HXH-3D-11.038' then begin grating = 2 Zloc = Yloc * SIN(angles(grating)*!DTOR) Yloc = Yloc * COS(angles(grating)*!DTOR) end if trw_id EQ 'D-HXH-EA-13.012' then begin ; Both MEG and HEG first orders are visible grating = 2 end if trw_id EQ 'D-LXH-EA-13.011' then begin ; The Y location is slightly off resulting in a split ; spectral peak... adjust for aligned 3rd order: Yloc = 0.05 ; mm ; (the neg first order is on HSI gap it appears) end ; Use location to help uncover grating being used ; if ALL,ALL if grating LT 0 then begin if cmdb_fields(this_meas(0),c_grating) EQ 'HETG' then begin if ABS(Yloc(0)) GT 0.5 then begin if TOTAL(Zloc/Yloc GT 0.05) GT 0 then begin ; MEG grating = 1 end if TOTAL(Zloc/Yloc LT -0.05) GT 0 then begin ; HEG grating = 2 end end end end ; print, location(s) for iloc = 0, n_elements(Xloc)-1 do begin print, ' Y,Z,X : '+STRING(Yloc(iloc)) + $ STRING(Zloc(iloc))+STRING(Xloc(iloc)) end if grating LT 0 then begin print, 'TRW ID '+trw_id+' not well-defined grating type' RETURN end else begin print, ' grating : '+names(grating) ; Calculate the order based on the location, energy, and ; grating disp_first = rowlands(grating)*(hc/src_energy)/periods(grating) sign = FIX(1.001*Yloc(0)/(ABS(Yloc(0)) > 1.0) ) ; gives -1,0,+1 order = sign * FIX((ABS(Yloc(0))+disp_first/2.0 )/disp_first) ; if the order is 0 then use 1 (probably a high-energy ; line so image was zero-centered...) if order EQ 0 then order = 1 print, ' order : '+ STRING(order,FORMAT='(I3)') end ; Special order selection if trw_id EQ 'E-HXH-3D-10.008a' then order = 1 if trw_id EQ 'D-LXH-3D-9.006' then order = 5 ; The iteration is generally 0 but in case it is not... iter = 0 ; Special cases... if trw_id EQ 'D-LXH-dF-16.001' then iter = 2 if trw_id EQ 'D-HXH-dF-16.002' then iter = 2 if trw_id EQ 'D-HXH-dF-16.003' then iter = 2 print, ' iter : '+ STRING(iter,FORMAT='(I3)') ; - - - - - - - - - - - - - - - - - - - - - - ; Done with preliminaries... if KEYWORD_SET(NO_ANAL) then RETURN ; Do only one iteration it_id = STRCOMPRESS(STRING(iter),/REMOVE_ALL) ; Read the data file file_name = data_prefix+date_code+'/hsi'+run_id+'i'+it_id+'.fits' ; Do QE depression correction when reading in the HSI data: qed_value = hsi_qe_depress(run_id_numerical) hsi_data = xrcf_hsi_read(file_name,/FAC,EXPOSURE=hsi_time, $ RAW_HSI= raw_hsi, $ QE_DEPRESS=qed_value, QE_INFO_STR=qe_info_str) ; EXPOSURE keyword is the actual exposure time in seconds from ; the FITS header print, 'hsi_time (sec) = ', hsi_time ; Create "global" coordinates w.r.t. zero-order glob_Y = hsi_data.Y + 1000.*Yloc(iter) glob_Z = hsi_data.Z + 1000.*Zloc(iter) ; Create wavelengths of each photon disp_dist = (ABS(glob_Y)/1000.0)/COS(!DTOR*angles(grating)) wavelengths = ((periods(grating)/ABS(FLOAT(order)))* $ (disp_dist/rowlands(grating))) ; Plot to .gif set_plot, 'X' window, 0, xsi=500, ysi=700 dd_load_ct, /OTHER ; Make the hsi "tile" plots here ; use the raw HSI data plot_title = trw_id + $ ' ('+date_code+'/hsi'+run_id+'i'+it_id+'.fits) ' hsi_tile_plots, raw_hsi, plot_title, QE_INFO_STR = qe_info_str if KEYWORD_SET(LPR) then begin ; Finish the .gif output wshow,0,iconic=0 image = tvrd() tvlct, red, green, blue, /GET gif_file = save_em_dir+'/'+trw_id+'_tile.gif' write_gif, gif_file, image, red, green, blue end else begin wait, 5. end ; Plot to .gif set_plot, 'X' window, 0, xsi=500, ysi=700 dd_load_ct, /OTHER !p.multi = [0,1,3] char_size = 1.5 ; plot the events plot, glob_Y, glob_Z, PSYM=3, BACK=0, COLOR=1, $ XRANGE = 1000.0*Yloc(iter)+[-9000.,9000.], $ YRANGE = 1000.0*Zloc(iter)+[-9000.,9000.], $ TITLE = trw_id+': '+source_str+', '+names(grating)+ $ '['+STRCOMPRESS(STRING(order),/REMOVE_ALL)+'] ('+ $ date_code+'/hsi'+run_id+'i'+ $ STRCOMPRESS(STRING(iter),/REMOVE_ALL)+'.fits)', $ XTITLE = 'Facility Y (um)', YTITLE = 'Facility Z (um)', $ CHARSIZE = char_size ; over plot extraction region for spectrum and background ys = 1000.0*Yloc(iter)+[-9000.,9000.] zs = ys*TAN(angles(grating)*!DTOR) oplot, ys, zs - 1000.0*spec_width/2.0, COLOR=10 oplot, ys, zs - 1000.0*spec_width, COLOR=10 oplot, ys, zs + 1000.0*spec_width/2.0, COLOR=10 oplot, ys, zs + 1000.0*spec_width, COLOR=10 ; OK, jump in here and do the HEAE_ANAL stuff if requested... if KEYWORD_SET(HEAE_ANAL) then begin ; read in the heae_whcounts.rdb structure heae_dir = !DDHETGCAL+'/cmp/eae/' heae = rdb_read(heae_dir+'heae_whcounts.rdb', /SILENT) ; find this trw_id in the heae this_one = where(heae.trw_id EQ trw_id, n_found) if n_found EQ 1 then begin this_one = this_one(0) ; load in mst_date, run_id, and iteration heae(this_one).mst_date = date_code heae(this_one).run_id = run_id heae(this_one).iteration = STRCOMPRESS(STRING(iter),/REMOVE_ALL) ; load in focal_time, focal_analysis=ROI,All heae(this_one).focal_time = hsi_time heae(this_one).focal_analysis = "ROI,All" ; load in focal_counts, focal_err_plus, focal_err_minus ; Find the number of events in circle of diameter heae.fp_aperture microns. ; about hsiY,Z = 0,0: ap_radius = heae(this_one).fp_aperture/2.0 in_circ = where ( SQRT(hsi_data.Y^2 + hsi_data.Z^2) LT $ ap_radius , n_found) heae(this_one).focal_counts = FLOAT(n_found) heae(this_one).focal_err_plus = SQRT(heae(this_one).focal_counts) heae(this_one).focal_err_minus = heae(this_one).focal_err_plus rdb_write, heae_dir+'heae_whcounts.rdb', heae, /SILENT ; overplot the HEAE extraction region circle facet_circle, 0.0, 360.0, circX, circY oplot, 1000.*Yloc(iter)+circX*ap_radius, $ 1000.*Zloc(iter)+circY*ap_radius, COLOR=100 ; and label the plot xyouts, 1000.*Yloc(iter), 1000.*Zloc(iter) + 1.15*5000.0, $ STRCOMPRESS(n_found)+' counts in ' + $ STRING(hsi_time,FORMAT='(F6.1)')+' s, ' + $ STRCOMPRESS(FIX(heae(this_one).fp_aperture)) + ' micron diameter', $ ALIGNMENT=0.5, CHARSIZE=1.2, color=100 ; all done, write out the heae_whcounts.rdb file end else begin print, '' print, ' * Did not find '+trw_id+' in heae_wcounts.rdb file ...' print, '' end end ; select spectrum and background photons expected_z = glob_Y * TAN(angles(grating)*!DTOR) spec_photons = where(glob_Z GT expected_z - 1000.0*spec_width/2.0 AND $ glob_Z LT expected_z + 1000.0*spec_width/2.0) bkg_photons = where( (glob_Z GT expected_z - 1000.0*2.*spec_width/2.0 AND $ glob_Z LT expected_z - 1000.0*spec_width/2.0) OR $ (glob_Z GT expected_z + 1000.0*spec_width/2.0 AND $ glob_Z LT expected_z + 1000.0*2.*spec_width/2.0) ) ; plot the selected events plot, glob_Y(spec_photons), glob_Z(spec_photons), PSYM=3, BACK=0, COLOR=1, $ XRANGE = 1000.0*Yloc(iter)+[-9000.,9000.], $ TITLE = 'Events in Spectral-extraction Region,' + $ STRING(spec_width,FORMAT='(F5.2)')+' mm wide', $ XTITLE = 'Facility Y (um)', YTITLE = 'Facility Z (um)', $ CHARSIZE = char_size oplot, ys, zs - 1000.0*spec_width/2.0, COLOR=10 oplot, ys, zs + 1000.0*spec_width/2.0, COLOR=10 ; plot the background events plot, glob_Y(bkg_photons), glob_Z(bkg_photons), PSYM=3, BACK=0, COLOR=1,$ XRANGE = 1000.0*Yloc(iter)+[-9000.,9000.], $ TITLE = 'Events in Background-extraction Regions', $ XTITLE = 'Facility Y (um)', YTITLE = 'Facility Z (um)', $ CHARSIZE = char_size oplot, ys, zs - 1000.0*spec_width/2.0, COLOR=10 oplot, ys, zs - 1000.0*spec_width, COLOR=10 oplot, ys, zs + 1000.0*spec_width/2.0, COLOR=10 oplot, ys, zs + 1000.0*spec_width, COLOR=10 ; Set the wavelength scale... disp_ave = ABS(Yloc(iter)) disp_max = disp_ave + 9.0 disp_min = (disp_ave - 9.0) > 1.0 lam_max = ((periods(grating)/ABS(FLOAT(order)))* $ (disp_max/rowlands(grating))) lam_min = ( ((periods(grating)/ABS(FLOAT(order)))* $ (disp_min/rowlands(grating))) > 1.24 ) dlam = (lam_max - lam_min)/500.0 lam_spec_hist = histogram(wavelengths(spec_photons), $ BINSIZE=dlam, MIN=lam_min, MAX=lam_max) lam_bkg_hist = histogram(wavelengths(bkg_photons), $ BINSIZE=dlam, MIN=lam_min, MAX=lam_max) lams = lam_min + (lam_max-lam_min)*(INDGEN(n_elements(lam_spec_hist))+0.5)/ $ FLOAT(n_elements(lam_spec_hist)) print, 'Lambdas: MIN, MAX: ',MIN(lams), MAX(lams) if KEYWORD_SET(LPR) then begin ; Finish the .gif output wshow,0,iconic=0 image = tvrd() tvlct, red, green, blue, /GET gif_file = save_em_dir+'/'+trw_id+'_image.gif' write_gif, gif_file, image, red, green, blue end else begin wait, 5. end ; Plot to .gif set_plot, 'X' window, 0, xsi=500, ysi=700 dd_load_ct, /OTHER !p.multi = [0,1,2] char_size = 1.3 plot_io, lams, lam_spec_hist, PSYM=10, BACK=0, COLOR=1, $ TITLE = trw_id+': '+source_str+', '+names(grating)+ $ '['+STRCOMPRESS(STRING(order),/REMOVE_ALL)+']', $ XTITLE = 'Wavelength (A)', YTITLE = 'counts/bin '+ $ '(bin size = '+STRING(dlam,FORMAT='(F7.4)')+' A)', $ YRANGE=[1.,MAX(lam_spec_hist(where(lams GE 0.5)))], $ XRANGE = [lam_min,lam_max], $ CHARSIZE = char_size oplot, lams, lam_bkg_hist, PSYM=1 ; Add indication of CMDB assumed energy oplot, hc/src_energy*[1.,1.], [1.,1.E6], LINESTYLE=2, COLOR=100 ; Make a plot vs Energy spec_es = hc/lams ; Convert the spectrum to units of counts/keV spec_per_kev = ((lams^2)/hc) * (lam_spec_hist - lam_bkg_hist)/dlam plot_oo, spec_es, spec_per_kev, PSYM=10, BACK=0, COLOR=1, $ TITLE = '('+ $ date_code+'/hsi'+run_id+'i'+ $ STRCOMPRESS(STRING(iter),/REMOVE_ALL)+'.fits)', $ XTITLE = 'Energy (keV)', YTITLE = 'counts/keV '+ $ '(bin size = '+STRING(dlam,FORMAT='(F7.4)')+' A)', $ YRANGE=[( MIN(spec_per_kev(where(spec_per_kev GT 0.))) > 1.), $ MAX(spec_per_kev(where(spec_es LT 10.0)))], $ XRANGE = [0.1,10.], $ CHARSIZE = char_size ; Add indication of CMDB assumed energy oplot, src_energy*[1.,1.], [1.E-3,1.E9], LINESTYLE=2, COLOR=100 if KEYWORD_SET(LPR) then begin ; Finish the .gif output wshow,0,iconic=0 image = tvrd() tvlct, red, green, blue, /GET gif_file = save_em_dir+'/'+trw_id+'_spec.gif' write_gif, gif_file, image, red, green, blue end if KEYWORD_SET(COMPARE) OR KEYWORD_SET(PS_COMPARE) then begin ; OK, now convert to photons/s cm^2 keV and compare with ; xss_sim model... ; HRMA area: cip_hrma - energy, shell1, shell3, shell4, shell6 ; HSI QE: cip_hsi - energy qe ; grating effic: cip_heg, cip_meg, cip_leg - energy, om1, oz, op1, etc. ; ; Expected "area" (cm^2 counts/photon) ; Use grating energy grid ; Calculate it on the spec_es energy grid CASE grating OF 0: begin ; LEG ; y = INTERPOL(Y X x) hrma_area = INTERPOL(cip_hrma.ea_hrma, cip_hrma.energy, spec_es) grating_effic = INTERPOL(cip_leg.(((n_tags(cip_leg))/2)+order), $ cip_leg.energy, spec_es) end 1: begin ; MEG hrma_area = INTERPOL((cip_hrma.ea_1 + cip_hrma.ea_3), $ cip_hrma.energy, spec_es) grating_effic = INTERPOL(cip_meg.(((n_tags(cip_meg))/2)+order), $ cip_meg.energy, spec_es) end 2: begin ; HEG hrma_area = INTERPOL((cip_hrma.ea_4 + cip_hrma.ea_6), $ cip_hrma.energy, spec_es) grating_effic = INTERPOL(cip_heg.(((n_tags(cip_heg))/2)+order), $ cip_heg.energy, spec_es) end endCASE hsi_qe = INTERPOL(cip_hsi.qe, cip_hsi.energy, spec_es) spec_sys_area = hrma_area * grating_effic * hsi_qe phot_spec = (spec_per_kev / (spec_sys_area * hsi_time)) > 0.0 ; Simulate the spectrum xss_sim, this_meas(0) + 1, /NO_BND, /NO_PLOT ; This creates trw_id.spec and trw_id.spec.ps files in the current ; directory. if KEYWORD_SET(LPR) then begin SPAWN, 'cp '+trw_id+'.spec '+save_em_dir+'/'+trw_id+'.spec' end SPAWN, 'rm '+trw_id+'.spec*' if KEYWORD_SET(PS_COMPARE) then begin ; Plot to .ps, e.g. for SPIE paper! ; Parameters for desired aspect ratio SET_PLOT, 'PS' device, /portrait, font_size = 12, XSIZE=17.5, YSIZE=7.0, $ YOFFSET=8.4, XOFFSET=2.0 device,/COLOR dd_load_ct, /OTHER back_color = 1 ; white plot_color = 0 ; black data_color = 0 data_style = 2 cmdb_color = 0 model_color = 0 ; Fit two plots on top of each other... !p.multi = [0,1,2] char_size = 0.7 end else begin ; Plot to .gif set_plot, 'X' window, 0, xsi=500, ysi=700 dd_load_ct, /OTHER back_color = 0 ; black plot_color = 1 ; white data_color = 10 ; red data_style = 0 cmdb_color = 100 ; yellow model_color = 250 ; green !p.multi = [0,1,2] char_size = 1.1 end ; Limit the range for the BOTH plots based on ; i) having data (lam_spec_hist > 2) min_e_meas = MIN(spec_es(where(lam_spec_hist GT 2))) max_e_meas = MAX(spec_es(where(lam_spec_hist GT 2))) ; ii) not being too far from the X-ray line min_e_meas = min_e_meas > 0.6*src_energy ; 0.5*src_energy max_e_meas = max_e_meas < 1.7*src_energy ; 2.*src_energy plot_io, spec_es, phot_spec, PSYM=10, BACK=back_color, $ COLOR=plot_color, $ TITLE = trw_id+': '+source_str+', '+names(grating)+ $ '['+STRCOMPRESS(STRING(order),/REMOVE_ALL)+']', $ XTITLE = 'Energy (keV) '+ $ '[bin size = '+STRING(dlam,FORMAT='(F7.4)')+' A]', $ YTITLE = 'photons / cm^2 s keV ', $ YRANGE=[MIN(phot_spec(where(spec_per_kev GT 0.))), $ MAX(phot_spec(where(spec_es LT 10.0)))], $ XRANGE = [min_e_meas, max_e_meas], XSTYLE=1, $ CHARSIZE = char_size, /NODATA oplot, spec_es, phot_spec, PSYM=10, COLOR=data_color ; Add indication of CMDB assumed energy oplot, src_energy*[1.,1.], [1.E-3,1.E9], LINESTYLE=2, COLOR=cmdb_color ; and model spectrum oplot, lx_es, lx_tubespec, COLOR=model_color ; photons/s cm^2 keV ; indicate the limits oplot, min_e_meas*[1.,1.], [0.001,1.E6],LINESTYLE=1, COLOR=cmdb_color oplot, max_e_meas*[1.,1.], [0.001,1.E6],LINESTYLE=1, COLOR=cmdb_color ; - - - - - - ; Make and Plot normalized cumulative plots plot, spec_es, phot_spec, BACK=back_color, COLOR=plot_color, $ TITLE = 'Normalized cummulative plots ('+ $ date_code+'/hsi'+run_id+'i'+ $ STRCOMPRESS(STRING(iter),/REMOVE_ALL)+'.fits)', $ XTITLE = 'Energy (keV)', YTITLE = 'Cummulative Fraction', $ XRANGE = [min_e_meas, max_e_meas], XSTYLE=1, $ YRANGE=[0.,1.0], $ CHARSIZE = char_size, /NODATA ; Add indication of CMDB assumed energy oplot, src_energy*[1.,1.], [1.E-3,1.E9], LINESTYLE=2, COLOR=cmdb_color ; indicate the limits oplot, min_e_meas*[1.,1.], [0.001,1.E9],LINESTYLE=1, COLOR=cmdb_color oplot, max_e_meas*[1.,1.], [0.001,1.E9],LINESTYLE=1, COLOR=cmdb_color meas_indices = where(spec_es GT min_e_meas AND $ spec_es LT max_e_meas, nwhere) if nwhere GT 2 then begin cum_meas = phot_spec(meas_indices) cum_meas_es = spec_es(meas_indices) cum_meas(0) = cum_meas(0)*ABS(cum_meas_es(1)-cum_meas_es(0)) for ib=1, n_elements(cum_meas)-1 do begin cum_meas(ib) = cum_meas(ib-1) + $ cum_meas(ib)*ABS(cum_meas_es(ib)-cum_meas_es(ib-1)) end cum_meas = 1.0 - cum_meas/MAX(cum_meas) oplot, cum_meas_es, cum_meas, PSYM=10, COLOR=data_color, $ LINESTYLE=data_style end model_indices = where(lx_es GT min_e_meas AND $ lx_es LT max_e_meas, nwhere) if nwhere GT 2 then begin cum_model = lx_tubespec(model_indices) cum_es = lx_es(model_indices) cum_model(0) = cum_model(0)*ABS(cum_es(1)-cum_es(0)) for ib=1, n_elements(cum_model)-1 do begin cum_model(ib) = cum_model(ib-1) + $ cum_model(ib)*ABS(cum_es(ib)-cum_es(ib-1)) end cum_model = cum_model/MAX(cum_model) oplot, cum_es, cum_model, PSYM=10, COLOR=model_color end if KEYWORD_SET(PS_COMPARE) then begin ; finish the .ps plot device, /close set_plot, 'X' SPAWN, 'cp idl.ps '+save_em_dir+'/'+trw_id+'_compare.ps' end if KEYWORD_SET(LPR) then begin ; Finish the .gif output wshow,0,iconic=0 image = tvrd() tvlct, red, green, blue, /GET gif_file = save_em_dir+'/'+trw_id+'_compare.gif' write_gif, gif_file, image, red, green, blue end endIF else begin ; COMPARE phot_spec = 0.0 * spec_per_kev end ;- - - - - - - - - - - - - - ; Output a file! ; RDB format TAB = STRING(9B) OPENW, out_unit, save_em_dir+'/'+trw_id+'_spec.rdb',/GET_LUN lam_fmt = '(F9.4)' cnt_fmt = '(I6)' e_fmt = '(F10.6)' flux_fmt = '(F10.3)' printf, out_unit, '#filename:'+TAB+trw_id+'_spec.rdb' printf, out_unit, '#fileCreationDate:'+TAB+SYSTIME() printf, out_unit, '#swVersion:'+TAB+'hsi_source_spectra.pro' printf, out_unit, '#trw_id:'+TAB+trw_id printf, out_unit, '#date_code:'+TAB+date_code printf, out_unit, '#run_id:'+TAB+run_id printf, out_unit, '#exposure:'+TAB+STRCOMPRESS(STRING(hsi_time),/REMOVE_ALL) printf, out_unit, '#iteration:'+TAB+STRCOMPRESS(STRING(iter),/REMOVE_ALL) printf, out_unit, '#source_string:'+TAB+source_str printf, out_unit, '#energy:'+TAB+STRING(src_energy,FORMAT='(F10.6)') printf, out_unit, '#X_hsi:'+TAB+STRING(Xloc(iter),FORMAT='(F10.4)') printf, out_unit, '#Y_hsi:'+TAB+STRING(Yloc(iter),FORMAT='(F10.4)') printf, out_unit, '#Z_hsi:'+TAB+STRING(Zloc(iter),FORMAT='(F10.4)') printf, out_unit, '#grating:'+TAB+names(grating) printf, out_unit, '#order:'+TAB+STRING(order,FORMAT='(I3)') printf, out_unit, '#spec_width:'+TAB+STRING(spec_width,FORMAT='(F10.6)') printf, out_unit, '#' printf, out_unit, '# Column definitions for humans:' printf, out_unit, '# wavelen wavelength of bin in Angstroms' printf, out_unit, '# (bins are all the same wavelength width)' printf, out_unit, '# counts counts in bin in extraction region' printf, out_unit, '# bkgcounts counts in bin in background regions' printf, out_unit, '# energy energy of the bin in keV' printf, out_unit, '# photspec photons / s cm^2 keV based on area, time, QE' printf, out_unit, '#' printf, out_unit, 'wavelen'+TAB+'counts'+TAB+'bkgcounts'+TAB+'energy'+ $ TAB+'photspec' printf, out_unit, 'N'+TAB+'N'+TAB+'N'+TAB+'N'+TAB+'N' for io=0,n_elements(lams)-1 do begin printf, out_unit, STRING(lams(io),FORMAT=lam_fmt) + TAB + $ STRING(FIX(lam_spec_hist(io)),FORMAT=cnt_fmt)+ TAB + $ STRING(FIX(lam_bkg_hist(io)),FORMAT=cnt_fmt)+ TAB + $ STRING(spec_es(io),FORMAT=e_fmt)+ TAB + $ STRING(phot_spec(io),FORMAT=flux_fmt) end close, out_unit free_lun, out_unit !p.multi = 0 RETURN END