Commit 55485070 authored by Médéric Boquien's avatar Médéric Boquien

More bits necessary to fit pseudo-magnitudes.

parent 9f46ac6d
......@@ -155,6 +155,47 @@ def _compute_scaling(models, obs, corr_dz, wz):
return num/denom
def _compute_scaling_mag(models, obs, corr_dz, wz):
"""Compute the scaling factor to be applied to the model fluxes to best fit
the observations. Note that we look over the bands to avoid the creation of
an array of the same size as the model_fluxes array. Because we loop on the
bands and not on the models, the impact on the performance should be small.
Parameters
----------
models: ModelsManagers class instance
Contains the models (fluxes, intensive, and extensive properties).
obs: Observation class instance
Contains the fluxes, intensive properties, extensive properties and
their errors, for a sigle observation.
corr_dz: float
Correction factor to scale the extensive properties to the right
distance
wz: slice
Selection of the models at the redshift of the observation or all the
redshifts in photometric-redshift mode.
Returns
-------
scaling: array
Scaling factors minimising the χ²
"""
_ = list(models.flux.keys())[0]
num = np.zeros_like(models.flux[_][wz])
denom = np.zeros_like(models.flux[_][wz])
for band, flux in obs.flux.items():
# Multiplications are faster than divisions, so we directly use the
# inverse error
inv_err2 = 1. / obs.flux_err[band] ** 2.
model = models.flux[band][wz]
num += (model - flux) * inv_err2
denom += inv_err2
return num/denom
def _correct_scaling_ul(scaling, mod, obs, wz):
"""Correct the scaling factor when one or more fluxes and/or properties are
upper limits.
......@@ -240,7 +281,7 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
scaling of the models to obtain the minimum χ²
"""
limits = lim_flag and (len(obs.flux_ul) > 0 or len(obs.extprop_ul) > 0)
scaling = _compute_scaling(models, obs, corr_dz, wz)
scaling = _compute_scaling_mag(models, obs, corr_dz, wz)
# Some observations may not have flux values in some filter(s), but
# they can have upper limit(s).
......@@ -256,7 +297,7 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
# inverse error
inv_flux_err = 1. / obs.flux_err[band]
model = models.flux[band][wz]
chi2 += ((flux - model * scaling) * inv_flux_err) ** 2.
chi2 += ((flux - model + scaling) * inv_flux_err) ** 2.
# Computation of the χ² from intensive properties
for name, prop in obs.intprop.items():
......
......@@ -185,17 +185,17 @@ def analysis(idx, obs):
else:
values = gbl_models.extprop[prop][wz]
_ = lambda x: x
mean, std = weighted_param(_(values[wlikely] * scaling_l * corr_dz),
mean, std = weighted_param(_(values[wlikely] * 10**(-.4*scaling_l) * corr_dz),
likelihood)
gbl_results.bayes.extmean[prop][idx] = mean
gbl_results.bayes.exterror[prop][idx] = std
if gbl_models.conf['analysis_params']['save_chi2'] is True:
save_chi2(obs, prop, gbl_models, chi2,
values * 10**(.4*scaling) * corr_dz)
values * 10**(-.4*scaling) * corr_dz)
for band in gbl_results.bayes.fluxmean:
values = gbl_models.flux[band][wz]
mean, std = weighted_param(values[wlikely] * scaling_l,
mean, std = weighted_param(values[wlikely] * 10**(-.4*scaling_l),
likelihood)
gbl_results.bayes.fluxmean[band][idx] = mean
gbl_results.bayes.fluxerror[band][idx] = std
......@@ -257,7 +257,7 @@ def bestfit(oidx, obs):
scaling = gbl_results.best.scaling[oidx] * corr_scaling
for band in gbl_results.best.flux:
gbl_results.best.flux[band][oidx] = sed.compute_fnu(band) * scaling
gbl_results.best.flux[band][oidx] = sed.compute_fnu(band) * 10**(-.4*scaling)
# If the distance is user defined, the redshift-based luminosity distance
# of the model is probably incorrect so we replace it
......@@ -266,7 +266,7 @@ def bestfit(oidx, obs):
gbl_results.best.intprop[prop][oidx] = sed.info[prop]
for prop in gbl_results.best.extprop:
gbl_results.best.extprop[prop][oidx] = sed.info[prop] * 10**(.4*scaling) \
gbl_results.best.extprop[prop][oidx] = sed.info[prop] * 10**(-.4*scaling) \
* corr_dz
if gbl_conf['analysis_params']["save_best_sed"]:
......
......@@ -35,7 +35,7 @@ class ObservationsManagerPassbands(object):
at each iteration.
"""
def __init__(self, config, params, defaulterror=0.1, modelerror=0.1,
def __init__(self, config, params, defaulterror=0.1, modelerror=0.,
threshold=-9990.):
self.conf = config
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment