PRO bp_examine, data_set_name, grat, RMFLISTYES=rmflistyes ; ; Look at the bp_extract fits files from ISIS Gaussian fits ; to data and simulations... ; where are the ISIS output files? ; bp_extract is run in here: ;;isis_out_dir = '/nfs/cxc/h3/dd/isis/v0.97/' ; results copied to here: isis_out_dir = '/nfs/spectra/d2/GtoHAK/ISIS_out/' ; ; Select the data set: ; ; SIMULATIONS: ;- ;;data_set_name = 'Capella_wbkg' ;;data_set_name = 'fil_CapellaSim' ;;data_set_name = 'cxc_CapellaSim' ;- ;;data_set_name = 'Narrow33_100ks' ;;data_set_name = 'fil_Narrow33' ;;data_set_name = 'cxc_Narrow33' ; ; REAL DATA: ;- ;;data_set_name = 'obsid_1103' ;;data_set_name = 'fil_1103' ;;data_set_name = 'cxc_1103' ; ; And grating type: ; ;;grat = 'meg' ; ; Get the "bump" count data: ; From ISIS bp_extract.i ... ; ; Minus-side data datm = mrdfits(isis_out_dir+ $ data_set_name + "_bp_"+grat+"m.fits",1) ; ; Plus-side data datp = mrdfits(isis_out_dir+ $ data_set_name + "_bp_"+grat+"p.fits",1) ; Select the valid fits sel = where(datm.ignore EQ 0 and datm.afit GT 0.0) ; X axis range... lam_max = 17.0 ; HEG if grat EQ 'meg' then lam_max = 26.0 ; ; Look at Absolute Wavelengths... ; ppmm = 1.E6 * (datm(sel).lfit - datm(sel).lam)/datm(sel).lam ppmp = 1.E6 * (datp(sel).lfit - datp(sel).lam)/datp(sel).lam ppmave = TOTAL([ppmm,ppmp])/n_elements([ppmm,ppmp]) pre_print_rect ; Compare wavelengths... plot, datm(sel).lam, ppmm, PSYM=4, $ yrange =[-2000.0, 2000.0], XRANGE=[0.0, lam_max], XSTYLE=1, YSTYLE=1, $ XTITLE="Wavelength (A)", $ YTITLE="lambda_fit/lambda_model - 1 (ppm)", $ TITLE=data_set_name+", "+grat+" : Wavelength Differences" oplot, datp(sel).lam, ppmp, PSYM=2 ; plot the overall average of plus-minus deviations... oplot, [0.0, lam_max], [1.0,1.0]*ppmave, LINESTYLE=2 ; and a key oplot, [1.0, 3.0], 1500.0*[1.0,1.0], LINESTYLE=2 xyouts, 3.4, 1500.0 - 50.0, "Ave. is " + string(ppmave,FORMAT='(F7.1)')+' ppm' device,/close spawn,"cp idl.ps "+data_set_name+"_"+grat+"_bp_abslam.ps" set_plot, 'X' ; ; Look at E/dE plots ; edem = datm(sel).lfit/datm(sel).wfit edep = datp(sel).lfit/datp(sel).wfit pre_print_rect plot_oo, !DDHC/datm(sel).lam, edem, PSYM=4, $ yrange =[70.0, 2000.0], XRANGE=[0.3, 10.0], XSTYLE=1, YSTYLE=1, $ XTITLE="Energy (keV)", $ YTITLE="Resolving Power, E/dE", $ TITLE=data_set_name+", "+grat+" : Resolving Power" oplot, !DDHC/datp(sel).lam, edep, PSYM=2 ; Overplot the conservative/optimistic E/dE expectations for ; this grating type ; - - - copied from useful/hist_lines.pro - - - if grat NE 'leg' then begin if STRPOS(!DDLOCATION, 'HAK') LT 0 then begin ; CALDB location of the files cipsdir = !DDHETGCAL+'/cip' end else begin ; HAK Stand-alone location of the files cipsdir = !DDHAKDATA end res_opt = rdb_read(cipsdir+'/hetg'+grat+'D1996-11-01res_optN0002.rdb', $ /SILENT) res_con = rdb_read(cipsdir+'/hetg'+grat+'D1996-11-01res_conN0002.rdb', $ /SILENT) end else begin if STRPOS(!DDLOCATION, 'HAK') LT 0 then begin ; CALDB location of the files cipsdir = !DDLETGCAL+'/cip' end else begin ; HAK Stand-alone location of the files cipsdir = !DDHAKDATA end res_opt = rdb_read(cipsdir+'/letgD1996-11-01res_optN0002.rdb', /SILENT) res_con = rdb_read(cipsdir+'/letgD1996-11-01res_conN0002.rdb', /SILENT) end oplot, res_opt.energies, res_opt.respower, LINESTYLE=2 oplot, res_con.energies, res_con.respower, LINESTYLE=1 ; - - - if KEYWORD_SET(RMFLISTYES) then begin ; Overplot the E/dE from the rmf list file... ; - - - rmf_list_file = '/nfs/bleda/d1/dsd/AXAF/RMF/test/'+grat+'1_list.out' readcol, rmf_list_file, SKIP=1, rmfl1, rmfl2, rmfl3 ; oplot, !DDHC/rmfl1, rmfl1/rmfl3, PSYM=-1 ; - - - end ; Overplot the E/dE from the polynomial Core FWHM approximation ; - - - if grat EQ 'meg' then begin c0 = 0.018851868 c1 = -0.00033946575 c2 = 2.9781558e-05 c3 = -4.1097023e-07 end else begin c0 = 0.010019662 c1 = -0.00014101982 c2 = 1.4591413e-05 c3 = -1.9632994e-07 end fit_lams = lam_max * FINDGEN(101)/100.0 fit_values = c0 + c1*fit_lams + c2*fit_lams^2 + c3*fit_lams^3 ; oplot, !DDHC/fit_lams, fit_lams/fit_values, THICK=3 ; - - - device,/close spawn,"cp idl.ps "+data_set_name+"_"+grat+"_bp_ede.ps" set_plot, 'X' ; ; Look at FWHM vs lambda plots ; fwhm_max = 0.015 if grat EQ 'meg' then fwhm_max = 0.030 pre_print_rect plot, datm(sel).lam, datm(sel).wfit, PSYM=4, $ yrange =[0.0, fwhm_max], XRANGE=[0.0, lam_max], XSTYLE=1, YSTYLE=1, $ XTITLE="Wavelength (A)", $ YTITLE="Fit FWHM (A)", $ TITLE=data_set_name+", "+grat+" : Gaussian FWHM" plot_errors, datm(sel).lam, datm(sel).lam, datm(sel).lam, $ datm(sel).wclo, datm(sel).wchi oplot, datp(sel).lam, datp(sel).wfit, PSYM=2 plot_errors, datp(sel).lam, datp(sel).lam, datp(sel).lam, $ datp(sel).wclo, datp(sel).wchi ; Do this to get the polynomial coefficients used above... if 0 GT 0 then begin ; Fit this with a polynomial fit_lams = [datm(sel).lam, datp(sel).lam] isort = sort(fit_lams) fit_lams = fit_lams(isort) fit_fwhms = [datm(sel).wfit, datp(sel).wfit] fit_fwhms = fit_fwhms(isort) fit_degree = 3 Result = POLY_FIT(fit_lams, fit_fwhms, fit_degree, fit_values) print, ' *** poly_fit output for FWHM of lambda :' print, Result print, ' *** ' end ; Overplot the polynomial FWHM values: oplot, fit_lams, fit_values, THICK=3 if KEYWORD_SET(RMFLISTYES) then begin ; ; Add the rmf list file's FWHM value oplot, rmfl1, rmfl3, PSYM=-1 end ; And the conservative/optimistic resolving power values: oplot, !DDHC/res_opt.energies, $ (!DDHC/res_opt.energies)/res_opt.respower, LINESTYLE=2 oplot, !DDHC/res_con.energies, $ (!DDHC/res_con.energies)/res_con.respower, LINESTYLE=1 device,/close spawn,"cp idl.ps "+data_set_name+"_"+grat+"_bp_fwhm.ps" set_plot, 'X' ; ; Look at per-cent error (pce) vs lambda plots ; ; The product of chi-squared times 100.0/SQRT(number of points) ; gives some idea of the per-cent error between model shape and real shape. (?!) ; pce_max = 40.0 ; chi of the fit converted to a rough per cent error pcem = 100.0 * datm(sel).chifit * (1.0/SQRT(datm(sel).afit)) pcep = 100.0 * datp(sel).chifit * (1.0/SQRT(datp(sel).afit)) pre_print_rect plot, datm(sel).lam, pcem, PSYM=4, $ yrange =[0.0, pce_max], XRANGE=[0.0, lam_max], XSTYLE=1, YSTYLE=1, $ XTITLE="Wavelength (A)", $ YTITLE='100.0 * Chi-squared / SQRT(Area)', $ TITLE=data_set_name+", "+grat+" : Per-cent Error Estimate" oplot, datp(sel).lam, pcep, PSYM=2 device,/close spawn,"cp idl.ps "+data_set_name+"_"+grat+"_bp_pce.ps" set_plot, 'X' RETURN END