Commit 3dc4cf8b authored by Mac Denis's avatar Mac Denis
Browse files

Merge branch 'develop' of gitlab.oamp.dev:yannick/pcigale into develop

parents 1b15860b ee8cdca4
# IRAS1
#photon
# IRAS 12um
70000.0 0.0000000000
75000.0 0.0080000000
80000.0 0.5350000000
85000.0 0.6890000000
90000.0 0.7350000000
95000.0 0.8150000000
100000.0 0.9000000000
105000.0 0.9040000000
110000.0 0.8340000000
115000.0 0.8160000000
120000.0 0.7930000000
125000.0 0.8540000000
130000.0 0.9380000000
135000.0 0.9910000000
140000.0 1.0000000000
145000.0 0.9340000000
150000.0 0.3880000000
155000.0 0.0000000000
# IRAS2
#photon
# IRAS 25um
155000.0 0.0000000000
160000.0 0.0070000000
165000.0 0.1010000000
170000.0 0.2880000000
175000.0 0.3880000000
180000.0 0.4520000000
185000.0 0.5210000000
190000.0 0.5620000000
195000.0 0.6260000000
200000.0 0.6830000000
205000.0 0.7290000000
210000.0 0.7780000000
215000.0 0.8320000000
220000.0 0.9120000000
225000.0 0.9140000000
230000.0 0.9380000000
235000.0 0.9330000000
240000.0 0.8750000000
245000.0 0.9100000000
250000.0 1.0000000000
255000.0 0.9110000000
260000.0 0.8400000000
265000.0 0.7630000000
270000.0 0.7490000000
275000.0 0.8290000000
280000.0 0.9140000000
285000.0 0.7900000000
290000.0 0.8770000000
295000.0 0.5580000000
300000.0 0.2740000000
305000.0 0.0690000000
310000.0 0.0120000000
315000.0 0.0000000000
# IRAS3
#photon
# IRAS 60um
270000.0 0.0000000000
300000.0 0.0100000000
330000.0 0.0360000000
360000.0 0.0680000000
390000.0 0.1740000000
420000.0 0.3150000000
450000.0 0.4830000000
480000.0 0.5850000000
510000.0 0.6580000000
540000.0 0.7160000000
570000.0 0.8240000000
600000.0 0.9150000000
630000.0 0.9870000000
660000.0 0.9900000000
690000.0 1.0000000000
720000.0 0.9460000000
750000.0 0.7130000000
780000.0 0.5310000000
810000.0 0.1740000000
840000.0 0.0470000000
870000.0 0.0000000000
# IRAS4
#photon
# IRAS 100um
650000.0 0.0000000000
700000.0 0.0100000000
750000.0 0.1130000000
800000.0 0.3060000000
850000.0 0.5050000000
900000.0 0.6950000000
950000.0 0.8240000000
1000000.0 0.9470000000
1050000.0 0.9390000000
1100000.0 1.0000000000
1150000.0 0.6310000000
1200000.0 0.3190000000
1250000.0 0.1950000000
1300000.0 0.1060000000
1350000.0 0.0530000000
1400000.0 0.0100000000
1450000.0 0.0000000000
# SCUBA450
#photon
# SCUBA 450um
375000.00 1.0E-4
380035.47 2.0E-4
390015.60 4.0E-4
400000.00 0.0012
410004.10 0.0046
420050.41 0.0178
430045.87 0.1134
435034.80 0.2508
440011.73 0.3843
445037.83 0.4042
450045.00 0.3963
455028.06 0.3803
460052.14 0.2996
465044.18 0.1942
470072.08 0.1072
475059.38 0.0482
480000.00 0.0212
485044.46 0.0117
490035.94 0.0105
500000.00 0.0074
510030.60 0.0047
520020.80 0.0029
540054.00 9.0E-4
560014.93 3.0E-4
580046.40 2.0E-4
600000.00 1.0E-4
# SCUBA850
#photon
# SCUBA 850um
7400098.67 0.0
7500000.00 0.001
7600709.40 0.004
7700205.34 0.008
7800312.01 0.016
7900974.45 0.03
8000000.00 0.059
8101539.29 0.124
8201202.84 0.24
8301051.47 0.365
8401008.12 0.45
8500991.78 0.452
8600917.43 0.446
8700696.06 0.443
8800234.67 0.418
8902077.15 0.352
9000900.09 0.271
9101941.75 0.17
9202453.99 0.104
9302325.58 0.062
9401441.55 0.04
9502692.43 0.024
9600000.00 0.016
9702457.96 0.01
9800718.72 0.007
9900990.10 0.005
1.000000E7 0.003
1.025290E7 0.001
1.050053E7 0.0
# VLA_C
#photon
# C-band VLA 4.85GHz - 6cm
219000000 0.0
319000000 1.0
619000000 1.0
919000000 1.0
1019000000 0.0
# VLA_L
#photon
# L-band VLA 1.49GHz - 20cm
1413000000 0.0
1513000000 1.0
2013000000 1.0
2513000000 1.0
2613000000 0.0
\ No newline at end of file
......@@ -192,7 +192,7 @@ def sed(idx):
format(n_computed, gbl_params.size,
np.around(t_elapsed, decimals=1),
np.around(n_computed/t_elapsed, decimals=1)),
end="\r")
end="\n" if n_computed == gbl_params.size else "\r")
def analysis(idx, obs):
......@@ -401,5 +401,5 @@ def analysis(idx, obs):
print("{}/{} objects analysed in {} seconds ({} objects/s)".
format(n_computed, gbl_n_obs, np.around(t_elapsed, decimals=1),
np.around(n_computed/t_elapsed, decimals=2)),
end="\r")
end="\n" if idx == gbl_n_obs-1 else "\r")
......@@ -100,4 +100,4 @@ def fluxes(idx):
format(n_computed, gbl_params.size,
np.around(t_elapsed, decimals=1),
np.around(n_computed/t_elapsed, decimals=1)),
end="\r")
end="\n" if n_computed == gbl_params.size else "\r")
......@@ -163,10 +163,10 @@ class M2005(CreationModule):
sed.add_info('stellar.mass_black_hole',
old_masses[4] + young_masses[4], True)
sed.add_contribution("ssp_old",
sed.add_contribution("stellar.old",
ssp.wavelength_grid,
old_spectrum)
sed.add_contribution("ssp_young",
sed.add_contribution("stellar.young",
ssp.wavelength_grid,
young_spectrum)
......
......@@ -18,6 +18,9 @@ import matplotlib.pyplot as plt
import multiprocessing as mp
import numpy as np
import os
import pkg_resources
from scipy.constants import c, parsec
from pcigale.sed.cosmology import cosmology
from pcigale.data import Database
from pcigale.utils import read_table
from pcigale.session.configuration import Configuration
......@@ -28,8 +31,8 @@ BEST_MODEL_FILE = "best_models.txt"
# Directory where the output files are stored
OUT_DIR = "out/"
# Wavelength limits (restframe) when plotting the best SED.
PLOT_L_MIN = 91.
PLOT_L_MAX = 1e6
PLOT_L_MIN = 0.1
PLOT_L_MAX = 2e6
def _chi2_worker(obj_name, var_name):
......@@ -88,7 +91,7 @@ def _pdf_worker(obj_name, var_name):
print("No PDF found for {}. No plot created.".format(obj_name))
def _sed_worker(obs, mod, filters):
def _sed_worker(obs, mod, filters, sed_type, nologo):
"""Plot the best SED with the associated fluxes in bands
Parameters
......@@ -99,61 +102,181 @@ def _sed_worker(obs, mod, filters):
Data from the best model of one object.
filters: ordered dictionary of Filter objects
The observed fluxes in each filter.
sed_type: string
Type of SED to plot. It can either be "mJy" (flux in mJy and observed
frame) or "lum" (luminosity in W and rest frame)
nologo: boolean
Do not add the logo when set to true.
"""
if os.path.isfile(OUT_DIR + "{}_best_model.xml".format(obs['id'])):
sed = Table.read(OUT_DIR + "{}_best_model.xml".format(obs['id']),
table_id="Fnu")
table_id="Flambda")
filters_wl = np.array([filt.effective_wavelength
for filt in filters.values()])
wavelength_spec = sed['wavelength']
obs_fluxes = np.array([obs[filt] for filt in filters.keys()])
obs_fluxes_err = np.array([obs[filt+'_err']
for filt in filters.keys()])
mod_fluxes = np.array([mod[filt] for filt in filters.keys()])
xmin = PLOT_L_MIN * (1. + obs['redshift'])
xmax = PLOT_L_MAX * (1. + obs['redshift'])
wsed = np.where((sed['wavelength'] > xmin)&(sed['wavelength'] < xmax))
z = obs['redshift']
DL = max(10, cosmology.luminosity_distance(z).value * 1e6) * parsec
if sed_type == 'lum':
xmin = PLOT_L_MIN
xmax = PLOT_L_MAX
k_corr_SED = 1e-29 * (4.*np.pi*DL*DL) * c / (filters_wl*1e-9)
obs_fluxes *= k_corr_SED
obs_fluxes_err *= k_corr_SED
mod_fluxes *= k_corr_SED
for cname in sed.colnames[1:]:
sed[cname] *= wavelength_spec
filters_wl /= 1. + z
wavelength_spec /= 1. + z
elif sed_type == 'mJy':
xmin = PLOT_L_MIN * (1. + z)
xmax = PLOT_L_MAX * (1. + z)
for cname in sed.colnames[1:]:
sed[cname] *= (wavelength_spec * 1e29 /
(c / (wavelength_spec * 1e-9)) /
(4. * np.pi * DL * DL))
else:
print("Unknown plot type")
filters_wl /= 1000.
wavelength_spec /= 1000.
wsed = np.where((wavelength_spec > xmin) & (wavelength_spec < xmax))
figure = plt.figure()
gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1])
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
ax1.loglog(sed['wavelength'][wsed], sed['F_nu'][wsed],
label="Model spectrum", color='k', nonposy='clip')
ax1.set_autoscale_on(False)
ax1.scatter(filters_wl, mod_fluxes, marker='o', color='r',
label="Model fluxes")
ax1.errorbar(filters_wl, obs_fluxes, yerr=obs_fluxes_err*3, ls='',
marker='_', label='Observed fluxes', color='b',
capsize=0.)
mask = np.where(obs_fluxes > 0.)
ax2.errorbar(filters_wl[mask],
(obs_fluxes[mask]-mod_fluxes[mask])/obs_fluxes[mask],
yerr=obs_fluxes_err[mask]/obs_fluxes[mask]*3, marker='_',
label="(Obs-Mod)/Obs", color='k', capsize=0.)
ax2.plot([xmin, xmax], [0., 0.], ls='--', color='k')
ax2.set_xscale('log')
ax2.minorticks_on()
figure.subplots_adjust(hspace=0.,wspace=0.)
ax1.set_xlim(xmin, xmax)
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(-1.0, 1.0)
ax2.set_xlabel("Observed wavelength [nm]")
ax1.set_ylabel("Flux [mJy]")
ax2.set_ylabel("Relative residual flux")
ax1.legend(loc='upper left', fontsize=10)
ax2.legend(loc='upper right', fontsize=10)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax1.get_yticklabels()[1], visible=False)
figure.suptitle("Best model for {} at z = {}. Reduced $\chi^2$={}".
format(
obs['id'], np.round(obs['redshift'], decimals=3),
np.round(mod['reduced_chi_square'], decimals=2))
)
figure.savefig(OUT_DIR + "{}_best_model.pdf".format(obs['id']))
plt.close(figure)
if (sed.columns[1][wsed] > 0.).any():
ax1 = plt.subplot(gs[0])
ax2 = plt.subplot(gs[1])
# Stellar emission
ax1.loglog(wavelength_spec[wsed],
(sed['stellar.young'][wsed] +
sed['attenuation.stellar.young'][wsed] +
sed['stellar.old'][wsed] +
sed['attenuation.stellar.old'][wsed]),
label="Stellar attenuated ", color='orange',
marker=None, nonposy='clip', linestyle='-',
linewidth=0.5)
ax1.loglog(wavelength_spec[wsed],
(sed['stellar.old'][wsed] +
sed['stellar.young'][wsed]),
label="Stellar unattenuated", color='b', marker=None,
nonposy='clip', linestyle='--', linewidth=0.5)
# Dust emission Draine & Li
if 'dust.Umin_Umin' in sed.columns:
ax1.loglog(wavelength_spec[wsed],
(sed['dust.Umin_Umin'][wsed] +
sed['dust.Umin_Umax'][wsed]),
label="Dust emission", color='r', marker=None,
nonposy='clip', linestyle='-', linewidth=0.5)
# Dust emission Dale
if 'dust' in sed.columns:
ax1.loglog(wavelength_spec[wsed], sed['dust'][wsed],
label="Dust emission", color='r', marker=None,
nonposy='clip', linestyle='-', linewidth=0.5)
# AGN emission Fritz
if 'agn_fritz2006_therm' in sed.columns:
ax1.loglog(wavelength_spec[wsed],
(sed['agn_fritz2006_therm'][wsed] +
sed['agn_fritz2006_scatt'][wsed] +
sed['agn_fritz2006_agn'][wsed]),
label="AGN emission", color='g', marker=None,
nonposy='clip', linestyle='-', linewidth=0.5)
# Radio emission
if 'radio_nonthermal' in sed.columns:
ax1.loglog(wavelength_spec[wsed],
sed['radio_nonthermal'][wsed],
label="Radio nonthermal", color='brown',
marker=None, nonposy='clip', linestyle='-',
linewidth=0.5)
ax1.loglog(wavelength_spec[wsed], sed['F_lambda_total'][wsed],
label="Model spectrum", color='k', nonposy='clip',
linestyle='-', linewidth=2.0)
ax1.set_autoscale_on(False)
ax1.scatter(filters_wl, mod_fluxes, marker='o', color='r', s=8,
label="Model fluxes")
mask_ok = np.logical_and(obs_fluxes > 0., obs_fluxes_err > 0.)
ax1.errorbar(filters_wl[mask_ok], obs_fluxes[mask_ok],
yerr=obs_fluxes_err[mask_ok]*3, ls='', marker='s',
label='Observed fluxes', markerfacecolor='None',
markersize=6, markeredgecolor='b', capsize=0.)
mask_uplim = np.logical_and(np.logical_and(obs_fluxes > 0.,
obs_fluxes_err < 0.),
obs_fluxes_err > -9990.)
if not mask_uplim.any() == False:
ax1.errorbar(filters_wl[mask_uplim], obs_fluxes[mask_uplim],
uplims=obs_fluxes_err[mask_uplim], ls='',
marker='v', label='Observed upper limits',
markerfacecolor='None', markersize=6,
markeredgecolor='b', capsize=0.)
mask_noerr = np.logical_and(obs_fluxes > 0.,
obs_fluxes_err < -9990.)
if not mask_noerr.any() == False:
ax1.errorbar(filters_wl[mask_noerr], obs_fluxes[mask_noerr],
ls='', marker='s', markerfacecolor='None',
markersize=6, markeredgecolor='g',
label='Observed fluxes, no errors', capsize=0.)
mask = np.where(obs_fluxes > 0.)
ax2.errorbar(filters_wl[mask],
(obs_fluxes[mask]-mod_fluxes[mask])/obs_fluxes[mask],
yerr=obs_fluxes_err[mask]/obs_fluxes[mask]*3,
marker='_', label="(Obs-Mod)/Obs", color='k',
capsize=0.)
ax2.plot([xmin, xmax], [0., 0.], ls='--', color='k')
ax2.set_xscale('log')
ax2.minorticks_on()
figure.subplots_adjust(hspace=0., wspace=0.)
ax1.set_xlim(xmin, xmax)
ymin = min(np.min(obs_fluxes[mask_ok]),
np.min(mod_fluxes[mask_ok]))
ymax = max(np.max(obs_fluxes[mask_ok]),
np.max(mod_fluxes[mask_ok]))
ax1.set_ylim(1e-2*ymin, 1e2*ymax)
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(-1.0, 1.0)
if sed_type == 'lum':
ax2.set_xlabel("Rest-frame wavelength [$\mu$m]")
ax1.set_ylabel("Luminosity [W]")
ax2.set_ylabel("Relative residual luminosity")
else:
ax2.set_xlabel("Observed wavelength [$\mu$m]")
ax1.set_ylabel("Flux [mJy]")
ax2.set_ylabel("Relative residual flux")
ax1.legend(fontsize=6, loc='best', fancybox=True, framealpha=0.5)
ax2.legend(fontsize=6, loc='best', fancybox=True, framealpha=0.5)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax1.get_yticklabels()[1], visible=False)
figure.suptitle("Best model for {} at z = {}. Reduced $\chi^2$={}".
format(obs['id'], np.round(obs['redshift'],
decimals=3),
np.round(mod['reduced_chi_square'],
decimals=2)))
if nologo is False:
image = plt.imread(pkg_resources.resource_filename(__name__,
"data/CIGALE.png"))
figure.figimage(image, 75, 330, origin='upper', zorder=10,
alpha=1)
figure.savefig(OUT_DIR + "{}_best_model.pdf".format(obs['id']))
plt.close(figure)
else:
print("No valid best SED found for {}. No plot created.".
format(obs['id']))
else:
print("No SED found for {}. No plot created.".format(obs['id']))
......@@ -172,20 +295,6 @@ def chi2(config):
pool.join()
def chi2(config):
"""Plot the χ² values of analysed variables.
"""
input_data = read_table(config.configuration['data_file'])
chi2_vars = (config.configuration['analysis_method_params']
['analysed_variables'])
with mp.Pool(processes=config.configuration['cores']) as pool:
items = product(input_data['id'], chi2_vars)
pool.starmap(_chi2_worker, items)
pool.close()
pool.join()
def pdf(config):
"""Plot the PDF of analysed variables.
"""
......@@ -200,18 +309,20 @@ def pdf(config):
pool.join()
def sed(config):
def sed(config, sed_type, nologo):
"""Plot the best SED with associated observed and modelled fluxes.
"""
obs = read_table(config.configuration['data_file'])
mod = Table.read(OUT_DIR + BEST_MODEL_FILE, format='ascii')
with Database() as base:
filters = OrderedDict([(name, base.get_filter(name))
for name in config.configuration['column_list']
if not name.endswith('_err')])
with mp.Pool(processes=config.configuration['cores']) as pool:
pool.starmap(_sed_worker, zip(obs, mod, repeat(filters)))
pool.starmap(_sed_worker, zip(obs, mod, repeat(filters),
repeat(sed_type), repeat(nologo)))
pool.close()
pool.join()
......@@ -232,6 +343,8 @@ def main():
chi2_parser.set_defaults(parser='chi2')
sed_parser = subparsers.add_parser('sed', help=sed.__doc__)
sed_parser.add_argument('--type', default='mJy')
sed_parser.add_argument('--nologo', action="store_true")
sed_parser.set_defaults(parser='sed')
args = parser.parse_args()
......@@ -246,4 +359,4 @@ def main():
elif args.parser == 'pdf':
pdf(config)
elif args.parser == 'sed':
sed(config)
sed(config, args.type, args.nologo)
......@@ -32,7 +32,7 @@ setup(
include_package_data=True,
cmdclass={"build": custom_build},
package_data={'': ['*.db']},
package_data={'pcigale': ['data/data.db'], 'pcigale-plots':['data/CIGALE.png']},
author="Yannick Roehlly",
author_email="yannick.roehlly@oamp.fr",
......
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