Add the possibility the indicate the distance of the object in Mpc in the input file.

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ## Unreleased
 ### Added
+- It is now possible to optionally indicate the distance in Mpc in the input file. If present it will be used in lieu of the distance computed from the redshift. This is especially useful in the nearby universe where the redshift is a very poor indicator of the actual distance. (Médéric Boquien)
+
 ### Changed
 ### Fixed
 ### Optimised

pcigale/analysis_modules/pdf_analysis/utils.py

@@ -11,6 +11,7 @@ from astropy import log
 from astropy.cosmology import WMAP7 as cosmo
 import numpy as np
 from scipy import optimize
+from scipy.constants import parsec
 from scipy.special import erf
 
 log.setLevel('ERROR')
@@ -29,7 +30,7 @@ def save_chi2(obs, variable, models, chi2, values):
 
 
 @lru_cache(maxsize=None)
-def compute_corr_dz(model_z, obs_z):
+def compute_corr_dz(model_z, obs_dist):
     """The mass-dependent physical properties are computed assuming the
     redshift of the model. However because we round the observed redshifts to
     two decimals, there can be a difference of 0.005 in redshift between the
@@ -42,16 +43,13 @@ def compute_corr_dz(model_z, obs_z):
     ----------
     model_z: float
         Redshift of the model.
-    obs_z: float
-        Redshift of the observed object.
+    obs_dist: float
+        Luminosity distance of the observed object.
 
     """
-    if model_z == obs_z:
-        return 1.
-    if model_z > 0.:
-        return (cosmo.luminosity_distance(obs_z).value /
-                cosmo.luminosity_distance(model_z).value)**2.
-    return (cosmo.luminosity_distance(obs_z).value * 1e5)**2.
+    if model_z == 0.:
+        return (obs_dist / (10. * parsec))**2.
+    return (obs_dist / cosmo.luminosity_distance(model_z).value)**2.
 
 
 def dchi2_over_ds2(s, obs_values, obs_errors, mod_values):

pcigale/analysis_modules/pdf_analysis/workers.py

@@ -173,7 +173,7 @@ def analysis(idx, obs):
         z = np.array(
             gbl_models.conf['sed_modules_params']['redshifting']['redshift'])
         wz = slice(np.abs(obs.redshift-z).argmin(), None, z.size)
-        corr_dz = compute_corr_dz(z[wz.start], obs.redshift)
+        corr_dz = compute_corr_dz(z[wz.start], obs.distance)
     else:  # We do not know the redshift so we use the full grid
         wz = slice(0, None, 1)
         corr_dz = 1.
@@ -263,12 +263,13 @@ def bestfit(oidx, obs):
     _, scaling = compute_chi2(fluxes[:, None], intprops[:, None],
                               extprops[:, None],  obs,
                               gbl_conf['analysis_params']['lim_flag'])
+    corr_dz = compute_corr_dz(obs.redshift, obs.distance)
 
     gbl_results.best.properties[oidx, :] = [sed.info[k] for k in
                                             gbl_results.best.propertiesnames]
     iprop = [i for i, k in enumerate(gbl_results.best.propertiesnames)
              if k in gbl_results.best.massproportional]
-    gbl_results.best.properties[oidx, iprop] *= scaling
+    gbl_results.best.properties[oidx, iprop] *= scaling * corr_dz
     gbl_results.best.fluxes[oidx, :] = fluxes * scaling
 
     if gbl_conf['analysis_params']["save_best_sed"]:

pcigale/managers/observations.py

@@ -3,8 +3,10 @@
 # Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
 # Author: Médéric Boquien
 
+from astropy.cosmology import WMAP7 as cosmo
 from astropy.table import Column
 import numpy as np
+from scipy.constants import parsec
 
 from ..utils import read_table
 from .utils import get_info
@@ -33,6 +35,15 @@ class Observation(object):
     def __init__(self, row, cls):
         self.redshift = row['redshift']
         self.id = row['id']
+        if 'distance' in row.colnames and np.isfinite(row['distance']):
+            self.distance = row['distance'] * parsec * 1e6
+        else:
+            if self.redshift == 0.:
+                self.distance = 10. * parsec
+            elif self.redshift > 0:
+                self.distance = cosmo.luminosity_distance(self.redshift).value
+            else:
+                self.distance = np.nan
         self.fluxes = np.array([row[band] for band in cls.bands])
         self.fluxes_err = np.array([row[band + '_err'] for band in cls.bands])
         self.intprops = np.array([row[prop] for prop in cls.intprops])
@@ -82,6 +93,10 @@ class ObservationsManagerPassbands(object):
                             threshold)
         self._add_model_error(modelerror)
 
+        # Rebuild the quantities to fit after vetting them
+        self.tofit = self.bands + self.intprops + self.extprops
+        self.tofit_err = self.bands_err + self.intprops_err + self.extprops_err
+
         self.observations = list([Observation(row, self) for row in self.table])
 
     def __len__(self):
@@ -116,8 +131,8 @@ class ObservationsManagerPassbands(object):
                                 "in the observation table.".format(item))
 
         for item in self.table.colnames:
-            if (item != 'id' and item != 'redshift' and item not in self.tofit +
-                self.tofit_err):
+            if (item != 'id' and item != 'redshift' and item != 'distance' and
+                item not in self.tofit + self.tofit_err):
                 self.table.remove_column(item)
                 print("Warning: {} in the input file but not to be taken into"
                       " account in the fit.".format(item))

pcigale/session/configuration.py

@@ -68,10 +68,12 @@ class Configuration(object):
         self.config.comments['data_file'] = wrap(
             "File containing the input data. The columns are 'id' (name of the"
             " object), 'redshift' (if 0 the distance is assumed to be 10 pc), "
-            "the filter names for the fluxes, and the filter names with the "
-            "'_err' suffix for the uncertainties. The fluxes and the "
-            "uncertainties must be in mJy. This file is optional to generate "
-            "the configuration file, in particular for the savefluxes module.")
+            "'distance' (Mpc, optional, if present it will be used in lieu "
+            "of the distance computed from the redshift), the filter names for"
+            " the fluxes, and the filter names with the '_err' suffix for the "
+            "uncertainties. The fluxes and the uncertainties must be in mJy. "
+            "This file is optional to generate the configuration file, in "
+            "particular for the savefluxes module.")
         self.spec['data_file'] = "string()"
 
         self.config['parameters_file'] = ""