From dc51142412dfc11255740cd0f8e1b6adffc14d74 Mon Sep 17 00:00:00 2001 From: Corentin Schreiber Date: Thu, 3 Aug 2023 09:09:41 +0100 Subject: [PATCH] Added SPEC_LSF_FILE --- CHANGELOG | 1 + README.md | 5 ++++ example/cosmos-20115/fast.param | 23 ++++++++++----- example/hdfn_fs99/fast.param | 23 ++++++++++----- src/fast++-gridder-custom.cpp | 6 ++-- src/fast++-gridder-ised.cpp | 4 +-- src/fast++-gridder.cpp | 52 ++++++++++++++++++++++++++++----- src/fast++-read_input.cpp | 49 +++++++++++++++++++++++++++++-- src/fast++.hpp | 13 +++++++-- 9 files changed, 145 insertions(+), 31 deletions(-) diff --git a/CHANGELOG b/CHANGELOG index 14c6f3b..54705aa 100644 --- a/CHANGELOG +++ b/CHANGELOG @@ -5,6 +5,7 @@ Features: - Added option TEMP_ERR_SPEC_FILE to control the template error function for the spectrum - The spectrum no longer inherits the same template error function as the photometry - AUTO_SCALE=1 no longer automatically disables the template error for the spectrum + - Added option SPEC_LSF_FILE to apply a line spread function to the models Bug fixes: - Fixed crash when saving extra outputs and using an absolute path in CATALOG diff --git a/README.md b/README.md index 30003d1..66ec19c 100644 --- a/README.md +++ b/README.md @@ -34,6 +34,7 @@ - [Using priors on the infrared luminosity](#using-priors-on-the-infrared-luminosity) - [Better treatment of spectra](#better-treatment-of-spectra) - [Velocity broadening](#velocity-broadening) + - [Line spread function](#line-spread-function) - [Continuum indices](#continuum-indices) - [Additional documentation](#additional-documentation) - [Adding new filters](#adding-new-filters) @@ -487,6 +488,10 @@ The solution implemented in FAST++ is to broaden all the spectral templates with NB: applying the velocity dispersion needs to be done for each SED library that is used in the fit. This incurs a small performance penalty at the beginning of the fit (and whenever the library is changed, e.g., when switching to a new metallicity). If your grid is small, this step can actually dominate the total computation time. +## Line spread function +```SPEC_LSF_FILE``` can be used to specify a file holding the line spread function of the instrument that was used to acquire the spectra. This file must contain the "sigma" of the Gaussian profile of the line spread function, expressed in Angstroms (observer frame), as a function of observed wavelength (also in Angstroms). It will be linearly interpolated at each observed wavelength of the models used in the fit, and the resulting value will be used to broaden the model by the specified value around that observed wavelength. + + ## Continuum indices Traditionally, absorption line strengths are measured as "equivalent widths" (in wavelength unit). These are part of a large set of continuum features call "spectral indices", see for example [Balogh et al. (1999)](http://adsabs.harvard.edu/abs/1999ApJ...527...54B), which includes other commonly used quantities such as the Dn4000 index. These features have been used extensively in the past to characterize the SEDs of galaxies and infer their star formation histories and/or metal abundances. diff --git a/example/cosmos-20115/fast.param b/example/cosmos-20115/fast.param index 18efa7e..43c59b3 100644 --- a/example/cosmos-20115/fast.param +++ b/example/cosmos-20115/fast.param @@ -106,9 +106,6 @@ MAX_QUEUED_FITS = 10000 # o TEMP_ERR_FILE: Template error function. The photometric errors are # in rest-frame multiplied by this function. # -# o TEMP_ERR_SPEC_FILE: Template error function for spectrum. The -# spectroscopic errors are in rest-frame multiplied by this function. -# # o NAME_ZPHOT: Header name of the column in your [CATALOG].zout file # that you want to use for your photometric redshifts. If not defined, # FAST will look for 'z_phot' @@ -148,7 +145,6 @@ AB_ZEROPOINT = 23.9 FILTERS_RES = '../../share/FILTER.RES.latest' FILTER_FORMAT = 1 TEMP_ERR_FILE = '../../share/TEMPLATE_ERROR.fast.v0.2' -TEMP_ERR_SPEC_FILE = '' NAME_ZPHOT = '' FORCE_ZPHOT = 0 # 0 / 1 BEST_AT_ZPHOT = 1 # 0 / 1 @@ -170,6 +166,9 @@ USE_LIR = 1 # 0 / 1 # - id[.]: should be an ID in the flux catalog # - Missing values can be signaled with negative errors or "NaN". # +# o TEMP_ERR_SPEC_FILE: Template error function for spectrum. The +# spectroscopic errors are in rest-frame multiplied by this function. +# # o AUTO_SCALE: This option automatically adjusts the spectrum to the # model separately from the rest of the photometry. The spectrum thus # does not participate in determining the best amplitude of the @@ -185,11 +184,21 @@ USE_LIR = 1 # 0 / 1 # when fitting spectra, and has an impact on performances (a couple # seconds of extra delay), so it is disabled by default. # +# o SPEC_LSF_FILE: Line spread function for spectrum. Must be given as the +# dispersion (in Angstrom) as a function of observed wavelength +# (in Angstrom). If specified, will broaden the templates with the +# provided dispersion. This is only relevant when fitting spectra, and +# has an large impact on performances. Only use if the spectral +# resolution of is finer than the LSF, and if binning the spectrum is +# not acceptable. +# #----------------------------------------------------------------------- -SPECTRUM = 'spectrum' -AUTO_SCALE = 1 # 0 / 1 -APPLY_VDISP = 350 # km/s +SPECTRUM = 'spectrum' +TEMP_ERR_SPEC_FILE = '' +AUTO_SCALE = 1 # 0 / 1 +APPLY_VDISP = 350 # km/s +SPEC_LSF_FILE = '' #--- OUTPUT INFORMATION ----------------------------------------------- diff --git a/example/hdfn_fs99/fast.param b/example/hdfn_fs99/fast.param index 8a4fad4..7258102 100644 --- a/example/hdfn_fs99/fast.param +++ b/example/hdfn_fs99/fast.param @@ -106,9 +106,6 @@ MAX_QUEUED_FITS = 1000 # o TEMP_ERR_FILE: Template error function. The photometric errors are # in rest-frame multiplied by this function. # -# o TEMP_ERR_SPEC_FILE: Template error function for spectrum. The -# spectroscopic errors are in rest-frame multiplied by this function. -# # o NAME_ZPHOT: Header name of the column in your [CATALOG].zout file # that you want to use for your photometric redshifts. If not defined, # FAST will look for 'z_phot' @@ -148,7 +145,6 @@ AB_ZEROPOINT = 25. FILTERS_RES = '../../share/FILTER.RES.latest' FILTER_FORMAT = 1 TEMP_ERR_FILE = '../../share/TEMPLATE_ERROR.fast.v0.2' -TEMP_ERR_SPEC_FILE = '' NAME_ZPHOT = 'z_m2' FORCE_ZPHOT = 0 # 0 / 1 BEST_AT_ZPHOT = 1 # 0 / 1 @@ -170,6 +166,9 @@ USE_LIR = 0 # 0 / 1 # - id[.]: should be an ID in the flux catalog # - Missing values can be signaled with negative errors or "NaN". # +# o TEMP_ERR_SPEC_FILE: Template error function for spectrum. The +# spectroscopic errors are in rest-frame multiplied by this function. +# # o AUTO_SCALE: This option automatically adjusts the spectrum to the # model separately from the rest of the photometry. The spectrum thus # does not participate in determining the best amplitude of the @@ -185,11 +184,21 @@ USE_LIR = 0 # 0 / 1 # when fitting spectra, and has an impact on performances (a couple # seconds of extra delay), so it is disabled by default. # +# o SPEC_LSF_FILE: Line spread function for spectrum. Must be given as the +# dispersion (in Angstrom) as a function of observed wavelength +# (in Angstrom). If specified, will broaden the templates with the +# provided dispersion. This is only relevant when fitting spectra, and +# has an large impact on performances. Only use if the spectral +# resolution of is finer than the LSF, and if binning the spectrum is +# not acceptable. +# #----------------------------------------------------------------------- -SPECTRUM = '' -AUTO_SCALE = 0 # 0 / 1 -APPLY_VDISP = 0 # km/s +SPECTRUM = '' +TEMP_ERR_SPEC_FILE = '' +AUTO_SCALE = 0 # 0 / 1 +APPLY_VDISP = 0 # km/s +SPEC_LSF_FILE = '' #--- OUTPUT INFORMATION ----------------------------------------------- diff --git a/src/fast++-gridder-custom.cpp b/src/fast++-gridder-custom.cpp index 4020251..40fe6cf 100644 --- a/src/fast++-gridder-custom.cpp +++ b/src/fast++-gridder-custom.cpp @@ -119,7 +119,7 @@ bool gridder_t::build_and_send_custom(fitter_t& fitter) { // Apply velocity dispersion if (is_finite(opts.apply_vdisp)) { - ssp.sed = convolve_vdisp(ssp.lambda, ssp.sed, opts.apply_vdisp); + ssp.sed = convolve_rest(ssp.lambda, ssp.sed); } // Pre-compute dust law & IGM absorption (they don't change with SFH) @@ -286,7 +286,7 @@ bool gridder_t::build_template_custom(uint_t iflat, vec1f& lam, vec1f& flux) con // Apply velocity dispersion if (is_finite(opts.apply_vdisp)) { - ssp->sed = convolve_vdisp(ssp->lambda, ssp->sed, opts.apply_vdisp); + ssp->sed = convolve_rest(ssp->lambda, ssp->sed); } } @@ -340,7 +340,7 @@ bool gridder_t::build_template_custom(uint_t iflat, vec1f& lam, vec1f& flux_youn // Apply velocity dispersion if (is_finite(opts.apply_vdisp)) { - ssp->sed = convolve_vdisp(ssp->lambda, ssp->sed, opts.apply_vdisp); + ssp->sed = convolve_rest(ssp->lambda, ssp->sed); } } diff --git a/src/fast++-gridder-ised.cpp b/src/fast++-gridder-ised.cpp index 4e01ad8..881a88c 100644 --- a/src/fast++-gridder-ised.cpp +++ b/src/fast++-gridder-ised.cpp @@ -308,7 +308,7 @@ bool gridder_t::build_and_send_ised(fitter_t& fitter) { // Apply velocity dispersion if (is_finite(opts.apply_vdisp)) { - ised.fluxes = convolve_vdisp(ised.lambda, ised.fluxes, opts.apply_vdisp); + ised.fluxes = convolve_rest(ised.lambda, ised.fluxes); } // Pre-compute dust law & IGM absorption (they don't change with SFH) @@ -434,7 +434,7 @@ bool gridder_t::build_template_ised(uint_t iflat, vec1f& lam, vec1f& flux) const // Apply velocity dispersion if (is_finite(opts.apply_vdisp)) { - ised->fluxes = convolve_vdisp(ised->lambda, ised->fluxes, opts.apply_vdisp); + ised->fluxes = convolve_rest(ised->lambda, ised->fluxes); } } } diff --git a/src/fast++-gridder.cpp b/src/fast++-gridder.cpp index 87ef838..6e606f4 100644 --- a/src/fast++-gridder.cpp +++ b/src/fast++-gridder.cpp @@ -702,6 +702,11 @@ void gridder_t::redshift_and_send(fitter_t& fitter, progress_t& pg, tpl_att_z_lam.safe[il] *= (1.0 + output_z.safe[iz]); } + if (!input.tpllsf_sigma.empty()) { + // Apply velocity dispersion + tpl_att_z_flux = convolve_obs(tpl_att_z_lam, tpl_att_z_flux); + } + // Integrate for (uint_t il : range(input.lambda)) { model.flux.safe[il] = astro::sed2flux( @@ -716,6 +721,8 @@ void gridder_t::redshift_and_send(fitter_t& fitter, progress_t& pg, } // See if we want to use this model or not + // NB we still want to compute the model fluxes for the cache, this is just to + // decide whether the model is used in the fit bool nofit = false; if (!opts.no_max_age && lage > auniv[iz]) { // Age greater than the age of the universe @@ -1014,6 +1021,11 @@ bool gridder_t::build_template_impl(uint_t iflat, bool nodust, lam.safe[il] *= (1.0 + z); } + if (!input.tpllsf_sigma.empty()) { + // Apply velocity dispersion + flux = convolve_obs(lam, flux); + } + // Integrate iflux.resize(input.lambda.size()); for (uint_t il : range(input.lambda)) { @@ -1046,8 +1058,10 @@ bool gridder_t::get_sfh(uint_t igrid, const vec1d& t, vec1d& sfh) const { } } -vec2d gridder_t::convolve_vdisp(const vec1d& lam, const vec2d& osed, double vdisp) const { - vec2d sed(osed.dims); +vec2d gridder_t::convolve_function(const vec1d& lam, const vec2d& osed, uint_t n_thread, + std::function sigma_fun) const { + + vec2d sed = replicate(0.0, osed.dims); const uint_t nlam = lam.size(); const double max_sigma = 5.0; @@ -1055,7 +1069,7 @@ vec2d gridder_t::convolve_vdisp(const vec1d& lam, const vec2d& osed, double vdis auto do_convolve = [&](uint_t lbegin, uint_t lend) { for (uint_t l : range(lbegin, lend)) { // Get sigma in wavelength units - double sigma = lam.safe[l]*(vdisp/2.99792e5); + double sigma = sigma_fun(lam.safe[l]); // Find bounds uint_t i0 = l, i1 = l; @@ -1077,16 +1091,16 @@ vec2d gridder_t::convolve_vdisp(const vec1d& lam, const vec2d& osed, double vdis } }; - if (opts.n_thread <= 1) { + if (n_thread <= 1) { // Single-threaded do_convolve(0, lam.size()); } else { // Multi-threaded - auto tp = thread::pool(opts.n_thread); - uint_t dl = lam.size()/opts.n_thread + 1; + auto tp = thread::pool(n_thread); + uint_t dl = lam.size()/n_thread + 1; uint_t l0 = 0; uint_t l1 = dl; - for (uint_t it : range(opts.n_thread)) { + for (uint_t it : range(n_thread)) { tp[it].start(do_convolve, l0, l1); l0 += dl; l1 += dl; @@ -1095,7 +1109,7 @@ vec2d gridder_t::convolve_vdisp(const vec1d& lam, const vec2d& osed, double vdis } } - for (uint_t it : range(opts.n_thread)) { + for (uint_t it : range(n_thread)) { tp[it].join(); } } @@ -1103,6 +1117,28 @@ vec2d gridder_t::convolve_vdisp(const vec1d& lam, const vec2d& osed, double vdis return sed; } +vec1d gridder_t::convolve_obs(const vec1d& lam, const vec1d& osed) const { + // We can't use multi-threading if 'generators' is chosen as parallel mode, since + // we will already be generating multiple models in different threads. + uint_t nthread = (opts.parallel == parallel_choice::generators ? 1 : opts.n_thread); + + return flatten(convolve_function(lam, reform(osed, 1, osed.size()), nthread, [&](double lobs) { + if (lobs < input.tpllsf_lam.front()) { + return input.tpllsf_sigma.front(); + } else if (lobs > input.tpllsf_lam.back()) { + return input.tpllsf_sigma.back(); + } else { + return interpolate(input.tpllsf_sigma, input.tpllsf_lam, lobs); + } + })); +} + +vec2d gridder_t::convolve_rest(const vec1d& lam, const vec2d& osed) const { + return convolve_function(lam, osed, opts.n_thread, [&](double lrest) { + return lrest*opts.apply_vdisp/2.99792e5; + }); +} + struct sed_id_pair { std::string id; uint_t igrid; diff --git a/src/fast++-read_input.cpp b/src/fast++-read_input.cpp index 2555be6..90bd43f 100644 --- a/src/fast++-read_input.cpp +++ b/src/fast++-read_input.cpp @@ -96,6 +96,7 @@ bool read_params(options_t& opts, input_state_t& state, const std::string& filen PARSE_OPTION_RENAME(filters_format, "filter_format") PARSE_OPTION(temp_err_file) PARSE_OPTION(temp_err_spec_file) + PARSE_OPTION(spec_lsf_file) PARSE_OPTION(name_zphot) PARSE_OPTION(spectrum) PARSE_OPTION(auto_scale) @@ -1703,7 +1704,7 @@ bool read_template_error(const options_t& opts, input_state_t& state) { } if (opts.verbose) { - note("apply template library error function to photometry"); + note("apply template library error function to photometry '", opts.temp_err_file, "'"); } if (!file::exists(opts.temp_err_file)) { @@ -1718,6 +1719,11 @@ bool read_template_error(const options_t& opts, input_state_t& state) { return false; } + if (!is_sorted(state.tplerr_lam)) { + error("template error function must be sorted by increasing wavelength"); + return false; + } + return true; } @@ -1727,7 +1733,7 @@ bool read_template_error_spec(const options_t& opts, input_state_t& state) { } if (opts.verbose) { - note("apply template library error function to spectrum"); + note("apply template library error function to spectrum '", opts.temp_err_spec_file, "'"); } if (!file::exists(opts.temp_err_spec_file)) { @@ -1742,6 +1748,40 @@ bool read_template_error_spec(const options_t& opts, input_state_t& state) { return false; } + if (!is_sorted(state.tplerr_spec_lam)) { + error("template spectrum error function must be sorted by increasing wavelength"); + return false; + } + + return true; +} + +bool read_lsf_spec(const options_t& opts, input_state_t& state) { + if (opts.spec_lsf_file.empty()) { + return true; + } + + if (opts.verbose) { + note("apply line spread function to models '", opts.spec_lsf_file, "'"); + } + + if (!file::exists(opts.spec_lsf_file)) { + error("could not open line spread function file '", opts.spec_lsf_file, "'"); + return false; + } + + ascii::read_table(opts.spec_lsf_file, state.tpllsf_lam, state.tpllsf_sigma); + + if (state.tpllsf_sigma.empty()) { + error("line spread function is empty, something must be wrong in the file"); + return false; + } + + if (!is_sorted(state.tpllsf_lam)) { + error("line spread function must be sorted by increasing wavelength"); + return false; + } + return true; } @@ -1946,6 +1986,11 @@ bool read_input(options_t& opts, input_state_t& state, const std::string& filena return false; } + // Read the spectrum line spread function, if any + if (!read_lsf_spec(opts, state)) { + return false; + } + // Read the continuum indices definitions, if any if (!read_continuum_indices(opts, state)) { return false; diff --git a/src/fast++.hpp b/src/fast++.hpp index f05bb42..b9ed026 100644 --- a/src/fast++.hpp +++ b/src/fast++.hpp @@ -49,7 +49,6 @@ struct options_t { // Templates parameter std::string temp_err_file; - std::string temp_err_spec_file; std::string library_dir; std::string library; std::string dust_law; @@ -99,6 +98,10 @@ struct options_t { // NB: parameters below were not in original FAST // ---------------------------------------------- + // Spectra + std::string temp_err_spec_file; + std::string spec_lsf_file; + // z-phot behavior bool force_zphot = false; bool best_at_zphot = false; @@ -229,6 +232,9 @@ struct input_state_t { vec1f tplerr_lam, tplerr_err; vec1f tplerr_spec_lam, tplerr_spec_err; + // Line spread function + vec1f tpllsf_lam, tpllsf_sigma; + // Continuum indices definitions vec<1,absorption_line_t> abs_lines; vec<1,continuum_ratio_t> cont_ratios; @@ -411,7 +417,10 @@ private : bool get_age_bounds(const vec1f& ised_age, float nage, std::array& p, double& x) const; void evaluate_sfh_custom(const vec1u& idm, const vec1d& t, vec1d& sfh) const; - vec2d convolve_vdisp(const vec1d& lam, const vec2d& osed, double vdisp) const; + vec2d convolve_function(const vec1d& lam, const vec2d& osed, uint_t n_thread, + std::function sigma_fun) const; + vec1d convolve_obs(const vec1d& lam, const vec1d& osed) const; + vec2d convolve_rest(const vec1d& lam, const vec2d& osed) const; mutable tinyexpr_wrapper sfh_expr; mutable tinyexpr_wrapper exclude_expr;