Merge branch 'automatic-fitting' into XSPEC-conversion

* automatic-fitting: Pass through additional useful keywords to the spectrum and srm FITS files Fix error in sci energy channel config file (i4Ds#96) Added automatic check for energy shift (i4Ds#101) Added a few functionalities to stx_science_data_lightcurve (i4Ds#99) stx_read_elut - fixed issue when filename keyword is not used # Conflicts: # stix/idl/processing/spectrogram/stx_convert_pixel_data.pro # stix/idl/processing/spectrogram/stx_convert_science_data2ospex.pro # stix/idl/processing/spectrogram/stx_convert_spectrogram.pro # stix/idl/processing/spectrogram/stx_fsw_sd_spectrogram2ospex.pro # stix/idl/processing/spectrogram/stx_write_ospex_fits.pro
grazwegian · Aug 8, 2022 · 044f7eb · 044f7eb
2 parents 04e9bcd + 370954d
commit 044f7eb
Show file tree

Hide file tree

Showing 9 changed files with 224 additions and 126 deletions.
diff --git a/stix/dbase/detector/ScienceEnergyChannels_1000.csv b/stix/dbase/detector/ScienceEnergyChannels_1000.csv
@@ -35,7 +35,7 @@ Channel Number,Channel Edge,Energy Edge ,Elower,Eupper ,BinWidth,dE/E,QL channel
 10,10,13,13,14,1,0.074,1
 11,11,14,14,15,1,0.069,1
 12,12,15,15,16,1,0.065,2
-13,13,16,16,17,1,0.061,2
+13,13,16,16,18,2,0.061,2
 14,14,18,18,20,2,0.105,2
 15,15,20,20,22,2,0.095,2
 16,16,22,22,25,3,0.128,2

diff --git a/stix/idl/processing/drm/stx_read_elut.pro b/stix/idl/processing/drm/stx_read_elut.pro
@@ -53,25 +53,32 @@
 ;    EKEV_ACTUAL   - channel edges with true energies based on full adc 4096 bins
 ;    SCALE1024     - 1 or 0, default 1.  If set, gain and offset are returned
 ;    for the 1024 ADC calibration data. EKEV_ACTUAL is the same regardless
+;    path          - Path to the folder in which to search for elut_table csv files, by default .../ssw/so/stix/dbase/detector is used. 
 ;
 ;
 ;
 ;
 ; :Author: rschwartz70@gmail.com, 2-jul-2019
 ; :History: 29-aug-2019, improved the file search for the elut file
 ; :RAS, 11-jun-2020, added ekev_actual, corrected the action of scale1024==0
-; :ECMD, 23-Feb-2022, fixed issue with finding file if current folder is not /stix 
+; :ECMD, 23-Feb-2022, fixed issue with finding file if current folder is not /stix
+; :ECMD, 03-Aug-2022, fixed issue when filename keyword is not used  
+;
 ;-
 pro stx_read_elut, gain, offset, adc4096_str, elut_filename = elut_filename, scale1024 = scale1024, $
-  ekev_actual = ekev_actual
+  ekev_actual = ekev_actual, path = path
 
   default, scale1024, 1
 
-  case 1 of
-    file_exist( elut_filename ) : elut_file = elut_filename
-    file_exist( concat_dir( concat_dir('SSW_STIX','dbase'),'detector') + get_delim() + elut_filename) : $
-      elut_file =  concat_dir( concat_dir('SSW_STIX','dbase'),'detector') + get_delim() + elut_filename
-  else : begin
+  if getenv('SSW_STIX') eq '' and ~keyword_set(path) then message, 'Path to ELUT table files not found.'
+
+  if keyword_set(elut_filename) then begin
+    case 1 of
+      file_exist( elut_filename ) : elut_file = elut_filename
+      file_exist( concat_dir( concat_dir('SSW_STIX','dbase'),'detector') + get_delim() + elut_filename) : $
+        elut_file =  concat_dir( concat_dir('SSW_STIX','dbase'),'detector') + get_delim() + elut_filename
+    end
+  endif else begin
 
     default, path, [concat_dir( concat_dir('SSW_STIX','dbase'),'detector')]
 
@@ -100,9 +107,8 @@ pro stx_read_elut, gain, offset, adc4096_str, elut_filename = elut_filename, sca
       ord = reverse( sort( anytim( date )))
       elut_file = elut_file[ ord[0] ]
     endif
-    end
-  endcase
-
+  end
+
   elut_str = reform_struct( read_csv( elut_file, n_table_header=3 ))
   ;Get the previous offset and gain suitable for the ADC1024 cal spectra
 

diff --git a/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro b/stix/idl/processing/spectrogram/stx_convert_pixel_data.pro
@@ -48,8 +48,8 @@
 ;                     N.B. WILL ONLY WORK WITH FULL TIME RESOLUTION DATA WHICH IS OFTEN NOT THE CASE FOR PIXEL DATA.
 ;
 ;     plot : in, type="boolean", default="1"
-;                     If set open OSPEX GUI and plot lightcurve in standard quicklook energy bands where there is data present 
-;                                  
+;                     If set open OSPEX GUI and plot lightcurve in standard quicklook energy bands where there is data present
+;
 ;      ospex_obj : out, type="OSPEX object"
 ;
 ;
@@ -62,12 +62,16 @@
 ;    18-Jun-2021 - ECMD (Graz), initial release
 ;    19-Jan-2022 - Andrea (FHNW), added keywords for selecting a subset of pixels and detectors for OSPEX
 ;    22-Feb-2022 - ECMD (Graz), documented, added default warnings, elut is determined by stx_date2elut_file, improved error calculation
-;    04-Jul-2022 - ECMD (Graz), added plot keyword 
+;    04-Jul-2022 - ECMD (Graz), added plot keyword
 ;    20-Jul-2022 - ECMD (Graz), distance factor now calculated in stx_convert_science_data2ospex
 ;
 ;-
-pro  stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fits_path_bk, time_shift = time_shift, energy_shift = energy_shift, distance = distance, $
-  flare_location= flare_location, xspec = xspec, plot = plot, ospex_obj = ospex_obj, det_ind = det_ind, pix_ind = pix_ind, shift_duration = shift_duration, no_attenuation=no_attenuation
+pro  stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fits_path_bk, $
+  time_shift = time_shift, energy_shift = energy_shift, distance = distance, flare_location= flare_location, $
+  det_ind = det_ind, pix_ind = pix_ind, $
+  shift_duration = shift_duration, no_attenuation=no_attenuation, sys_uncert = sys_uncert, $
+  generate_fits = generate_fits, specfile = specfile, srmfile = srmfile, xspec = xspec$
+  plot = plot, ospex_obj = ospex_obj
 
 
   if n_elements(time_shift) eq 0 then begin
@@ -76,14 +80,14 @@ pro  stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit
     print, 'using stx_get_header_corrections.pro.'
     time_shift = 0.
   endif
-
 
   default, flare_location, [0.,0.]
   default, shift_duration, 0
   default, xspec, 0
   default, plot, 1
 
 
+
   g10=[3,20,22]-1
   g09=[16,14,32]-1
   g08=[21,26,4]-1
@@ -113,6 +117,11 @@ pro  stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit
   start_time = atime(stx_time2any((t_axis.time_start)[0]))
 
   elut_filename = stx_date2elut_file(start_time)
+  uid = control_str.request_id
+
+  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, $
+    generate_fits = generate_fits, specfile = specfile, srmfile = srmfile, elut_file = elut_filename)
 
   counts_in = data_str.counts
 
@@ -214,9 +223,9 @@ pro  stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit
   ; ******************** 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 
+  ; 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.
@@ -252,7 +261,7 @@ pro  stx_convert_pixel_data, fits_path_data = fits_path_data, fits_path_bk = fit
   specpar = { sp_atten_state :  {time:ut_rcr[index], state:state} }
 
   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,$
-    distance = distance, fits_path_data = fits_path_data, flare_location= flare_location, eff_ewidth = eff_ewidth, xspec = xspec, plot = plot, ospex_obj = ospex_obj
+    distance = distance, fits_path_data = fits_path_data, flare_location = flare_location, eff_ewidth = eff_ewidth, xspec = xspec, sys_uncert = sys_uncert, plot = plot, fits_info_params = fits_info_params, ospex_obj = ospex_obj
 
 end
 
diff --git a/stix/idl/processing/spectrogram/stx_convert_science_data2ospex.pro b/stix/idl/processing/spectrogram/stx_convert_science_data2ospex.pro
@@ -14,21 +14,21 @@
 ;
 ;
 ; :keywords:
-; 
+;
 ;    spectrogram : in, type="stx_fsw_sd_spectrogram", default="1.0"
-;               the spectrogram structure containing the data 
-;               
+;               the spectrogram structure containing the data
+;
 ;    specpar : in, type="float", default="1.0"
-;               spectrum control and information parameters for OSPEX 
+;               spectrum control and information parameters for OSPEX
 ;
 ;    time_shift : in, type="float", default="1.0"
-;               Applied light travel time correction 
+;               Applied light travel time correction
 ;
 ;    data_level : in, type="int", default="1.0"
-;               The science data x-ray compaction level used to create the spectrogram 
+;               The science data x-ray compaction level used to create the spectrogram
 ;
 ;    data_dims : in, type="float", default="1.0"
-;               The dimensions of the spectrogram count data 
+;               The dimensions of the spectrogram count data
 ;
 ;    fits_path_bk : in, optional, type="string"
 ;              The path to file containing the background observation this should be in pixel data format i.e. sci-xray-cpd (or sci-xray-l1)
@@ -43,10 +43,10 @@
 ;               an output float value
 ;
 ;    plot : in, type="boolean", default="1"
-;                     If set open OSPEX GUI and plot lightcurve in standard quicklook energy bands where there is data present 
-;                                  
+;                     If set open OSPEX GUI and plot lightcurve in standard quicklook energy bands where there is data present
+;
 ;      ospex_obj : out, type="OSPEX object",
-;               the output OSPEX object containing the data 
+;               the output OSPEX object containing the data
 ;
 ;
 ; :history:
@@ -57,24 +57,25 @@
 ;
 ;-
 pro 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,$
-  distance = distance, fits_path_data = fits_path_data, flare_location= flare_location, eff_ewidth = eff_ewidth, xspec = xspec, plot = plot, ospex_obj = ospex_obj
+  distance = distance, fits_path_data = fits_path_data, xspec = xspec,fits_info_params = fits_info_params,  flare_location = flare_location, eff_ewidth = eff_ewidth, sys_uncert = sys_uncert,  plot = plot, generate_fits = fits,ospex_obj = ospex_obj
 
   default, plot, 0
 
   ;if distance is not set use the average value from the fits header
   stx_get_header_corrections, fits_path_data, distance = header_distance
   default, distance, header_distance
-  print, 'Using Solar Orbiter distance of :' + distance +  ' AU'
+  print, 'Using Solar Orbiter distance of :' + strtrim(distance,2) +  ' AU'
 
   dist_factor = 1./(distance^2.)
 
+
   n_energies = data_dims[0]
   n_detectors = data_dims[1]
   n_pixels = data_dims[2]
   n_times = data_dims[3]
 
   counts_spec  = spectrogram.counts
-  livetime_frac =  stx_spectrogram_livetime(  spectrogram, corrected_counts = corrected_counts, corrected_error = corrected_error, level = data_level )
+  livetime_frac =  stx_spectrogram_livetime( spectrogram, corrected_counts = corrected_counts, corrected_error = corrected_error, level = data_level )
 
   corrected_counts = total(reform(corrected_counts, [n_energies, n_detectors, n_times ]),2)
 
@@ -154,23 +155,24 @@ pro stx_convert_science_data2ospex, spectrogram = spectrogram, specpar = specpar
     error_bk =  reform(error_bk,[n_elements(energy_bins), n_times])
 
 
-    spec_in_corr = corrected_counts - corrected_counts_bk 
+    spec_in_corr = corrected_counts - corrected_counts_bk
 
-    spec_in_uncorr = counts_spec - spec_in_bk 
+    spec_in_uncorr = counts_spec - spec_in_bk
 
     total_error = sqrt(corrected_error^2. + error_bk^2. )
 
 
   endif else begin
 
-    spec_in_corr = corrected_counts 
+    spec_in_corr = corrected_counts
 
-    spec_in_uncorr = counts_spec 
+    spec_in_uncorr = counts_spec
 
     total_error = corrected_error
 
   endelse
 
+
   eff_livetime_fraction = f_div(total(counts_spec,1) , total(corrected_counts,1) , default = 1 )
   eff_livetime_fraction_expanded = transpose(rebin([eff_livetime_fraction],n_elements(eff_livetime_fraction),n_energies))
   spec_in_corr *= eff_livetime_fraction_expanded
@@ -201,30 +203,41 @@ pro stx_convert_science_data2ospex, spectrogram = spectrogram, specpar = specpar
     rcr           : spectrogram.rcr,$
     error         : total_error}
 
-  fstart_time = time2fid(atime(stx_time2any((spectrogram.time_axis.time_start)[0])),/full,/time)
+  uid = fits_info_params.uid
 
-  srmfilename = 'stx_spectrum_srm_' + fstart_time + '.fits'
-  specfilename =  'stx_spectrum_'   + fstart_time + '.fits'
+;  fstart_time = time2fid(atime(stx_time2any((spectrogram.time_axis.time_start)[0])),/full,/time)
+
+  default, specfilename, 'stx_spectrum_' + strtrim(uid,2) + '.fits'
+  default, srmfilename,  'stx_srm_'      + strtrim(uid,2) + '.fits'
+
+;  outfile = dialog_pickfile(file = specfilename, path=curdir(), filter='*.fits', $
+;    title = 'Select output file name')
+
+   fits_info_params.distance = distance
+   fits_info_params.specfile = specfilename
+   fits_info_params.srmfile = srmfilename
 
   transmission = read_csv(loc_file( 'stix_trans_by_component.csv', path = getenv('STX_GRID')))
 
   phe = transmission.field9
   phe = phe[where(phe gt emin-1 and phe lt 2*emax)]
   edge_products, phe, mean = mean_phe, width = w_phe
-  ph_in = [mean_phe[0]- w_phe[0], mean_phe]
+  ph_in = [mean_phe[0] - w_phe[0], mean_phe]
+
+  ospex_obj =   stx_fsw_sd_spectrogram2ospex( spectrogram, specpar = specpar, time_shift= time_shift, ph_energy_edges = ph_in, $
+    /include_damage, generate_fits = generate_fits, xspec = xspec , /tail, livetime_fraction = eff_livetime_fraction, $
+    dist_factor = dist_factor, flare_location= flare_location, sys_uncert = sys_uncert, fits_info_params = fits_info_params )
 
-  ospex_obj =   stx_fsw_sd_spectrogram2ospex( spectrogram, specpar = specpar, time_shift= time_shift, ph_energy_edges = ph_in, /include_damage, /fits , /tail, livetime_fraction = eff_livetime_fraction, $
-    dist_factor = dist_factor, flare_location= flare_location, xspec = xspec )
-
   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
+  ospex_obj->set, spex_error_use_expected = 0
 
   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
+    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