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/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/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_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/io/stx_get_calibration_file.pro b/stix/idl/io/stx_get_calibration_file.pro new file mode 100644 index 00000000..64d7144e --- /dev/null +++ b/stix/idl/io/stx_get_calibration_file.pro @@ -0,0 +1,211 @@ +;+ +; +; 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 +; 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 + + 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='/') + + ;; Concatenate files from the input date and from the day before/after + type_path = '/fits/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) + + 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) + + 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) + + path_after = get_fid(this_start_time_after,this_end_time_after,/full,delim='/') + + 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 n_elements(found_files_array) eq 0 then begin + + message, $ + [" ", " ", "No STIX calibration file was found within 1 day before or after " +date_path[0]+". Please, download an appropriate calibration file manually.", " ", " "] + + endif else begin + + start_time_file = [] + end_time_file = [] + elut_check = [] + + ;; Extract file names + 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' + + 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) + + 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) + 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)] + + elut_check = [elut_check, (elut_input_time eq stx_date2elut_file(this_start_time))] + + endfor + + ;; Compute mid point of the time interval of the calibration file + mid_time_file = (start_time_file + end_time_file) / 2. + + ;; Compute duration file + duration = end_time_file - start_time_file + + ;; ELUT check + idx_elut = where(elut_check, n_elut_check) + + 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] + end_time_file = end_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] + 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 + if n_duration gt 1 then begin + + ;; 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] + + selected_file = found_files[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 = start_time_file + selected_file_time_end = end_time_file + + endelse + + endelse + + endelse + + endelse + + ;; Download file + 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 + + + +end 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/stx_calibration_data.pro b/stix/idl/processing/calib_data/stx_calibration_data.pro new file mode 100644 index 00000000..20d1dc08 --- /dev/null +++ b/stix/idl/processing/calib_data/stx_calibration_data.pro @@ -0,0 +1,43 @@ +;+ +; +; NAME: +; +; stx_calibration_data +; +; PURPOSE: +; +; Defines a STIX calibration data structure +; +; CALLING SEQUENCE: +; +; cal_data = stx_calibration_data() +; +; OUTPUTS: +; +; 'stx_calibration_data' structure containing data from L2 calibration fits file +; +; +; HISTORY: March 2026, Massa P. (FHNW), created +; +; CONTACT: +; paolo.massa@fhnw.ch +;- + +function stx_calibration_data + + struct = { $ + type : 'stx_calibration_data', $ + time_start : stx_time(), $ + time_end : stx_time(), $ + t_mean : stx_time(), $ + energy_bin_low : dblarr(32, 12, 32), $ + energy_bin_high : dblarr(32, 12, 32), $ + gain : fltarr(12, 32), $ + offset : fltarr(12, 32), $ + live_time : double(0), $ + elut_name : '' $ + } + + return, struct + +end diff --git a/stix/idl/processing/imaging/stx_em.pro b/stix/idl/processing/imaging/stx_em.pro index 86483601..af057c9c 100644 --- a/stix/idl/processing/imaging/stx_em.pro +++ b/stix/idl/processing/imaging/stx_em.pro @@ -13,17 +13,23 @@ ; pixel data structure containing photon counts per time, energy, detector, and pixel ; (the counts registered are summed as if they were recorded by 4 virtual pixels per ; detector). For details on the summation see the header of 'stx_pixel_data_sum.pro'. +; +; aux_data: auxiliary data structure corresponding to the selected time range +; ; KEYWORDS: -; DET_USED: array containing the indices of detector used -; (default is 0-31, 8 and 9 excluded) ; IMSIZE: output map size in pixels (default is [129, 129]) ; PIXEL: pixel size in arcsec (default is [1., 1.]) +; MAPCENTER: two-element array containing the coordinates of the center of the map to reconstruct (STIX coordinate frame, arcsec). +; Default, (0.,0.) +; 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. Default, (0,0) +; SUBC_INDEX: array containing the indices of the selected imginging detectors. Default, indices of +; the detectors from 10 to 3 ; MAXITER: max number of iterations (default is 5000) ; TOLERANCE: parameter for the stopping rule (default is 0.01) ; SILENT: if not set, plots the STD (variable to test convergence) and the ; C-statistic every 25 iterations ; MAKEMAP: if set, returns the map structure. Otherwise returns the 2D matrix -; XYOFFSET: array containing the map center coordinates. ; ; RETURNS: ; an image (2D matrix) or an image map in the structure format provided by the @@ -40,12 +46,13 @@ ; in the EM iterative scheme ; July 2023, Massa P., made it compatible with the new definition of (u,v)-points (see stx_uv_points) ; March 2026, Massa P., new subcollimator transmission is implemented +; April 2026, Massa P., made it compatible with the fact that bkg subtraction is performed in 'stx_construct_pixel_data' ; ;CONTACT: paolo.massa@fhnw.ch FUNCTION stx_em, pixel_data_summed, aux_data, imsize=imsize, pixel=pixel, $ mapcenter=mapcenter, xy_flare=xy_flare, subc_index=subc_index, $ - maxiter=maxiter, tolerance=tolerance, silent=silent, makemap=makemap + maxiter=maxiter, tolerance=tolerance, silent=silent, makemap=makemap, _extra=extra default, subc_index, stx_label2ind(['3a','3b','3c','4a','4b','4c','5a','5b','5c','6a','6b','6c',$ '7a','7b','7c','8a','8b','8c','9a','9b','9c','10a','10b','10c']) @@ -108,7 +115,7 @@ phase_corr *= !dtor ;;**************** 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 @@ -125,9 +132,6 @@ n_det_used = n_elements(subc_index) countrates = pixel_data_summed.COUNT_RATES[subc_index,*] y = reform(countrates, n_det_used*4) -countrates_bkg = pixel_data_summed.COUNT_RATES_BKG[subc_index,*] -b = reform(countrates_bkg, n_det_used*4) - ;;**************** EXPECTATION MAXIMIZATION ALGORITHM ;Initialization @@ -141,17 +145,17 @@ if ~keyword_set(silent) then print, 'EM iterations: ' & print, 'N. Iter: ST ; Loop of the algorithm for iter = 1, maxiter do begin Hx = H # x - z = f_div(y , Hx + b) + z = f_div(y , Hx) Hz = H ## z x = x * transpose(f_div(Hz, Ht1)) - cstat = 2. / n_elements(y[y_index]) * total(y[y_index] * alog(f_div(y[y_index],Hx[y_index] + b[y_index])) + Hx[y_index] + b[y_index] - y[y_index]) + cstat = 2. / n_elements(y[y_index]) * total(y[y_index] * alog(f_div(y[y_index],Hx[y_index])) + Hx[y_index] - y[y_index]) ; Stopping rule if iter gt 10 and (iter mod 25) eq 0 then begin emp_back_res = total((x * (Ht1 - Hz))^2) - std_back_res = total(x^2 * (f_div(1.0, Hx + b) # H2)) + std_back_res = total(x^2 * (f_div(1.0, Hx) # H2)) std_index = f_div(emp_back_res, std_back_res) if ~keyword_set(silent) then print, iter, std_index, cstat @@ -165,8 +169,8 @@ x_im = reform(x, imsize[0],imsize[1]) ;;**************** Create visibility structure to be used in 'stx_make_map' -vis = stx_pixel_data_summed2visibility(pixel_data_summed,subc_index=subc_index, mapcenter=mapcenter) -vis = stx_calibrate_visibility(vis, xy_flare=xy_flare) +vis = stx_pixel_data_summed2visibility(pixel_data_summed,subc_index=subc_index) +vis = stx_calibrate_visibility(vis, mapcenter=mapcenter, xy_flare=xy_flare) em_map = stx_make_map(x_im, aux_data, pixel, "EM", vis) diff --git a/stix/idl/processing/livetime/stx_cpd_livetime.pro b/stix/idl/processing/livetime/stx_cpd_livetime.pro new file mode 100644 index 00000000..fcf046f9 --- /dev/null +++ b/stix/idl/processing/livetime/stx_cpd_livetime.pro @@ -0,0 +1,59 @@ +;+ +; +; NAME: +; +; stx_cpd_livetime +; +; PURPOSE: +; +; Calculate the live time and associated uncertainty for a Compressed Pixel Data (CPD) file. +; +; CALLING SEQUENCE: +; +; live_time_data = stx_cpd_livetime(triggers, triggers_err, t_axis) +; +; INPUTS: +; +; triggers: array (dimensions: 16 x num. time bins) containing the trigger values +; triggers_err: array (dimensions: 16 x num. time bins) containing the uncertainty on the trigger values +; t_axis: 'stx_time_axis' structure returned by 'stx_read_pixel_data_fits_file' +; +; OUTPUTS: +; +; Structure containing: +; - LIVE_TIME_BINS: array (dimensions: 32 × number of time bins) containing the live time value +; for each time bin in the CPD file +; - LIVE_TIME_BINS_ERR: array (dimensions: 32 × number of time bins) containing the uncertainty on +; the live time value for each time bin in the CPD file +; - LIVETIME_FRACTION: array (dimensions: 32 × number of time bins) containing the value of +; the live time fraction for each time bin in the CPD file +; - LIVETIME_FRACTION_ERR: array (dimensions: 32 × number of time bins) containing the uncertainty on +; the live time fraction for each time bin in the CPD file +; +; HISTORY: +; May 2025, Massa P. (FHNW), first release +; +; CONTACT: +; paolo.massa@fhnw.ch +;- + +function stx_cpd_livetime, triggers, triggers_err, t_axis + + triggergram = stx_triggergram(triggers, triggers_err, t_axis) + + livetime_fraction_data = stx_livetime_fraction(triggergram) + livetime_fraction = livetime_fraction_data.livetime_fraction + livetime_fraction_err = livetime_fraction_data.livetime_fraction_err + + 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_bins_err = livetime_fraction_err * duration_time_bins + + return, {live_time_bins: live_time_bins,$ + live_time_bins_err: live_time_bins_err,$ + livetime_fraction: livetime_fraction, $ + livetime_fraction_err: livetime_fraction_err} + +end diff --git a/stix/idl/processing/livetime/stx_livetime_fraction.pro b/stix/idl/processing/livetime/stx_livetime_fraction.pro index e1bab61e..436a8f20 100755 --- a/stix/idl/processing/livetime/stx_livetime_fraction.pro +++ b/stix/idl/processing/livetime/stx_livetime_fraction.pro @@ -1,10 +1,10 @@ ;+ ; :Categories: ; STIX imaging and spectroscopy -; +; ; :Name: ; stx_livetime_fraction -; +; ; :Examples: ; livetime_fraction = stx_livetime_fraction( triggergram, det_select, tau_array = tau_array ) ; @@ -21,7 +21,7 @@ ; ADG_IDX INT Array[16] ;accumulator id's 1-16 ; T_AXIS STRUCT -> 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_spectrogram_livetime.pro b/stix/idl/processing/livetime/stx_spectrogram_livetime.pro index 95d17424..9b7a16b5 100644 --- a/stix/idl/processing/livetime/stx_spectrogram_livetime.pro +++ b/stix/idl/processing/livetime/stx_spectrogram_livetime.pro @@ -37,7 +37,9 @@ ; :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 +; 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 @@ -48,55 +50,42 @@ 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] - 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 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. + livetime = time_bin_duration * livetime_fraction + livetime_err = time_bin_duration * livetime_fraction_err + + corrected_counts = f_div(spectrogram.counts,livetime) - 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_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/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.pro b/stix/idl/processing/livetime/stx_triggergram.pro old mode 100644 new mode 100755 index bbb3d283..f0292e9b --- 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 +end 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/pixel_data/stx_construct_pixel_data.pro b/stix/idl/processing/pixel_data/stx_construct_pixel_data.pro index 44d0871a..c40fcc36 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,40 @@ 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])) + elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) + + ;; 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 +298,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 +483,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 +497,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_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro index b5d93f89..162bdf29 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,33 @@ ; 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. (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, $ - 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, silent, 0 default, xspec, 0 + default, sys_uncert, 0.05 + default, delta_time_min, 20. + default, tailing, 1 + default, include_damage, 1 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,32 +139,82 @@ 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 + ;; 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' - 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) + + ;; 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 + elut_filename = stx_date2elut_file(stx_time2any(t_axis.TIME_START[0])) + 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) - 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 +234,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,140 +242,282 @@ 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]) + elut_filename_bkg = stx_date2elut_file(stx_time2any(t_axis_bkg.TIME_START)) - counts_spec = spec_in[energy_bins,*, *] + ;; Compare ELUT tables + elut_comp = STRCMP(elut_filename, elut_filename_bkg) - counts_spec = reform(counts_spec,[n_energies, n_detectors, 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 = reform(data_str.counts_err,[dim_counts[0:2], n_times]) + 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 = sqrt(total(reform(counts_err[*,pixels_used,detectors_used,*]^2.,[32,n_pixels,n_detectors,n_times]),2)) + endif else begin - counts_err = reform(counts_err,[dim_counts[0],n_detectors, n_times]) + counts_bkg = dblarr(size(counts, /dim)) + counts_error_bkg = dblarr(size(counts, /dim)) - counts_err = counts_err[energy_bins,*, *] - - counts_err = reform(counts_err,[n_energies, n_detectors, n_times]) + endelse - triggers = transpose(reform(data_str.triggers,[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 - triggers_err = transpose(reform(data_str.triggers_err,[16, n_times])) + 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]) - rcr = data_str.rcr - - if keyword_set(elut_correction) then begin + 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 - 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) + endif else begin - edge_products, ave_edge, width = ewidth + live_time_bkg = dblarr(32) + 1. + live_time_error_bkg = dblarr(32) - eff_ewidth = (e_axis.width)/ewidth + 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]) - counts_spec = counts_spec * reform(reproduce(eff_ewidth, n_detectors*n_times),n_energies, n_detectors, n_times) + ;; 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. ) - counts_spec = reform(counts_spec,[n_energies, n_detectors, n_times]) + 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 - counts_err = counts_err * reform(reproduce(eff_ewidth, n_detectors*n_times),n_energies, n_detectors, n_times) + 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] ) - 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 - - ;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(counts_spec,2))[2,*] - shift((total(counts_spec,2))[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 - closest_jumps = value_closest(jumps_use, index) - index = jumps_use[closest_jumps] + ;; 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 + + if n_elements(pixels_used) gt 1 then begin + + spectrum = total(rebinned_count_rates[*,pixels_used,*], 2) + + endif else begin + + spectrum = reform(rebinned_count_rates[*,pixels_used,*]) + + endelse + + if n_elements(detectors_used) gt 1 then begin + + spectrum = total(spectrum[*,detectors_used], 2) + + endif else begin + + 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, 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 + + endfor + + endfor + + count_rates = count_rates_elut + count_rates_error = count_rates_error_elut 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 SPECTROGRAM - 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 + if n_elements(pixels_used) gt 1 then begin -end + 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 + + + ;; 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} + + + ;;------------------------------------------------------ + + 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..4d53507e 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,107 @@ ; 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, 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 +194,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 +211,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 +225,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 +238,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_elut_correction.pro b/stix/idl/processing/spectrogram/stx_elut_correction.pro new file mode 100644 index 00000000..7ff5b600 --- /dev/null +++ b/stix/idl/processing/spectrogram/stx_elut_correction.pro @@ -0,0 +1,274 @@ +;+ +; +; 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. +; - SP_INDEX: array containing the value of the spectral indices in the different energy bins +; +; 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: +; March 2026, Massa P., first release +; +; 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 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..b3705531 --- /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 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 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..4bd88a8e 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 +; February 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) }