Added SPEC_LSF_FILE

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

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.
 

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  -----------------------------------------------

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  -----------------------------------------------

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);
         }
     }
 

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);
             }
         }
     }

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,16 +1058,18 @@ 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<double(double)> sigma_fun) const {
+
+    vec2d sed = replicate(0.0, osed.dims);
 
     const uint_t nlam = lam.size();
     const double max_sigma = 5.0;
 
     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,14 +1109,36 @@ 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();
         }
     }
 
     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;

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;