;--------------------------------------------------------------------------- ; Document name: spex_read_fit_results.pro ; Created by: Sandhia Bansal ; ; Time-stamp: ;--------------------------------------------------------------------------- ; ;+ ; PROJECT: ; HESSI ; ; NAME: ; SPEX_READ_FIT_RESULTS ; ; PURPOSE: ; Reads the FITS or geny file that contains fit results and copies them to ; a ospex_summ structure. ; ; CATEGORY: ; gen/fits ; ; CONSTRUCTION: ; fit_results = spex_read_fit_results(filename) ; ; SEE ALSO ; HISTORY: ; Sep-2004: created - Sandhia Bansal ; 25-Sep-2004, Kim. If filename not passed in, open dialog_pickfile ; 29-Sep-2004, Kim. For geny files, add any tags that weren't defined when ; file was created, and fill with zeros. For FITS files, this is already taken ; care of since we call spex_summ_define_str which gets current struct def, and ; then fills in what it gets from the FITS file. ; 7-Apr-2006, Kim. Added ncomp to spex_summ_define_str args and added spectrum, ; model to fields to retrieve from file. Also spex_summ_abun_table ; SEE NOTE IN SPEX__FITSWRITE header) ; 20-Dec-2007, Kim. handle files with several '.' by using /last keyword ; 16-Apr-2014, Kim. Call spex_convert_results if version lt 1.2 now (was 1.1) ; 15-May-2014, Kim. Use is_fits to check if FITS format (previously just checked extension='.fits', ; and changed case statement to if statements. ;-------------------------------------------------------------------- function spex_read_fit_results, filename if size(filename, /tname) ne 'STRING' then begin filename = dialog_pickfile (path=curdir(), filter=['*.fits','*.fts','*.geny'], $ file=filename, $ title = 'Select geny or FITS file to restore fit results from', $ group=group, $ get_path=path) endif if filename eq '' then begin message, /info, 'No file selected.' return, 0 endif if not file_exist(filename) then begin message, /info, 'Cannot find file: ' + filename + ' Returning.' return, 0 endif break_file, filename, dum,dum, dum, ext, /last is_geny = ext eq '.geny' is_fits = is_fits(filename) if ~is_geny and ~is_fits then begin message, /info, 'File type not recognized: ' + filename + ' File must be FITS or GENY format. No fit results read.' return, 0 endif print,'Reading fit results from file ', filename if is_geny then begin restgenx, file=filename, fit_file ; any geny file must be Version 1. since don't allow genys anymore. spex_convert_results_geny, fit_file ; ensure that all fields are in structure (in case we've added new information ; to fit_results structure since this file was written) ; spex_summ_define_str will define a structure with all fields as currently ; defined. struct_assign will merge the fields in fit_file into ; fit_results. Any fields in fit_results that weren't in fit_file will be ; left at their default state. nintervals = n_elements( fit_file.spex_summ_ct_rate[0,*] ) nenergies = n_elements( fit_file.spex_summ_ct_rate[*,0] ) nparams = n_elements( fit_file.spex_summ_params[*,0] ) ncomp = fit_function_query (fit_file.spex_summ_fit_function, /ncomp) fit_results = spex_summ_define_str(nintervals, nenergies, nparams, ncomp) struct_assign, fit_file, fit_results, /nozero endif if is_fits then begin ; Read Primary extension ; If file was not written by OSPEX, then consider it invalid. a = mrdfits(filename, 0, primhdr) author = fxpar(primhdr, 'AUTHOR') if (strupcase(trim(author)) ne 'OSPEX') then begin message, 'Invalid OSPEX fit results file: <' + filename + '>..., aborting', /info return, 0 endif ; Read RATE extension extno = get_fits_extno(filename, 'RATE', message=message) if message ne '' then message, message, /info if extno eq -1 then begin message, "Invalid file, aborting. No 'RATE' extension found in file " + filename, /info return, 0 endif rates_table = mrdfits(filename, extno, rates_hdr) ; Read EBOUNDS extension extno = get_fits_extno(filename, 'ENEBAND', message=message) if message ne '' then message, message, /info if extno eq -1 then begin message, "Invalid file, aborting. No 'ENEBAND' extension found in file " + filename, /info return, 0 endif ebounds_table = mrdfits(filename, extno, ebounds_hdr) version = fxpar(primhdr, 'VERSION') if float(version) lt 1.2 then spex_convert_results,primhdr, $ rates_table, rates_hdr, $ ebounds_table, ebounds_hdr ; Create a ospex_summ structure and copy fits results to it ntimes = n_elements(rates_table) nchan = n_elements(ebounds_table) nparams = n_elements(rates_table.params) / ntimes fit_func = trim(fxpar(rates_hdr, 'FITFUNC')) ncomp = fit_function_query (fit_func, /ncomp) fit_results = spex_summ_define_str(ntimes, nchan, nparams, ncomp) fit_results.SPEX_SUMM_AREA = fxpar(rates_hdr, 'GEOAREA') fit_results.SPEX_SUMM_FIT_FUNCTION = fit_func fit_results.SPEX_SUMM_CHISQ = rates_table.chisq fit_results.SPEX_SUMM_CONV = rates_table.convfac fit_results.SPEX_SUMM_CT_RATE = rates_table.rate fit_results.SPEX_SUMM_CT_ERROR = rates_table.stat_err fit_results.SPEX_SUMM_BK_RATE = rates_table.bk_rate fit_results.SPEX_SUMM_BK_ERROR = rates_table.bk_error fit_results.SPEX_SUMM_EMASK = rates_table.emask fit_results.SPEX_SUMM_FILTER = rates_table.filter fit_results.SPEX_SUMM_FIT_DONE = rates_table.fit_done fit_results.SPEX_SUMM_FREE_MASK = rates_table.freemask timezero = fxpar(rates_hdr, 'TIMEZERO') mjdref = fxpar(rates_hdr, 'MJDREF') fit_results.spex_summ_time_interval[0,*] = mjd2any(timezero+mjdref) + rates_table.time - rates_table.timedel/2. fit_results.spex_summ_time_interval[1,*] = fit_results.spex_summ_time_interval[0,*] + rates_table.timedel fit_results.SPEX_SUMM_PH_MODEL = rates_table.phmodel fit_results.SPEX_SUMM_FUNC_SPECTRUM = rates_table.func_spectrum fit_results.SPEX_SUMM_FUNC_MODEL = rates_table.func_model fit_results.SPEX_SUMM_CHIANTI_VERSION = rates_table[0].chianti_version fit_results.SPEX_SUMM_ABUN_TABLE = rates_table[0].abun_table fit_results.SPEX_SUMM_SOURCE_XY = rates_table[0].source_xy fit_results.SPEX_SUMM_SOURCE_ANGLE = rates_table[0].source_angle fit_results.SPEX_SUMM_PARAMS = rates_table.params fit_results.SPEX_SUMM_STARTING_PARAMS = rates_table.startpar fit_results.SPEX_SUMM_MINIMA = rates_table.minima fit_results.SPEX_SUMM_MAXIMA = rates_table.maxima fit_results.SPEX_SUMM_MAXITER = rates_table.maxiter fit_results.SPEX_SUMM_NITER = rates_table.niter fit_results.SPEX_SUMM_RESID = rates_table.residual fit_results.SPEX_SUMM_SIGMAS = rates_table.sigmas fit_results.SPEX_SUMM_UNCERT = rates_table.uncert fit_results.SPEX_SUMM_STOP_MSG = trim(rates_table.stopmsg) fit_results.spex_summ_energy[0,*] = ebounds_table.e_min fit_results.spex_summ_energy[1,*] = ebounds_table.e_max endif return, fit_results end