Added TEMP_ERR_SPEC_FILE

CHANGELOG

@@ -1,6 +1,11 @@
-v1.3.3 (14/07/2023)
+v1.4.0 (02/08/2023)
 ======
 
+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
+
 Bug fixes:
  - Fixed crash when saving extra outputs and using an absolute path in CATALOG
 

example/cosmos-20115/fast.param

@@ -106,6 +106,9 @@ 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'
@@ -140,16 +143,17 @@ MAX_QUEUED_FITS = 10000
 #
 #-----------------------------------------------------------------------
 
-CATALOG        = 'photometry'
-AB_ZEROPOINT   = 23.9
-FILTERS_RES    = '../../share/FILTER.RES.latest'
-FILTER_FORMAT  = 1
-TEMP_ERR_FILE  = '../../share/TEMPLATE_ERROR.fast.v0.2'
-NAME_ZPHOT     = ''
-FORCE_ZPHOT    = 0          # 0 / 1
-BEST_AT_ZPHOT  = 1          # 0 / 1
-ZPHOT_CONF     = 68         # 68 / 95 / 99
-USE_LIR        = 1          # 0 / 1
+CATALOG             = 'photometry'
+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
+ZPHOT_CONF          = 68         # 68 / 95 / 99
+USE_LIR             = 1          # 0 / 1
 
 
 #--- SPECTROSCOPIC INFORMATION -----------------------------------------

example/hdfn_fs99/fast.param

@@ -106,6 +106,9 @@ 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'
@@ -140,16 +143,17 @@ MAX_QUEUED_FITS = 1000
 #
 #-----------------------------------------------------------------------
 
-CATALOG        = 'hdfn_fs99'
-AB_ZEROPOINT   = 25.
-FILTERS_RES    = '../../share/FILTER.RES.latest'
-FILTER_FORMAT  = 1
-TEMP_ERR_FILE  = '../../share/TEMPLATE_ERROR.fast.v0.2'
-NAME_ZPHOT     = 'z_m2'
-FORCE_ZPHOT    = 0          # 0 / 1
-BEST_AT_ZPHOT  = 1          # 0 / 1
-ZPHOT_CONF     = 68         # 68 / 95 / 99
-USE_LIR        = 0          # 0 / 1
+CATALOG             = 'hdfn_fs99'
+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
+ZPHOT_CONF          = 68         # 68 / 95 / 99
+USE_LIR             = 0          # 0 / 1
 
 
 #--- SPECTROSCOPIC INFORMATION -----------------------------------------

src/fast++-fitter.cpp

@@ -7,27 +7,38 @@ fitter_t::fitter_t(const options_t& opt, const input_state_t& inp, const gridder
     vec1f& output_z = output.grid[grid_id::z];
 
     // Pre-compute template error
-    if (!opts.temp_err_file.empty()) {
+    if (!opts.temp_err_file.empty() || !opts.temp_err_spec_file.empty()) {
         if (opts.verbose) note("initializing template error function...");
-        double l0 = median(input.tplerr_lam);
+        tpl_err = replicate(0.0f, output_z.size(), input.lambda.size());
+    }
+
+    auto precompute_tpl_err = [&](uint_t is, uint_t nlambda, const vec1f& lam, const vec1f& err) {
+        double l0 = median(lam);
 
-        tpl_err.resize(output_z.size(), input.lambda.size());
         for (uint_t iz : range(output_z))
-        for (uint_t il : range(input.lambda)) {
-            tpl_err.safe(iz,il) = sqr(astro::sed2flux(
+        for (uint_t il : range(is, nlambda)) {
+            tpl_err.safe(iz,il) = astro::sed2flux(
                 input.filters[il].wl, input.filters[il].tr,
-                input.tplerr_lam*(1.0 + output_z[iz]), input.tplerr_err
-            ));
+                lam*(1.0 + output_z[iz]), err
+            );
 
             if (!is_finite(tpl_err.safe(iz,il))) {
                 // The filter goes out of the template error function, extrapolate
                 if (input.lambda.safe[il] < l0) {
-                    tpl_err.safe(iz,il) = sqr(input.tplerr_err.front());
+                    tpl_err.safe(iz,il) = err.front();
                 } else {
-                    tpl_err.safe(iz,il) = sqr(input.tplerr_err.back());
+                    tpl_err.safe(iz,il) = err.back();
                 }
             }
         }
+    };
+
+    if (!opts.temp_err_file.empty()) {
+        precompute_tpl_err(input.phot_start, input.phot_end, input.tplerr_lam, input.tplerr_err);
+    }
+
+    if (!opts.temp_err_spec_file.empty()) {
+        precompute_tpl_err(input.spec_start, input.spec_end, input.tplerr_spec_lam, input.tplerr_spec_err);
     }
 
     // Pre-identify zgrid for galaxies with zspec or zphot
@@ -522,18 +533,24 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
             ++iflx;
         }
 
-        for (uint_t il : range(iflx, nscale)) {
-            wsp.weight[il] = (opts.temp_err_file.empty() ?
-                1.0/input.eflux.safe(is,il-iflx) :
-                1.0/sqrt((sqr(input.eflux.safe(is,il-iflx)) +
-                    tpl_err.safe(iz,il-iflx)*sqr(input.flux.safe(is,il-iflx))))
+        bool use_template_error = !tpl_err.empty();
+
+        auto fill_workspace = [&](uint_t iw, uint_t il) {
+            wsp.weight[iw] = (!use_template_error ?
+                1.0/input.eflux.safe(is,il) :
+                1.0/sqrt((sqr(input.eflux.safe(is,il)) +
+                    sqr(tpl_err.safe(iz,il)*input.flux.safe(is,il))))
             );
 
-            wsp.wflux[il] = input.flux.safe(is,il-iflx)*wsp.weight[il];
-            wsp.wmodel[il] = model.flux.safe[il-iflx]*wsp.weight[il];
+            wsp.wflux[iw] = input.flux.safe(is,il)*wsp.weight[iw];
+            wsp.wmodel[iw] = model.flux.safe[il]*wsp.weight[iw];
 
-            wfm += wsp.wmodel[il]*wsp.wflux[il];
-            wmm += sqr(wsp.wmodel[il]);
+        };
+
+        for (uint_t iw : range(iflx, nscale)) {
+            fill_workspace(iw, iw-iflx);
+            wfm += wsp.wmodel[iw]*wsp.wflux[iw];
+            wmm += sqr(wsp.wmodel[iw]);
         }
 
         double scale = wfm/wmm;
@@ -546,13 +563,10 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
         // the shape of the spectrum.
         double swfm = 0, swmm = 0;
         if (spec_auto_scale) {
-            for (uint_t il : range(nscale, ndata)) {
-                wsp.weight[il] = 1.0/input.eflux.safe(is,il-iflx);
-                wsp.wflux[il] = input.flux.safe(is,il-iflx)*wsp.weight[il];
-                wsp.wmodel[il] = model.flux.safe[il-iflx]*wsp.weight[il];
-
-                swfm += wsp.wmodel[il]*wsp.wflux[il];
-                swmm += sqr(wsp.wmodel[il]);
+            for (uint_t iw : range(nscale, ndata)) {
+                fill_workspace(iw, iw-iflx);
+                swfm += wsp.wmodel[iw]*wsp.wflux[iw];
+                swmm += sqr(wsp.wmodel[iw]);
             }
 
             spec_scale = swfm/swmm;
@@ -565,11 +579,11 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
 
         // Compute chi2
         double tchi2 = 0;
-        for (uint_t il : range(nscale)) {
-            tchi2 += sqr(wsp.wflux[il] - scale*wsp.wmodel[il]);
+        for (uint_t iw : range(nscale)) {
+            tchi2 += sqr(wsp.wflux[iw] - scale*wsp.wmodel[iw]);
         }
-        for (uint_t il : range(nscale, ndata)) {
-            tchi2 += sqr(wsp.wflux[il] - spec_scale*wsp.wmodel[il]);
+        for (uint_t iw : range(nscale, ndata)) {
+            tchi2 += sqr(wsp.wflux[iw] - spec_scale*wsp.wmodel[iw]);
         }
 
         // Add LIR as a contribution to chi2 (if given in log units)
@@ -621,15 +635,15 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
                 // all models (and each galaxy) will use the same random numbers
                 wfm = 0;
                 swfm = 0;
-                for (uint_t il : range(nscale)) {
+                for (uint_t iw : range(nscale)) {
                     // In weighted units, the random perturbations have a sigma of unity
-                    wsp.rflux[il] = wsp.wflux[il] + sim_rnd.safe(im,il);
-                    wfm += wsp.wmodel[il]*wsp.rflux[il];
+                    wsp.rflux[iw] = wsp.wflux[iw] + sim_rnd.safe(im,iw);
+                    wfm += wsp.wmodel[iw]*wsp.rflux[iw];
                 }
-                for (uint_t il : range(nscale, ndata)) {
+                for (uint_t iw : range(nscale, ndata)) {
                     // When auto scale is ON, spectral data doesn't participate in scale factor
-                    wsp.rflux[il] = wsp.wflux[il] + sim_rnd.safe(im,il);
-                    swfm += wsp.wmodel[il]*wsp.rflux[il];
+                    wsp.rflux[iw] = wsp.wflux[iw] + sim_rnd.safe(im,iw);
+                    swfm += wsp.wmodel[iw]*wsp.rflux[iw];
                 }
 
                 scale = wfm/wmm;
@@ -641,11 +655,11 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
 
                 // Compute chi2
                 tchi2 = 0;
-                for (uint_t il : range(nscale)) {
-                    tchi2 += sqr(wsp.rflux[il] - scale*wsp.wmodel[il]);
+                for (uint_t iw : range(nscale)) {
+                    tchi2 += sqr(wsp.rflux[iw] - scale*wsp.wmodel[iw]);
                 }
-                for (uint_t il : range(nscale, ndata)) {
-                    tchi2 += sqr(wsp.rflux[il] - spec_scale*wsp.wmodel[il]);
+                for (uint_t iw : range(nscale, ndata)) {
+                    tchi2 += sqr(wsp.rflux[iw] - spec_scale*wsp.wmodel[iw]);
                 }
 
                 // Add LIR as a contribution to chi2 (if given in log units)

src/fast++-read_input.cpp

@@ -95,6 +95,7 @@ bool read_params(options_t& opts, input_state_t& state, const std::string& filen
         PARSE_OPTION(filters_res)
         PARSE_OPTION_RENAME(filters_format, "filter_format")
         PARSE_OPTION(temp_err_file)
+        PARSE_OPTION(temp_err_spec_file)
         PARSE_OPTION(name_zphot)
         PARSE_OPTION(spectrum)
         PARSE_OPTION(auto_scale)
@@ -1062,6 +1063,9 @@ bool read_fluxes(const options_t& opts, input_state_t& state) {
         state.eflux(i,idbb) = 1e19*abzp*astro::uJy2cgs(state.lambda*1e-4, state.eflux(i,idbb));
     }
 
+    state.phot_start = 0;
+    state.phot_end = state.flux.dims[1];
+
     if (opts.verbose) {
         note("fitting ", state.flux.dims[0], " source", (state.flux.dims[0] > 1 ? "s" : ""),
             " with ", state.flux.dims[1], " fluxes", (state.flux.dims[0] > 1 ? " each" : ""));
@@ -1717,6 +1721,30 @@ bool read_template_error(const options_t& opts, input_state_t& state) {
     return true;
 }
 
+bool read_template_error_spec(const options_t& opts, input_state_t& state) {
+    if (opts.temp_err_spec_file.empty()) {
+        return true;
+    }
+
+    if (opts.verbose) {
+        note("apply template library error function to spectrum");
+    }
+
+    if (!file::exists(opts.temp_err_spec_file)) {
+        error("could not open template spectrum error function file '", opts.temp_err_spec_file, "'");
+        return false;
+    }
+
+    ascii::read_table(opts.temp_err_spec_file, state.tplerr_spec_lam, state.tplerr_spec_err);
+
+    if (state.tplerr_spec_lam.empty()) {
+        error("template spectrum error function is empty, something must be wrong in the file");
+        return false;
+    }
+
+    return true;
+}
+
 bool check_input(options_t& opts, input_state_t& state) {
     if (!state.zspec.empty()) {
         // Check that zspecs are covered by the redshift grid
@@ -1913,6 +1941,11 @@ bool read_input(options_t& opts, input_state_t& state, const std::string& filena
         return false;
     }
 
+    // Read the spectrum template error function, if any
+    if (!read_template_error_spec(opts, state)) {
+        return false;
+    }
+
     // Read the continuum indices definitions, if any
     if (!read_continuum_indices(opts, state)) {
         return false;

src/fast++-write_output.cpp

@@ -69,6 +69,9 @@ void write_catalog(const options_t& opts, const input_state_t& input, const grid
     if (!opts.temp_err_file.empty()) {
     fout << "# Template error function: " << opts.temp_err_file << std::endl;
     }
+    if (!opts.temp_err_spec_file.empty()) {
+    fout << "# Template spectrum error function: " << opts.temp_err_spec_file << std::endl;
+    }
     fout << "# AB ZP:       " << opts.ab_zeropoint << std::endl;
     fout << "# Library:     " << pretty_library(opts.library) << std::endl;
     switch (opts.sfh) {

src/fast++.hpp

@@ -49,6 +49,7 @@ 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;
@@ -197,6 +198,8 @@ struct sfh_quantity_t {
 struct input_state_t {
     // List of filter ID used in the photometric catalog
     vec1u no_filt;                      // [nfilt]
+    // First and last ID of photometric measurements in the flux array
+    uint_t phot_start = npos, phot_end = npos; // [nphot]
     // First and last ID of spectral measurements in the flux array
     uint_t spec_start = npos, spec_end = npos; // [nspec]
     // Central wavelength of the filters
@@ -222,8 +225,9 @@ struct input_state_t {
     vec<1,fast_filter_t> filters; // [nfilt+nspec]
     vec<1,fast_filter_t> rf_filters; // [nrffilt]
 
-    // Template error function
+    // Template error functions
     vec1f tplerr_lam, tplerr_err;
+    vec1f tplerr_spec_lam, tplerr_spec_err;
 
     // Continuum indices definitions
     vec<1,absorption_line_t> abs_lines;