Commit 5f54be4d authored by Guang's avatar Guang

Change alpha_ox to det_alpha_ox, to use the empirical alpha_ox-L_2500 relation

parent 0cd0d6c7
......@@ -254,6 +254,15 @@ def compute_chi2(models, obs, corr_dz, wz, lim_flag):
model = models.flux[band][wz]
chi2 += ((flux - model * scaling) * inv_flux_err) ** 2.
# Penalize det_alpha_ox out of the user-set range
if 'xray' in models.params.modules:
# Calculate expected alpha_ox from Lnu_2500 (Just et al. 2007)
exp_alpha_ox = -0.137*np.log10(models.extprop['agn.agn_intrin_Lnu_2500A'][wz]*1e7*scaling)+2.638
# Calculate det_alpha_ox = alpha_ox - alpha_ox(Lnu_2500)
det_alpha_ox = models.intprop['xray.alpha_ox'][wz] - exp_alpha_ox
# If det_alpha_ox out of range, set corresponding chi2 to nan
chi2[np.abs(det_alpha_ox) > models.intprop['xray.max_dev_alpha_ox'][wz]] = np.nan
# Computation of the χ² from intensive properties
for name, prop in obs.intprop.items():
model = models.intprop[name][wz]
......
......@@ -33,6 +33,12 @@ class ModelsManager(object):
props_nolog = set([prop[:-4] if prop.endswith('log') else prop
for prop in conf['analysis_params']['variables']])
# If X-ray module used, add some necessary properties for internal use
if 'xray' in self.params.modules:
if 'xray.alpha_ox' not in props_nolog: props_nolog.add('xray.alpha_ox')
if 'xray.max_dev_alpha_ox' not in props_nolog: props_nolog.add('xray.max_dev_alpha_ox')
if 'agn.agn_intrin_Lnu_2500A' not in props_nolog: props_nolog.add('agn.agn_intrin_Lnu_2500A')
self.intpropnames = (self.allintpropnames & set(obs.intprops) |
self.allintpropnames & props_nolog)
self.extpropnames = (self.allextpropnames & set(obs.extprops) |
......
......@@ -42,6 +42,11 @@ def complete_parameters(given_parameters, parameter_list):
if (key not in given_parameters) and (
parameter_list[key][2] is not None):
given_parameters[key] = parameter_list[key][2]
# alpha_ox is added internally to given_parameters
# Now, add it to parameter_list for consistency
if 'alpha_ox' in given_parameters:
parameter_list['alpha_ox'] = ('cigale_list()', 'power-law slope connecting L_nu at rest-frame 2500 A and 2 keV', \
given_parameters['alpha_ox'])
# Check parameter consistency between the parameter list and the given
# parameters.
if not set(given_parameters) == set(parameter_list):
......
......@@ -33,11 +33,15 @@ class Xray(SedModule):
"The photon index (Gamma) of intrinsic X-ray spectrum.",
1.6
)),
("alpha_ox", (
("max_dev_alpha_ox", (
"cigale_list()",
"Power-law slope connecting L_nu at rest-frame 2500 A and 2 keV. "
"Defined as alpha_ox=0.3838*log(Lnu_2keV/Lnu_2500A)",
-1.5
"Maximum deviation of alpha_ox from the empirical alpha_ox-Lnu_2500A relation (Just et al. 2007), "
"i.e. |alpha_ox-alpha_ox(Lnu_2500A)| <= max_dev_alpha_ox. "
"alpha_ox is the power-law slope connecting L_nu at rest-frame 2500 A and 2 keV, "
"defined as alpha_ox = 0.3838*log(Lnu_2keV/Lnu_2500A), "
"which is often modeled as a function of Lnu_2500A."
"The alpha_ox-Lnu_2500A relation has a 1-sigma scatter of ~ 0.1.",
0.3
))
])
......@@ -45,7 +49,7 @@ class Xray(SedModule):
"""Build the model for a given set of parameters."""
self.gam = float(self.parameters["gam"])
self.alpha_ox = float(self.parameters["alpha_ox"])
self.max_dev_alpha_ox = float(self.parameters["max_dev_alpha_ox"])
# We define various constants necessary to compute the model. For
# consistency, we define speed of light in units of nm s¯¹
......@@ -102,6 +106,8 @@ class Xray(SedModule):
logT = np.log10( sed.info['stellar.age_m_star']/1e3 )
# metallicity, units: none
Z = sed.info['stellar.metallicity']
# alpha_ox, added internally
alpha_ox = self.parameters['alpha_ox']
# AGN 2500A intrinsic luminosity
if 'agn.agn_intrin_Lnu_2500A' not in sed.info:
sed.add_info('agn.agn_intrin_Lnu_2500A', 1., True)
......@@ -110,7 +116,8 @@ class Xray(SedModule):
# Add the configuration for X-ray module
sed.add_module(self.name, self.parameters)
sed.add_info("xray.gam", self.gam)
sed.add_info("xray.alpha_ox", self.alpha_ox)
sed.add_info("xray.max_dev_alpha_ox", self.max_dev_alpha_ox)
sed.add_info("xray.alpha_ox", self.parameters['alpha_ox'])
# Calculate 0.5-2 keV hot-gas luminosities
# Mezcua et al. 2018, Eq. 5
......@@ -125,7 +132,7 @@ class Xray(SedModule):
10**(33.276 - 1.503*logT - 0.423*logT**2 + 0.425*logT**3 + 0.136*logT**4)
# Calculate L_lam_2keV from Lnu_2500A
Lnu_2keV = 10**(self.alpha_ox/0.3838) * Lnu_2500A
Lnu_2keV = 10**(alpha_ox/0.3838) * Lnu_2500A
L_lam_2keV = Lnu_2keV * self.nu_2keV**2/self.c
# Calculate total AGN corona X-ray luminosity
l_agn_xray_total = np.trapz(L_lam_2keV * self.lumin_corona,
......
......@@ -267,6 +267,7 @@ class Configuration(object):
sys.exit(1)
self.complete_redshifts()
self.complete_xray()
self.complete_analysed_parameters()
vdt = validate.Validator(validation.functions)
......@@ -359,6 +360,12 @@ class Configuration(object):
raise Exception("No flux file and no redshift indicated. "
"The spectra cannot be computed. Aborting.")
def complete_xray(self):
"""Add alpha_ox for internal usage.
"""
self.config['sed_modules_params']['xray']['alpha_ox'] = \
[-2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8]
def complete_analysed_parameters(self):
"""Complete the configuration when the variables are missing from the
configuration file and must be extracted from a dummy run."""
......
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