From 52b7f73b06f173bf5a0af36dc6b0a5ebbb846ace Mon Sep 17 00:00:00 2001 From: paolomassa Date: Mon, 23 Mar 2026 13:17:02 +0100 Subject: [PATCH 01/13] First release of new ELUT correction based on daily calibration fits files and propagation of livetime uncertainty. --- stix/idl/io/stx_read_calibration_data.pro | 99 +++++ stix/idl/processing/calib_data/.DS_Store | Bin 0 -> 6148 bytes .../calib_data/stx_calibration_data.pro | 43 +++ stix/idl/processing/imaging/stx_em.pro | 30 +- .../processing/livetime/stx_cpd_livetime.pro | 59 +++ .../livetime/stx_livetime_fraction.pro | 36 +- .../processing/livetime/stx_triggergram.pro | 53 +-- .../pixel_data/stx_construct_pixel_data.pro | 353 +++++++++++------- .../stx_construct_pixel_data_summed.pro | 9 +- .../pixel_data/stx_sum_pixel_data.pro | 86 +---- .../spectrogram/stx_elut_correction.pro | 272 ++++++++++++++ .../stx_estimate_spectral_index.pro | 119 ++++++ .../vis/stx_calibrate_visibility.pro | 24 +- .../stx_construct_calibrated_visibility.pro | 24 +- .../vis/stx_construct_visibility.pro | 25 +- .../vis/stx_pixel_data_summed2visibility.pro | 31 +- stix/idl/structures/stx_pixel_data.pro | 13 +- stix/idl/structures/stx_pixel_data_summed.pro | 8 +- 18 files changed, 941 insertions(+), 343 deletions(-) create mode 100644 stix/idl/io/stx_read_calibration_data.pro create mode 100644 stix/idl/processing/calib_data/.DS_Store create mode 100644 stix/idl/processing/calib_data/stx_calibration_data.pro create mode 100644 stix/idl/processing/livetime/stx_cpd_livetime.pro mode change 100644 => 100755 stix/idl/processing/livetime/stx_triggergram.pro create mode 100644 stix/idl/processing/spectrogram/stx_elut_correction.pro create mode 100644 stix/idl/processing/spectrogram/stx_estimate_spectral_index.pro diff --git a/stix/idl/io/stx_read_calibration_data.pro b/stix/idl/io/stx_read_calibration_data.pro new file mode 100644 index 00000000..9de17fcf --- /dev/null +++ b/stix/idl/io/stx_read_calibration_data.pro @@ -0,0 +1,99 @@ +;+ +; +; name: +; stx_read_calibration_data +; +; :description: +; Read the values contained in an 'solo_CAL_stix-cal-energy' fits file and returns a corresponding +; 'stx_calibration_data' structure +; +; +; :params: +; fits_path : in, required, type="string" +; the path of the FITS file to be read. Passed through to mrdfits. +; +; +; :keywords: +; +; primary_header : out, type="string array" +; The header of the primary HDU of the pixel data files +; +; data_header : out, type="string array" +; The header of the data extension of the auxiliary file +; +; data_str : out, type="structure" +; The contents of the data extension of the auxiliary file +; +; control_header : out, type="string array" +; The header of the control extension of the auxiliary file +; +; control_str : out, type="structure" +; The contents of the control extension of the auxiliary file +; +; idb_version_header : out, type="string array" +; The header of the idb version extension of the auxiliary file +; +; idb_version_str : out, type="structure" +; The contents of the idb version extension of the auxiliary file +; +; :returns: +; +; a 'stx_aux_data' structure containing the values read from an auxiliary fits file +; +; :examples: +; +; stx_read_aux_fits, fits_path, aux_data=aux_data +; +; :history: +; +; March 2026: Massa P. (FHNW), created +; +;- + +pro stx_read_calibration_data, fits_path, calib_data=calib_data, primary_header = primary_header, $ + data_str = data, data_header = data_header, control_str = control, $ + control_header= control_header, idb_version_str = idb_version, $ + idb_version_header = idb_version_header + + + ;; Read fits file + !null = stx_read_fits(fits_path, 0, primary_header, mversion_full = mversion_full, /silent) + control = stx_read_fits(fits_path, 'control', control_header, mversion_full = mversion_full, /silent) + data = stx_read_fits(fits_path, 'data', data_header, mversion_full = mversion_full, /silent) + idb_version = stx_read_fits(fits_path, 'idb_versions', idb_version_header, mversion_full = mversion_full, /silent) + + ;; Create time object + start_time = sxpar(primary_header, 'date-beg') + stx_time_obj = stx_time() + stx_time_obj.value = anytim(start_time , /mjd) + + time_bin_center = data.TIME ;; in cs + duration = data.TIMEDEL ;; in cs + + t_start = stx_time_add( stx_time_obj, seconds = [ time_bin_center/100 - duration/200 ] ) + t_end = stx_time_add( stx_time_obj, seconds = [ time_bin_center/100 + duration/200 ] ) + t_mean = stx_time_add( stx_time_obj, seconds = [ time_bin_center/100 ] ) + + ;; Define 'calib_data' structure + calib_data = stx_calibration_data() + + calib_data.time_start = t_start + calib_data.time_end = t_end + calib_data.t_mean = t_mean + + ;; Actual energy bins in keV + e_edges_actual = data.E_EDGES_ACTUAL + calib_data.energy_bin_low = reform(e_edges_actual[0:31,*,*]) + calib_data.energy_bin_high = reform(e_edges_actual[1:32,*,*]) + + ;; Gain + calib_data.gain = 1./data.GAIN + ;; Offset + calib_data.offset = data.OFFSET + ;; Live time + calib_data.live_time = data.LIVE_TIME + + ;; ELUT name + calib_data.elut_name = control.OB_ELUT_NAME + +end diff --git a/stix/idl/processing/calib_data/.DS_Store b/stix/idl/processing/calib_data/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5008ddfcf53c02e82d7eee2e57c38e5672ef89f6 GIT binary patch literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0 Array[1]; stx time_axis structure ; Det_select - sub-collimator numbers, 1-32 -; +; ; :Description: ; 28-Apr-13 gh ; @@ -120,20 +120,21 @@ ; :File_comments: ; Uses a default Tau, deadtime per event, of 9.6 microseconds. This may need to ; change based on tests of the Caliste detectors. -; +; ; :Author: ; richard.schwartz@nasa.gov -; +; ; :History: ; 29-april-2013, created ; 18-april-2015, richard.schwartz@nasa.gov, major revision ; 03-dec-2018, ECMD (Graz), change of calculation including eta and tau ; 03-mar-2022, ECMD (Graz), update of eta and tau values and definition tau is now only readout time -; 21-apr-2022, ECMD (Graz), added pileup correction parameter -; 17-oct-2023, ECMD (Graz), update of default eta value using empirical high trigger rate data +; 21-apr-2022, ECMD (Graz), added pileup correction parameter +; 17-oct-2023, ECMD (Graz), update of default eta value using empirical high trigger rate data +; 13-may-2024, P. Massa (FHNW), modified to return statistical error on live time fraction ; ;- -function stx_livetime_fraction, triggergram, det_select, tau_array = tau_array, eta_array=eta_array, error=error +function stx_livetime_fraction, triggergram, det_select, tau_array = tau_array, eta_array=eta_array, error=error error = 1 adg_sc = stx_adg_sc_table() @@ -143,7 +144,7 @@ function stx_livetime_fraction, triggergram, det_select, tau_array = tau_array, default, tau_array, 10.1e-6 + fltarr(ntrig) ;10.1 microseconds readout time per event default, eta_array, 1.10e-6 + fltarr(ntrig) ;1.1 microseconds latency time per event (best fit May 2023) - beta = stx_pileup_corr_parameter() ; get estimate of pileup correction parameter + beta = stx_pileup_corr_parameter() ; get estimate of pileup correction parameter idx_select = ( adg_sc[ where_arr( adg_sc.sc, det_select ) ] ).adg_idx ;these are the agd id needed (1-16) test_triggers = where_arr( triggergram.adg_idx, idx_select, /notequal, test_forzero ) ;which triggers to use @@ -151,14 +152,25 @@ function stx_livetime_fraction, triggergram, det_select, tau_array = tau_array, ;triggergram must be sorted into adg_idx order! ix_fordet = value_locate( triggergram.adg_idx, idx_select ) - ndt = n_elements( triggergram.t_axis.duration ) + ndt = n_elements(triggergram.t_axis.duration ) duration = transpose( rebin( triggergram.t_axis.duration, ndt, ntrig )) tau_rate = rebin( tau_array, ntrig, ndt ) / duration eta_rate = rebin( eta_array, ntrig, ndt ) / duration nin = triggergram.triggerdata / (1. - triggergram.triggerdata *(tau_rate+eta_rate)) + nin_err = 1. / (1. - triggergram.triggerdata *(tau_rate+eta_rate))^2. * triggergram.triggerdata_err + ;triggergram.triggerdata / (1. - triggergram.triggerdata *(tau_rate+eta_rate)) livetime_fraction = exp( -1.*beta*eta_rate*nin) /(1. + (tau_rate+eta_rate)* nin) - result = livetime_fraction[ ix_fordet[sort((where_arr( adg_sc.sc, det_select,/map ))[where_arr( adg_sc.sc, det_select)])], * ] + + ;; Compute statistical error on live time fraction + livetime_fraction_err = exp(-1.*beta*eta_rate*nin) * $ + abs(-1.*beta*eta_rate - 1.*beta*eta_rate * (tau_rate+eta_rate) * nin - (tau_rate+eta_rate)) / $ + (1. + (tau_rate+eta_rate)* nin)^2. * nin_err + livetime_fraction = livetime_fraction[ ix_fordet[sort((where_arr( adg_sc.sc, det_select,/map ))[where_arr( adg_sc.sc, det_select)])], * ] + livetime_fraction_err = livetime_fraction_err[ ix_fordet[sort((where_arr( adg_sc.sc, det_select,/map ))[where_arr( adg_sc.sc, det_select)])], * ] error = 0 + result = {livetime_fraction:livetime_fraction,$ + livetime_fraction_err: livetime_fraction_err} + return, result - + end diff --git a/stix/idl/processing/livetime/stx_triggergram.pro b/stix/idl/processing/livetime/stx_triggergram.pro old mode 100644 new mode 100755 index bbb3d283..db6b7546 --- a/stix/idl/processing/livetime/stx_triggergram.pro +++ b/stix/idl/processing/livetime/stx_triggergram.pro @@ -13,14 +13,14 @@ ; Data exchange structure for triggergrams ; ; CATEGORY: -; +; ; ; CALLING SEQUENCE: ; structure = stx_triggergram(triggerdata, t_axis) ; ; HISTORY: ; 17-apr-2015, based on stx_spectrogram -; +; 13-May-2024, P. Massa (FHNW), modified to return 'triggerdata_err' ;- ;+ @@ -32,9 +32,9 @@ ; triggerdata is the 2D-data array containing triggers over time, either all 16 triggers by time or 1 trigger id by time ; t_axis is the 1D-data array containing the time axis ; adg_idx - Keyword, if passed it is the adg_idx (1-16) of the single trigger accumulator data passed -; +; ;- -function stx_triggergram, triggerdata, t_axis, adg_idx = adg_idx +function stx_triggergram, triggerdata, triggerdata_err, t_axis, adg_idx = adg_idx error = 0 catch, error if (error ne 0)then begin @@ -43,33 +43,42 @@ function stx_triggergram, triggerdata, t_axis, adg_idx = adg_idx message, err return, 0 endif - + t_axis_dim = size(t_axis.duration, /dimension) - if t_axis_dim[0] eq 1 then triggerdata = reform(triggerdata,n_elements(triggerdata),1) - + if t_axis_dim[0] eq 1 then begin + triggerdata = reform(triggerdata,n_elements(triggerdata),1) + triggerdata_err = reform(triggerdata_err,n_elements(triggerdata_err),1) + endif + sizetriggerdata = size(triggerdata, /str) + sizetriggerdata_err = size(triggerdata_err, /str) ; Do some parameter checking if ~(isarray(triggerdata) and sizetriggerdata.dimensions[0] eq 16 or sizetriggerdata.dimensions[0] eq 1) then $ - message, "Parameter 'triggerdata' must be a 16 x M or 1 x M int array" + message, "Parameter 'triggerdata' must be a 16 x M or 1 x M int array" + if ~(isarray(triggerdata_err) and sizetriggerdata_err.dimensions[0] eq 16 or sizetriggerdata_err.dimensions[0] eq 1) then $ + message, "Parameter 'triggerdata' must be a 16 x M or 1 x M int array" if ~(ppl_typeof(t_axis,compareto='stx_time_axis')) then message, "Parameter 't_axis' must be a stx_time_axis structure" - if sizetriggerdata.dimensions[0] eq 16 then adg_idx = indgen(16) + 1 - - + if sizetriggerdata_err.dimensions[0] eq 16 then adg_idx = indgen(16) + 1 + + ; Do some final checking on dimenstions (do they agree) triggerdata_dim = sizetriggerdata.dimensions - - -; if (triggerdata_dim[1] ne t_axis_dim[0] or t_axis_dim[0] eq 1) then $ - if (triggerdata_dim[1] ne t_axis_dim[0]) then $ - message, "'t_axis' dimensions do not agree with 'triggerdata' dimensions" - - + triggerdata_err_dim = sizetriggerdata_err.dimensions + + ; if (triggerdata_dim[1] ne t_axis_dim[0] or t_axis_dim[0] eq 1) then $ + if (triggerdata_dim[1] ne t_axis_dim[0]) then $ + message, "'t_axis' dimensions do not agree with 'triggerdata' dimensions" + if (triggerdata_err_dim[1] ne t_axis_dim[0]) then $ + message, "'t_axis' dimensions do not agree with 'triggerdata_err' dimensions" + + ; If all goes well, put the hsp_triggergram structure together and return it to caller trigstr = { type : 'stx_triggergram', $ - triggerdata : ULONG64(triggerdata), $ - adg_idx : adg_idx, $ - t_axis : t_axis } ; $ + triggerdata : ULONG64(triggerdata), $ + triggerdata_err : DOUBLE(triggerdata_err), $ + adg_idx : adg_idx, $ + t_axis : t_axis } ; $ + - return, trigstr end \ No newline at end of file diff --git a/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro b/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro index 44d0871a..53b5be54 100644 --- a/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro +++ b/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro @@ -29,27 +29,25 @@ ; - COUNTS_ERROR: array 32x12 containing the errors (statistics + compression) associated with the number of counts ; recorded by the detector pixels ; - LIVE_TIME: 32-element array containing the live time of each detector in the considered time interval +; - LIVE_TIME_ERROR: 32-element array containing the uncertainty on the live time of each detector in the considered time interval ; - TIME_RANGE: two-element 'stx_time' array containing the lower and upper edge of the selected time interval ; (the time bins containing the start and the end time provided as input are included in the ; selected interval) ; - ENERGY_RANGE: array containing the lower and the upper edge of the selected energy interval ; (the energy bins containing the lower and the upper edges provided as input are included in the interval) -; - COUNTS_BKG: array 32x12 containing an estimate of the number of background counts recorded by the detector pixels in the selected -; time and energy intervals. If no bakground measurement is provided, it is filled with zeros -; - COUNTS_ERROR_BKG: array 32x12 containing the errors (statistics + compression) associated with the estimate of the -; number of background counts recorded by the detector pixels in the selected time and energy intervals -; - LIVE_TIME_BKG: 32-element array containing the live time of each detector for the background measurement ; - RCR: rate control regime status in the selcted time interval. If the RCR changes in that interval, an error is thrown ; - PIXEL_MASKS: matrix containing info on the pixels used in the selected time interval ; - DETECTOR_MASKS: matrix containing information on the detectors used in the selected time interval +; - TOT_COUNTS: total number of counts recorded by the imaging subcollimators selected by means of 'subc_index' +; - TOT_COUNTS_BKG: estimate of the total number of background counts recorded by the imaging subcollimators (selected with 'subc_index') +; in the considered time and energy interval ; ; KEYWORDS: ; ; path_bkg_file: path of a background L1 fits file. If provided, the fields 'COUNTS_BKG', 'COUNTS_ERROR_BKG' and 'LIVE_TIME_BKG' ; of the pixel_data structure are filled with the values read from the background measurement file ; -; elut_corr: if set, a correction based on a ELUT table is applied to the measured counts -; +; calib_data: if a 'stx_calibration_data' structure is passed as input, then the ELUT correction is applied ; ; subc_index: array containing the indices of the selected imaging detectors. Used only for plotting the lightcurve by means of ; 'stx_plot_selected_time_range'. Default, indices of @@ -74,17 +72,17 @@ ; October 2023, Massa P., fixed bug in the selection of the energy bin indices ; November 2023, Massa P., use simplified version of the subcollimator transmission (temporary solution) ; January 2026, Massa P., removed 'xy_flare' keyword and grid transmission correction +; March 2026, Massa P., made it compatible with new ELUT correction based on daily calibration data ; ; CONTACT: ; paolo.massa@fhnw.ch ;- -function stx_construct_pixel_data, path_sci_file, time_range, energy_range, elut_corr=elut_corr, $ +function stx_construct_pixel_data, path_sci_file, time_range, energy_range, calib_data=calib_data, $ path_bkg_file=path_bkg_file, subc_index=subc_index, $ sumcase=sumcase, silent=silent, no_small=no_small, no_rcr_check=no_rcr_check, $ shift_duration=shift_duration, _extra=extra -default, elut_corr, 1 default, silent, 0 default, subc_index, stx_label2ind(['10a','10b','10c','9a','9b','9c','8a','8b','8c','7a','7b','7c',$ '6a','6b','6c','5a','5b','5c','4a','4b','4c','3a','3b','3c']) @@ -143,16 +141,9 @@ if (energy_range[0] lt min(energy_low)) then $ if (energy_range[1] gt max(energy_high)) then $ message, 'The upper edge of the selected energy interval is outside the science energy interval (' + $ num2str(fix(min(energy_low))) + ' - ' + num2str(fix(max(energy_high))) + ' keV)' - -energy_ind_min = where(energy_low le energy_range[0]) -energy_ind_min = energy_ind_min[-1] -energy_ind_max = where(energy_high ge energy_range[1]) -energy_ind_max = energy_ind_max[0] - -energy_ind = [energy_ind_min:energy_ind_max] -if n_elements(energy_ind) eq 1 then energy_ind = energy_ind[0] -this_energy_range = [energy_low[energy_ind_min], energy_high[energy_ind_max]] +energy_min = min(energy_low) +energy_max = max(energy_high) if keyword_set(path_bkg_file) then begin @@ -162,49 +153,81 @@ if keyword_set(path_bkg_file) then begin energy_low_bkg = e_axis_bkg.LOW energy_high_bkg = e_axis_bkg.HIGH - + if (energy_range[0] lt min(energy_low_bkg)) then $ message, 'The lower edge of the selected energy interval is outside the background energy interval (' + $ num2str(fix(min(energy_low_bkg))) + ' - ' + num2str(fix(max(energy_high_bkg))) + ' keV)' - + if (energy_range[1] gt max(energy_high_bkg)) then $ message, 'The upper edge of the selected energy interval is outside the background energy interval (' + $ num2str(fix(min(energy_low_bkg))) + ' - ' + num2str(fix(max(energy_high_bkg))) + ' keV)' - + + ;; Extract energy range in common between science and background file + energy_min = max([energy_min,min(energy_low_bkg)]) + energy_max = min([energy_max,max(energy_high_bkg)]) + idx_energy_bkg = where((energy_low_bkg ge energy_min) and (energy_high_bkg le energy_max)) + + energy_bin_idx_bkg = energy_bin_idx_bkg[idx_energy_bkg] + energy_low_bkg = energy_low_bkg[idx_energy_bkg] + energy_high_bkg = energy_high_bkg[idx_energy_bkg] + energy_ind_min_bkg = where(energy_low_bkg le energy_range[0]) energy_ind_min_bkg = energy_ind_min_bkg[-1] energy_ind_max_bkg = where(energy_high_bkg ge energy_range[1]) energy_ind_max_bkg = energy_ind_max_bkg[0] - + energy_ind_bkg = [energy_ind_min_bkg:energy_ind_max_bkg] if n_elements(energy_ind_bkg) eq 1 then energy_ind_bkg = energy_ind_bkg[0] - + this_energy_range_bkg = [energy_low_bkg[energy_ind_min_bkg], energy_high_bkg[energy_ind_max_bkg]] endif +;; Extract energy range in common between science and background file +;; (in the case a background measurement file is passed as input) +idx_energy = where((energy_low ge energy_min) and (energy_high le energy_max)) + +energy_bin_idx = energy_bin_idx[idx_energy] +energy_low = energy_low[idx_energy] +energy_high = energy_high[idx_energy] + +energy_ind_min = where(energy_low le energy_range[0]) +energy_ind_min = energy_ind_min[-1] +energy_ind_max = where(energy_high ge energy_range[1]) +energy_ind_max = energy_ind_max[0] + +energy_ind = [energy_ind_min:energy_ind_max] +if n_elements(energy_ind) eq 1 then energy_ind = energy_ind[0] + +this_energy_range = [energy_low[energy_ind_min], energy_high[energy_ind_max]] + ;;************** Compute livetime -triggergram = stx_triggergram(data.TRIGGERS, t_axis) -livetime_fraction = stx_livetime_fraction(triggergram) -duration_time_bins = t_axis.DURATION -duration_time_bins = transpose(cmreplicate(duration_time_bins, 32)) -live_time_bins = livetime_fraction*duration_time_bins -live_time = n_elements(time_ind) eq 1? reform(live_time_bins[*,time_ind]) : $ - total(live_time_bins[*,time_ind],2) +live_time_data = stx_cpd_livetime(data.TRIGGERS, data.TRIGGERS_ERR, t_axis) +live_time_bins = live_time_data.LIVE_TIME_BINS +live_time_bins_err = live_time_data.LIVE_TIME_BINS_ERR + +live_time = n_elements(time_ind) eq 1? reform(live_time_bins[*,time_ind]) : total(live_time_bins[*,time_ind],2) +live_time_error = n_elements(time_ind) eq 1? reform(live_time_bins_err[*,time_ind]) : $ + sqrt(total(live_time_bins_err[*,time_ind]^2.,2)) if keyword_set(path_bkg_file) then begin - - triggergram_bkg = stx_triggergram(data_bkg.TRIGGERS, t_axis_bkg) - livetime_fraction_bkg = stx_livetime_fraction(triggergram_bkg) - duration_time_bins_bkg = t_axis_bkg.DURATION - live_time_bkg = duration_time_bins_bkg[0]*livetime_fraction_bkg - -endif + + live_time_bkg_data = stx_cpd_livetime(data_bkg.TRIGGERS, data_bkg.TRIGGERS_ERR, t_axis_bkg) + live_time_bkg = live_time_bkg_data.LIVE_TIME_BINS + live_time_error_bkg = live_time_bkg_data.LIVE_TIME_BINS_ERR + +endif else begin + + live_time_bkg = dblarr(32) + 1. + live_time_error_bkg = dblarr(32) + +endelse + ;;************** Define count matrix and bkg count matrix -;; WE ASSUME THAT THE SCIENCE DATA FILE AND THE BKG DATA FILE CONTAIN MORE THAN 1 ENERGY BIN +;; WE ASSUME THAT THE SCIENCE DATA FILE AND THE BKG DATA FILE CONTAIN MORE THAN 1 ENERGY BIN ;; (I.E., THE NUMBER OF ELEMENTS IN energy_bin_idx AND IN energy_bin_idx_bkg IS ASSUMED TO BE ;; LARGER THAN 1) @@ -223,24 +246,43 @@ counts_error = counts_error[energy_bin_idx,*,*,*] if keyword_set(path_bkg_file) then begin + ;; Check if science and background files are reconrded with the same ELUT + + elut_filename = stx_date2elut_file(stx_time2any(this_time_range[0])) + stx_read_elut, elut_gain, elut_offset, adc4096_str, elut_filename = elut_filename, ekev_actual = ekev_actual_elut + + elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) + stx_read_elut, ekev_actual = ekev_actual_bkg, elut_filename = elut_filename_bkg + + ;; Compare ELUT tables + elut_comp = STRCMP(elut_filename, elut_filename_bkg) + + if not elut_comp then $ + message, 'The background file must be recorded when the same ELUT as the science file was uploaded. Please choose a different background file that is closer in time to the science file.' + counts_bkg = data_bkg.COUNTS counts_error_bkg = data_bkg.COUNTS_ERR counts_bkg = counts_bkg[energy_bin_idx_bkg,*,*] counts_error_bkg = counts_error_bkg[energy_bin_idx_bkg,*,*] -endif +endif else begin + + counts_bkg = dblarr(size(counts, /dim)) + counts_error_bkg = dblarr(size(counts, /dim)) + +endelse ;;************** Plot lightcurve (if ~silent) if ~silent and (n_elements(size(live_time_bins, /dim)) gt 1) then begin - + if keyword_set(path_bkg_file) then begin stx_plot_selected_time_range, stx_time2any(t_axis.MEAN), energy_ind, time_ind, counts, live_time_bins, subc_index, $ - sumcase, this_energy_range, this_time_range, counts_bkg=counts_bkg, $ - live_time_bkg=live_time_bkg, energy_ind_bkg=energy_ind_bkg + sumcase, this_energy_range, this_time_range, counts_bkg=counts_bkg, $ + live_time_bkg=live_time_bkg, energy_ind_bkg=energy_ind_bkg endif else begin stx_plot_selected_time_range, stx_time2any(t_axis.MEAN), energy_ind, time_ind, counts, live_time_bins, subc_index, $ - sumcase, this_energy_range, this_time_range + sumcase, this_energy_range, this_time_range endelse endif @@ -259,111 +301,144 @@ endif else begin endelse -;;************** Sum counts (and related errors) in energy +;;************** Background subtraction -;; Compute elut correction (if 'elut_corr' is set) - Correct just the first and the last energy bins (flat spectrum is assumed) -if elut_corr then begin - - elut_filename = stx_date2elut_file(stx_time2any(this_time_range[0])) - stx_read_elut, ekev_actual = ekev_actual, elut_filename = elut_filename - - ;; We assume that the first energy bin starts from 0. The right end of the last energy bin is set to NaN - ;; We assume that the energy interval is always partitioned into 32 energy bins - energy_bin_low = fltarr(32,12,32) - energy_bin_low[1:31,*,*] = ekev_actual - - energy_bin_high = fltarr(32,12,32) - energy_bin_high[0:30,*,*] = ekev_actual - energy_bin_high[31,*,*] = !VALUES.F_NaN - - energy_bin_low = energy_bin_low[energy_bin_idx,*,*] - energy_bin_high = energy_bin_high[energy_bin_idx,*,*] - - energy_bin_size = energy_bin_high - energy_bin_low - - if keyword_set(path_bkg_file) then begin - - ;; The bkg elut is potentially different from the one used for correcting the science data - elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) - stx_read_elut, ekev_actual = ekev_actual_bkg, elut_filename = elut_filename_bkg - - ;; We assume that the first energy bin starts from 0. The right end of the last energy bin is set to NaN - ;; We assume that the energy interval is always partitioned into 32 energy bins - energy_bin_low_bkg = fltarr(32,12,32) - energy_bin_low_bkg[1:31,*,*] = ekev_actual_bkg - - energy_bin_high_bkg = fltarr(32,12,32) - energy_bin_high_bkg[0:30,*,*] = ekev_actual_bkg - energy_bin_high_bkg[31,*,*] = !VALUES.F_NaN - - energy_bin_low_bkg = energy_bin_low_bkg[energy_bin_idx_bkg,*,*] - energy_bin_high_bkg = energy_bin_high_bkg[energy_bin_idx_bkg,*,*] - - energy_bin_size_bkg = energy_bin_high_bkg - energy_bin_low_bkg - - endif +live_time_rep = transpose(cmreplicate(live_time, [n_elements(energy_bin_idx),12]), [1,2,0]) +live_time_error_rep = transpose(cmreplicate(live_time_error, [n_elements(energy_bin_idx),12]), [1,2,0]) +live_time_bkg_rep = keyword_set(path_bkg_file)? transpose(cmreplicate(live_time_bkg, [n_elements(energy_bin_idx_bkg),12]), [1,2,0]) : fltarr(size(counts_bkg, /dim)) + 1. +live_time_error_bkg_rep = keyword_set(path_bkg_file)? transpose(cmreplicate(live_time_error_bkg, [n_elements(energy_bin_idx_bkg),12]), [1,2,0]) : fltarr(size(counts_bkg, /dim)) + +counts_bkg_estimate = f_div( live_time_rep * counts_bkg, live_time_bkg_rep ) +error_numerator = abs(live_time_rep * counts_bkg) * sqrt( f_div(live_time_error_rep,live_time_rep)^2. + $ + f_div(counts_error_bkg,counts_bkg)^2.) +counts_bkg_estimate_error = abs(counts_bkg_estimate) * sqrt( f_div(error_numerator,live_time_rep * counts_bkg)^2. + $ + f_div(live_time_error_bkg_rep,live_time_bkg_rep)^2.) + +if keyword_set(calib_data) and ~silent then begin + + ;; To be used for plot of the bkg subtracted spectrum + spectrum_with_bkg = total(total(counts[*,0:7,subc_index], 3), 2) / (energy_high - energy_low) + spectrum_bkg = keyword_set(path_bkg_file)? total(total(counts_bkg_estimate[*,0:7,subc_index], 3), 2) / (energy_high_bkg - energy_low_bkg) : total(total(counts_bkg[*,0:7,subc_index], 3), 2) endif -if n_elements(energy_ind) eq 1 then begin - - energy_corr_factor = elut_corr ? (this_energy_range[1] - this_energy_range[0]) / reform(energy_bin_size[energy_ind,*,*]) : $ - dblarr(12,32)+1 +;; Determine indices of the selected pixels. It will be used for estimating total number of BKG counts and for estimating the amount of ELUT correction +case sumcase of + + 'TOP': begin + pixel_ind = [0] + end + + 'BOT': begin + pixel_ind = [1] + end + + 'TOP+BOT': begin + pixel_ind = [0,1] + end + + 'ALL': begin + pixel_ind = [0,1,2] + end + + 'SMALL': begin + pixel_ind = [2] + end +end + +;; Compute total number of counts (to be used for comparison with bkg counts) +counts_reshaped = reform(counts,n_elements(energy_bin_idx), 4, 3, 32) +tot_counts = total(counts_reshaped[energy_ind,*,pixel_ind,subc_index]) + +;; Apply BKG subtraction +counts = counts - counts_bkg_estimate +counts_error = sqrt(counts_error^2. + counts_bkg_estimate_error^2.) + +if keyword_set(path_bkg_file) then begin + + ;; Compute total number of background counts. Select only imaging detectors + counts_bkg = reform(counts_bkg_estimate, n_elements(energy_bin_idx_bkg), 4, 3, 32) + tot_counts_bkg = total(counts_bkg[energy_ind_bkg,*,pixel_ind,subc_index]) - counts = reform(counts[energy_ind,*,*]) * energy_corr_factor - counts_error = reform(counts_error[energy_ind,*,*]) * energy_corr_factor - - endif else begin - - energy_corr_factor_low = elut_corr ? (reform(energy_bin_high[energy_ind[0],*,*]) - this_energy_range[0]) / $ - reform(energy_bin_size[energy_ind[0],*,*]) : dblarr(12,32)+1 - - energy_corr_factor_high = elut_corr ? (this_energy_range[1] - reform(energy_bin_low[energy_ind[-1],*,*])) / $ - reform(energy_bin_size[energy_ind[-1],*,*]) : dblarr(12,32)+1 - - counts[energy_ind[0],*,*] *= energy_corr_factor_low - counts[energy_ind[-1],*,*] *= energy_corr_factor_high - counts_error[energy_ind[0],*,*] *= energy_corr_factor_low - counts_error[energy_ind[-1],*,*] *= energy_corr_factor_high - counts = total(counts[energy_ind,*,*],1) - counts_error = sqrt(total(counts_error[energy_ind,*,*]^2.,1)) - + tot_counts_bkg = 0. + endelse +;;************** Print total number of counts -if keyword_set(path_bkg_file) then begin - - if n_elements(energy_ind_bkg) eq 1 then begin +if ~silent then begin - ;; Background - energy_corr_factor_bkg = elut_corr ? (this_energy_range[1] - this_energy_range[0]) / $ - reform(energy_bin_size_bkg[energy_ind_bkg,*,*]) : dblarr(12,32)+1 - - counts_bkg = reform(counts_bkg[energy_ind_bkg,*,*]) * energy_corr_factor_bkg - counts_error_bkg = reform(counts_error_bkg[energy_ind_bkg,*,*]) * energy_corr_factor_bkg - - endif else begin + print + print + print,'***********************************************************************' + print,'Total number of counts in image: '+strtrim(tot_counts) + print,'Background counts: '+strtrim(tot_counts_bkg) + print,'Counts above background: '+strtrim(tot_counts-tot_counts_bkg) + print,'Total to background: '+strtrim(tot_counts/tot_counts_bkg) + print,'***********************************************************************' + print + print + +endif + +;;************** Sum counts (and related errors) in energy + +;; Compute elut correction (if 'calib_data' is provided) - Correct just the first and the last energy bins +if keyword_set(calib_data) then begin - ;; Background - energy_corr_factor_low_bkg = elut_corr ? (reform(energy_bin_high_bkg[energy_ind_bkg[0],*,*]) - this_energy_range[0]) / $ - reform(energy_bin_size_bkg[energy_ind_bkg[0],*,*]) : dblarr(12,32)+1 + energy_bin_low = calib_data.ENERGY_BIN_LOW + energy_bin_high = calib_data.ENERGY_BIN_HIGH - energy_corr_factor_high_bkg = elut_corr ? (this_energy_range[1] - reform(energy_bin_low_bkg[energy_ind_bkg[-1],*,*])) / $ - reform(energy_bin_size_bkg[energy_ind_bkg[-1],*,*]) : dblarr(12,32)+1 + energy_bin_low = energy_bin_low[energy_bin_idx,*,*] + energy_bin_high = energy_bin_high[energy_bin_idx,*,*] - counts_bkg[energy_ind_bkg[0],*,*] *= energy_corr_factor_low_bkg - counts_bkg[energy_ind_bkg[-1],*,*] *= energy_corr_factor_high_bkg - counts_error_bkg[energy_ind_bkg[0],*,*] *= energy_corr_factor_low_bkg - counts_error_bkg[energy_ind_bkg[-1],*,*] *= energy_corr_factor_high_bkg + ;; Apply ELUT correction + spectrum = total(total(counts[*,0:7,subc_index], 3), 2) / (energy_high - energy_low) + elut_data = stx_elut_correction(counts, counts_error, $ + energy_bin_idx, energy_bin_low, energy_bin_high, energy_high, energy_low, energy_ind, this_energy_range, $ + spectrum, pixel_ind, subc_index, $ + spectrum_with_bkg=spectrum_with_bkg, spectrum_bkg=spectrum_bkg,silent=silent) + + counts = elut_data.COUNTS + counts_error = elut_data.COUNTS_ERROR + counts_no_elut = elut_data.COUNTS_NO_ELUT + counts_elut = elut_data.COUNTS_ELUT + + if ~silent then begin + + ;; Print minimun and maximum ELUT correction in the different pixels + elut_corr_perc = f_div(abs(counts_no_elut - counts_elut),counts_no_elut) * 100 ;; in percentage + counts_error_reshaped = reform(counts_error, 4, 3, 32) + counts_error_elut=reform(counts_error_reshaped[*,pixel_ind,subc_index]) + rel_error=abs(counts_error_elut/counts_elut)*100. + print + print + print,'***********************************************************************' + print,'Min/Max ELUT correction over all pixels: '+num2str(min(elut_corr_perc), format='(f7.2)')+'% - '+num2str(max(elut_corr_perc), format='(f7.2)')+'% ' + print,'Average ELUT correction over all pixels: '+num2str(average(elut_corr_perc), format='(f7.2)')+'%' + print,'Standard deviation of ELUT correction over all pixels: '+num2str(stdev(elut_corr_perc), format='(f7.2)')+'%' + print,'Average count error over all pixels: '+num2str(average(rel_error), format='(f7.2)')+'%' + print,'***********************************************************************' + print + print + + endif +endif else begin - counts_bkg = total(counts_bkg[energy_ind_bkg,*,*],1) - counts_error_bkg = sqrt(total(counts_error_bkg[energy_ind_bkg,*,*]^2.,1)) + if n_elements(energy_ind) eq 1 then begin + + counts = reform(counts[energy_ind,*,*]) + counts_error = reform(counts_error[energy_ind,*,*]) + + endif else begin + + counts = total(counts[energy_ind,*,*],1) + counts_error = sqrt(total(counts_error[energy_ind,*,*]^2.,1)) - endelse - -endif + endelse +endelse ;;************** RCR @@ -411,15 +486,11 @@ pixel_data = stx_pixel_data() pixel_data.TIME_RANGE = this_time_range pixel_data.ENERGY_RANGE = this_energy_range pixel_data.LIVE_TIME = live_time +pixel_data.LIVE_TIME_ERROR = live_time_error pixel_data.COUNTS = transpose(counts) pixel_data.COUNTS_ERROR = transpose(counts_error) -if keyword_set(path_bkg_file) then begin - - pixel_data.LIVE_TIME_BKG = live_time_bkg - pixel_data.COUNTS_BKG = transpose(counts_bkg) - pixel_data.COUNTS_ERROR_BKG = transpose(counts_error_bkg) - -endif +pixel_data.TOT_COUNTS = tot_counts +if keyword_set(path_bkg_file) then pixel_data.TOT_COUNTS_BKG = tot_counts_bkg pixel_data.RCR = rcr[0] @@ -429,12 +500,12 @@ if alpha then begin endif else begin pixel_data.PIXEL_MASKS = pixel_masks[*,0] endelse - + pixel_data.DETECTOR_MASKS = detector_masks[*,0] if ~silent then stx_plot_moire_pattern, pixel_data, no_small=no_small -return, pixel_data +return, pixel_data end diff --git a/stix/idl/processing/pixel_data/stx_construct_pixel_data_summed.pro b/stix/idl/processing/pixel_data/stx_construct_pixel_data_summed.pro index 99b488e6..281d9f71 100644 --- a/stix/idl/processing/pixel_data/stx_construct_pixel_data_summed.pro +++ b/stix/idl/processing/pixel_data/stx_construct_pixel_data_summed.pro @@ -30,9 +30,7 @@ ; ; path_bkg_file: path of a background L1 fits file. If provided, the fields 'COUNT_RATES_BKG', 'COUNT_RATES_ERROR_BKG' ; and 'LIVE_TIME_BKG' of the output 'stx_pixel_data_summed' structure are filled with the values read from the -; background measurement file -; -; elut_corr: if set, a correction based on a ELUT table is applied to the measured counts +; background measurement file ; ; subc_index: array containing the indices of the selected imginging detectors. Used only for plotting the lightcurve by means of ; 'stx_plot_selected_time_range' and for computing the total number of counts in the image. Default, indices of @@ -50,18 +48,19 @@ ; HISTORY: July 2022, Massa P., created ; January 2026, Massa P., removed 'xy_flare' keyword. Grid transmission correction is not performed ; at this stage +; March 2026, Massa P., removed 'elut_corr' keyword as not necessary ; ; CONTACT: ; paolo.massa@fhnw.ch ;- function stx_construct_pixel_data_summed, path_sci_file, time_range, energy_range, path_bkg_file=path_bkg_file, $ - elut_corr=elut_corr, subc_index=subc_index, sumcase=sumcase, silent=silent, no_small=no_small,$ + subc_index=subc_index, sumcase=sumcase, silent=silent, no_small=no_small,$ no_rcr_check=no_rcr_check, _extra=extra ;;************** Construct pixel data -pixel_data = stx_construct_pixel_data(path_sci_file, time_range, energy_range, elut_corr=elut_corr, $ +pixel_data = stx_construct_pixel_data(path_sci_file, time_range, energy_range, $ path_bkg_file=path_bkg_file, subc_index=subc_index, sumcase=sumcase, $ silent=silent, no_small=no_small, no_rcr_check=no_rcr_check, _extra=extra) ;;************** Sum pixel data diff --git a/stix/idl/processing/pixel_data/stx_sum_pixel_data.pro b/stix/idl/processing/pixel_data/stx_sum_pixel_data.pro index cdd27a1f..6c1689d4 100644 --- a/stix/idl/processing/pixel_data/stx_sum_pixel_data.pro +++ b/stix/idl/processing/pixel_data/stx_sum_pixel_data.pro @@ -46,13 +46,7 @@ ; Optionally, they are corrected for the grid internal shadowing and transmission ; - COUNTS_RATES_ERROR: array 32x4 containing the errors associated with the countrates A,B,C,D recorded by each detector ; (compression errors and statistical errors are taken into account, no systematic errors are added) -; - TOT_COUNTS: total number of counts recorded by the imaging subcollimators selected by means of 'subc_index' -; - LIVE_TIME_BKG: 32-element array containing the live time of each detector during the background measurement -; - COUNT_RATES_BKG: array 32x4 containing the background countrates A,B,C,D recorded by each subcollimator. Pixel measurements -; are summed (from 12 to 4) and countrates are normalized by live time, incident area, length of the considered energy -; interval. Optionally, they are corrected for the grid internal shadowing and transmission -; - COUNT_RATES_ERROR_BKG: array 32x4 containing the errors associated with the background countrates A,B,C,D recorded by each detector -; (compression errors and statistical errors are taken into account, no systematic errors are added) +; - TOT_COUNTS: total number of counts recorded by the imaging subcollimators selected by means of 'subc_index' ; - TOT_COUNTS_BKG: estimate of the total number of background counts recorded by the imaging subcollimators (selected with 'subc_index') ; in the considered time and energy interval ; - RCR: Rate Control Regime status during the flare measurement @@ -63,9 +57,10 @@ ; HISTORY: August 2022, Massa P., created ; October 2023, Massa P., fixed bug in the estimation of the total number of bkg counts ; January 2026, Massa P., removed 'xy_flare' entry as grid transmission correction is not applied anymore to the raw counts +; March 2026, Massa P., made it compatible with new ELUT correction ; ; CONTACT: -; paolo.massa@wku.edu +; paolo.massa@fhnw.ch ;- function stx_sum_pixel_data, pixel_data, subc_index=subc_index, sumcase=sumcase, silent=silent @@ -105,66 +100,23 @@ pixel_masks = reform(pixel_data.PIXEL_MASKS,4,3) if total(pixel_masks[*,pixel_ind]) lt 4.*n_elements(pixel_ind) then message, "Change 'sumcase': one of the selected pixels is not available" -count_rates = reform(pixel_data.COUNTS, 32, 4, 3) +counts = reform(pixel_data.COUNTS, 32, 4, 3) ;; Compute total counts: saved in the pixel data structure -tot_counts = total(count_rates[subc_index,*,pixel_ind]) - -count_rates = n_elements( pixel_ind ) eq 1 ? reform(count_rates[*, *, pixel_ind]) : $ - total( count_rates[*, *, pixel_ind], 3 ) +tot_counts = total(counts[subc_index,*,pixel_ind]) -counts_rates_error = reform(pixel_data.COUNTS_ERROR, 32, 4, 3) -counts_rates_error = n_elements( pixel_ind ) eq 1 ? reform(counts_rates_error[*, *, pixel_ind]) : $ - sqrt(total( counts_rates_error[*, *, pixel_ind]^2, 3 )) +counts = n_elements( pixel_ind ) eq 1 ? reform(counts[*, *, pixel_ind]) : $ + total( counts[*, *, pixel_ind], 3 ) -count_rates_bkg = reform(pixel_data.COUNTS_BKG, 32, 4, 3) -;; Compute total background counts: saved in the pixel data structure -tot_counts_bkg = n_elements(pixel_ind) eq 1 ? total(reform(count_rates_bkg[subc_index,*,pixel_ind]),2) : $ - total(total(count_rates_bkg[subc_index,*,pixel_ind],2),2) - -count_rates_bkg = n_elements( pixel_ind ) eq 1 ? reform(count_rates_bkg[*, *, pixel_ind]) : $ - total( count_rates_bkg[*, *, pixel_ind], 3 ) - -count_rates_error_bkg = reform(pixel_data.COUNTS_ERROR_BKG, 32, 4, 3) -count_rates_error_bkg = n_elements( pixel_ind ) eq 1 ? reform(count_rates_error_bkg[*, *, pixel_ind]) : $ - sqrt(total( count_rates_error_bkg[*, *, pixel_ind]^2, 3 )) +counts_error = reform(pixel_data.COUNTS_ERROR, 32, 4, 3) +counts_error = n_elements( pixel_ind ) eq 1 ? reform(counts_error[*, *, pixel_ind]) : $ + sqrt(total( counts_error[*, *, pixel_ind]^2, 3 )) ;;************** Livetime correction: units are counts s^-1 live_time = cmreplicate(pixel_data.LIVE_TIME, 4) -count_rates = f_div(count_rates,live_time) -counts_rates_error = f_div(counts_rates_error,live_time) - -;; Compute live time fraction of each detector. It is used for computing the total number of bkg counts -time_range = stx_time2any(pixel_data.TIME_RANGE) -live_time_fraction = pixel_data.LIVE_TIME/(time_range[1]-time_range[0]) - -live_time_bkg = pixel_data.LIVE_TIME_BKG -;; Estimate of the background counts in the image: only a fraction proportional to the time range -;; Multiply total number of bkg counts by the live time fraction of the science data. -;; In this way, we keep into account that the live time fraction can be different between science and bkg data -tot_counts_bkg = total(f_div(tot_counts_bkg*live_time_fraction[subc_index],live_time_bkg[subc_index])) $ - * (time_range[1]-time_range[0]) - -live_time_bkg = cmreplicate(pixel_data.LIVE_TIME_BKG, 4) -count_rates_bkg = f_div(count_rates_bkg,live_time_bkg) -count_rates_error_bkg = f_div(count_rates_error_bkg,live_time_bkg) - -;;************** Print total number of counts - -if ~silent then begin - - print - print - print,'***********************************************************************' - print,'Total number of counts in image: '+strtrim(tot_counts) - print,'Background counts: '+strtrim(tot_counts_bkg) - print,'Counts above background: '+strtrim(tot_counts-tot_counts_bkg) - print,'Total to background: '+strtrim(tot_counts/tot_counts_bkg) - print,'***********************************************************************' - print - print - -endif +live_time_error = cmreplicate(pixel_data.LIVE_TIME_ERROR, 4) +count_rates = f_div(counts,live_time) +counts_rates_error = abs(count_rates) * sqrt( f_div(counts_error,counts)^2. + f_div(live_time_error,live_time)^2. ) ;;************** Normalization per keV: units are counts s^-1 keV^-1 @@ -172,9 +124,6 @@ energy_range = pixel_data.ENERGY_RANGE count_rates = count_rates/(energy_range[1]-energy_range[0]) counts_rates_error = counts_rates_error/(energy_range[1]-energy_range[0]) -count_rates_bkg = count_rates_bkg/(energy_range[1]-energy_range[0]) -count_rates_error_bkg = count_rates_error_bkg/(energy_range[1]-energy_range[0]) - ;;************** Normalization for effective area: units are counts s^-1 keV^-1 cm^-2 subc_str = stx_construct_subcollimator() @@ -184,8 +133,6 @@ eff_area = n_elements( pixel_ind ) eq 1 ? reform(eff_area[*, *, pixel_ind]) : to count_rates = count_rates/eff_area counts_rates_error = counts_rates_error/eff_area -count_rates_bkg = count_rates_bkg/eff_area -count_rates_error_bkg = count_rates_error_bkg/eff_area ;;************** Fill in calibrated pixel data structure @@ -196,11 +143,8 @@ pixel_data_summed.TIME_RANGE = pixel_data.TIME_RANGE pixel_data_summed.ENERGY_RANGE = pixel_data.ENERGY_RANGE pixel_data_summed.COUNT_RATES = count_rates pixel_data_summed.COUNTS_RATES_ERROR = counts_rates_error -pixel_data_summed.TOT_COUNTS = tot_counts -pixel_data_summed.LIVE_TIME_BKG = pixel_data.LIVE_TIME_BKG -pixel_data_summed.COUNT_RATES_BKG = count_rates_bkg -pixel_data_summed.COUNT_RATES_ERROR_BKG = count_rates_error_bkg -pixel_data_summed.TOT_COUNTS_BKG = tot_counts_bkg +pixel_data_summed.TOT_COUNTS = pixel_data.tot_counts +pixel_data_summed.TOT_COUNTS_BKG = pixel_data.tot_counts_bkg pixel_data_summed.RCR = pixel_data.RCR pixel_data_summed.SUMCASE = sumcase pixel_data_summed.DETECTOR_MASKS = pixel_data.DETECTOR_MASKS diff --git a/stix/idl/processing/spectrogram/stx_elut_correction.pro b/stix/idl/processing/spectrogram/stx_elut_correction.pro new file mode 100644 index 00000000..070c5312 --- /dev/null +++ b/stix/idl/processing/spectrogram/stx_elut_correction.pro @@ -0,0 +1,272 @@ +;+ +; +; NAME: +; +; stx_elut_correction +; +; PURPOSE: +; +; Applies ELUT correction to pixel data counts before summing different energy bins. +; +; +; CALLING SEQUENCE: +; +; elut_data = stx_elut_correction(counts, counts_error, energy_bin_idx, $ +; energy_bin_low, energy_bin_high, energy_high, energy_low, energy_ind, this_energy_range, $ +; spectrum, spectrum_with_bkg, spectrum_bkg, $ +; pixel_ind, subc_index, silent=silent) +; +; INPUTS: +; +; COUNTS_IN: float array (dimension: number of energy bins x number of pixels x number of detectors x number of time bins) containing the BKG subtracted counts +; COUNTS_IN_ERROR: float array (dimension: number of energy bins x number of pixels x number of detectors ) containing the uncertainties of BKG subtracted counts +; ENERGY_BIN_IDX: array containing the indices of the energy bins contained in the fits file. +; ENERGY_BIN_LOW: float array of dimension 32 x 12 x 32 (number of energy bins x number of pixels x number of detectors) +; containing the low energy edges of the daily ELUT. +; ENERGY_BIN_HIGH: float array of dimension 32 x 12 x 32 (number of energy bins x number of pixels x number of detectors) +; containing the high energy edges of the daily ELUT. +; ENERGY_HIGH: float array containing the high energy edges of the nominal energy bins. +; ENERGY_LOW: float array containing the low energy edges of the nominal energy bins. +; ENERGY_IND: array containing the indices of the energy bins selected for integration along the energy dimension. +; ENERGY_RANGE: array containing the low and high edge of the energy range selected for integration (e.g., [10,15]) +; SPECTRUM: array containing the BKG-subtracted spectrum. It is used for deriving the spectral index for each energy bin. +; PIXEL_IND: array containing the indices of the considered pixels (e.g., [0] for top row, [0,1] for top and bottom row, etc.) +; SUBC_INDEX: array containing the indices of the considered sub-collimators. +; +; OUTPUTS: +; Structure containing: +; - COUNTS: array of dimension 12 x 32 x number of time bins containing the ELUT corrected counts integrated within the selected energy interval. +; - COUNTS_ERROR: array of dimension 12 x 32 x number of time bins containing the uncertainty associated with the ELUT corrected counts integrated within the selected energy interval. +; - COUNTS_NO_ELUT: array of dimension 4 (A,B,C,D) x number of rows (TOP, BOT, SMALL) x number of subcollimators x number of time bins containing the number of counts integrated in energy without ELUT correction. +; It is used to compute the amount of ELUT correction for all pixels. +; - COUNTS_ELUT: array of dimension 4 (A,B,C,D) x number of rows (TOP, BOT, SMALL) x number of subcollimators x number of time bins containing the number of counts integrated in energy after ELUT correction. +; It is used to compute the amount of ELUT correction for all pixels. +; +; KEYWORDS: +; +; SPECTRUM_WITH_BKG: array containing the observed spectrum (including BKG). It is used for display (if silent is not set to 0) +; SPECTRUM_BKG: array containing the BKG spectrum. It is used for display (if silent is not set to 0) +; SILENT: if set to 1, the spectrum plot is not displayed. +; +; HISTORY: +; +; CONTACT: +; paolo.massa@fhnw.ch +;- + +function stx_elut_correction, counts_in, counts_in_error, $ + energy_bin_idx, energy_bin_low, energy_bin_high, energy_high, energy_low, energy_ind, energy_range, $ + spectrum, pixel_ind, subc_index, $ + spectrum_with_bkg=spectrum_with_bkg, spectrum_bkg=spectrum_bkg, silent=silent + +default, gain_offset_version, 'median' +default, silent, 0 +default, spectrum_with_bkg, fltarr(spectrum.DIM) +default, spectrum_bkg, fltarr(spectrum.DIM) + +;;**************** Compute spectral index to be used for ELUT correction + +index_data = stx_estimate_spectral_index(energy_low, energy_high, spectrum) + +sp_index = index_data.index_final +idx_peak = index_data.idx_peak + +if ~silent then begin + + charsize = 1.8 + + loadct,5 + device, Window_State=win_state + if not win_state[10] then window,10,xsize=1000,ysize=500 + wset,10 + cleanplot + + !p.multi = [0,2,1] + + energy_axis = (energy_low+energy_high)/2. + + plot, energy_axis, spectrum_with_bkg, psym=10, /xst, /yst, /ylog, /xlog, yrange=[max([1.,min(spectrum)]),max(spectrum)*10.],charsize=charsize, $ + title='STIX spectrum',xtitle='Energy [keV]', ytitle = 'STIX spectrum [counts s!U-1!n keV!U-1!n]' + oplot, [energy_range[0],energy_range[0]], [max([1.,min(spectrum)]),max(spectrum)*10.], linestyle=1 + oplot, [energy_range[1],energy_range[1]], [max([1.,min(spectrum)]),max(spectrum)*10.], linestyle=1 + oplot, energy_axis, spectrum_bkg, psym=10, linestyle=2 + oplot, energy_axis, spectrum, psym=10, color=122 + leg_text = ['Observed spectrum', 'Background', 'BKG-subtracted'] + leg_color = [255,255,122] + leg_style = [0, 2, 0] + ssw_legend, leg_text, color=leg_color, linest=leg_style, box=0, charsize=1.5, thick=1.5, /right + + + plot,energy_axis, sp_index, psym=10, /xst, /yst, /xlog, charsize=charsize, $ + title='Estimate of the spectral index', xtitle='Energy [keV]', ytitle = 'Spectral index' + oplot, [energy_range[0],energy_range[0]], [-20,20], linestyle=1 + oplot, [energy_range[1],energy_range[1]], [-20,20], linestyle=1 + oplot, [4,150], [0,0], linestyle=2 + +endif + +;; Determine whether the time axis is present +if n_elements(counts_in.dim) eq 4 then begin + + dimensions = counts_in.dim + n_times = dimensions[3] + +endif else begin + + n_times = 1 + +endelse + +;; We assume the spectrum has a powerlaw distribution E^-sp_index at any energy bin. +; +; The number of counts in the energy range [a,b] is +; +; C = \int_a^b E^-sp_index dE = (b^(-sp_index+1) - a^(-sp_index+1)) / (-sp_index+1) +; +; We distinguish 3 cases: +; +; 1. The considered energy range consists of a single energy bin which is not at the peak +; (i.e., different from 6-7 keV; example: 9-10 keV). The correction factor is +; +; corr = (10^(-sp_index+1) - 9^(-sp_index+1)) / (10.03^(-sp_index+1) - 8.97^(-sp_index+1)), +; +; where 10.03 and 8.97 are derived from daily ELUT +; +; 2. The considered energy range consists of a single energy bin which is at the peak of the spectrum +; (i.e., 6-7 keV when attenuator is not inserted). We consider the middle point of the energy range (e.g., 6.5 keV) +; and we assume that the spectrum follows a powerlaw distribution below and above the middle point, viz. +; +; A E^-alpha if E < 6.5 keV +; F(E) = +; B E^-beta if E >= 6.5 keV +; +; where alpha < 0 and beta > 0. Assuming that the spectrum is continuous at 6.5 keV, +; we obtain that A and B are related to the following equation: +; +; A = B 6.5^(alpha - beta) +; +; Assuming that the ELUT edges of the considered energy bin are 5.98 and 7.05 keV, we obtain +; +; corr = (\int_6^7 F(E) dE) / (\int_5.98^7.05 F(E) dE) , +; +; which leads to +; +; corr = (6.5^(alpha - beta) * (6.5^(-alpha+1) - 6^(-alpha+1)) / (-alpha + 1) ) + (7^(-beta+1) - 6.5^(-beta+1)) / (-beta+1) ) / +; (6.5^(alpha - beta) * (6.5^(-alpha+1) - 5.98^(-alpha+1)) / (-alpha + 1) ) + (7.05^(-beta+1) - 6.5^(-beta+1)) / (-beta+1) ) +; +; 3. The considered energy range consists of multiple energy bins (e.g, 5-10 keV). Therefore, we apply a correction factor +; only to the lowest and to the highest bin (corr_low and corr_high). Assuming that the energy edges of these bins contained +; in the ELUT table are 5.01, 5.98 keV and 8.97, 10.02 keV, we have +; +; corr_low = (5.98^(-alpha+1) - 5^(-alpha+1)) / (5.98^(-alpha+1) - 5.01^(-alpha+1)) +; +; and +; +; corr_high = (10^(-beta+1) - 8.97^(-beta+1)) / (10.02^(-beta+1) - 8.97^(-beta+1)) +; +; where alpha and beta are the powerlaw indices in the first and the last energy bins, respectively. + +if n_elements(energy_ind) eq 1 then begin + + if energy_ind eq idx_peak then begin + + ;; energy_bin_idx contains at least two indices (see control above) + case energy_ind of + + 0: begin + + sp_index_low = sp_index[energy_ind+1] + sp_index_high = sp_index[energy_ind+1] + + end + + n_elements(energy_bin_idx)-1: begin + + sp_index_low = sp_index[energy_ind-1] + sp_index_high = sp_index[energy_ind-1] + + end + + else: begin + + sp_index_low = sp_index[energy_ind-1] + sp_index_high = sp_index[energy_ind+1] + + end + + endcase + + energy_bin_mid = (energy_high[energy_ind] + energy_low[energy_ind]) / 2. + + corr_factor_low_ELUT = (energy_bin_mid^(-sp_index_low+1.) - reform(energy_bin_low[energy_ind,*,*])^(-sp_index_low+1.)) / $ + (-sp_index_low+1.) + corr_factor_high_ELUT = (reform(energy_bin_high[energy_ind,*,*])^(-sp_index_high+1.) - energy_bin_mid^(-sp_index_high+1.)) / $ + (-sp_index_high+1.) + + corr_factor_low_SCI = (energy_bin_mid^(-sp_index_low+1.) - energy_low[energy_ind]^(-sp_index_low+1.)) / (-sp_index_low+1.) + corr_factor_high_SCI = (energy_high[energy_ind]^(-sp_index_high+1.) - energy_bin_mid^(-sp_index_high+1.)) / (-sp_index_high+1.) + + norm_factor = energy_bin_mid^(-sp_index_high + sp_index_low) + energy_corr_factor = (norm_factor * corr_factor_low_SCI + corr_factor_high_SCI) / (norm_factor * corr_factor_low_ELUT + corr_factor_high_ELUT) + + endif else begin + + energy_corr_factor = (energy_high[energy_ind]^(-sp_index[energy_ind]+1.) - energy_low[energy_ind]^(-sp_index[energy_ind]+1.)) / $ + (reform(energy_bin_high[energy_ind,*,*])^(-sp_index[energy_ind]+1.) - reform(energy_bin_low[energy_ind,*,*])^(-sp_index[energy_ind]+1.)) + + endelse + + energy_corr_factor = reform(cmreplicate(energy_corr_factor, n_times)) + + ;; Compute total number of counts BEFORE ELUT correction is applied. It is used for estimating the amount of the ELUT correction + counts_reshaped = reform(counts_in, n_elements(energy_bin_idx), 4, 3, 32, n_times) + counts_no_elut = reform(total(counts_reshaped[energy_ind,*,pixel_ind,subc_index,*], 1)) + + counts = reform(counts_in[energy_ind,*,*,*]) * energy_corr_factor + counts_error = reform(counts_in_error[energy_ind,*,*,*]) * energy_corr_factor + + ;; Compute total number of counts AFTER ELUT correction is applied. It is used for estimating the amount of the ELUT correction + counts_reshaped = reform(counts, 4, 3, 32, n_times) + counts_elut = reform(counts_reshaped[*,pixel_ind,subc_index,*]) + + +endif else begin + + energy_corr_factor_low = (reform(energy_bin_high[energy_ind[0],*,*])^(-sp_index[energy_ind[0]]+1.) - energy_low[energy_ind[0]]^(-sp_index[energy_ind[0]]+1.)) / $ + (reform(energy_bin_high[energy_ind[0],*,*])^(-sp_index[energy_ind[0]]+1.) - reform(energy_bin_low[energy_ind[0],*,*])^(-sp_index[energy_ind[0]]+1.)) + + energy_corr_factor_high = (energy_high[energy_ind[-1]]^(-sp_index[energy_ind[-1]]+1.) - reform(energy_bin_low[energy_ind[-1],*,*])^(-sp_index[energy_ind[-1]]+1.)) / $ + (reform(energy_bin_high[energy_ind[-1],*,*])^(-sp_index[energy_ind[-1]]+1.) - reform(energy_bin_low[energy_ind[-1],*,*])^(-sp_index[energy_ind[-1]]+1.)) + + ;; Compute total number of counts BEFORE ELUT correction is applied. It is used for estimating the amount of the ELUT correction + counts_reshaped = reform(counts_in, n_elements(energy_bin_idx), 4, 3, 32, n_times) + counts_no_elut = reform(total(counts_reshaped[energy_ind,*,pixel_ind,subc_index,*], 1)) + + energy_corr_factor_low = reform(cmreplicate(energy_corr_factor_low, n_times)) + energy_corr_factor_high = reform(cmreplicate(energy_corr_factor_high, n_times)) + + counts = counts_in + counts_error = counts_in_error + + counts[energy_ind[0],*,*,*] *= energy_corr_factor_low + counts[energy_ind[-1],*,*,*] *= energy_corr_factor_high + counts_error[energy_ind[0],*,*,*] *= energy_corr_factor_low + counts_error[energy_ind[-1],*,*,*] *= energy_corr_factor_high + + counts = total(counts[energy_ind,*,*,*],1) + counts_error = sqrt(total(counts_error[energy_ind,*,*,*]^2.,1)) + + ;; Compute total number of counts AFTER ELUT correction is applied. It is used for estimating the amount of the ELUT correction + counts_reshaped = reform(counts, 4, 3, 32, n_times) + counts_elut = reform(counts_reshaped[*,pixel_ind,subc_index,*]) + +endelse + +return, {counts: counts, $ + counts_error: counts_error, $ + counts_no_elut: counts_no_elut, $ + counts_elut: counts_elut, $ + sp_index: sp_index} + + +end \ No newline at end of file diff --git a/stix/idl/processing/spectrogram/stx_estimate_spectral_index.pro b/stix/idl/processing/spectrogram/stx_estimate_spectral_index.pro new file mode 100644 index 00000000..b92d741c --- /dev/null +++ b/stix/idl/processing/spectrogram/stx_estimate_spectral_index.pro @@ -0,0 +1,119 @@ +;+ +; +; NAME: +; +; stx_estimate_spectral_index +; +; PURPOSE: +; +; Computes the spectral index in the different energy bins of a STIX spectrum, assuming a powerlaw dostribution in each energy bin. +; +; CALLING SEQUENCE: +; +; results = stx_estimate_spectral_index(e_low, e_high, spectrum) +; +; INPUTS: +; +; e_low: array containing the lower edges of the energy bins +; +; e_high: array containing the higher edges of the energy bins +; +; spectrum: array containing the values of the STIX spectrum in the different energy bins +; +; OUTPUTS: +; +; Structure containing the following fields: +; +; index_final: array containing the final estimate of the spectral indices in the different energy bins +; +; idx_peak: int, index corresponding to the peak of the spectrum. The spectral index at the peak is 0 by default. +; +; +; HISTORY: April 2025, Krucker S. and Massa P., first release +; +; CONTACT: +; samuel.krucker@fhnw.ch +; paolo.massa@fhnw.ch +; +;- +function stx_estimate_spectral_index, e_low, e_high, spectrum + + e_mean=(e_low+e_high)/2. + ;weighted average calculated below + e_weighted=e_mean*0. + de=(e_high-e_low) + + + ;;**** Calculate power law using the center of the energy bin + + ;value using center of bin + index_tmp=fltarr(n_elements(spectrum)) + ;value using weighted average of energy bin + index_final=fltarr(n_elements(spectrum)) + + idx_peak = where(spectrum eq max(spectrum)) + idx_peak = idx_peak[0] + + if idx_peak gt 0 then begin + + for i=0,idx_peak-1 do begin + this_x=alog(e_mean(i+1))-alog(e_mean(i)) + this_y=alog(spectrum(i+1))-alog(spectrum(i)) + index_tmp(i)=this_y/this_x + ;use initial value to calculate weighted center + h0=(e_high(i)^(index_tmp(i)+2)-e_low(i)^(index_tmp(i)+2))/(index_tmp(i)+2) + h1=(e_high(i)^(index_tmp(i)+1)-e_low(i)^(index_tmp(i)+1))/(index_tmp(i)+1) + e_weighted(i)=h0/h1 + endfor + + endif + + ;above peak (negative index) + for i=idx_peak+1,n_elements(spectrum)-1 do begin + this_x=alog(e_mean(i))-alog(e_mean(i-1)) + this_y=alog(spectrum(i))-alog(spectrum(i-1)) + index_tmp(i)=this_y/this_x + ;use initial value to calculate weighted center + h0=(e_high(i)^(index_tmp(i)+2)-e_low(i)^(index_tmp(i)+2))/(index_tmp(i)+2) + h1=(e_high(i)^(index_tmp(i)+1)-e_low(i)^(index_tmp(i)+1))/(index_tmp(i)+1) + e_weighted(i)=h0/h1 + endfor + + ;for peak, no value is calculated + ;set e_weighted to center + e_weighted(idx_peak)=e_mean(idx_peak) + + ;recalculate slopes with weighted center + for i=idx_peak+1,n_elements(spectrum)-1 do begin + this_x=alog(e_weighted(i))-alog(e_weighted(i-1)) + this_y=alog(spectrum(i))-alog(spectrum(i-1)) + index_final(i)=this_y/this_x + endfor + + if idx_peak gt 0 then begin + + for i=0,idx_peak-1 do begin + this_x=alog(e_weighted(i+1))-alog(e_weighted(i)) + this_y=alog(spectrum(i+1))-alog(spectrum(i)) + index_final(i)=this_y/this_x + endfor + + endif + + ;; Change sign to be consistent with definition of powerlaw: E^-idx + index_final = -index_final + + ;; Set NAN or INF values equal to 0 + idx_nan_inf = where(~finite(index_final), /null) + index_final[idx_nan_inf] = 0. + + ;; Set indices >= 8 or <= -8 equal to 8 and -8, respectively. + idx_large_abs = where(abs(index_final) ge 8.) + index_final[idx_large_abs] = signum(index_final[idx_large_abs]) * 8. + + results = {index_final: index_final, $ + idx_peak: idx_peak} + + return, results + +end \ No newline at end of file diff --git a/stix/idl/processing/vis/stx_calibrate_visibility.pro b/stix/idl/processing/vis/stx_calibrate_visibility.pro index 5b5d36a4..64347237 100644 --- a/stix/idl/processing/vis/stx_calibrate_visibility.pro +++ b/stix/idl/processing/vis/stx_calibrate_visibility.pro @@ -19,6 +19,11 @@ ; ; KEYWORDS: ; +; mapcenter: two-element array containing the coordinates of the center of the map to reconstruct +; from the visibility values (STIX coordinate frame, arcsec). It is used during the visibility phase calibration +; process. A phase factor is added to the visibilities so that coordinates saved in 'mapcenter' become the +; center of the reconstructed map. +; ; xy_flare: two-element array containing the coordinates of the estimated flare location (STIX coordinate frame, arcsec). ; It is used for computing the grid transmission correction within the visibility amplitude calibration. Default, (0,0) ; @@ -35,12 +40,6 @@ ; ; syserr_sigamp: float, percentage of systematic error to be added to the visibility amplitude errors. ; Default, 5% -; -; f2r_sep: separation between the front and the rear grid (mm). Default, 550 mm. It is used for computing the default -; projection correction factors -; -; r2d_sep: separation between the rear grid and the detector (mm, used for the phase projection correction). -; Default, 47 mm. It is used for computing the default ; ; OUTPUTS: ; @@ -50,16 +49,16 @@ ; July 2023, Massa P., removed visibility phase 'projection correction' since the new definition of ; (u,v)-points is adopted (see stx_uv_points). ; February 2026, Massa P., new visibility amplitude calibration is implemented +; March 2026, Massa P., 'mapcenter' keyword is added here. Also 'f2r_sep' and 'r2d_sep' keyword are removed as not needed ; ; CONTACT: -; paolo.massa@wku.edu +; paolo.massa@fhnw.ch ;- -function stx_calibrate_visibility, vis, xy_flare=xy_flare, phase_calib_factors=phase_calib_factors, amp_calib_factors=amp_calib_factors, $ - syserr_sigamp = syserr_sigamp, r2d_sep=r2d_sep, f2r_sep=f2r_sep +function stx_calibrate_visibility, vis, mapcenter=mapcenter, xy_flare=xy_flare, phase_calib_factors=phase_calib_factors, amp_calib_factors=amp_calib_factors, $ + syserr_sigamp = syserr_sigamp -default, f2r_sep, 545.30 -default, r2d_sep, 47.78 +default, mapcenter, [0.,0.] default, xy_flare, [0.,0.] n_vis = n_elements(vis) @@ -82,7 +81,7 @@ tmp = read_csv(loc_file( 'PhaseCorrFactors.csv', path = getenv('STX_VIS_PHASE')) ad_hoc_phase_corr = tmp.field2[vis.ISC - 1] ;; Mapcenter correction -phase_mapcenter_corr = -2 * !pi * (vis.XYOFFSET[0] * vis.U + vis.XYOFFSET[1] * vis.V ) * !radeg +phase_mapcenter_corr = -2 * !pi * (mapcenter[0] * vis.U + mapcenter[1] * vis.V ) * !radeg default, amp_calib_factors, fltarr(n_vis) + modulation_efficiency default, phase_calib_factors, grid_phase_corr + ad_hoc_phase_corr + phase_mapcenter_corr @@ -117,6 +116,7 @@ calibrated_vis.sigamp = calibrated_sigamp calibrated_vis.CALIBRATED = 1 calibrated_vis.XY_FLARE = xy_flare +calibrated_vis.XYOFFSET = mapcenter return, calibrated_vis diff --git a/stix/idl/processing/vis/stx_construct_calibrated_visibility.pro b/stix/idl/processing/vis/stx_construct_calibrated_visibility.pro index 8bb6a256..5ebb4507 100644 --- a/stix/idl/processing/vis/stx_construct_calibrated_visibility.pro +++ b/stix/idl/processing/vis/stx_construct_calibrated_visibility.pro @@ -29,20 +29,12 @@ ; path_bkg_file: if provided, the background counts are subtracted before computing the visibility values and the ; field 'TOT_COUNTS_BKG' of the visibility structure is filled in ; -; elut_corr: if set, a correction based on a ELUT table is applied to the measured counts -; ; xy_flare: two-element array containing the coordinates of the estimated flare location (STIX coordinate frame, arcsec). ; It is used for computing the grid transmission correction within the visibility amplitude calibration. ; If it is not provided, then the flare location is set equal to (0,0) ; ; sumcase: string containing information on the pixels to be considered for computing the visibilities ; (see the header of 'stx_sum_pixel_data' for more details). -; -; f2r_sep: distance between the front and the rear grid (mm, used for computing the values of the (u,v) frequencies). -; Default, 550 mm -; -; r2d_sep: distance between the rear grid and the detectors (mm, used for computing the projection correction factors). -; Default, 47 mm ; ; silent: if set, no message is printed ; @@ -68,14 +60,16 @@ ; amplitudes and phases. For information on the fields of the visibility structure see the header of 'stx_construct_visibility' ; ; HISTORY: August 2022, Massa P., created +; March 2026, Massa P., removed 'mapcenter' keyword from 'stx_construct_visibility'. Also, 'r2d_sep', 'f2r_sep=f2r_sep', and 'elut_corr' +; keywords are removed as not necessary ; ; CONTACT: ; paolo.massa@fhnw.ch ;- function stx_construct_calibrated_visibility, path_sci_file, time_range, energy_range, mapcenter, $ - path_bkg_file=path_bkg_file, elut_corr=elut_corr, xy_flare=xy_flare, $ - sumcase=sumcase, f2r_sep=f2r_sep, r2d_sep=r2d_sep, silent=silent, $ + path_bkg_file=path_bkg_file, xy_flare=xy_flare, $ + sumcase=sumcase, silent=silent, $ subc_index=subc_index, phase_calib_factors=phase_calib_factors, $ amp_calib_factors=amp_calib_factors, syserr_sigamp = syserr_sigamp, $ no_small=no_small, no_rcr_check=no_rcr_check, _extra=extra @@ -84,15 +78,15 @@ function stx_construct_calibrated_visibility, path_sci_file, time_range, energy_ ;;*********** Create visibility structure -vis = stx_construct_visibility(path_sci_file, time_range, energy_range, mapcenter, path_bkg_file=path_bkg_file, $ - elut_corr=elut_corr, sumcase=sumcase, f2r_sep=f2r_sep, silent=silent, $ - subc_index=subc_index, no_small=no_small, no_rcr_check=no_rcr_check, _extra=extra) +vis = stx_construct_visibility(path_sci_file, time_range, energy_range, path_bkg_file=path_bkg_file, $ + sumcase=sumcase, silent=silent, subc_index=subc_index, $ + no_small=no_small, no_rcr_check=no_rcr_check, _extra=extra) ;;*********** Calibrate visibility -calibrated_vis = stx_calibrate_visibility(vis, xy_flare=xy_flare, phase_calib_factors=phase_calib_factors, $ +calibrated_vis = stx_calibrate_visibility(vis, mapcenter=mapcenter, xy_flare=xy_flare, phase_calib_factors=phase_calib_factors, $ amp_calib_factors=amp_calib_factors, $ - syserr_sigamp = syserr_sigamp, r2d_sep=r2d_sep, f2r_sep=f2r_sep) + syserr_sigamp = syserr_sigamp, _extra=extra) ;;********** Plot visibility amplitudes vs resolution (if ~silent) diff --git a/stix/idl/processing/vis/stx_construct_visibility.pro b/stix/idl/processing/vis/stx_construct_visibility.pro index d2dae6e9..7d5f3492 100644 --- a/stix/idl/processing/vis/stx_construct_visibility.pro +++ b/stix/idl/processing/vis/stx_construct_visibility.pro @@ -20,25 +20,15 @@ ; time_range: string array containing the start and the end of the time interval to consider ; ; energy_range: array containing the values of the lower and upper edge of the energy interval to consider -; -; mapcenter: two-element array containing the coordinates of the center of the map to reconstruct -; from the visibility values (STIX coordinate frame, arcsec). It is used during the visibility phase calibration -; process. A phase factor is added to the visibilities so that coordinates saved in 'mapcenter' appear in the -; center of the reconstructed map. Usually corresponds to an estimate of the flare location ; ; KEYWORDS: ; ; path_bkg_file: if provided, the background counts are subtracted to the countrates before computing the visibility values. ; Further, the field 'TOT_COUNTS_BKG' of the visibility structure is filled in with an estimate of the total ; number of background counts that are measured in the selected time and energy intervals -; -; elut_corr: if set, a correction based on a ELUT table is applied to the measured counts ; ; sumcase: string indicating which pixels are summed for computing the visibilities (see the header of ; 'stx_calibrate_pixel_data' for more information). Default, 'TOP+BOT' -; -; f2r_sep: distance between the front and the rear grid (mm, used for computing the values of the (u,v) frequencies). -; Default, 550 mm ; ; silent: if set, no message is printed ; @@ -67,7 +57,7 @@ ; - V: v coordinate of the frequencies sampled by the sub-collimators ; - PHASE_SENSE: array containing the sense of the phase measured by the sub-collimator (-1 or 1 values) ; - XYOFFSET: two-element array containing the coordinates of the center of the map to renconstruct from the -; visibility values (it is contains the values passed in 'mapcenter') +; visibility values (Default, (0,0)) ; - XY_FLARE: two-element array containing the coordinates of the estimated flare location (STIX coordinate frame, arcsec), ; which is used for computing the grid transmission correction. It is initialized with NaN values. ; If no correction is applied later (see stx_calibrate_visibility.pro), it remains filled with NaNs. @@ -76,26 +66,25 @@ ; ; HISTORY: August 2022, Massa P., created ; January 2026, Massa P., removed 'xy_flare' keyword. Grid transmission correction is used only for visibility amplitude calibration +; March 2026, Massa P., removed 'mapcenter', 'elut_corr', and 'f2r_sep' keywords as not necessary ; ; CONTACT: ; paolo.massa@fhnw.ch ;- -function stx_construct_visibility, path_sci_file, time_range, energy_range, mapcenter, path_bkg_file=path_bkg_file, $ - elut_corr=elut_corr, sumcase=sumcase, f2r_sep=f2r_sep, silent=silent, $ - subc_index=subc_index, no_small=no_small, no_rcr_check=no_rcr_check, _extra=extra +function stx_construct_visibility, path_sci_file, time_range, energy_range, path_bkg_file=path_bkg_file, $ + sumcase=sumcase, silent=silent, subc_index=subc_index, no_small=no_small, $ + no_rcr_check=no_rcr_check, _extra=extra ;;************** Construct a 'stx_pixel_data_summed' structure pixel_data_summed = stx_construct_pixel_data_summed(path_sci_file, time_range, energy_range, $ - path_bkg_file=path_bkg_file, $ - elut_corr=elut_corr, subc_index=subc_index, $ + path_bkg_file=path_bkg_file, subc_index=subc_index, $ sumcase=sumcase, silent=silent,no_small=no_small, $ no_rcr_check=no_rcr_check, _extra=extra) ;;************** Create an uncalibrated 'stx_visibility' structure from a the correpsonding 'stx_pixel_data_summed' structure -vis = stx_pixel_data_summed2visibility(pixel_data_summed, subc_index=subc_index, $ - f2r_sep=f2r_sep, mapcenter=mapcenter) +vis = stx_pixel_data_summed2visibility(pixel_data_summed, subc_index=subc_index, _extra=extra) return, vis diff --git a/stix/idl/processing/vis/stx_pixel_data_summed2visibility.pro b/stix/idl/processing/vis/stx_pixel_data_summed2visibility.pro index 3bc67a63..70ecbe21 100644 --- a/stix/idl/processing/vis/stx_pixel_data_summed2visibility.pro +++ b/stix/idl/processing/vis/stx_pixel_data_summed2visibility.pro @@ -19,13 +19,7 @@ ; KEYWORDS: ; ; subc_index: array containing the indexes of the subcollimators to be considered for computing the visibility -; values -; -; mapcenter: two-element array containing the coordinates of the center of the map to reconstruct -; from the visibility values (STIX coordinate frame, arcsec) -; -; f2r_sep: distance between the front and the rear grid (mm, used for computing the values of the (u,v) frequencies) -; +; values ; ; OUTPUTS: ; @@ -54,16 +48,17 @@ ; ; HISTORY: August 2022, Massa P., created ; July 2023, Massa P., made it compatible with the new definition of (u,v)-points (see stx_uv_points) +; March 2026, Massa P., removed 'mapcenter' keyword as it is always (0,0) at this stage. +; A different value for 'mapcenter' can be set in 'stx_calibrate_pixel_data.pro'. +; Further, removed 'f2r_sep' keyword and added '_extra=extra'. Finally, removed background +; subtraction step as performed already in 'stx_construct_pixel_data.pro' ; ; CONTACT: -; paolo.massa@wku.edu +; paolo.massa@fhnw.ch ;- -function stx_pixel_data_summed2visibility, pixel_data_summed, subc_index=subc_index, mapcenter=mapcenter, $ - f2r_sep=f2r_sep +function stx_pixel_data_summed2visibility, pixel_data_summed, subc_index=subc_index, _extra=extra -default, mapcenter, [0.,0.] -default, f2r_sep, 545.30 default, subc_index, stx_label2ind(['10a','10b','10c','9a','9b','9c','8a','8b','8c','7a','7b','7c',$ '6a','6b','6c','5a','5b','5c','4a','4b','4c','3a','3b','3c']) @@ -79,7 +74,7 @@ subc_str = subc_str[subc_index] ;;************** Define (u,v) points -uv_data = stx_uv_points(subc_index) +uv_data = stx_uv_points(subc_index, _extra=extra) u = uv_data.u v = uv_data.v @@ -90,15 +85,6 @@ count_rates = count_rates[subc_index,*] counts_rates_error = pixel_data_summed.COUNTS_RATES_ERROR counts_rates_error = counts_rates_error[subc_index,*] -count_rates_bkg = pixel_data_summed.COUNT_RATES_BKG -count_rates_bkg = count_rates_bkg[subc_index,*] -count_rates_error_bkg = pixel_data_summed.COUNT_RATES_ERROR_BKG -count_rates_error_bkg = count_rates_error_bkg[subc_index,*] - -;; Background subtraction -count_rates = count_rates - count_rates_bkg -counts_rates_error = sqrt(counts_rates_error^2 + count_rates_error_bkg^2) - ;; A,B,C,D A = count_rates[*,0] B = count_rates[*,1] @@ -170,7 +156,6 @@ vis.LIVE_TIME = pixel_data_summed.LIVE_TIME[subc_index] vis.ENERGY_RANGE = pixel_data_summed.ENERGY_RANGE vis.TIME_RANGE = pixel_data_summed.TIME_RANGE vis.PHASE_SENSE = subc_str.PHASE -vis.XYOFFSET = mapcenter return, vis diff --git a/stix/idl/structures/stx_pixel_data.pro b/stix/idl/structures/stx_pixel_data.pro index b85c27c0..5c9bfc6e 100644 --- a/stix/idl/structures/stx_pixel_data.pro +++ b/stix/idl/structures/stx_pixel_data.pro @@ -18,7 +18,10 @@ ; ; ; HISTORY: July 2022, Massa P., created -; January 2026, Massa P., removed 'xy_flare' entry as grid transmission correction is not applied anymore to the raw counts +; April 2025, Massa P., made it compatible with new ELUT correction. Bkg-subtraction is applied to the counts. +; Therefore, information on bkg is not stored anymore +; Fubruary 2026, Massa P., Removed 'xy_flare' keyword as grid transmission correction is not applied anymore to raw counts +; March 2026, Massa P., updated to make it compatible with new ELUT correction. ; ; CONTACT: ; paolo.massa@fhnw.ch @@ -29,13 +32,13 @@ function stx_pixel_data return, { $ type : 'stx_pixel_data', $ live_time : fltarr(32), $ ; Live time of the 32 detectors + live_time_error : fltarr(32), $ ; Live time error of the 32 detectors time_range : replicate(stx_time(),2), $ ; Selected time range (edges) energy_range : fltarr(2), $ ; Selected energy range (edges) counts : dblarr(32,12), $ ; Counts recorded by the detector pixels (summed in time and energy) - counts_error : dblarr(32,12), $ ; Errors associated with the measured counts (statistics + compression) - live_time_bkg : fltarr(32), $ ; Live time of the 32 detectors during the background measurement - counts_bkg : dblarr(32,12), $ ; Counts recorded by the detector pixels during the background measurement (summed in time and energy) - counts_error_bkg : dblarr(32,12), $ ; Errors associated with the measured background counts (statistics + compression) + counts_error : dblarr(32,12), $ ; Errors associated with the measured counts (statistics + compression) + tot_counts : double(0), $ ; Estimate of the total number of flare counts recorded during the flaring event + tot_counts_bkg : double(0), $ ; Estimate of the total number of background counts recorded during the flaring event rcr : byte(0), $ ; Rate Control Regime (RCR) status pixel_masks : bytarr(12), $ ; Matrix containing information on the pixels used for the measurement detector_masks : bytarr(32) $ ; Matrix containing information on the detectors used for the measurement diff --git a/stix/idl/structures/stx_pixel_data_summed.pro b/stix/idl/structures/stx_pixel_data_summed.pro index 59cf70d4..39dfdf89 100644 --- a/stix/idl/structures/stx_pixel_data_summed.pro +++ b/stix/idl/structures/stx_pixel_data_summed.pro @@ -19,6 +19,7 @@ ; ; HISTORY: July 2022, Massa P., created ; January 2026, Massa P., removed 'xy_flare' entry as grid transmission correction is not applied anymore to the raw counts +; March 2026, Massa P., made it compatible with new ELUT correction ; ; CONTACT: ; paolo.massa@fhnw.ch @@ -37,15 +38,10 @@ function stx_pixel_data_summed counts_rates_error : dblarr(32,4), $ ; Errors associated with the measured counts rates (statistics + compression, ; no systematic errors are included) tot_counts : double(0), $ ; Total number of counts recorded during the event - live_time_bkg : fltarr(32), $ ; Live time of the 32 detectors during the background measurement - count_rates_bkg : dblarr(32,4), $ ; Counts rates recored by the the sub-collimators (summed in time and energy) during the - ; background measurement - count_rates_error_bkg : dblarr(32,4), $ ; Errors associated with the measured background counts (statistics + compression, - ; no systematic errors are included) tot_counts_bkg : double(0), $ ; Estimate of the total number of background counts recorded during the flaring event rcr : byte(0), $ ; Rate Control Regime (RCR) status sumcase : string(""), $ ; Which pixels are summed: 'TOP' (top row), 'BOT' (bottom row), 'SMALL' (small pixels) - ; 'TOP+BOT' (top and bottom row), 'ALL', (all pixels) + ; 'TOP+BOT' (top and bottom row), 'ALL', (all pixels) detector_masks : bytarr(32) $ ; array containing information on the detectors used for the measurement ; (1 if the detector is used, 0 otherwise) } From 7e6321c7bfefe4586f17d2eefa7cd72da8e12545 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Mon, 23 Mar 2026 17:43:05 +0100 Subject: [PATCH 02/13] New calibration files for spectroscopy and minor fixes New calibration files are now used for spectroscopy (only for CPD files, not for spectrograms for now) --- .../dbase/grid/real_bkg_grid_transmission.txt | 14 + .../livetime/stx_spectrogram_livetime.pro | 43 +- .../processing/livetime/stx_test_livetime.pro | 12 +- .../stx_triggergram_construction_4test.pro | 7 +- .../spectrogram/stx_convert_pixel_data.pro | 554 +++++++++++++++--- .../spectrogram/stx_elut_correction.pro | 2 + 6 files changed, 505 insertions(+), 127 deletions(-) create mode 100644 stix/dbase/grid/real_bkg_grid_transmission.txt diff --git a/stix/dbase/grid/real_bkg_grid_transmission.txt b/stix/dbase/grid/real_bkg_grid_transmission.txt new file mode 100644 index 00000000..eb8d57b3 --- /dev/null +++ b/stix/dbase/grid/real_bkg_grid_transmission.txt @@ -0,0 +1,14 @@ +;Real transmission factors for on-axis grid transmission in BGK detector pixels MZS - 18-Jun-24 +;N.B. tranmission for fully covered pixels is currently set to 0 +0.1341096 +0.0010084 +0.0115397 +0.0000000 +0.0000000 +0.0113317 +0.0009876 +0.1393076 +0.0000000 +0.0000000 +0.0000000 +0.0000000 diff --git a/stix/idl/processing/livetime/stx_spectrogram_livetime.pro b/stix/idl/processing/livetime/stx_spectrogram_livetime.pro index 95d17424..d5425e29 100644 --- a/stix/idl/processing/livetime/stx_spectrogram_livetime.pro +++ b/stix/idl/processing/livetime/stx_spectrogram_livetime.pro @@ -37,6 +37,7 @@ ; :history: ; 24-Feb-2021 - ECMD (Graz), initial release ; 22-Feb-2022 - ECMD (Graz), documented, added livetime error estimation +; 23-Mar-2026 - Massa P. (FHNW), made it compatible with live time uncertainty propagation ; ;- function stx_spectrogram_livetime, spectrogram, corrected_counts =corrected_counts, corrected_error= corrected_error, level = level @@ -52,36 +53,21 @@ default, level, 1 case level of 1: begin dim_counts = [nenergies, ndet, ntimes] - triggergram = stx_triggergram(transpose( spectrogram.trigger ), spectrogram.time_axis) - livetime_fraction = stx_livetime_fraction(triggergram, det_used) - livetime_fraction = transpose( rebin(reform(livetime_fraction, ndet, ntimes),[dim_counts[1:2],dim_counts[0]]),[2,0,1]) - - triggergram_lower = stx_triggergram(transpose( spectrogram.trigger - spectrogram.trigger_err > 0), spectrogram.time_axis) - livetime_fraction_lower = stx_livetime_fraction(triggergram_lower, det_used) - livetime_fraction_lower = transpose( rebin(reform(livetime_fraction_lower, ndet, ntimes),[dim_counts[1:2],dim_counts[0]]),[2,0,1]) - - triggergram_upper = stx_triggergram(transpose( spectrogram.trigger + spectrogram.trigger_err), spectrogram.time_axis) - livetime_fraction_upper = stx_livetime_fraction(triggergram_upper, det_used) - livetime_fraction_upper = transpose( rebin(reform(livetime_fraction_upper, ndet, ntimes),[dim_counts[1:2],dim_counts[0]]),[2,0,1]) + triggergram = stx_triggergram(transpose( spectrogram.trigger ), transpose( spectrogram.trigger_err ), spectrogram.time_axis) + livetime_fraction_data = stx_livetime_fraction(triggergram, det_used) + livetime_fraction = transpose( rebin(reform(livetime_fraction_data.livetime_fraction, ndet, ntimes),[dim_counts[1:2],dim_counts[0]]),[2,0,1]) + livetime_fraction_err = transpose( rebin(reform(livetime_fraction_data.livetime_fraction_err, ndet, ntimes),[dim_counts[1:2],dim_counts[0]]),[2,0,1]) end 4:begin dim_counts = [nenergies, ntimes] trig = (fltarr(16)+1./16.)##transpose(spectrogram.trigger) - triggergram = stx_triggergram(transpose(trig), spectrogram.time_axis) - livetime_fraction = stx_livetime_fraction(triggergram, det_used) - livetime_fraction = transpose( rebin(reform(livetime_fraction[0,*]),[dim_counts[1],dim_counts[0]])) - - trig_lower = (fltarr(16)+1./16.)##transpose(spectrogram.trigger - spectrogram.trigger_err > 0) - triggergram_lower = stx_triggergram(transpose(trig_lower), spectrogram.time_axis) - livetime_fraction_lower = stx_livetime_fraction(triggergram_lower, det_used) - livetime_fraction_lower = transpose( rebin(reform(livetime_fraction_lower[0,*]),[dim_counts[1],dim_counts[0]])) - - trig_upper = (fltarr(16)+1./16.)##transpose(spectrogram.trigger + spectrogram.trigger_err) - triggergram_upper = stx_triggergram(transpose(trig_upper), spectrogram.time_axis) - livetime_fraction_upper = stx_livetime_fraction(triggergram_upper, det_used) - livetime_fraction_upper = transpose( rebin(reform(livetime_fraction_upper[0,*]),[dim_counts[1],dim_counts[0]])) + trig_err = (fltarr(16)+1./16.)##transpose(spectrogram.trigger_err) + triggergram = stx_triggergram(transpose(trig), transpose(trig_err), spectrogram.time_axis) + livetime_fraction_data = stx_livetime_fraction(triggergram, det_used) + livetime_fraction = transpose( rebin(reform(livetime_fraction_data.livetime_fraction[0,*]),[dim_counts[1],dim_counts[0]])) + livetime_fraction_err = transpose( rebin(reform(livetime_fraction_data.livetime_fraction_err[0,*]),[dim_counts[1],dim_counts[0]])) end @@ -89,13 +75,10 @@ default, level, 1 else: message, 'Currently supported compaction levels are 1 (pixel data) and 4 (spectrogram)' endcase - corrected_counts = spectrogram.counts/livetime_fraction - - corrected_counts_upper = spectrogram.counts/livetime_fraction_upper - corrected_counts_lower = spectrogram.counts/livetime_fraction_lower - error_from_livetime = (corrected_counts_upper - corrected_counts_lower)/2. + corrected_counts = f_div(spectrogram.counts,livetime_fraction) - corrected_error = sqrt( (spectrogram.error/livetime_fraction)^2. + error_from_livetime^2.) + corrected_error = abs(corrected_counts) * sqrt( f_div(spectrogram.error,spectrogram.counts)^2. + $ + f_div(livetime_fraction_err,livetime_fraction)^2. ) return, livetime_fraction diff --git a/stix/idl/processing/livetime/stx_test_livetime.pro b/stix/idl/processing/livetime/stx_test_livetime.pro index e3178f22..b0035ece 100644 --- a/stix/idl/processing/livetime/stx_test_livetime.pro +++ b/stix/idl/processing/livetime/stx_test_livetime.pro @@ -1,7 +1,9 @@ ;stx_test_livetime, main program to build a testbest for the livetime computaton ; triggergram = stx_triggergram_construction_4test( rate, duration = duration ) -livetime_fraction = stx_livetime_fraction( triggergram, det_select, tau_array = tau_array ) +livetime_fraction_data = stx_livetime_fraction( triggergram, det_select, tau_array = tau_array ) +livetime_fraction = livetime_fraction_data.livetime_fraction +livetime_fraction_err = livetime_fraction_data.livetime_fraction_err help, rate, duration help, triggergram help, triggergram, /st @@ -16,13 +18,17 @@ pmm, rate, livetime_fraction rate2 = 10.0 * rate triggergram = stx_triggergram_construction_4test( rate2, duration = duration ) -livetime_fraction = stx_livetime_fraction( triggergram, det_select, tau_array = tau_array ) +livetime_fraction_data = stx_livetime_fraction( triggergram, det_select, tau_array = tau_array ) +livetime_fraction = livetime_fraction_data.livetime_fraction +livetime_fraction_err = livetime_fraction_data.livetime_fraction_err print, 'Min and Max rate and Livetime_fraction pmm, rate2, livetime_fraction rate3 = rate2 * 10.0 triggergram = stx_triggergram_construction_4test( rate3, duration = duration ) -livetime_fraction = stx_livetime_fraction( triggergram, det_select, tau_array = tau_array ) +livetime_fraction_data = stx_livetime_fraction( triggergram, det_select, tau_array = tau_array ) +livetime_fraction = livetime_fraction_data.livetime_fraction +livetime_fraction_err = livetime_fraction_data.livetime_fraction_err print, 'Min and Max rate and Livetime_fraction pmm, rate3, livetime_fraction diff --git a/stix/idl/processing/livetime/stx_triggergram_construction_4test.pro b/stix/idl/processing/livetime/stx_triggergram_construction_4test.pro index f253e168..0fd744a9 100644 --- a/stix/idl/processing/livetime/stx_triggergram_construction_4test.pro +++ b/stix/idl/processing/livetime/stx_triggergram_construction_4test.pro @@ -14,6 +14,7 @@ ; :History: ; ; 18-april-2015, richard.schwartz@nasa.gov +; 17-april-2025, Massa P. (FHNW), make it compatible with new version of stx_triggergram.pro ;- function stx_triggergram_construction_4test, rate, duration = duration @@ -25,7 +26,9 @@ function stx_triggergram_construction_4test, rate, duration = duration ntbin = n_elements( duration ) ratem = rate+fltarr(32) - det_select = (stx_vis()).isc ;imaging subcollimators + subc_str = stx_construct_subcollimator() + idx_imaging_subc = where((subc_str.LABEL ne 'cfl') and (subc_str.LABEL ne 'bkg')) + det_select = subc_str[idx_imaging_subc].DET_N ;(stx_vis()).isc ;imaging subcollimators all_det = indgen( 32 ) + 1 nofourier = all_det [ where_arr( all_det, det_select, /notequal) ] ratem[ nofourier-1 ] = 0 ;very small rates in the background and flare location detectors @@ -39,7 +42,7 @@ function stx_triggergram_construction_4test, rate, duration = duration ixpair = stx_ltpair_assignment( /pairs ) -1 ;det_index (0:31) of components of each trigger for i = 0, 15 do trigger[i,0] = transpose( total( events[ ixpair[*,i], * ], 1) ) - triggergram = stx_triggergram( trigger, time_axis ) + triggergram = stx_triggergram( trigger, sqrt(trigger), time_axis ) ; we estimate trigger errors as the square root of the triggers (i.e., we assume no compression error) return, triggergram end \ No newline at end of file diff --git a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro index b5d93f89..63e30903 100644 --- a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro +++ b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro @@ -34,12 +34,6 @@ ; Shift all energies by this value in keV. Rarely needed only for cases where there is a significant shift ; in calibration before a new ELUT can be uploaded. ; -; flare_location_hpc : in, type="2 element float array" -; the location of the flare (X,Y) in Helioprojective Cartesian coordinates as seen from Solar Orbiter [arcsec] -; -; aux_fits_file : in, required if flare_location_hpc is passed in, type="string" -; the path of the auxiliary ephemeris FITS file to be read." -; ; det_ind : in, type="int array", default="all detectors present in observation" ; indices of detectors to sum when making spectrogram ; @@ -63,12 +57,9 @@ ; ; srmfile : in, type="string", default="'stx_srm_'+ UID + '.fits'" ; File name to use when saving the srm FITS file for OSPEX input. -; +; ; silent : in, type="int", default="0" ; If set prevents informational messages being displayed. -; -; background_data : out, type="stx_background_data structure" -; Structure containing the subtracted background for external plotting. ; ; plot : in, type="boolean", default="1" ; If set open OSPEX GUI and plot lightcurve in standard quicklook energy bands @@ -79,6 +70,11 @@ ; ; ospex_obj : out, type="OSPEX object" ; +; delta_time_min: in, type="float". Pixel data counts are rebinned in time in such a way that the +; time resolution is at least equal to delta_time_min (defined in seconds). +; Rebinned data are used to construct spectra needed for the ELUT correction. +; +; calib_data: if a 'stx_calibration_data' structure is passed as input, then the ELUT correction is applied ; ; :examples: ; fits_path_data = 'solo_L1A_stix-sci-xray-l1-2104170007_20210417T153019-20210417T171825_009610_V01.fits' @@ -99,28 +95,34 @@ ; 16-Jun-2023 - ECMD (Graz), for a source location dependent response estimate, the location in HPC and the auxiliary ephemeris file must be provided. ; 06-Dec-2023 - ECMD (Graz), added silent keyword, more information is now printed if not set ; 2024-07-12, F. Schuller (AIP): added optional keyword xspec +; 07-May-2025 - Stiefel M. and Massa P., re-implemented to take into account new ELUT correction, new subcollimator transmission and live time normalization ; ;- pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fits_path_bk, $ time_shift = time_shift, energy_shift = energy_shift, distance = distance, $ - aux_fits_file = aux_fits_file, flare_location_hpc = flare_location_hpc, flare_location_stx = flare_location_stx, $ - det_ind = det_ind, pix_ind = pix_ind, elut_correction = elut_correction, shift_duration = shift_duration, $ + flare_location_stx = flare_location_stx, det_ind = det_ind, pix_ind = pix_ind, calib_data = calib_data, shift_duration = shift_duration, $ no_attenuation = no_attenuation, sys_uncert = sys_uncert, generate_fits = generate_fits, specfile = specfile, $ - srmfile = srmfile, silent = silent, background_data = background_data, plot = plot, xspec=xspec, ospex_obj = ospex_obj + srmfile = srmfile, silent = silent, plot = plot, xspec=xspec, ospex_obj = ospex_obj, $ + delta_time_min = delta_time_min, tailing=tailing, include_damage=include_damage, _extra=extra default, shift_duration, 0 default, plot, 1 default, det_ind, 'top24' - default, elut_correction, 1 + default, elut_correction, 1 default, silent, 0 default, xspec, 0 + default, sys_uncert, 0.05 + default, delta_time_min, 20. + default, tailing, 1 + default, include_damage, 1 + default, flare_location_stx, [0.,0.] if n_elements(time_shift) eq 0 then begin - if ~keyword_set(silent) then begin - message, 'Time shift value not set, using default value of 0 [s].', /info - print, 'File averaged values can be obtained from the FITS file header' - print, 'using stx_get_header_corrections.pro.' - endif + if ~keyword_set(silent) then begin + message, 'Time shift value not set, using default value of 0 [s].', /info + print, 'File averaged values can be obtained from the FITS file header' + print, 'using stx_get_header_corrections.pro.' + endif time_shift = 0. endif @@ -138,21 +140,60 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit mask_use_pixels = intarr(12) if n_elements(pix_ind) eq 0 then mask_use_pixels[*] = 1 else mask_use_pixels[pix_ind] = 1 + ;;***************** READ SCIENCE AND BKG DATA stx_read_pixel_data_fits_file, fits_path_data, time_shift, primary_header = primary_header, data_str = data_str, data_header = data_header, control_str = control_str, $ control_header= control_header, energy_str = energy_str, energy_header = energy_header, t_axis = t_axis, energy_shift = energy_shift, e_axis = e_axis , use_discriminators = 0, $ shift_duration = shift_duration, silent=silent - data_level = 1 + ;; Select indices of the energy bins (among the 32) that are actually present in the pixel data science file + energy_bin_mask = data_str.energy_bin_mask + energy_bin_idx = where(energy_bin_mask eq 1) - start_time = atime(stx_time2any((t_axis.time_start)[0])) + energy_low = e_axis.LOW + energy_high = e_axis.HIGH + + energy_min = min(energy_low) + energy_max = max(energy_high) + + ;; Read BKG data + if keyword_set(fits_path_bk) then begin + + stx_read_pixel_data_fits_file, fits_path_bk, data_str = data_bkg, t_axis = t_axis_bkg, e_axis = e_axis_bkg, _extra=extra + + if n_elements(t_axis_bkg.DURATION) gt 1 then message, 'The chosen file does not contain a background measurement' + + ;; Select indices of the energy bins (among the 32) that are actually present in the pixel data bkg file + energy_bin_mask_bkg = data_bkg.energy_bin_mask + energy_bin_idx_bkg = where(energy_bin_mask_bkg eq 1) + + energy_low_bkg = e_axis_bkg.LOW + energy_high_bkg = e_axis_bkg.HIGH + + ;; Extract energy range in common between science and background file + energy_min = max([energy_min,min(energy_low_bkg)]) + energy_max = min([energy_max,max(energy_high_bkg)]) + idx_energy_bkg = where((energy_low_bkg ge energy_min) and (energy_high_bkg le energy_max)) + + energy_bin_idx_bkg = energy_bin_idx_bkg[idx_energy_bkg] + energy_low_bkg = energy_low_bkg[idx_energy_bkg] + energy_high_bkg = energy_high_bkg[idx_energy_bkg] - elut_filename = stx_date2elut_file(start_time) - - if ~keyword_set(silent) then begin - print, 'Using ELUT file ' + elut_filename endif - + + ;; Extract energy range in common between science and background file + idx_energy = where((energy_low ge energy_min) and (energy_high le energy_max)) + + energy_bin_idx = energy_bin_idx[idx_energy] + energy_low = energy_low[idx_energy] + energy_high = energy_high[idx_energy] + + ct_edges = get_uniq( [energy_low,energy_high],epsilon=0.0001) + + ;;***************** READ FITS INFO PARAMS + + data_level = 1 + uid = control_str.request_id if n_elements(distance) ne 0 then fits_distance = distance @@ -161,9 +202,12 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit distance = fits_distance, time_shift = time_shift, fits_path_bk = fits_path_bk, uid = uid, $ generate_fits = generate_fits, specfile = specfile, srmfile = srmfile, elut_file = elut_filename, silent = silent) - counts_in = data_str.counts + ;;***************** CHECK FOR POTENTIAL PIXEL SHADOWING + + counts = data_str.COUNTS + counts_error = data_str.COUNTS_ERR - dim_counts = counts_in.dim + dim_counts = counts.dim n_times = n_elements(dim_counts) gt 3 ? dim_counts[3] : 1 @@ -183,7 +227,6 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit endelse - pixel_mask_used = intarr(12) pixel_mask_used[pixels_used] = 1 n_pixels = total(pixel_mask_used) @@ -192,94 +235,280 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit detector_mask_used[detectors_used] = 1 n_detectors = total(detector_mask_used) + ;; Check if the top or bottom row of pixels is not fully-illuminated if ~keyword_set(silent) then begin - if total(pixel_mask_used[0:3]) eq total(pixel_mask_used[4:7]) then begin - count_ratio_threshold = 1.05 - counts_top = total(counts_in[1:25,0:3,detectors_used,*]) - counts_bottom = total(counts_in[1:25,4:7,detectors_used,*]) - case 1 of - f_div(counts_top, counts_bottom, default = 2) gt count_ratio_threshold : message, 'Top pixel total 5% higher than bottom row. Possible pixel shadowing. Recommend using only top pixels for analysis.',/info - f_div(counts_bottom, counts_top, default = 2) gt count_ratio_threshold : message, 'Bottom pixel total 5% higher than top row. Possible pixel shadowing. Recommend using only bottom pixels for analysis.',/info - else: - endcase + if total(pixel_mask_used[0:3]) eq total(pixel_mask_used[4:7]) then begin + count_ratio_threshold = 1.05 + counts_top = total(counts[1:25,0:3,detectors_used,*]) + counts_bottom = total(counts[1:25,4:7,detectors_used,*]) + case 1 of + f_div(counts_top, counts_bottom, default = 2) gt count_ratio_threshold : message, 'Top pixel total 5% higher than bottom row. Possible pixel shadowing. Recommend using only top pixels for analysis.',/info + f_div(counts_bottom, counts_top, default = 2) gt count_ratio_threshold : message, 'Bottom pixel total 5% higher than top row. Possible pixel shadowing. Recommend using only bottom pixels for analysis.',/info + else: + endcase + endif endif - endif + ;;***************** Compute count rates (normalize by livetime) - counts_in = reform(counts_in,[dim_counts[0:2], n_times]) + counts = counts[energy_bin_idx,*,*,*] + counts_error = counts_error[energy_bin_idx,*,*,*] - spec_in = total(reform(counts_in[*,pixels_used,detectors_used,*],[32,n_pixels,n_detectors,n_times]),2) + if keyword_set(fits_path_bk) then begin - spec_in = reform(spec_in,[dim_counts[0],n_detectors, n_times]) + ;; Check if science and background files are reconrded with the same ELUT + elut_filename = stx_date2elut_file(stx_time2any(t_axis.TIME_START[0])) + stx_read_elut, elut_gain, elut_offset, adc4096_str, elut_filename = elut_filename - counts_spec = spec_in[energy_bins,*, *] + elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) + stx_read_elut, ekev_actual = ekev_actual_bkg, elut_filename = elut_filename_bkg - counts_spec = reform(counts_spec,[n_energies, n_detectors, n_times]) + ;; Compare ELUT tables + elut_comp = STRCMP(elut_filename, elut_filename_bkg) - counts_err = reform(data_str.counts_err,[dim_counts[0:2], n_times]) + if not elut_comp then $ + message, 'The background file must be recorded when the same ELUT as the science file was uploaded. Please choose a different background file that is closer in time to the science file.' - counts_err = sqrt(total(reform(counts_err[*,pixels_used,detectors_used,*]^2.,[32,n_pixels,n_detectors,n_times]),2)) + counts_bkg = data_bkg.COUNTS + counts_error_bkg = data_bkg.COUNTS_ERR + counts_bkg = counts_bkg[energy_bin_idx_bkg,*,*] + counts_error_bkg = counts_error_bkg[energy_bin_idx_bkg,*,*] - counts_err = reform(counts_err,[dim_counts[0],n_detectors, n_times]) + endif else begin - counts_err = counts_err[energy_bins,*, *] - - counts_err = reform(counts_err,[n_energies, n_detectors, n_times]) + counts_bkg = dblarr(size(counts, /dim)) + counts_error_bkg = dblarr(size(counts, /dim)) - triggers = transpose(reform(data_str.triggers,[16, n_times])) + endelse - triggers_err = transpose(reform(data_str.triggers_err,[16, n_times])) + ;; Compute live time + live_time_data = stx_cpd_livetime(data_str.TRIGGERS, data_str.TRIGGERS_ERR, t_axis) + live_time_bins = live_time_data.LIVE_TIME_BINS + live_time_bins_error = live_time_data.LIVE_TIME_BINS_ERR + live_time_fraction_bins = live_time_data.LIVETIME_FRACTION - rcr = data_str.rcr - - if keyword_set(elut_correction) then begin + live_time_bins_rep = transpose(cmreplicate(live_time_bins, [n_elements(energy_bin_idx),12]), [2,3,0,1]) + live_time_bins_error_rep = transpose(cmreplicate(live_time_bins_error, [n_elements(energy_bin_idx),12]), [2,3,0,1]) + + if keyword_set(fits_path_bk) then begin + + live_time_bkg_data = stx_cpd_livetime(data_bkg.TRIGGERS, data_bkg.TRIGGERS_ERR, t_axis_bkg) + live_time_bkg = live_time_bkg_data.LIVE_TIME_BINS + live_time_error_bkg = live_time_bkg_data.LIVE_TIME_BINS_ERR + + endif else begin + + live_time_bkg = dblarr(32) + 1. + live_time_error_bkg = dblarr(32) + + endelse + + live_time_bkg_rep = transpose(cmreplicate(live_time_bkg, [n_elements(energy_bin_idx),12]), [1,2,0]) + live_time_error_bkg_rep = transpose(cmreplicate(live_time_error_bkg, [n_elements(energy_bin_idx),12]), [1,2,0]) + + ;; Normalize by livetime + count_rates = f_div( counts, live_time_bins_rep ) + count_rates_error = count_rates * sqrt( f_div( counts_error, counts )^2. + f_div( live_time_bins_error_rep, live_time_bins_rep )^2. ) + + count_rates_bkg = f_div( counts_bkg, live_time_bkg_rep ) + count_rates_bkg_error = count_rates_bkg * sqrt( f_div( counts_error_bkg, counts_bkg )^2. + f_div( live_time_error_bkg_rep, live_time_bkg_rep )^2. ) + + if n_times gt 1 then begin + + count_rates_bkg = cmreplicate( count_rates_bkg, n_times);[1,n_times] ) + count_rates_bkg_error = cmreplicate( count_rates_bkg_error, n_times);[1,n_times] ) + + endif + + ;; Apply BKG subtraction + count_rates = count_rates - count_rates_bkg + count_rates_error = sqrt( count_rates_error^2. + count_rates_bkg_error^2. ) + + ;;***************** APPLY ELUT CORRECTION + + if (total(pixel_mask_used[0:3]) gt 0.) and (total(pixel_mask_used[4:7]) gt 0.) $ + and (total(pixel_mask_used[8:11]) gt 0.) then sumcase = 'ALL' + + if (total(pixel_mask_used[0:3]) gt 0.) and (total(pixel_mask_used[4:7]) gt 0.) $ + and (total(pixel_mask_used[8:11]) eq 0.) then sumcase = 'TOP+BOT' + + if (total(pixel_mask_used[0:3]) gt 0.) and (total(pixel_mask_used[4:7]) eq 0.) $ + and (total(pixel_mask_used[8:11]) eq 0.) then sumcase = 'TOP' + + if (total(pixel_mask_used[0:3]) eq 0.) and (total(pixel_mask_used[4:7]) gt 0.) $ + and (total(pixel_mask_used[8:11]) eq 0.) then sumcase = 'BOT' + + if (total(pixel_mask_used[0:3]) eq 0.) and (total(pixel_mask_used[4:7]) eq 0.) $ + and (total(pixel_mask_used[8:11]) gt 0.) then sumcase = 'SMALL' + + case sumcase of + + 'TOP': begin + pixel_ind = [0] + end + + 'BOT': begin + pixel_ind = [1] + end + + 'TOP+BOT': begin + pixel_ind = [0,1] + end + + 'ALL': begin + pixel_ind = [0,1,2] + end + + 'SMALL': begin + pixel_ind = [2] + end + end + + + + if keyword_set(calib_data) then begin + + ;; Create daily ELUT + energy_bin_low = calib_data.ENERGY_BIN_LOW + energy_bin_high = calib_data.ENERGY_BIN_HIGH + + energy_bin_low = energy_bin_low[energy_bin_idx,*,*] + energy_bin_high = energy_bin_high[energy_bin_idx,*,*] + + ; Rebin counts in time. From stx_science_data_lightcurve: + ; determine time bins with minimum duration - keep adding consecutive bins until the minimum + ; value is at least reached + + duration = t_axis.DURATION + + i=0 + j=0 + total_time=0 + iall=[] + + while (i lt n_elements(duration)-1) do begin + while (total_time lt delta_time_min) and (i+j le n_elements(duration)-1) do begin + total_time = total(duration[i:i+j]) + j++ + endwhile + iall = [iall,i] + i = i+j + j = 0 + total_time = 0 + endwhile + + idx_time_min = iall[0:-2] + idx_time_max = iall[1:-1]-1 + + + + count_rates_elut = fltarr(count_rates.dim) + count_rates_error_elut = fltarr(count_rates_error.dim) + + n_times_rebinned = n_elements(idx_time_min) + + for t_bin = 0,n_times_rebinned-1 do begin + + this_count_rates = reform(count_rates[*,*,*,idx_time_min[t_bin]:idx_time_max[t_bin]]) + this_count_rates_error = reform(count_rates_error[*,*,*,idx_time_min[t_bin]:idx_time_max[t_bin]]) + + if idx_time_max[t_bin]-idx_time_min[t_bin] ge 1 then begin + + rebinned_count_rates = average(this_count_rates, 4) + + endif else begin + + rebinned_count_rates = this_count_rates + endelse - stx_read_elut, ekev_actual = ekev_actual, elut_filename = elut_filename + if n_elements(pixels_used) gt 1 then begin - ave_edge = mean(reform(ekev_actual[energy_edges_used-1, pixels_used, detectors_used, 0 ],n_energy_edges, n_pixels, n_detectors), dim= 2) - ave_edge = mean(reform(ave_edge,n_energy_edges, n_detectors), dim= 2) + spectrum = total(rebinned_count_rates[*,pixels_used,*], 2) - edge_products, ave_edge, width = ewidth + endif else begin - eff_ewidth = (e_axis.width)/ewidth + spectrum = reform(rebinned_count_rates[*,pixels_used,*]) + endelse - counts_spec = counts_spec * reform(reproduce(eff_ewidth, n_detectors*n_times),n_energies, n_detectors, n_times) + if n_elements(detectors_used) gt 1 then begin - counts_spec = reform(counts_spec,[n_energies, n_detectors, n_times]) + spectrum = total(spectrum[*,detectors_used], 2) + endif else begin - counts_err = counts_err * reform(reproduce(eff_ewidth, n_detectors*n_times),n_energies, n_detectors, n_times) + spectrum = reform(spectrum[*,detectors_used]) + + endelse + + spectrum = spectrum / (energy_high - energy_low) + + ;; Apply ELUT correction + for e_bin=0,n_elements(idx_energy)-1 do begin + + this_energy_range = [energy_low[e_bin], energy_high[e_bin]] + + elut_data = stx_elut_correction(this_count_rates, this_count_rates_error, $ + energy_bin_idx, energy_bin_low, energy_bin_high, energy_high, energy_low, e_bin, this_energy_range, $ + spectrum, pixel_ind, det_ind, /silent) + + count_rates_elut[e_bin,*,*,idx_time_min[t_bin]:idx_time_max[t_bin]] = elut_data.COUNTS + count_rates_error_elut[e_bin,*,*,idx_time_min[t_bin]:idx_time_max[t_bin]] = elut_data.COUNTS_ERROR + + endfor + + endfor + + count_rates = count_rates_elut + count_rates_error = count_rates_error_elut - counts_err = reform(counts_err,[n_energies, n_detectors, n_times]) - endif - ;insert the information from the telemetry file into the expected stx_fsw_sd_spectrogram structure - spectrogram = { $ - type : "stx_fsw_sd_spectrogram", $ - counts : counts_spec, $ - trigger : triggers, $ - trigger_err : triggers_err, $ - time_axis : t_axis , $ - energy_axis : e_axis, $ - pixel_mask : pixel_mask_used , $ - detector_mask : detector_mask_used, $ - rcr : rcr, $ - error : counts_err} - - data_dims = lonarr(4) - data_dims[0] = n_energies - data_dims[1] = n_detectors - data_dims[2] = n_pixels - data_dims[3] = n_times + ;;***************** CREATE SPECTROGRAM + + if n_elements(pixels_used) gt 1 then begin + + spec = total(count_rates[*,pixels_used,*,*], 2) + spec_error = sqrt(total(count_rates_error[*,pixels_used,*,*]^2., 2)) + + endif else begin + + spec = reform(count_rates[*,pixels_used,*,*]) + spec_error = reform(count_rates_error[*,pixels_used,*,*]) + + endelse + + if n_elements(detectors_used) gt 1 then begin + + spec = total(spec[*,detectors_used,*], 2) + spec_error = sqrt(total(spec_error[*,detectors_used,*]^2., 2)) + + endif else begin + + spec = reform(spec[*,detectors_used,*]) + spec_error = reform(spec_error[*,detectors_used,*]) + + endelse + + ;; Compute average livetime and livetime fraction + avg_live_time_bins = average(live_time_bins[detectors_used,*], 1) + avg_live_time_fraction_bins = average(live_time_fraction_bins[detectors_used,*], 1) + + ;; Multiply by average live time. The units of the spectrogram are counts + avg_live_time_bins_rep = transpose(cmreplicate(avg_live_time_bins, n_elements(energy_bin_idx))) + + spec *= avg_live_time_bins_rep + spec_error *= avg_live_time_bins_rep + + ;;***************** CHECK FOR RCR CHANGES + + rcr = data_str.rcr ;get the rcr states and the times of rcr changes from the ql_lightcurves structure ut_rcr = stx_time2any(t_axis.time_end) find_changes, rcr, index, state, count=count - ; ************************************************************ ; ******************** TEMPORARY FIX ************************* ; ***** Andrea: 2022-April-05 @@ -303,8 +532,8 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit ; large change in counts in the counts of the 5 - 6 keV energy bin. ; find all time intervals where the difference between adjacent bins is large if max(rcr) gt 0 then begin; skip if in the standard state of RCR0 for the full time range - - jumps = where(abs((total(counts_spec,2))[2,*] - shift((total(counts_spec,2))[2,*],-1)) gt 1e4) + jumps = where(abs(total(spec,1) - shift(total(spec,1),-1)) * 24. / n_elements(detectors_used) gt 1e4) + ;jumps = where(abs(total(spec,1) - shift(total(spec,1),-1)) gt 1e4) ; include the starting bin jumps = [0, jumps] ; as the attenuator motion can be present in two consecutive bins select only the first @@ -316,16 +545,157 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit index = jumps_use[closest_jumps] endif - ; ************************************************************ ;add the rcr information to a specpar structure so it can be included in the spectrum FITS file specpar = { sp_atten_state : {time:ut_rcr[index], state:state}, flare_xyoffset : fltarr(2), use_flare_xyoffset:0 } - stx_convert_science_data2ospex, spectrogram = spectrogram, specpar=specpar, time_shift = time_shift, $ - data_level = data_level, data_dims = data_dims, fits_path_bk = fits_path_bk, fits_path_data = fits_path_data,$ - aux_fits_file = aux_fits_file, flare_location_hpc = flare_location_hpc, flare_location_stx = flare_location_stx, $ - eff_ewidth = eff_ewidth, sys_uncert = sys_uncert, plot = plot, background_data = background_data, silent = silent, $ - elut_correction = elut_correction, fits_info_params = fits_info_params, xspec=xspec, ospex_obj = ospex_obj + ;;***************** CREATE SRM + + pixel_mask =detector_mask_used ## pixel_mask_used + + transmission = read_csv(loc_file( 'stix_transmission_highres_20251110.csv', path = getenv('STX_GRID'))) + + emin = 4 + emax = 150 + phe = transmission.(0) + phe = phe[where(phe gt emin-1 and phe lt 3.5*emax)] + edge_products, phe, mean = mean_phe, width = w_phe + ph_edges = [mean_phe[0] - w_phe[0], mean_phe] + + distance = fits_info_params.distance + dist_factor = 1./(distance^2.) + + ;make the srm for the appropriate pixel mask and energy edges + ;srm = stx_build_pixel_drm(ct_edges, pixel_mask, ph_energy_edges = ph_edges, dist_factor = dist_factor, tailing = tailing, include_damage = include_damage, _extra=extra) + ph_edges = get_uniq( [ph_edges,ct_edges],epsilon=0.0001) + + ;; Compute subc. transmission + edge_products,ph_edges, mean=ph_in + subc_transmission = stx_subc_transmission(flare_location_stx, ph_in, _extra=extra) + + if n_elements(detectors_used) eq 1 then begin + grid_factor = reform(subc_transmission[*,detectors_used]) + endif else begin + grid_factor = total(subc_transmission[*,detectors_used],2) / n_elements(detectors_used) + endelse + + + ;; Check if BKG or CFL detectors are used + + idx_cfl = where(detectors_used eq 8, n_cfl) + idx_bkg = where(detectors_used eq 9, n_bkg) + + if n_cfl eq 1 then message, "CFL detector can not be selected for spectral fitting." + if (n_elements(detectors_used) gt 1) and (n_bkg eq 1) then message, "BKG detector can not be selected together with imaging detectors for spectral fitting." + + if n_elements(detectors_used) eq 1 then if detectors_used eq 9 then begin + grid_transmission_file = concat_dir(getenv('STX_GRID'), 'real_bkg_grid_transmission.txt') + readcol, grid_transmission_file, bk_grid_factors, format = 'f', skip = 2, silent = silent + + grid_factor = average(bk_grid_factors[pixels_used]) + + endif + + ;; Creates appropriate SRM for different attenuator states + rcr_states = specpar.sp_atten_state.state + rcr_states = rcr_states[uniq(rcr_states, sort(rcr_states))] + nrcr_states = n_elements(rcr_states) + + srm_atten = replicate( {rcr:0, srm:fltarr(n_elements(ct_edges)-1,n_elements(ph_edges)-1)},nrcr_states ) + + for i =0, nrcr_states-1 do begin + ;make the srm for the appropriate pixel mask and energy edges + rcr = rcr_states[i] + + srm = stx_build_pixel_drm(ct_edges, pixel_mask, rcr = rcr, grid_factor= grid_factor, ph_energy_edges = ph_edges, dist_factor = dist_factor, tailing = tailing, include_damage = include_damage, _extra=extra) + srm_atten[i].srm = srm.smatrix + srm_atten[i].rcr = rcr + + endfor + + ;;***************** SAVE FITS + detector_label = stx_det_mask2label(detector_mask_used) + pixel_label = stx_pix_mask2label(pixel_mask_used) + + ospex_obj = ospex(/no) + + ;if the fits keyword is set write the spectrogram and srm data to fits files and then read them in to the ospex object + if fits_info_params.generate_fits eq 1 then begin + utime = transpose([stx_time2any( t_axis.time_start )]) + + ;spectrogram structure for passing to fits writer routine + spectrum_in = { type : 'stx_spectrogram', $ + data : spec, $ + t_axis : t_axis, $ + e_axis : e_axis, $ + ltime : avg_live_time_fraction_bins, $ + attenuator_state : data_str.rcr , $ + error : spec_error } + + specfilename = fits_info_params.specfile + srmfilename = fits_info_params.srmfile + + fits_info_params.grid_factor.add, grid_factor + fits_info_params.detused = detector_label + ', Pixels: ' + pixel_label + + if keyword_set(xspec) then begin + ;xspec in general works with energy depandent systematic errors + e_axis = spectrum_in.e_axis + n_energies = n_elements(e_axis.mean) + sys_err = fltarr(n_energies) + + idx_below10kev = where(e_axis.mean lt 10, cb10) + sys_err[*] = 0.03 + if cb10 gt 0 then sys_err[idx_below10kev] = 0.05 + idx_below7kev = where(e_axis.mean lt 7, cb7) + if cb7 gt 0 then sys_err[idx_below7kev] = 0.07 + + sys_err = rebin(sys_err, n_energies,n_times) + endif + + stx_write_ospex_fits, spectrum = spectrum_in, srmdata = srm, specpar = specpar, time_shift = time_shift, $ + srm_atten = srm_atten, specfilename = specfilename, srmfilename = srmfilename, ph_edges = ph_edges, $ + fits_info_params = fits_info_params, xspec = xspec, silent = silent + + ospex_obj->set, spex_file_reader = 'stx_read_sp' + ospex_obj->set, spex_specfile = specfilename ; name of your spectrum file + ospex_obj->set, spex_drmfile = srmfilename + + endif else begin + ;if the generate_fits keyword is not set use the spex_user_data strategy to pass in the data directly to the ospex object + + energy_edges = e_axis.edges_2 + Edge_Products, ph_edges, edges_2 = ph_edges2 + + utime2 = transpose(stx_time2any( [[t_axis.time_start], [t_axis.time_end]] )) + + ospex_obj->set, spex_data_source = 'spex_user_data' + ospex_obj->set, spectrum = float(spec), $ + spex_ct_edges = energy_edges, $ + spex_ut_edges = utime2, $ + livetime = avg_live_time_bins_rep, $ + errors = spec_error + srm = rep_tag_name(srm,'smatrix','drm') + ospex_obj->set, spex_respinfo = srm + ospex_obj->set, spex_area = srm.area + ospex_obj->set, spex_detectors = 'STIX' + ospex_obj->set, spex_drm_ct_edges = energy_edges + ospex_obj->set, spex_drm_ph_edges = ph_edges2 + endelse + + ospex_obj->set, spex_uncert = sys_uncert + ospex_obj->set, spex_error_use_expected = 0 + + counts_str = ospex_obj->getdata(spex_units='counts') + origunits = ospex_obj->get(/spex_data_origunits) + origunits.data_name = 'STIX' + ospex_obj->set, spex_data_origunits = origunits + + if keyword_set(plot) then begin + ospex_obj ->gui + ospex_obj ->set, spex_eband = get_edges([4.,10.,15.,25, 50, 84.], /edges_2) + ospex_obj ->plot_time, spex_units='flux', /show_err, obj = plotman_object + endif end diff --git a/stix/idl/processing/spectrogram/stx_elut_correction.pro b/stix/idl/processing/spectrogram/stx_elut_correction.pro index 070c5312..7384e38e 100644 --- a/stix/idl/processing/spectrogram/stx_elut_correction.pro +++ b/stix/idl/processing/spectrogram/stx_elut_correction.pro @@ -41,6 +41,7 @@ ; It is used to compute the amount of ELUT correction for all pixels. ; - COUNTS_ELUT: array of dimension 4 (A,B,C,D) x number of rows (TOP, BOT, SMALL) x number of subcollimators x number of time bins containing the number of counts integrated in energy after ELUT correction. ; It is used to compute the amount of ELUT correction for all pixels. +; - SP_INDEX: array containing the value of the spectral indices in the different energy bins ; ; KEYWORDS: ; @@ -49,6 +50,7 @@ ; SILENT: if set to 1, the spectrum plot is not displayed. ; ; HISTORY: +; March 2026, Massa P., first release ; ; CONTACT: ; paolo.massa@fhnw.ch From 42585b82c19b1cc0d49e7f0e789793f58ed60eac Mon Sep 17 00:00:00 2001 From: paolomassa Date: Tue, 24 Mar 2026 18:01:58 +0100 Subject: [PATCH 03/13] First version of the function for automatic download of the calibration file. --- stix/idl/demo/stx_imaging_demo.pro | 32 +++- stix/idl/io/stx_get_calibration_file.pro | 215 +++++++++++++++++++++++ 2 files changed, 244 insertions(+), 3 deletions(-) create mode 100644 stix/idl/io/stx_get_calibration_file.pro diff --git a/stix/idl/demo/stx_imaging_demo.pro b/stix/idl/demo/stx_imaging_demo.pro index 15018d96..4d75ab56 100644 --- a/stix/idl/demo/stx_imaging_demo.pro +++ b/stix/idl/demo/stx_imaging_demo.pro @@ -73,6 +73,31 @@ aux_data = stx_create_auxiliary_data(aux_fits_file, time_range) stop +;************************************* DOWNLOAD CALIBRATION FITS FILE ********************************* + +; This file is used for creating a structure containing daily calibration data to be used for ELUT +; correction (see next section) +path_calib_file = stx_get_calibration_file(time_range[0], time_range[1], out_dir=out_dir) + +stop + +;;********************************** CONSTRUCT CALIBRATION DATA STRUCTURE ***************************** + +; Create a structure containing daily calibration data to use for ELUT correction: +; - TIME_START: start time of the calibration file which is used to derive gain/offset values +; - TIME_END: end time of the calibration file which is used to derive gain/offset values +; - T_MEAN: mid point of the time range of the calibration file which is used to derive gain/offset values +; - ENERGY_BIN_LOW: actual lower energy edges (in keV) of the science channels for each pixel +; - ENERGY_BIN_HIGH: actual higher energy edges (in keV) of the science channels for each pixel +; - GAIN: gain value for each pixel +; - OFFSET: offset value for each pixel +; - LIVE_TIME: live time of the calibration file +; - ELUT_NAME: name of the ELUT that was utilized onboard when the calibration file was recorded + +stx_read_calibration_data, path_calib_file, calib_data=calib_data + +stop + ;*************************************** ESTIMATE FLARE LOCATION ************************************** ; Returns the coordinates of the estimated flare location (arcsec, Helioprojective Cartesian coordinates @@ -80,7 +105,7 @@ stop ; center of the maps to be reconstructed stx_estimate_flare_location, path_sci_file, time_range, aux_data, flare_loc=flare_loc, $ - path_bkg_file=path_bkg_file + path_bkg_file=path_bkg_file, calib_data=calib_data stop @@ -102,7 +127,7 @@ xy_flare = mapcenter ; '6a','6b','6c','5a','5b','5c','4a','4b','4c','3a','3b','3c']) vis=stx_construct_calibrated_visibility(path_sci_file, time_range, energy_range, mapcenter, subc_index=subc_index, $ - path_bkg_file=path_bkg_file, xy_flare=xy_flare) + path_bkg_file=path_bkg_file, xy_flare=xy_flare, calib_data=calib_data) stop @@ -178,7 +203,8 @@ stop ; Expectation Maximization algorithm from STIX counts (see Massa P. et al (2019) for details). Takes as input a ; summed pixel data structure -pixel_data_summed = stx_construct_pixel_data_summed(path_sci_file, time_range, energy_range, path_bkg_file=path_bkg_file, /silent) +pixel_data_summed = stx_construct_pixel_data_summed(path_sci_file, time_range, energy_range, path_bkg_file=path_bkg_file, $ + calib_data=calib_data, /silent) em_map = stx_em(pixel_data_summed, aux_data, imsize=imsize, pixel=pixel, xy_flare=xy_flare, mapcenter=mapcenter) diff --git a/stix/idl/io/stx_get_calibration_file.pro b/stix/idl/io/stx_get_calibration_file.pro new file mode 100644 index 00000000..ddc41277 --- /dev/null +++ b/stix/idl/io/stx_get_calibration_file.pro @@ -0,0 +1,215 @@ +;+ +; +; name: +; stx_get_calibration_file +; +; :description: +; This procedure checks the STIX data archive for a given observation time and finds any +; STIX CAL calibration file. If no file is present for the input date, the procedure searches for +; the calibration file which is closest in time. If multiple files are present for that date, +; this procedure downloads the one with largest duration. +; +; :categories: +; template, example +; +; :params: +; start_time : in, required, type="string" +; the start time of the observation +; end_time : in, required, type="string" +; the end time of the observation +; +; :keywords +; out_dir: path of the folder where the STIX calibration FITS files are saved. Default is the current directory +; +; clobber: 0 or 1. If set to 0, the code does not download the file again if it is already present in 'out_dir'. +; +; :returns: +; Path of the downloaded STIX calibration FITS file +; +; :examples: +; out_file = stx_get_calibration_file('09-May-23 06:14:37.094', '09-May-23 06:36:12.194') +; +; :history: +; 24-Mar-2026 - Massa P. (FHNW), first release +;- +function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobber=clobber + + cd, current=current + + default, out_dir, current + default, clobber, 0 + + day_in_s = 86400.d ;; Number of seconds in a day. To be used later to identify the most appropriate calibration file + + site = 'http://dataarchive.stix.i4ds.net' + date_path = get_fid(start_time,end_time,/full,delim='/') + + type_path = '/fits/CAL/' + path = type_path + date_path[0] +'/CAL' + filter = '*stix-cal-energy*.fits' + found_files=sock_find(site,filter,path=path,count=count) + + if count eq 0 then begin + message, "No STIX calibration files found for this day ("+date_path[0]+"). Searching for calibration file which is closest in time.", /continue + + ;; Forward in time (for 30 days max - arbitrary choice) + day_n_forward=1 + count_forward=0 + while (day_n_forward le 30) and (count_forward eq 0) do begin + + this_start_time_forward = anytim(anytim(start_time) + day_n_forward * day_in_s, /vms) + this_end_time_forward = anytim(anytim(end_time) + day_n_forward * day_in_s, /vms) + + date_path_forward = get_fid(this_start_time_forward,this_end_time_forward,/full,delim='/') + + path = type_path + date_path_forward[0] +'/CAL' + found_files_forward = sock_find(site,filter,path=path,count=count_forward) + + day_n_forward += 1 + + endwhile + + ;; Backward in time (for 30 days max - arbitrary choice) + day_n_backward=1 + count_backward=0 + while (day_n_backward le 30) and (count_backward eq 0) do begin + + this_start_time_backward = anytim(anytim(start_time) - day_n_backward * day_in_s, /vms) + this_end_time_backward = anytim(anytim(end_time) - day_n_backward * day_in_s, /vms) + + date_path_backward = get_fid(this_start_time_backward,this_end_time_backward,/full,delim='/') + + path = type_path + date_path_backward[0] +'/CAL' + found_files_backward = sock_find(site,filter,path=path,count=count_backward) + + day_n_backward += 1 + + endwhile + + if (count_forward eq 0) and (count_backward eq 0) then begin + + message, "No STIX calibration file was found in the month before or after " +date_path[0]+". Please, download a calibration file manually." + + endif else begin + + diff_days = day_n_forward - day_n_backward + ;; if the difference is negative, it means that the closest file is the one found forward in time + ;; if the difference is positive, it means that the closest file is the one found backward in time + ;; if the difference is zero, we arbitrarily select the file found forward in time + + if diff_days le 0 then begin + + count = count_forward + found_files = found_files_forward + date_path = date_path_forward + + endif else begin + + count = count_backward + found_files = found_files_backward + date_path = date_path_backward + + endelse + + endelse + + endif + + message, "Download STIX calibration file recorded at " +date_path[0], /continue + + len_path = STRLEN(site+path) + + if count eq 1 then begin + + selected_file = found_files[0] + + ;; Check that the ELUT that was onboard when the calibration file was recorded is the same as the on that was on board + ;; during the input time range. If not, raise an error + + ;; Find start time of the calibration file + len_full_path = STRLEN(found_files[0]) + filename = STRMID(found_files[0], len_path+1, len_full_path) + + string_dates = STRMID(filename, STRLEN('solo_CAL_stix-cal-energy_'), STRLEN(filename)-STRLEN('solo_CAL_stix-cal-energy_')-9) ;; 9 is the number of characters in the string '_V02.fits' + + mid_part = STRPOS(string_dates, '-') + + this_start_time_file = STRMID(string_dates, 0, mid_part) + year = STRMID(this_start_time_file, 0, 4) + month = STRMID(this_start_time_file, 4, 2) + day = STRMID(this_start_time_file, 6, 2) + hour = STRMID(this_start_time_file, 9, 2) + min = STRMID(this_start_time_file, 11, 2) + sec = STRMID(this_start_time_file, 13, 2) + start_time_file = year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec + + ;; Compare ELUT + elut_file = stx_date2elut_file(start_time_file) + elut_time_range = stx_date2elut_file(start_time) + if elut_file ne elut_time_range then message, $ + "The closest-in-time calibration file was acquired with an onboard ELUT that differs from the one used during the selected time range. Please, manually download a different calibration file." + + ;; Download file + sock_copy, selected_file, out_name, local_file=out_file, out_dir = out_dir, clobber=clobber + + endif else begin + + message, "More than 1 STIX calibration file found for this day ("+date_path[0]+").", /continue + + start_time_file = [] + end_time_file = [] + + ;; Extract file names + for i = 0,n_elements(found_files)-1 do begin + + len_full_path = STRLEN(found_files[i]) + filename = STRMID(found_files[i], len_path+1, len_full_path) + + string_dates = STRMID(filename, STRLEN('solo_CAL_stix-cal-energy_'), STRLEN(filename)-STRLEN('solo_CAL_stix-cal-energy_')-9) ;; 9 is the number of characters in the string '_V02.fits' + + mid_part = STRPOS(string_dates, '-') + + this_start_time_file = STRMID(string_dates, 0, mid_part) + year = STRMID(this_start_time_file, 0, 4) + month = STRMID(this_start_time_file, 4, 2) + day = STRMID(this_start_time_file, 6, 2) + hour = STRMID(this_start_time_file, 9, 2) + min = STRMID(this_start_time_file, 11, 2) + sec = STRMID(this_start_time_file, 13, 2) + start_time_file = [start_time_file, anytim(year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec)] + + this_end_time_file = STRMID(string_dates, mid_part+1, STRLEN(string_dates)) + year = STRMID(this_end_time_file, 0, 4) + month = STRMID(this_end_time_file, 4, 2) + day = STRMID(this_end_time_file, 6, 2) + hour = STRMID(this_end_time_file, 9, 2) + min = STRMID(this_end_time_file, 11, 2) + sec = STRMID(this_end_time_file, 13, 2) + end_time_file = [end_time_file, anytim(year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec)] + + endfor + + ;; Select file with largest duration + duration = end_time_file - start_time_file + + idx_file = where(duration eq max(duration)) + selected_file = found_files[idx_file] + + ;; Check that the ELUT that was onboard when the calibration file was recorded is the same as the on that was on board + ;; during the input time range. If not, raise an error + + ;; Compare ELUT + elut_file = stx_date2elut_file(start_time_file[idx_file]) + elut_time_range = stx_date2elut_file(start_time) + if elut_file ne elut_time_range then message, $ + "The closest-in-time calibration file was acquired with an onboard ELUT that differs from the one used during the selected time range. Please, manually download a different calibration file." + + sock_copy, selected_file, out_name, local_file=out_file, out_dir = out_dir, clobber=clobber + + endelse + + return, out_file + + + +end \ No newline at end of file From 58ed4ab0cbf754ba85a2b1b8b64995523f0bfeb8 Mon Sep 17 00:00:00 2001 From: Laszlo Etesi Date: Mon, 30 Mar 2026 12:45:40 +0200 Subject: [PATCH 04/13] Review cleanup: remove .DS_Store, fix typos, add trailing newlines MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Remove committed .DS_Store and add it to .gitignore - Fix date typos: "April 22026" → "April 2026", "Fubruary" → "February" - Add missing trailing newlines to five new files --- .gitignore | 3 +++ stix/idl/io/stx_get_calibration_file.pro | 2 +- stix/idl/processing/calib_data/.DS_Store | Bin 6148 -> 0 bytes .../calib_data/stx_calibration_data.pro | 2 +- stix/idl/processing/imaging/stx_em.pro | 2 +- stix/idl/processing/livetime/stx_triggergram.pro | 2 +- .../spectrogram/stx_elut_correction.pro | 2 +- .../spectrogram/stx_estimate_spectral_index.pro | 2 +- stix/idl/structures/stx_pixel_data.pro | 2 +- 9 files changed, 10 insertions(+), 7 deletions(-) delete mode 100644 stix/idl/processing/calib_data/.DS_Store diff --git a/.gitignore b/.gitignore index 8e502732..657ba3d7 100644 --- a/.gitignore +++ b/.gitignore @@ -40,3 +40,6 @@ stx_scenario_* .vscode +# macOS metadata +.DS_Store + diff --git a/stix/idl/io/stx_get_calibration_file.pro b/stix/idl/io/stx_get_calibration_file.pro index ddc41277..748f420c 100644 --- a/stix/idl/io/stx_get_calibration_file.pro +++ b/stix/idl/io/stx_get_calibration_file.pro @@ -212,4 +212,4 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe -end \ No newline at end of file +end diff --git a/stix/idl/processing/calib_data/.DS_Store b/stix/idl/processing/calib_data/.DS_Store deleted file mode 100644 index 5008ddfcf53c02e82d7eee2e57c38e5672ef89f6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0 Date: Mon, 13 Apr 2026 10:31:07 +0200 Subject: [PATCH 05/13] Update to search for calibration files within 2 days before or after the input date --- stix/idl/io/stx_get_calibration_file.pro | 63 ++++++++++++------------ 1 file changed, 32 insertions(+), 31 deletions(-) diff --git a/stix/idl/io/stx_get_calibration_file.pro b/stix/idl/io/stx_get_calibration_file.pro index 748f420c..a3cce056 100644 --- a/stix/idl/io/stx_get_calibration_file.pro +++ b/stix/idl/io/stx_get_calibration_file.pro @@ -31,6 +31,7 @@ ; ; :history: ; 24-Mar-2026 - Massa P. (FHNW), first release +; 13-Apr-2026 - Massa P. (FHNW), updated to search for calibration files within 2 days before or after the input date ;- function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobber=clobber @@ -40,6 +41,7 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe default, clobber, 0 day_in_s = 86400.d ;; Number of seconds in a day. To be used later to identify the most appropriate calibration file + day_limit = 2 ;; Maximum number of days before or after the input date to consider when searching for a calibration file (arbitrary choice) site = 'http://dataarchive.stix.i4ds.net' date_path = get_fid(start_time,end_time,/full,delim='/') @@ -50,12 +52,12 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe found_files=sock_find(site,filter,path=path,count=count) if count eq 0 then begin - message, "No STIX calibration files found for this day ("+date_path[0]+"). Searching for calibration file which is closest in time.", /continue + message, [" ", " ", "No STIX calibration files found for this day ("+date_path[0]+"). Searching for calibration file which is closest in time.", " ", " "], /continue - ;; Forward in time (for 30 days max - arbitrary choice) + ;; Forward in time (for 2 days max - arbitrary choice) day_n_forward=1 count_forward=0 - while (day_n_forward le 30) and (count_forward eq 0) do begin + while (day_n_forward le day_limit) and (count_forward eq 0) do begin this_start_time_forward = anytim(anytim(start_time) + day_n_forward * day_in_s, /vms) this_end_time_forward = anytim(anytim(end_time) + day_n_forward * day_in_s, /vms) @@ -69,10 +71,12 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe endwhile - ;; Backward in time (for 30 days max - arbitrary choice) + day_n_forward -= 1 ;; Remove 1 day to compensate that when the loop exits after finding a file, each counter is +1 from the actual day offset that produced the match. + + ;; Backward in time (for 2 days max - arbitrary choice) day_n_backward=1 count_backward=0 - while (day_n_backward le 30) and (count_backward eq 0) do begin + while (day_n_backward le day_limit) and (count_backward eq 0) do begin this_start_time_backward = anytim(anytim(start_time) - day_n_backward * day_in_s, /vms) this_end_time_backward = anytim(anytim(end_time) - day_n_backward * day_in_s, /vms) @@ -86,9 +90,11 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe endwhile + day_n_backward -= 1 ;; Remove 1 day to compensate that when the loop exits after finding a file, each counter is +1 from the actual day offset that produced the match. + if (count_forward eq 0) and (count_backward eq 0) then begin - message, "No STIX calibration file was found in the month before or after " +date_path[0]+". Please, download a calibration file manually." + message, [" ", " ", "No STIX calibration file was found within 2 days before or after " +date_path[0]+". Please, download a calibration file manually.", " ", " "] endif else begin @@ -115,7 +121,7 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe endif - message, "Download STIX calibration file recorded at " +date_path[0], /continue + message, [" ", " ", "Download STIX calibration file recorded at " +date_path[0], " ", " "], /continue len_path = STRLEN(site+path) @@ -123,10 +129,8 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe selected_file = found_files[0] - ;; Check that the ELUT that was onboard when the calibration file was recorded is the same as the on that was on board - ;; during the input time range. If not, raise an error + ;; Extract the calibration file start time from the filename. This is used below to verify that the ELUT onboard at the time of calibration matches the one used during the input time range. - ;; Find start time of the calibration file len_full_path = STRLEN(found_files[0]) filename = STRMID(found_files[0], len_path+1, len_full_path) @@ -142,19 +146,10 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe min = STRMID(this_start_time_file, 11, 2) sec = STRMID(this_start_time_file, 13, 2) start_time_file = year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec - - ;; Compare ELUT - elut_file = stx_date2elut_file(start_time_file) - elut_time_range = stx_date2elut_file(start_time) - if elut_file ne elut_time_range then message, $ - "The closest-in-time calibration file was acquired with an onboard ELUT that differs from the one used during the selected time range. Please, manually download a different calibration file." - - ;; Download file - sock_copy, selected_file, out_name, local_file=out_file, out_dir = out_dir, clobber=clobber endif else begin - message, "More than 1 STIX calibration file found for this day ("+date_path[0]+").", /continue + message, [" ", " ", "More than 1 STIX calibration file found for this day ("+date_path[0]+"). Download the one with largest duration.", " ", " "], /continue start_time_file = [] end_time_file = [] @@ -195,19 +190,25 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe idx_file = where(duration eq max(duration)) selected_file = found_files[idx_file] - ;; Check that the ELUT that was onboard when the calibration file was recorded is the same as the on that was on board - ;; during the input time range. If not, raise an error - - ;; Compare ELUT - elut_file = stx_date2elut_file(start_time_file[idx_file]) - elut_time_range = stx_date2elut_file(start_time) - if elut_file ne elut_time_range then message, $ - "The closest-in-time calibration file was acquired with an onboard ELUT that differs from the one used during the selected time range. Please, manually download a different calibration file." - - sock_copy, selected_file, out_name, local_file=out_file, out_dir = out_dir, clobber=clobber + start_time_file = start_time_file[idx_file] endelse - + + + ;; Check that the ELUT onboard when the calibration file was recorded matches the ELUT onboard during the input time range; otherwise, raise an error. + + ;; Compare ELUT + elut_file = stx_date2elut_file(start_time_file) + elut_time_range = stx_date2elut_file(start_time) + if elut_file ne elut_time_range then message, $ + [" ", " ", "The closest-in-time calibration file was acquired with an onboard ELUT that differs from the one used during the selected time range. Please, manually download a different calibration file.", " ", " "] + + sock_copy, selected_file, out_name, local_file=out_file, out_dir = out_dir, clobber=clobber + + filename = STRMID(out_file, STRLEN(out_dir), STRLEN(out_file)-STRLEN(out_dir)) + + message, [" ", " ", "Downloaded file: " + filename, " ", " "], /continue + return, out_file From 470e80e5ea578bb9248489ae286bb4fb74541817 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Mon, 27 Apr 2026 10:42:03 +0200 Subject: [PATCH 06/13] Update: select files with duration larger than half a day --- stix/idl/io/stx_get_calibration_file.pro | 228 +++++++++++------------ 1 file changed, 109 insertions(+), 119 deletions(-) diff --git a/stix/idl/io/stx_get_calibration_file.pro b/stix/idl/io/stx_get_calibration_file.pro index a3cce056..68729f66 100644 --- a/stix/idl/io/stx_get_calibration_file.pro +++ b/stix/idl/io/stx_get_calibration_file.pro @@ -31,7 +31,7 @@ ; ; :history: ; 24-Mar-2026 - Massa P. (FHNW), first release -; 13-Apr-2026 - Massa P. (FHNW), updated to search for calibration files within 2 days before or after the input date +; 13-Apr-2026 - Massa P. (FHNW), updated to search for calibration files the day before or after the input date (in case there are no files for the input date) ;- function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobber=clobber @@ -41,124 +41,69 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe default, clobber, 0 day_in_s = 86400.d ;; Number of seconds in a day. To be used later to identify the most appropriate calibration file - day_limit = 2 ;; Maximum number of days before or after the input date to consider when searching for a calibration file (arbitrary choice) site = 'http://dataarchive.stix.i4ds.net' date_path = get_fid(start_time,end_time,/full,delim='/') + ;; Concatenate files from the input date and from the day before/after type_path = '/fits/CAL/' - path = type_path + date_path[0] +'/CAL' filter = '*stix-cal-energy*.fits' + + ;; Load the ELUT file for the iput time (it is used later for checks) + elut_input_time = stx_date2elut_file(start_time) + + + found_files_array=[] + + ;; Input day + path = type_path + date_path[0] +'/CAL' found_files=sock_find(site,filter,path=path,count=count) + + if count gt 0 then found_files_array=[found_files_array, found_files] + + ;; Day before + this_start_time_before = anytim(anytim(start_time) - day_in_s, /vms) + this_end_time_before = anytim(anytim(end_time) - day_in_s, /vms) - if count eq 0 then begin - message, [" ", " ", "No STIX calibration files found for this day ("+date_path[0]+"). Searching for calibration file which is closest in time.", " ", " "], /continue - - ;; Forward in time (for 2 days max - arbitrary choice) - day_n_forward=1 - count_forward=0 - while (day_n_forward le day_limit) and (count_forward eq 0) do begin - - this_start_time_forward = anytim(anytim(start_time) + day_n_forward * day_in_s, /vms) - this_end_time_forward = anytim(anytim(end_time) + day_n_forward * day_in_s, /vms) - - date_path_forward = get_fid(this_start_time_forward,this_end_time_forward,/full,delim='/') - - path = type_path + date_path_forward[0] +'/CAL' - found_files_forward = sock_find(site,filter,path=path,count=count_forward) - - day_n_forward += 1 - - endwhile - - day_n_forward -= 1 ;; Remove 1 day to compensate that when the loop exits after finding a file, each counter is +1 from the actual day offset that produced the match. - - ;; Backward in time (for 2 days max - arbitrary choice) - day_n_backward=1 - count_backward=0 - while (day_n_backward le day_limit) and (count_backward eq 0) do begin - - this_start_time_backward = anytim(anytim(start_time) - day_n_backward * day_in_s, /vms) - this_end_time_backward = anytim(anytim(end_time) - day_n_backward * day_in_s, /vms) - - date_path_backward = get_fid(this_start_time_backward,this_end_time_backward,/full,delim='/') - - path = type_path + date_path_backward[0] +'/CAL' - found_files_backward = sock_find(site,filter,path=path,count=count_backward) + path_before = get_fid(this_start_time_before,this_end_time_before,/full,delim='/') + + path = type_path + path_before[0] +'/CAL' + found_files=sock_find(site,filter,path=path,count=count) - day_n_backward += 1 + if count gt 0 then found_files_array=[found_files_array, found_files] + + ;; Day after + this_start_time_after = anytim(anytim(start_time) + day_in_s, /vms) + this_end_time_after = anytim(anytim(end_time) + day_in_s, /vms) - endwhile - - day_n_backward -= 1 ;; Remove 1 day to compensate that when the loop exits after finding a file, each counter is +1 from the actual day offset that produced the match. - - if (count_forward eq 0) and (count_backward eq 0) then begin - - message, [" ", " ", "No STIX calibration file was found within 2 days before or after " +date_path[0]+". Please, download a calibration file manually.", " ", " "] - - endif else begin - - diff_days = day_n_forward - day_n_backward - ;; if the difference is negative, it means that the closest file is the one found forward in time - ;; if the difference is positive, it means that the closest file is the one found backward in time - ;; if the difference is zero, we arbitrarily select the file found forward in time - - if diff_days le 0 then begin - - count = count_forward - found_files = found_files_forward - date_path = date_path_forward - - endif else begin - - count = count_backward - found_files = found_files_backward - date_path = date_path_backward - - endelse - - endelse - - endif + path_after = get_fid(this_start_time_after,this_end_time_after,/full,delim='/') - message, [" ", " ", "Download STIX calibration file recorded at " +date_path[0], " ", " "], /continue + path = type_path + path_after[0] +'/CAL' + found_files=sock_find(site,filter,path=path,count=count) + if count gt 0 then found_files_array=[found_files_array, found_files] + + + ;;********* + len_path = STRLEN(site+path) - - if count eq 1 then begin - - selected_file = found_files[0] + + if n_elements(found_files_array) eq 0 then begin - ;; Extract the calibration file start time from the filename. This is used below to verify that the ELUT onboard at the time of calibration matches the one used during the input time range. + message, $ + [" ", " ", "No STIX calibration file was found within 1 day before or after " +date_path[0]+". Please, download an appropriate calibration file manually.", " ", " "] - len_full_path = STRLEN(found_files[0]) - filename = STRMID(found_files[0], len_path+1, len_full_path) - - string_dates = STRMID(filename, STRLEN('solo_CAL_stix-cal-energy_'), STRLEN(filename)-STRLEN('solo_CAL_stix-cal-energy_')-9) ;; 9 is the number of characters in the string '_V02.fits' - - mid_part = STRPOS(string_dates, '-') - - this_start_time_file = STRMID(string_dates, 0, mid_part) - year = STRMID(this_start_time_file, 0, 4) - month = STRMID(this_start_time_file, 4, 2) - day = STRMID(this_start_time_file, 6, 2) - hour = STRMID(this_start_time_file, 9, 2) - min = STRMID(this_start_time_file, 11, 2) - sec = STRMID(this_start_time_file, 13, 2) - start_time_file = year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec - endif else begin - - message, [" ", " ", "More than 1 STIX calibration file found for this day ("+date_path[0]+"). Download the one with largest duration.", " ", " "], /continue start_time_file = [] end_time_file = [] + elut_check = [] ;; Extract file names - for i = 0,n_elements(found_files)-1 do begin - - len_full_path = STRLEN(found_files[i]) - filename = STRMID(found_files[i], len_path+1, len_full_path) + for i = 0,n_elements(found_files_array)-1 do begin + + len_full_path = STRLEN(found_files_array[i]) + filename = STRMID(found_files_array[i], len_path+1, len_full_path) string_dates = STRMID(filename, STRLEN('solo_CAL_stix-cal-energy_'), STRLEN(filename)-STRLEN('solo_CAL_stix-cal-energy_')-9) ;; 9 is the number of characters in the string '_V02.fits' @@ -171,7 +116,9 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe hour = STRMID(this_start_time_file, 9, 2) min = STRMID(this_start_time_file, 11, 2) sec = STRMID(this_start_time_file, 13, 2) - start_time_file = [start_time_file, anytim(year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec)] + + this_start_time = anytim(year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec) + start_time_file = [start_time_file, this_start_time] this_end_time_file = STRMID(string_dates, mid_part+1, STRLEN(string_dates)) year = STRMID(this_end_time_file, 0, 4) @@ -182,35 +129,78 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe sec = STRMID(this_end_time_file, 13, 2) end_time_file = [end_time_file, anytim(year + '-' + month + '-' + day + 'T' + hour + ':' + min + ':' + sec)] + elut_check = [elut_check, (elut_input_time eq stx_date2elut_file(this_start_time))] + endfor - ;; Select file with largest duration - duration = end_time_file - start_time_file + ;; Compute mid point of the time interval of the calibration file + mid_time_file = (start_time_file + end_time_file) / 2. - idx_file = where(duration eq max(duration)) - selected_file = found_files[idx_file] + ;; Compute duration file + duration = end_time_file - start_time_file + + ;; ELUT check + idx_elut = where(elut_check, n_elut_check) - start_time_file = start_time_file[idx_file] + if n_elut_check eq 0 then begin + + message, [" ", " ", "Available STIX calibration files within 1 day before or after " +date_path[0]+" have been registered with an onboard ELUT different from that used during the selected time range and cannot be used. Please, download an appropriate calibration file manually.", " ", " "] + + endif else begin + + ;; Select calibration files recorded with the same ELUT as the one that was onboard during the input time range + found_files = found_files_array[idx_elut] + duration = duration[idx_elut] + start_time_file = start_time_file[idx_elut] + mid_time_file = mid_time_file[idx_elut] + + ;; Select calibration file with largest duration + idx_duration = where(duration ge day_in_s / 2., n_duration) + + if n_duration eq 0 then begin + + message, [" ", " ", "Available STIX calibration files within 1 day before or after " +date_path[0]+" have a duration shorter than half a day and should not be used. Please, download an appropriate calibration file manually.", " ", " "] + + endif else begin + + found_files = found_files[idx_duration] + duration = duration[idx_duration] + start_time_file = start_time_file[idx_duration] + mid_time_file = mid_time_file[idx_duration] + + ;; If multiple files have the largest duration, select the closest one to the input time range + if n_duration gt 1 then begin - endelse - - - ;; Check that the ELUT onboard when the calibration file was recorded matches the ELUT onboard during the input time range; otherwise, raise an error. + ;; Select the file closest in time to input time range + mid_time_input = (anytim(start_time) + anytim(end_time)) / 2. + time_diff = abs(mid_time_file - mid_time_input) + + idx_time_diff = sort(time_diff) + + idx_file = idx_time_diff[0] - ;; Compare ELUT - elut_file = stx_date2elut_file(start_time_file) - elut_time_range = stx_date2elut_file(start_time) - if elut_file ne elut_time_range then message, $ - [" ", " ", "The closest-in-time calibration file was acquired with an onboard ELUT that differs from the one used during the selected time range. Please, manually download a different calibration file.", " ", " "] + selected_file = found_files[idx_file] + selected_file_time = start_time_file[idx_file] + + endif else begin + + selected_file = found_files + selected_file_time = start_time_file + + endelse + + endelse + + endelse + + endelse + ;; Download file + message, [" ", " ", "Download STIX calibration file recorded at " +anytim(selected_file_time, /vms), " ", " "], /continue sock_copy, selected_file, out_name, local_file=out_file, out_dir = out_dir, clobber=clobber - - filename = STRMID(out_file, STRLEN(out_dir), STRLEN(out_file)-STRLEN(out_dir)) - - message, [" ", " ", "Downloaded file: " + filename, " ", " "], /continue - + return, out_file - + end From ee478d01fcc492a495f58b43bef393b81792021c Mon Sep 17 00:00:00 2001 From: paolomassa Date: Mon, 4 May 2026 10:39:34 +0200 Subject: [PATCH 07/13] Change to consider detectors and pixels that are actually present in the fits file --- stix/idl/processing/spectrogram/stx_convert_pixel_data.pro | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro index 63e30903..6bec9b82 100644 --- a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro +++ b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro @@ -451,7 +451,7 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit elut_data = stx_elut_correction(this_count_rates, this_count_rates_error, $ energy_bin_idx, energy_bin_low, energy_bin_high, energy_high, energy_low, e_bin, this_energy_range, $ - spectrum, pixel_ind, det_ind, /silent) + spectrum, pixels_used, detectors_used, /silent) count_rates_elut[e_bin,*,*,idx_time_min[t_bin]:idx_time_max[t_bin]] = elut_data.COUNTS count_rates_error_elut[e_bin,*,*,idx_time_min[t_bin]:idx_time_max[t_bin]] = elut_data.COUNTS_ERROR From 9a964cec10b0ebc094ee98fcab1f2004c9c06040 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Mon, 4 May 2026 11:10:54 +0200 Subject: [PATCH 08/13] Keyword elut_correction removed from stx_convert_pixel_data --- stix/idl/processing/spectrogram/stx_convert_pixel_data.pro | 1 - 1 file changed, 1 deletion(-) diff --git a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro index 6bec9b82..dceec249 100644 --- a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro +++ b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro @@ -108,7 +108,6 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit default, shift_duration, 0 default, plot, 1 default, det_ind, 'top24' - default, elut_correction, 1 default, silent, 0 default, xspec, 0 default, sys_uncert, 0.05 From c8560e90c6dc1351836fd6003cddb247385fa4a5 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Mon, 4 May 2026 18:56:27 +0200 Subject: [PATCH 09/13] Moved lines to provide elut_file as input to stx_fits_info_params. --- .../idl/processing/spectrogram/stx_convert_pixel_data.pro | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro index dceec249..4877f969 100644 --- a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro +++ b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro @@ -197,6 +197,10 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit if n_elements(distance) ne 0 then fits_distance = distance + ;; Check if science and background files are reconrded with the same ELUT + elut_filename = stx_date2elut_file(stx_time2any(t_axis.TIME_START[0])) + stx_read_elut, elut_gain, elut_offset, adc4096_str, elut_filename = elut_filename + fits_info_params = stx_fits_info_params( fits_path_data = fits_path_data, data_level = data_level, $ distance = fits_distance, time_shift = time_shift, fits_path_bk = fits_path_bk, uid = uid, $ generate_fits = generate_fits, specfile = specfile, srmfile = srmfile, elut_file = elut_filename, silent = silent) @@ -255,10 +259,6 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit if keyword_set(fits_path_bk) then begin - ;; Check if science and background files are reconrded with the same ELUT - elut_filename = stx_date2elut_file(stx_time2any(t_axis.TIME_START[0])) - stx_read_elut, elut_gain, elut_offset, adc4096_str, elut_filename = elut_filename - elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) stx_read_elut, ekev_actual = ekev_actual_bkg, elut_filename = elut_filename_bkg From 7b0d7a1cdd70203d39fb2e62a39748cf1b49b7b1 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Tue, 5 May 2026 15:24:47 +0200 Subject: [PATCH 10/13] Fixed bug which made 'stx_science_data_lightcurve' break --- stix/idl/processing/spectrogram/stx_convert_pixel_data.pro | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro index 4877f969..fb5c6323 100644 --- a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro +++ b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro @@ -145,6 +145,12 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit control_header= control_header, energy_str = energy_str, energy_header = energy_header, t_axis = t_axis, energy_shift = energy_shift, e_axis = e_axis , use_discriminators = 0, $ shift_duration = shift_duration, silent=silent + ;; Define default file name for spectrum and srm + uid = control_str.request_id + default, specfile, 'stx_spectrum_' + strtrim(uid,2) + '.fits' + if ~keyword_set(srmfile) then $ + srmfile = xspec ? 'stx_srm_' + strtrim(uid,2) + '_XSPEC.fits' : 'stx_srm_' + strtrim(uid,2) + '.fits' + ;; Select indices of the energy bins (among the 32) that are actually present in the pixel data science file energy_bin_mask = data_str.energy_bin_mask energy_bin_idx = where(energy_bin_mask eq 1) From 1a24f459ae35f5ea02bdadc3c26ed6c3f985f792 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Fri, 8 May 2026 10:42:50 +0200 Subject: [PATCH 11/13] Use new transmission calibration for spectroscopy The procedure 'stx_convert_spectrogram' now utilizes the new transmission calibration for the SRM calculation. The code has been refactored so that 'stx_convert_pixel_data' and 'stx_convert_spectrogram' utilize the same function ('stx_convert_spectrogram2ospex') to avoid code replication. Finally, 'stx_convert_spectrogram2ospex' prints a warning if the fine-resolution sub-collimators are considered, as the high-energy calibration is not completed yet. --- stix/idl/io/stx_fits_info_params.pro | 7 +- .../livetime/stx_spectrogram_livetime.pro | 20 +- .../pixel_data/stx_construct_pixel_data.pro | 3 - .../spectrogram/stx_convert_pixel_data.pro | 227 ++----------- .../spectrogram/stx_convert_spectrogram.pro | 267 ++++++++------- .../stx_convert_spectrogram2ospex.pro | 305 ++++++++++++++++++ .../stx_read_spectrogram_fits_file.pro | 24 +- 7 files changed, 511 insertions(+), 342 deletions(-) create mode 100644 stix/idl/processing/spectrogram/stx_convert_spectrogram2ospex.pro diff --git a/stix/idl/io/stx_fits_info_params.pro b/stix/idl/io/stx_fits_info_params.pro index 303c242d..8e959899 100644 --- a/stix/idl/io/stx_fits_info_params.pro +++ b/stix/idl/io/stx_fits_info_params.pro @@ -61,7 +61,8 @@ ; ; :history: ; 17-Aug-2022 - ECMD (Graz), initial release -; 05-Apr-2023 - ECMD (Graz), fix issue with taking header distance as default +; 05-Apr-2023 - ECMD (Graz), fix issue with taking header distance as default +; 07-May-2026 - Massa P. (FHNW), print warning if generate_fits = 0 (regardless of 'specfile' or 'srmfile' are defined) ; ;- function stx_fits_info_params, fits_path_data = fits_path_data, data_level = data_level, $ @@ -69,8 +70,8 @@ function stx_fits_info_params, fits_path_data = fits_path_data, data_level = dat generate_fits = generate_fits, specfile = specfile, srmfile = srmfile, elut_file = elut_file, silent = silent if n_elements(generate_fits) ne 0 then begin - if generate_fits eq 0 and keyword_set(specfile) || keyword_set(srmfile) then begin - message, 'FITS file generation has been set to 0 but an output filename has been specified.' + if generate_fits eq 0 then begin + message, [" ", " ", 'FITS file generation has been set to 0.', " ", " "], /continue endif endif diff --git a/stix/idl/processing/livetime/stx_spectrogram_livetime.pro b/stix/idl/processing/livetime/stx_spectrogram_livetime.pro index d5425e29..9b7a16b5 100644 --- a/stix/idl/processing/livetime/stx_spectrogram_livetime.pro +++ b/stix/idl/processing/livetime/stx_spectrogram_livetime.pro @@ -38,7 +38,8 @@ ; 24-Feb-2021 - ECMD (Graz), initial release ; 22-Feb-2022 - ECMD (Graz), documented, added livetime error estimation ; 23-Mar-2026 - Massa P. (FHNW), made it compatible with live time uncertainty propagation -; +; 07-May-2026 - Massa P. (FHNW), fixed bug. In the previous version, counts were normalized by livetime fraction instead of livetime. +; Now, both the livetime and the livetime fraction (and corresponding uncertainties) are returned ;- function stx_spectrogram_livetime, spectrogram, corrected_counts =corrected_counts, corrected_error= corrected_error, level = level ;convert the triggers to livetime @@ -49,7 +50,8 @@ default, level, 1 nenergies = (spectrogram.counts.dim)[0] det_used = where(spectrogram.detector_mask eq 1, ndet) + 1 ; stx_livetime_fraction expects detector number in the range 1 - 32 - + time_bin_duration = transpose(cmreplicate(spectrogram.time_axis.duration, nenergies)) + case level of 1: begin dim_counts = [nenergies, ndet, ntimes] @@ -57,7 +59,6 @@ default, level, 1 livetime_fraction_data = stx_livetime_fraction(triggergram, det_used) livetime_fraction = transpose( rebin(reform(livetime_fraction_data.livetime_fraction, ndet, ntimes),[dim_counts[1:2],dim_counts[0]]),[2,0,1]) livetime_fraction_err = transpose( rebin(reform(livetime_fraction_data.livetime_fraction_err, ndet, ntimes),[dim_counts[1:2],dim_counts[0]]),[2,0,1]) - end 4:begin @@ -69,17 +70,22 @@ default, level, 1 livetime_fraction = transpose( rebin(reform(livetime_fraction_data.livetime_fraction[0,*]),[dim_counts[1],dim_counts[0]])) livetime_fraction_err = transpose( rebin(reform(livetime_fraction_data.livetime_fraction_err[0,*]),[dim_counts[1],dim_counts[0]])) - end else: message, 'Currently supported compaction levels are 1 (pixel data) and 4 (spectrogram)' endcase + + livetime = time_bin_duration * livetime_fraction + livetime_err = time_bin_duration * livetime_fraction_err - corrected_counts = f_div(spectrogram.counts,livetime_fraction) + corrected_counts = f_div(spectrogram.counts,livetime) corrected_error = abs(corrected_counts) * sqrt( f_div(spectrogram.error,spectrogram.counts)^2. + $ - f_div(livetime_fraction_err,livetime_fraction)^2. ) + f_div(livetime_err,livetime)^2. ) - return, livetime_fraction + return, {livetime: livetime, $ + livetime_err: livetime_err, $ + livetime_fraction:livetime_fraction, $ + livetime_fraction_err: livetime_fraction_err} end \ No newline at end of file diff --git a/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro b/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro index 53b5be54..c40fcc36 100644 --- a/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro +++ b/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro @@ -249,10 +249,7 @@ if keyword_set(path_bkg_file) then begin ;; Check if science and background files are reconrded with the same ELUT elut_filename = stx_date2elut_file(stx_time2any(this_time_range[0])) - stx_read_elut, elut_gain, elut_offset, adc4096_str, elut_filename = elut_filename, ekev_actual = ekev_actual_elut - elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) - stx_read_elut, ekev_actual = ekev_actual_bkg, elut_filename = elut_filename_bkg ;; Compare ELUT tables elut_comp = STRCMP(elut_filename, elut_filename_bkg) diff --git a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro index fb5c6323..162bdf29 100644 --- a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro +++ b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro @@ -74,7 +74,7 @@ ; time resolution is at least equal to delta_time_min (defined in seconds). ; Rebinned data are used to construct spectra needed for the ELUT correction. ; -; calib_data: if a 'stx_calibration_data' structure is passed as input, then the ELUT correction is applied +; calib_data: if a 'stx_calibration_data' structure is passed as input, then the ELUT correction is applied ; ; :examples: ; fits_path_data = 'solo_L1A_stix-sci-xray-l1-2104170007_20210417T153019-20210417T171825_009610_V01.fits' @@ -95,10 +95,11 @@ ; 16-Jun-2023 - ECMD (Graz), for a source location dependent response estimate, the location in HPC and the auxiliary ephemeris file must be provided. ; 06-Dec-2023 - ECMD (Graz), added silent keyword, more information is now printed if not set ; 2024-07-12, F. Schuller (AIP): added optional keyword xspec -; 07-May-2025 - Stiefel M. and Massa P., re-implemented to take into account new ELUT correction, new subcollimator transmission and live time normalization +; 07-May-2025 - Stiefel M. and Massa P. (FHNW), re-implemented to take into account new ELUT correction, new subcollimator transmission and live time normalization +; 07-May-2026 - Massa P. (FHNW), this function now calls 'stx_convert_spectrogram2ospex' ; ;- -pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fits_path_bk, $ +pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fits_path_bk, $ time_shift = time_shift, energy_shift = energy_shift, distance = distance, $ flare_location_stx = flare_location_stx, det_ind = det_ind, pix_ind = pix_ind, calib_data = calib_data, shift_duration = shift_duration, $ no_attenuation = no_attenuation, sys_uncert = sys_uncert, generate_fits = generate_fits, specfile = specfile, $ @@ -114,7 +115,6 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit default, delta_time_min, 20. default, tailing, 1 default, include_damage, 1 - default, flare_location_stx, [0.,0.] if n_elements(time_shift) eq 0 then begin if ~keyword_set(silent) then begin @@ -203,9 +203,7 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit if n_elements(distance) ne 0 then fits_distance = distance - ;; Check if science and background files are reconrded with the same ELUT elut_filename = stx_date2elut_file(stx_time2any(t_axis.TIME_START[0])) - stx_read_elut, elut_gain, elut_offset, adc4096_str, elut_filename = elut_filename fits_info_params = stx_fits_info_params( fits_path_data = fits_path_data, data_level = data_level, $ distance = fits_distance, time_shift = time_shift, fits_path_bk = fits_path_bk, uid = uid, $ @@ -266,7 +264,6 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit if keyword_set(fits_path_bk) then begin elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) - stx_read_elut, ekev_actual = ekev_actual_bkg, elut_filename = elut_filename_bkg ;; Compare ELUT tables elut_comp = STRCMP(elut_filename, elut_filename_bkg) @@ -372,7 +369,7 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit if keyword_set(calib_data) then begin - + ;; Create daily ELUT energy_bin_low = calib_data.ENERGY_BIN_LOW energy_bin_high = calib_data.ENERGY_BIN_HIGH @@ -404,8 +401,8 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit idx_time_min = iall[0:-2] idx_time_max = iall[1:-1]-1 - - + + count_rates_elut = fltarr(count_rates.dim) count_rates_error_elut = fltarr(count_rates_error.dim) @@ -462,7 +459,7 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit count_rates_error_elut[e_bin,*,*,idx_time_min[t_bin]:idx_time_max[t_bin]] = elut_data.COUNTS_ERROR endfor - + endfor count_rates = count_rates_elut @@ -506,201 +503,21 @@ pro stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit spec *= avg_live_time_bins_rep spec_error *= avg_live_time_bins_rep - ;;***************** CHECK FOR RCR CHANGES - - rcr = data_str.rcr - - ;get the rcr states and the times of rcr changes from the ql_lightcurves structure - ut_rcr = stx_time2any(t_axis.time_end) - - find_changes, rcr, index, state, count=count - ; ************************************************************ - ; ******************** TEMPORARY FIX ************************* - ; ***** Andrea: 2022-April-05 - ; Temporarily creation of the no_attenuation keyword in order - ; to avoid attenuation of the fitted curve. This is useful for - ; obtaining thermal fit parameters with the BKG detector in the - ; case the attenuator is inserted. We tested it with the X - ; class flare on 2021-Oct-26 and it works nicely. - if keyword_set(no_attenuation) then begin - rcr = rcr*0. - index = 0 - state = 0 - endif - ; ************************************************************ - ; ************************************************************ - - ; ******************** TEMPORARY FIX ************************* - ; ***** ECMD: 2022-Jun-27 - ; As the reported time of the RCR status change can be inaccurate - ; up to several seconds correct this by finding the times where there is a - ; large change in counts in the counts of the 5 - 6 keV energy bin. - ; find all time intervals where the difference between adjacent bins is large - if max(rcr) gt 0 then begin; skip if in the standard state of RCR0 for the full time range - jumps = where(abs(total(spec,1) - shift(total(spec,1),-1)) * 24. / n_elements(detectors_used) gt 1e4) - ;jumps = where(abs(total(spec,1) - shift(total(spec,1),-1)) gt 1e4) - ; include the starting bin - jumps = [0, jumps] - ; as the attenuator motion can be present in two consecutive bins select only the first - idx_jumps = where(abs(jumps - shift(jumps, -1)) gt 2) - jumps_use= [jumps[idx_jumps]] - ; each transition should correspond close in time to a recorded transition in the FITS file - ; adjust the time indexes of these transitions to the closest jumps - closest_jumps = value_closest(jumps_use, index) - index = jumps_use[closest_jumps] - - endif - - ;add the rcr information to a specpar structure so it can be included in the spectrum FITS file - specpar = { sp_atten_state : {time:ut_rcr[index], state:state}, flare_xyoffset : fltarr(2), use_flare_xyoffset:0 } - - ;;***************** CREATE SRM - - pixel_mask =detector_mask_used ## pixel_mask_used - - transmission = read_csv(loc_file( 'stix_transmission_highres_20251110.csv', path = getenv('STX_GRID'))) - - emin = 4 - emax = 150 - phe = transmission.(0) - phe = phe[where(phe gt emin-1 and phe lt 3.5*emax)] - edge_products, phe, mean = mean_phe, width = w_phe - ph_edges = [mean_phe[0] - w_phe[0], mean_phe] - - distance = fits_info_params.distance - dist_factor = 1./(distance^2.) - ;make the srm for the appropriate pixel mask and energy edges - ;srm = stx_build_pixel_drm(ct_edges, pixel_mask, ph_energy_edges = ph_edges, dist_factor = dist_factor, tailing = tailing, include_damage = include_damage, _extra=extra) - ph_edges = get_uniq( [ph_edges,ct_edges],epsilon=0.0001) - - ;; Compute subc. transmission - edge_products,ph_edges, mean=ph_in - subc_transmission = stx_subc_transmission(flare_location_stx, ph_in, _extra=extra) - - if n_elements(detectors_used) eq 1 then begin - grid_factor = reform(subc_transmission[*,detectors_used]) - endif else begin - grid_factor = total(subc_transmission[*,detectors_used],2) / n_elements(detectors_used) - endelse - - - ;; Check if BKG or CFL detectors are used - - idx_cfl = where(detectors_used eq 8, n_cfl) - idx_bkg = where(detectors_used eq 9, n_bkg) - - if n_cfl eq 1 then message, "CFL detector can not be selected for spectral fitting." - if (n_elements(detectors_used) gt 1) and (n_bkg eq 1) then message, "BKG detector can not be selected together with imaging detectors for spectral fitting." - - if n_elements(detectors_used) eq 1 then if detectors_used eq 9 then begin - grid_transmission_file = concat_dir(getenv('STX_GRID'), 'real_bkg_grid_transmission.txt') - readcol, grid_transmission_file, bk_grid_factors, format = 'f', skip = 2, silent = silent - - grid_factor = average(bk_grid_factors[pixels_used]) - - endif + ;; Create spectrogram structure + spec_data = {type : 'stx_spectrogram', $ + data : spec, $ + t_axis : t_axis, $ + e_axis : e_axis, $ + ltime : avg_live_time_fraction_bins, $ + attenuator_state : data_str.RCR , $ + error : spec_error} - ;; Creates appropriate SRM for different attenuator states - rcr_states = specpar.sp_atten_state.state - rcr_states = rcr_states[uniq(rcr_states, sort(rcr_states))] - nrcr_states = n_elements(rcr_states) - - srm_atten = replicate( {rcr:0, srm:fltarr(n_elements(ct_edges)-1,n_elements(ph_edges)-1)},nrcr_states ) - - for i =0, nrcr_states-1 do begin - ;make the srm for the appropriate pixel mask and energy edges - rcr = rcr_states[i] - - srm = stx_build_pixel_drm(ct_edges, pixel_mask, rcr = rcr, grid_factor= grid_factor, ph_energy_edges = ph_edges, dist_factor = dist_factor, tailing = tailing, include_damage = include_damage, _extra=extra) - srm_atten[i].srm = srm.smatrix - srm_atten[i].rcr = rcr - - endfor - - ;;***************** SAVE FITS - detector_label = stx_det_mask2label(detector_mask_used) - pixel_label = stx_pix_mask2label(pixel_mask_used) - - ospex_obj = ospex(/no) - - ;if the fits keyword is set write the spectrogram and srm data to fits files and then read them in to the ospex object - if fits_info_params.generate_fits eq 1 then begin - utime = transpose([stx_time2any( t_axis.time_start )]) - - ;spectrogram structure for passing to fits writer routine - spectrum_in = { type : 'stx_spectrogram', $ - data : spec, $ - t_axis : t_axis, $ - e_axis : e_axis, $ - ltime : avg_live_time_fraction_bins, $ - attenuator_state : data_str.rcr , $ - error : spec_error } - - specfilename = fits_info_params.specfile - srmfilename = fits_info_params.srmfile - - fits_info_params.grid_factor.add, grid_factor - fits_info_params.detused = detector_label + ', Pixels: ' + pixel_label - - if keyword_set(xspec) then begin - ;xspec in general works with energy depandent systematic errors - e_axis = spectrum_in.e_axis - n_energies = n_elements(e_axis.mean) - sys_err = fltarr(n_energies) - - idx_below10kev = where(e_axis.mean lt 10, cb10) - sys_err[*] = 0.03 - if cb10 gt 0 then sys_err[idx_below10kev] = 0.05 - idx_below7kev = where(e_axis.mean lt 7, cb7) - if cb7 gt 0 then sys_err[idx_below7kev] = 0.07 - - sys_err = rebin(sys_err, n_energies,n_times) - endif - - stx_write_ospex_fits, spectrum = spectrum_in, srmdata = srm, specpar = specpar, time_shift = time_shift, $ - srm_atten = srm_atten, specfilename = specfilename, srmfilename = srmfilename, ph_edges = ph_edges, $ - fits_info_params = fits_info_params, xspec = xspec, silent = silent - - ospex_obj->set, spex_file_reader = 'stx_read_sp' - ospex_obj->set, spex_specfile = specfilename ; name of your spectrum file - ospex_obj->set, spex_drmfile = srmfilename - - endif else begin - ;if the generate_fits keyword is not set use the spex_user_data strategy to pass in the data directly to the ospex object - - energy_edges = e_axis.edges_2 - Edge_Products, ph_edges, edges_2 = ph_edges2 - - utime2 = transpose(stx_time2any( [[t_axis.time_start], [t_axis.time_end]] )) - - ospex_obj->set, spex_data_source = 'spex_user_data' - ospex_obj->set, spectrum = float(spec), $ - spex_ct_edges = energy_edges, $ - spex_ut_edges = utime2, $ - livetime = avg_live_time_bins_rep, $ - errors = spec_error - srm = rep_tag_name(srm,'smatrix','drm') - ospex_obj->set, spex_respinfo = srm - ospex_obj->set, spex_area = srm.area - ospex_obj->set, spex_detectors = 'STIX' - ospex_obj->set, spex_drm_ct_edges = energy_edges - ospex_obj->set, spex_drm_ph_edges = ph_edges2 - endelse - - ospex_obj->set, spex_uncert = sys_uncert - ospex_obj->set, spex_error_use_expected = 0 - - counts_str = ospex_obj->getdata(spex_units='counts') - origunits = ospex_obj->get(/spex_data_origunits) - origunits.data_name = 'STIX' - ospex_obj->set, spex_data_origunits = origunits - - if keyword_set(plot) then begin - ospex_obj ->gui - ospex_obj ->set, spex_eband = get_edges([4.,10.,15.,25, 50, 84.], /edges_2) - ospex_obj ->plot_time, spex_units='flux', /show_err, obj = plotman_object - endif + ;;------------------------------------------------------ + + stx_convert_spectrogram2ospex, spec_data, pixel_mask_used, detector_mask_used, fits_info_params, ct_edges, $ + no_attenuation=no_attenuation, flare_location_stx=flare_location_stx, time_shift = time_shift, $ + sys_uncert=sys_uncert, silent=silent, plot=plot, xspec=xspec, ospex_obj = ospex_obj, $ + tailing=tailing, include_damage=include_damage, _extra=extra end - diff --git a/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro b/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro index faf72921..dd840e3b 100644 --- a/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro +++ b/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro @@ -25,6 +25,9 @@ ; ; distance : in, optional, type="float", default= "1." ; The distance between Solar Orbiter and the Sun centre in Astronomical Units needed to correct flux. +; +; flare_location_stx : in, type="2 element float array" +; the location of the flare (X,Y) in the STIX imaging frame [arcsec] ; ; time_shift : in, optional, type="float", default="0." ; The difference in seconds in light travel time between the Sun and Earth and the Sun and Solar Orbiter @@ -34,12 +37,6 @@ ; Shift all energies by this value in keV. Rarely needed only for cases where there is a significant shift ; in calibration before a new ELUT can be uploaded. ; -; flare_location_hpc : in, type="2 element float array" -; the location of the flare (X,Y) in Helioprojective Cartesian coordinates as seen from Solar Orbiter [arcsec] -; -; aux_fits_file : in, required if flare_location_hpc is passed in, type="string" -; the path of the auxiliary ephemeris FITS file to be read." -; ; shift_duration : in, type="boolean", default="0" ; Shift all time bins by 1 to account for FSW time input discrepancy prior to 09-Dec-2021. ; N.B. WILL ONLY WORK WITH FULL TIME RESOLUTION DATA WHICH IS OFTEN NOT THE CASE FOR PIXEL DATA. @@ -58,16 +55,13 @@ ; srmfile : in, type="string", default="'stx_srm_'+ UID + '.fits'" ; File name to use when saving the srm FITS file for OSPEX input. ; -; background_data : out, type="stx_background_data structure" -; Structure containing the subtracted background for external plotting. -; ; plot : in, type="boolean", default="1" ; If set open OSPEX GUI and plot lightcurve in standard quicklook energy bands ; where there is data present ; ; xspec : in, type="boolean", default="0" ; If set, generate SRM file compatible with XSPEC rather than OSPEX. -; +; ; ospex_obj : out, type="OSPEX object" ; ; @@ -89,50 +83,108 @@ ; 16-Jun-2023 - ECMD (Graz), for a source location dependent response estimate, the location in HPC and the auxiliary ephemeris file must be provided. ; 06-Dec-2023 - ECMD (Graz), added silent keyword, more information is now printed if not set ; 2024-07-12, F. Schuller (AIP): added optional keyword xspec -; +; 07-May-2026 - Massa P., removed ELUT correction as it is not possible to apply it to spectrogram data. Furter, this routine now calls 'stx_convert_spectrogram2ospex' +; which is called also by 'stx_convert_pixel_data' ;- -pro stx_convert_spectrogram, fits_path_data = fits_path_data, fits_path_bk = fits_path_bk,$ - aux_fits_file = aux_fits_file, flare_location_hpc = flare_location_hpc, flare_location_stx = flare_location_stx, $ - time_shift = time_shift, energy_shift = energy_shift, distance = distance, replace_doubles = replace_doubles, $ - keep_short_bins = keep_short_bins, apply_time_shift = apply_time_shift, elut_correction = elut_correction, $ +pro stx_convert_spectrogram, fits_path_data = fits_path_data, fits_path_bk = fits_path_bk,$ + flare_location_stx = flare_location_stx, time_shift = time_shift, energy_shift = energy_shift, $ + distance = distance, replace_doubles = replace_doubles, keep_short_bins = keep_short_bins, $ shift_duration = shift_duration, no_attenuation = no_attenuation, sys_uncert = sys_uncert, $ generate_fits = generate_fits, specfile = specfile, srmfile = srmfile, silent = silent, $ - background_data = background_data, plot = plot, xspec=xspec, ospex_obj = ospex_obj + plot = plot, xspec=xspec, ospex_obj = ospex_obj, tailing=tailing, include_damage=include_damage, _extra=extra - default, plot, 1 - default, silent, 0 - default, xspec, 0 if n_elements(time_shift) eq 0 then begin if ~keyword_set(silent) then begin - message, 'Time shift value is not set. Using default value of 0 [s].', /info - print, 'File averaged values can be obtained from the FITS file header' - print, 'using stx_get_header_corrections.pro.' + message, 'Time shift value is not set. Using default value of 0 [s].', /info + print, 'File averaged values can be obtained from the FITS file header' + print, 'using stx_get_header_corrections.pro.' endif time_shift = 0. endif - + + default, shift_duration, 0 + default, sys_uncert, 0.05 + default, silent, 0 default, plot, 1 - default, elut_correction, 1 + default, xspec, 0 + default, tailing, 1 + default, include_damage, 1 + ;;------------------------------------------------------------ + stx_read_spectrogram_fits_file, fits_path_data, time_shift, primary_header = primary_header, data_str = data_str, data_header = data_header, control_str = control_str, $ control_header= control_header, energy_str = energy_str, energy_header = energy_header, t_axis = t_axis, energy_shift = energy_shift, e_axis = e_axis , use_discriminators = 0,$ replace_doubles = replace_doubles, keep_short_bins = keep_short_bins, shift_duration = shift_duration - data_level = 4 - start_time = atime(stx_time2any((t_axis.time_start)[0])) + ;; Select indices of the energy bins (among the 32) that are actually present in the pixel data science file + energy_bin_mask = data_str.energy_bin_mask + energy_bin_idx = where(energy_bin_mask eq 1) + + energy_low = e_axis.LOW + energy_high = e_axis.HIGH + + energy_min = min(energy_low) + energy_max = max(energy_high) + + if keyword_set(fits_path_bk) then begin + + stx_read_pixel_data_fits_file, fits_path_bk, data_str = data_bkg, t_axis = t_axis_bkg, e_axis = e_axis_bkg, _extra=extra + + if n_elements(t_axis_bkg.DURATION) gt 1 then message, 'The chosen file does not contain a background measurement' + + ;; Check if the spectrogram and the BKG file are recorded with the same ELUT + elut_filename = stx_date2elut_file(stx_time2any(t_axis.TIME_START[0])) + elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) + + elut_comp = STRCMP(elut_filename, elut_filename_bkg) ;; Compare ELUT tables + + if not elut_comp then $ + message, 'The background file must be recorded when the same ELUT as the science file was uploaded. Please choose a different background file that is closer in time to the science file.' + + + ;; Select indices of the energy bins (among the 32) that are actually present in the pixel data bkg file + energy_bin_mask_bkg = data_bkg.energy_bin_mask + energy_bin_idx_bkg = where(energy_bin_mask_bkg eq 1) + + energy_low_bkg = e_axis_bkg.LOW + energy_high_bkg = e_axis_bkg.HIGH + + ;; Extract energy range in common between science and background file + energy_min = max([energy_min,min(energy_low_bkg)]) + energy_max = min([energy_max,max(energy_high_bkg)]) + idx_energy_bkg = where((energy_low_bkg ge energy_min) and (energy_high_bkg le energy_max)) + + energy_bin_idx_bkg = energy_bin_idx_bkg[idx_energy_bkg] + energy_low_bkg = energy_low_bkg[idx_energy_bkg] + energy_high_bkg = energy_high_bkg[idx_energy_bkg] - elut_filename = stx_date2elut_file(start_time) - - if ~keyword_set(silent) then begin - print, 'Using ELUT file ' + elut_filename endif + ;; Extract energy range in common between science and background file + idx_energy = where((energy_low ge energy_min) and (energy_high le energy_max)) + + energy_bin_idx = energy_bin_idx[idx_energy] + energy_low = energy_low[idx_energy] + energy_high = energy_high[idx_energy] + + ct_edges = get_uniq( [energy_low,energy_high],epsilon=0.0001) + + ;;--------- Create spectrogram structure + + data_level = 4 uid = control_str.request_id + + ;; Define default file name for spectrum and srm + default, specfile, 'stx_spectrum_' + strtrim(uid,2) + '.fits' + if ~keyword_set(srmfile) then $ + srmfile = xspec ? 'stx_srm_' + strtrim(uid,2) + '_XSPEC.fits' : 'stx_srm_' + strtrim(uid,2) + '.fits' + + if n_elements(distance) ne 0 then fits_distance = distance fits_info_params = stx_fits_info_params( fits_path_data = fits_path_data, data_level = data_level, $ - distance = distance, time_shift = time_shift, fits_path_bk = fits_path_bk, uid = uid, $ + distance = fits_distance, time_shift = time_shift, fits_path_bk = fits_path_bk, uid = uid, $ generate_fits = generate_fits, specfile = specfile, srmfile = srmfile, elut_file = elut_filename, silent = silent) counts_in = data_str.counts @@ -143,10 +195,9 @@ pro stx_convert_spectrogram, fits_path_data = fits_path_data, fits_path_bk = fi pixels_used = where(data_str.pixel_masks eq 1) detectors_used = where(control_str.detector_masks eq 1) + - energy_edges_used = where(control_str.energy_bin_edge_mask eq 1, n_energy_edges) - energy_bin_mask = stx_energy_edge2bin_mask(control_str.energy_bin_edge_mask) - energy_bins = where(energy_bin_mask eq 1, n_energies) + n_energies = n_elements(energy_bin_idx) pixel_mask_used = intarr(12) pixel_mask_used[pixels_used] = 1 @@ -161,11 +212,11 @@ pro stx_convert_spectrogram, fits_path_data = fits_path_data, fits_path_bk = fi spec_in = counts_in - counts_spec = spec_in[energy_bins, *] + counts_spec = spec_in[energy_bin_idx, *] counts_spec = reform(counts_spec,[n_energies, n_times]) - counts_err = data_str.counts_err[energy_bins,*] + counts_err = data_str.counts_err[energy_bin_idx,*] counts_err = reform(counts_err,[n_energies, n_times]) @@ -175,30 +226,8 @@ pro stx_convert_spectrogram, fits_path_data = fits_path_data, fits_path_bk = fi rcr = data_str.rcr - -if keyword_set(elut_correction) then begin - stx_read_elut, ekev_actual = ekev_actual, elut_filename = elut_filename - - ave_edge = mean(reform(ekev_actual[energy_edges_used-1, pixels_used, detectors_used, 0 ], n_energy_edges, n_pixels, n_detectors), dim = 2) - ave_edge = mean(reform(ave_edge,n_energy_edges, n_detectors), dim = 2) - - edge_products, ave_edge, width = ewidth - - eff_ewidth = (e_axis.width)/ewidth + ;; Create spectrogram structure for live time correction - counts_spec = counts_spec * reproduce(eff_ewidth, n_times) - - counts_spec = reform(counts_spec,[n_energies, n_times]) - - counts_err = counts_err * reproduce(eff_ewidth, n_times) - - counts_err = reform(counts_err,[n_energies, n_times]) -endif - - - - - ;insert the information from the telemetry file into the expected stx_fsw_sd_spectrogram structure spectrogram = { $ type : "stx_fsw_sd_spectrogram", $ counts : counts_spec, $ @@ -210,62 +239,74 @@ endif detector_mask : detector_mask_used, $ rcr : rcr, $ error : counts_err} - - data_dims = lonarr(4) - data_dims[0] = n_energies - data_dims[1] = 1 - data_dims[2] = 1 - data_dims[3] = n_times - - ;get the rcr states and the times of rcr changes from the ql_lightcurves structure - ut_rcr = stx_time2any(t_axis.time_start) - find_changes, rcr, index, state, count=count - - ; ************************************************************ - ; ******************** TEMPORARY FIX ************************* - ; ***** Andrea: 2022-April-05 - ; Temporarily creation of the no_attenuation keyword in order - ; to avoid attenuation of the fitted curve. This is useful for - ; obtaining thermal fit parameters with the BKG detector in the - ; case the attenuator is inserted. We tested it with the X - ; class flare on 2021-Oct-26 and it works nicely. - if keyword_set(no_attenuation) then begin - rcr = rcr*0. - index = 0 - state = 0 - endif - ; ************************************************************ - ; ************************************************************ - - ; ******************** TEMPORARY FIX ************************* - ; ***** ECMD: 2022-Jun-27 - ; As the reported time of the RCR status change can be inaccurate - ; up to several seconds correct this by finding the times where there is a - ; large change in counts in the counts of the 5 - 6 keV energy bin. - ; find all time intervals where the difference between adjacent bins is large - if max(rcr) gt 0 then begin; skip if in the standard state of RCR0 for the full time range - - jumps = where(abs( counts_spec[2,*] - shift(counts_spec[2,*],-1) ) gt 1e4) - ; include the starting bin - jumps = [0, jumps] - ; as the attenuator motion can be present in two consecutive bins select only the first - idx_jumps = where(abs(jumps - shift(jumps, -1)) gt 2) - jumps_use= [jumps[idx_jumps]] - ; each transition should correspond close in time to a recorded transition in the FITS file - ; adjust the time indexes of these transitions to the closest jumps - index = jumps_use - + + + ;; Live time correction + livetime_data = stx_spectrogram_livetime( spectrogram, corrected_counts = counts_spec, corrected_error = counts_spec_error, level = data_level ) + livetime = livetime_data.livetime + livetime_err = livetime_data.livetime_err + + ;;----- BKG subtraction + + if keyword_set(fits_path_bk) then begin + + counts_bkg = data_bkg.COUNTS + counts_error_bkg = data_bkg.COUNTS_ERR + counts_bkg = counts_bkg[energy_bin_idx_bkg,*,*] + counts_error_bkg = counts_error_bkg[energy_bin_idx_bkg,*,*] + + live_time_bkg_data = stx_cpd_livetime(data_bkg.TRIGGERS, data_bkg.TRIGGERS_ERR, t_axis_bkg) + live_time_bkg = live_time_bkg_data.LIVE_TIME_BINS + live_time_error_bkg = live_time_bkg_data.LIVE_TIME_BINS_ERR + + live_time_bkg_rep = transpose(cmreplicate(live_time_bkg, [n_elements(energy_bin_idx),12]), [1,2,0]) + live_time_error_bkg_rep = transpose(cmreplicate(live_time_error_bkg, [n_elements(energy_bin_idx),12]), [1,2,0]) + + ;; Normalize by livetime + count_rates_bkg = f_div( counts_bkg, live_time_bkg_rep ) + count_rates_bkg_error = count_rates_bkg * sqrt( f_div( counts_error_bkg, counts_bkg )^2. + f_div( live_time_error_bkg_rep, live_time_bkg_rep )^2. ) + + ;; Sum over detectors and pixels + counts_spec_bkg = total(total(count_rates_bkg, 2) ,2) + counts_spec_bkg_error = sqrt(total(total(count_rates_bkg_error^2., 2) ,2)) + + if n_times gt 1 then begin + + counts_spec_bkg = cmreplicate( counts_spec_bkg, n_times ) + counts_spec_bkg_error = cmreplicate( counts_spec_bkg_error, n_times ) + + endif + + counts_spec = counts_spec - counts_spec_bkg + counts_spec_error = sqrt(counts_spec_error^2. + counts_spec_bkg_error^2.) + + endif - ; ************************************************************ - ;add the rcr information to a specpar structure so it can be included in the spectrum FITS file - specpar = { sp_atten_state : {time:ut_rcr[index], state:state}, flare_xyoffset : fltarr(2), use_flare_xyoffset:0 } - - stx_convert_science_data2ospex, spectrogram = spectrogram, specpar = specpar, time_shift = time_shift, data_level = data_level, $ - data_dims = data_dims, fits_path_bk = fits_path_bk, fits_path_data = fits_path_data, $ - elut_correction = elut_correction, eff_ewidth = eff_ewidth, fits_info_params = fits_info_params, sys_uncert = sys_uncert, $ - aux_fits_file = aux_fits_file, flare_location_hpc = flare_location_hpc, flare_location_stx = flare_location_stx, $ - silent = silent, background_data = background_data, plot = plot, generate_fits = generate_fits, xspec=xspec, ospex_obj = ospex_obj + ;;------ Multiply by live time. The units of the spectrogram are counts + spec = counts_spec * livetime + spec_error = abs(spec) * sqrt( f_div(counts_spec_error,counts_spec)^2. + f_div(livetime_err,livetime)^2. ) + + + livetime_fraction = livetime_data.livetime_fraction + livetime_fraction_dim = size(livetime_fraction, /dim) + if n_elements(livetime_fraction_dim) ge 2 then livetime_fraction = reform(livetime_fraction[0,*]) + + ;; Create spectrogram structure + spec_data = {type : 'stx_spectrogram', $ + data : spec, $ + t_axis : t_axis, $ + e_axis : e_axis, $ + ltime : livetime_fraction, $ + attenuator_state : data_str.RCR , $ + error : spec_error} + + ;;------------------------------------------------------ + + stx_convert_spectrogram2ospex, spec_data, pixel_mask_used, detector_mask_used, fits_info_params, ct_edges, $ + no_attenuation=no_attenuation, flare_location_stx=flare_location_stx, time_shift = time_shift, $ + sys_uncert=sys_uncert, silent=silent, plot=plot, xspec=xspec, ospex_obj = ospex_obj, $ + tailing=tailing, include_damage=include_damage, _extra=extra end diff --git a/stix/idl/processing/spectrogram/stx_convert_spectrogram2ospex.pro b/stix/idl/processing/spectrogram/stx_convert_spectrogram2ospex.pro new file mode 100644 index 00000000..e5c9d91b --- /dev/null +++ b/stix/idl/processing/spectrogram/stx_convert_spectrogram2ospex.pro @@ -0,0 +1,305 @@ +;--------------------------------------------------------------------------- +;+ +; :project: +; STIX +; +; :name: +; stx_convert_spectrogram2ospex +; +; :description: +; This procedure reads a spectrogram structure (which is created from either a spectrogram file or a compressed pixel data file) and converts it to a spectrogram +; file which can be read in by OSPEX. +; +; :categories: +; spectroscopy +; +; :keywords: +; +; spec_data: in, mandatory, type="stx_spectrogram" +; Spectrogram structure created either by 'stx_convert_pixel_data.pro' or 'stx_convert_spectrogram.pro' +; +; pixel_mask_used: in, mandatory, type = "int array" +; 12-element array with entries equal to 1 if the corresponding detector pixels are used (and 0 otherwise) +; +; detector_mask_used: in, mandatory, type = "int array" +; 32-element array with entries equal to 1 if the corresponding detector is used (and 0 otherwise) +; +; fits_info_params: in, mandatory, type="struct" +; Structure created by 'stx_fits_info_params.pro' +; +; ct_edges: in, mandatory, type="float array" +; Energy edges of the count energy bins of the spectrogram +; +; no_attenuation: in, optional, type="bool" +; If set, avoid attenuation of the fitted curve. This is useful for obtaining thermal fit parameters with the BKG detector (temporary fix by Andrea) +; +; flare_location_stx : in, optional, type="2 element float array" +; the location of the flare (X,Y) in the STIX imaging frame [arcsec]. Default, [0.,0.]. It is used to compute the grid transmission +; +; time_shift : in, optional, type="float", default="0." +; The difference in seconds in light travel time between the Sun and Earth and the Sun and Solar Orbiter +; i.e. Time(Sun to Earth) - Time(Sun to S/C) +; +; sys_uncert : in, type="float", default="0.05" +; The fractional systematic uncertainty to be added +; +; silent : in, type="int", default="0" +; If set prevents informational messages being displayed. +; +; plot : in, type="boolean", default="1" +; If set open OSPEX GUI and plot lightcurve in standard quicklook energy bands where there is data present +; +; xspec : in, type="boolean", default="0" +; If set, generate SRM file compatible with XSPEC rather than OSPEX. +; +; ospex_obj : out, type="OSPEX object" + +; +; :history: 07-May-2026 - Massa P. (FHNW), created +; +;- + +pro stx_convert_spectrogram2ospex, spec_data, pixel_mask_used, detector_mask_used, fits_info_params, ct_edges, $ + no_attenuation=no_attenuation, flare_location_stx=flare_location_stx, time_shift = time_shift, $ + sys_uncert=sys_uncert, silent=silent, plot=plot, xspec=xspec, ospex_obj = ospex_obj, $ + tailing=tailing, include_damage=include_damage, _extra=extra + + default, no_attenuation, 0 + default, sys_uncert, 0.05 + + if n_elements(time_shift) eq 0 then begin + if ~keyword_set(silent) then begin + message, 'Time shift value not set, using default value of 0 [s].', /info + print, 'File averaged values can be obtained from the FITS file header' + print, 'using stx_get_header_corrections.pro.' + endif + time_shift = 0. + endif + + ;;****************** Define parameters + + ut_rcr = stx_time2any(spec_data.t_axis.time_start) + rcr = spec_data.attenuator_state + spec = spec_data.data + + detectors_used = where(detector_mask_used eq 1, n_det) + pixels_used = where(pixel_mask_used eq 1, n_pix) + + ;; Print warning if fine-resolution sub-collimators are considered + subc_fine_label = ['1a','1b','1c','2a','2b','2c'] + subc_fine_index = stx_label2ind(subc_fine_label) + subc_fine_used = detector_mask_used[subc_fine_index] + + idx = where(subc_fine_used eq 1b, n_det) + + if n_det gt 0 then message, [" ", " ", "The high-energy calibration of sub-collimator " + subc_fine_label[idx] + " is not completed yet. This could affect the spectral fit results.", " ", " "], /continue + + ;;***************** CHECK FOR RCR CHANGES + + find_changes, rcr, index, state, count=count + ; ************************************************************ + ; ******************** TEMPORARY FIX ************************* + ; ***** Andrea: 2022-April-05 + ; Temporarily creation of the no_attenuation keyword in order + ; to avoid attenuation of the fitted curve. This is useful for + ; obtaining thermal fit parameters with the BKG detector in the + ; case the attenuator is inserted. We tested it with the X + ; class flare on 2021-Oct-26 and it works nicely. + if keyword_set(no_attenuation) then begin + rcr = rcr*0. + index = 0 + state = 0 + endif + ; ************************************************************ + ; ************************************************************ + + ; ******************** TEMPORARY FIX ************************* + ; ***** ECMD: 2022-Jun-27 + ; As the reported time of the RCR status change can be inaccurate + ; up to several seconds correct this by finding the times where there is a + ; large change in counts in the counts of the 5 - 6 keV energy bin. + ; find all time intervals where the difference between adjacent bins is large + if max(rcr) gt 0 then begin; skip if in the standard state of RCR0 for the full time range + jumps = where(abs(total(spec,1) - shift(total(spec,1),-1)) * 24. / n_elements(detectors_used) gt 1e4) + ;jumps = where(abs(total(spec,1) - shift(total(spec,1),-1)) gt 1e4) + ; include the starting bin + jumps = [0, jumps] + ; as the attenuator motion can be present in two consecutive bins select only the first + idx_jumps = where(abs(jumps - shift(jumps, -1)) gt 2) + jumps_use= [jumps[idx_jumps]] + ; each transition should correspond close in time to a recorded transition in the FITS file + ; adjust the time indexes of these transitions to the closest jumps + closest_jumps = value_closest(jumps_use, index) + index = jumps_use[closest_jumps] + + endif + + ;add the rcr information to a specpar structure so it can be included in the spectrum FITS file + if not keyword_set(flare_location_stx) then begin + flare_location_stx = [0.,0.] + use_flare_xyoffset = 0 + endif else begin + use_flare_xyoffset = 1 + endelse + + specpar = { sp_atten_state : {time:ut_rcr[index], state:state}, flare_xyoffset : flare_location_stx, use_flare_xyoffset:use_flare_xyoffset } + + ;;***************** CREATE SRM + + pixel_mask =detector_mask_used ## pixel_mask_used + + transmission = read_csv(loc_file( 'stix_transmission_highres_20251110.csv', path = getenv('STX_GRID'))) + + emin = 4 + emax = 150 + phe = transmission.(0) + phe = phe[where(phe gt emin-1 and phe lt 3.5*emax)] + edge_products, phe, mean = mean_phe, width = w_phe + ph_edges = [mean_phe[0] - w_phe[0], mean_phe] + + distance = fits_info_params.distance + dist_factor = 1./(distance^2.) + + ;make the srm for the appropriate pixel mask and energy edges + ;srm = stx_build_pixel_drm(ct_edges, pixel_mask, ph_energy_edges = ph_edges, dist_factor = dist_factor, tailing = tailing, include_damage = include_damage, _extra=extra) + ph_edges = get_uniq( [ph_edges,ct_edges],epsilon=0.0001) + + ;; Compute subc. transmission + edge_products,ph_edges, mean=ph_in + subc_transmission = stx_subc_transmission(flare_location_stx, ph_in, _extra=extra) + + if n_elements(detectors_used) eq 1 then begin + grid_factor = reform(subc_transmission[*,detectors_used]) + endif else begin + grid_factor = total(subc_transmission[*,detectors_used],2) / n_elements(detectors_used) + endelse + + + ;; Check if BKG or CFL detectors are used + + idx_cfl = where(detectors_used eq 8, n_cfl) + idx_bkg = where(detectors_used eq 9, n_bkg) + + if n_cfl eq 1 then message, "CFL detector can not be selected for spectral fitting." + if (n_elements(detectors_used) gt 1) and (n_bkg eq 1) then message, "BKG detector can not be selected together with imaging detectors for spectral fitting." + + if n_elements(detectors_used) eq 1 then if detectors_used eq 9 then begin + grid_transmission_file = concat_dir(getenv('STX_GRID'), 'real_bkg_grid_transmission.txt') + readcol, grid_transmission_file, bk_grid_factors, format = 'f', skip = 2, silent = silent + + grid_factor = average(bk_grid_factors[pixels_used]) + + endif + + ;; Creates appropriate SRM for different attenuator states + rcr_states = specpar.sp_atten_state.state + rcr_states = rcr_states[uniq(rcr_states, sort(rcr_states))] + nrcr_states = n_elements(rcr_states) + + srm_atten = replicate( {rcr:0, srm:fltarr(n_elements(ct_edges)-1,n_elements(ph_edges)-1)},nrcr_states ) + + for i =0, nrcr_states-1 do begin + ;make the srm for the appropriate pixel mask and energy edges + rcr = rcr_states[i] + + srm = stx_build_pixel_drm(ct_edges, pixel_mask, rcr = rcr, grid_factor= grid_factor, ph_energy_edges = ph_edges, dist_factor = dist_factor, tailing = tailing, include_damage = include_damage, xspec=xspec, _extra=extra) + srm_atten[i].srm = srm.smatrix + srm_atten[i].rcr = rcr + + endfor + + ;;***************** SAVE FITS + detector_label = stx_det_mask2label(detector_mask_used) + pixel_label = stx_pix_mask2label(pixel_mask_used) + + ospex_obj = ospex(/no) + + ;if the fits keyword is set write the spectrogram and srm data to fits files and then read them in to the ospex object + if fits_info_params.generate_fits eq 1 then begin + + utime = transpose([stx_time2any( spec_data.t_axis.time_start )]) + + specfilename = fits_info_params.specfile + srmfilename = fits_info_params.srmfile + + fits_info_params.grid_factor.add, grid_factor + fits_info_params.detused = detector_label + ', Pixels: ' + pixel_label + + if keyword_set(xspec) then begin + ;xspec in general works with energy depandent systematic errors + e_axis = spec_data.e_axis + n_energies = n_elements(e_axis.mean) + sys_err = fltarr(n_energies) + + n_times = n_elements(spec_data.t_axis.time_start) + + idx_below10kev = where(e_axis.mean lt 10, cb10) + sys_err[*] = 0.03 + if cb10 gt 0 then sys_err[idx_below10kev] = 0.05 + idx_below7kev = where(e_axis.mean lt 7, cb7) + if cb7 gt 0 then sys_err[idx_below7kev] = 0.07 + + sys_err = rebin(sys_err, n_energies,n_times) + endif + + stx_write_ospex_fits, spectrum = spec_data, srmdata = srm, specpar = specpar, time_shift = time_shift, $ + srm_atten = srm_atten, specfilename = specfilename, srmfilename = srmfilename, ph_edges = ph_edges, $ + fits_info_params = fits_info_params, xspec = xspec, silent = silent + + ospex_obj->set, spex_file_reader = 'stx_read_sp' + ospex_obj->set, spex_specfile = specfilename ; name of your spectrum file + ospex_obj->set, spex_drmfile = srmfilename + + endif else begin + ;if the generate_fits keyword is not set use the spex_user_data strategy to pass in the data directly to the ospex object + e_axis = spec_data.e_axis + t_axis = spec_data.t_axis + + energy_edges = e_axis.edges_2 + Edge_Products, ph_edges, edges_2 = ph_edges2 + + utime2 = transpose(stx_time2any( [[t_axis.time_start], [t_axis.time_end]] )) + + livetime_frac = spec_data.LTIME + nchan = n_elements( ct_edges ) - 1 + duration = spec_data.t_axis.duration + nrows = n_elements(duration) + duration_array=rebin(duration, nrows, nchan) + livetime_frac_array = rebin([livetime_frac], nrows, nchan) + data=f_div(float(spec),transpose(duration_array*livetime_frac_array)) + data_error = f_div(spec_data.ERROR,transpose(duration_array*livetime_frac_array)) + + ospex_obj->set, spex_data_source = 'spex_user_data' + ospex_obj->set, spectrum = data, $ + spex_ct_edges = energy_edges, $ + spex_ut_edges = utime2, $ + livetime = spec_data.LTIME * spec_data.t_axis.duration, $ + errors = data_error + + srm = rep_tag_name(srm,'smatrix','drm') + ospex_obj->set, spex_respinfo = srm + ospex_obj->set, spex_area = srm.area + ospex_obj->set, spex_detectors = 'STIX' + ospex_obj->set, spex_drm_ct_edges = energy_edges + ospex_obj->set, spex_drm_ph_edges = ph_edges2 + endelse + + ospex_obj->set, spex_uncert = sys_uncert + ospex_obj->set, spex_error_use_expected = 0 + + counts_str = ospex_obj->getdata(spex_units='counts') + origunits = ospex_obj->get(/spex_data_origunits) + origunits.data_name = 'STIX' + ospex_obj->set, spex_data_origunits = origunits + + if keyword_set(plot) then begin + ospex_obj ->gui + ospex_obj ->set, spex_eband = get_edges([4.,10.,15.,25, 50, 84.], /edges_2) + ospex_obj ->plot_time, spex_units='flux', /show_err, obj = plotman_object + endif + +end + + + + \ No newline at end of file diff --git a/stix/idl/processing/spectrogram/stx_read_spectrogram_fits_file.pro b/stix/idl/processing/spectrogram/stx_read_spectrogram_fits_file.pro index 70ef2a16..be8da780 100644 --- a/stix/idl/processing/spectrogram/stx_read_spectrogram_fits_file.pro +++ b/stix/idl/processing/spectrogram/stx_read_spectrogram_fits_file.pro @@ -80,7 +80,7 @@ ; 15-Mar-2023 - ECMD (Graz), updated to handle release version of L1 FITS files ; 27-Mar-2023 - ECMD (Graz), added check for duration shift already applied in FITS file ; 13-Sep-2023 - ECMD (Graz), respect shift_duration if file is "possibly summed on board". -; +; 07-May-2026 - Massa P. (FHNW), added 'energy_bin_mask' to the output 'data' structure ;- pro stx_read_spectrogram_fits_file, fits_path, time_shift, primary_header = primary_header, data_str = data, data_header = data_header, control_str = control, $ control_header= control_header, energy_str = energy, energy_header = energy_header, t_axis = t_axis, e_axis = e_axis, $ @@ -121,7 +121,7 @@ pro stx_read_spectrogram_fits_file, fits_path, time_shift, primary_header = prim edges_used = where( control.energy_bin_edge_mask eq 1, nedges) energy_bin_mask = stx_energy_edge2bin_mask(control.energy_bin_edge_mask) energies_used = where(energy_bin_mask eq 1) - + data.counts_comp_err = sqrt(data.counts_comp_err^2. + data.counts) data.triggers_comp_err = sqrt( data.triggers_comp_err^2. + data.triggers) @@ -135,7 +135,7 @@ pro stx_read_spectrogram_fits_file, fits_path, time_shift, primary_header = prim if ~keyword_set(keep_short_bins) and (anytim(hstart_time) lt anytim('2020-11-25T00:00:00') ) then $ message, 'Automatic short bin removal should not be attempted on observations before 25-Nov-20' - shift_duration = shift_duration && ~duration_shifted + shift_duration = shift_duration && ~duration_shifted if keyword_set(shift_duration) and (anytim(hstart_time) gt anytim('2021-12-09T00:00:00') ) then $ message, 'Shift of duration with respect to time bins is no longer needed after 09-Dec-21' @@ -278,6 +278,8 @@ pro stx_read_spectrogram_fits_file, fits_path, time_shift, primary_header = prim t_axis.time_end = t_end t_axis.duration = duration + energy_edges_2 = transpose([[energy.e_low], [energy.e_high]]) + data = {time: time_bin_center,$ timedel:duration , $ triggers:triggers, $ @@ -286,12 +288,11 @@ pro stx_read_spectrogram_fits_file, fits_path, time_shift, primary_header = prim counts_err: counts_err ,$ rcr: rcr,$ pixel_masks: pixel_masks,$ - control_index:control_index} - - energy_edges_2 = transpose([[energy.e_low], [energy.e_high]]) + control_index:control_index,$ + detector_masks:control.detector_masks} if control.energy_bin_edge_mask[0] and ~keyword_set(use_discriminators) then begin - + control.energy_bin_edge_mask[0] = 0 data.counts[0,*] = 0. data.counts_err[0,*] = 0. @@ -311,7 +312,7 @@ pro stx_read_spectrogram_fits_file, fits_path, time_shift, primary_header = prim edges_used = where( control.energy_bin_edge_mask eq 1, nedges) energy_bin_mask = stx_energy_edge2bin_mask(control.energy_bin_edge_mask) energies_used = where(energy_bin_mask eq 1) - + ; changed 2023-03-15: Now only the used energies are in table energy, therefore: edge_products, energy_edges_2, edges_1 = energy_edges_1 @@ -320,6 +321,7 @@ pro stx_read_spectrogram_fits_file, fits_path, time_shift, primary_header = prim e_axis = stx_construct_energy_axis(energy_edges = energy_edges_1 + energy_shift, select = use_edges ) e_axis.low_fsw_idx = energies_used e_axis.high_fsw_idx = energies_used - - -end \ No newline at end of file + + data = add_tag(data,energy_bin_mask,'energy_bin_mask') + +end From edd1eb8175ef003082405c69d958c2cf62f0f260 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Tue, 12 May 2026 15:34:15 +0200 Subject: [PATCH 12/13] Print start and time of the calibration file --- stix/idl/io/stx_get_calibration_file.pro | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/stix/idl/io/stx_get_calibration_file.pro b/stix/idl/io/stx_get_calibration_file.pro index 68729f66..64d7144e 100644 --- a/stix/idl/io/stx_get_calibration_file.pro +++ b/stix/idl/io/stx_get_calibration_file.pro @@ -152,6 +152,7 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe found_files = found_files_array[idx_elut] duration = duration[idx_elut] start_time_file = start_time_file[idx_elut] + end_time_file = end_time_file[idx_elut] mid_time_file = mid_time_file[idx_elut] ;; Select calibration file with largest duration @@ -166,6 +167,7 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe found_files = found_files[idx_duration] duration = duration[idx_duration] start_time_file = start_time_file[idx_duration] + end_time_file = end_time_file[idx_duration] mid_time_file = mid_time_file[idx_duration] ;; If multiple files have the largest duration, select the closest one to the input time range @@ -180,12 +182,14 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe idx_file = idx_time_diff[0] selected_file = found_files[idx_file] - selected_file_time = start_time_file[idx_file] + selected_file_time_start = start_time_file[idx_file] + selected_file_time_end = end_time_file[idx_file] endif else begin selected_file = found_files - selected_file_time = start_time_file + selected_file_time_start = start_time_file + selected_file_time_end = end_time_file endelse @@ -196,7 +200,8 @@ function stx_get_calibration_file, start_time, end_time, out_dir=out_dir, clobbe endelse ;; Download file - message, [" ", " ", "Download STIX calibration file recorded at " +anytim(selected_file_time, /vms), " ", " "], /continue + message, [" ", " ", "Download STIX calibration file recorded between " +anytim(selected_file_time_start, /vms)+$ + " and "+anytim(selected_file_time_end, /vms), " ", " "], /continue sock_copy, selected_file, out_name, local_file=out_file, out_dir = out_dir, clobber=clobber return, out_file From 79afc9f7b53704391c2faac016c081be59e77372 Mon Sep 17 00:00:00 2001 From: paolomassa Date: Tue, 19 May 2026 16:31:25 +0200 Subject: [PATCH 13/13] Fixed bug --- stix/idl/processing/spectrogram/stx_convert_spectrogram.pro | 1 - 1 file changed, 1 deletion(-) diff --git a/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro b/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro index dd840e3b..4d53507e 100644 --- a/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro +++ b/stix/idl/processing/spectrogram/stx_convert_spectrogram.pro @@ -103,7 +103,6 @@ pro stx_convert_spectrogram, fits_path_data = fits_path_data, fits_path_bk = fit time_shift = 0. endif - default, shift_duration, 0 default, sys_uncert, 0.05 default, silent, 0 default, plot, 1