Commit 5921476d authored by Yannick Roehlly's avatar Yannick Roehlly
Browse files

Allow multiple module use and change the way information are stored in the SED

The modules can now be named with any prefix using a dot (e.g.
dh2002_ir.1). That means that a 'pristine' module can't have a dot in
it's name.

The add_component SED method was removed. To add components to a sed,
add_info, add_module (NEW) and add_contribution (NEW) must be used.

The modules now add more information in the sed info dictionary, even
the parametres that can be used for the statistical analysis. In
particular, mh2005_sfh add the SFR and the average SFR.

The sfr variable was renamed to sfh.
parent 66963e35
......@@ -131,14 +131,27 @@ class SED(object):
raise KeyError("The information %s is yet present "
"in the SED. " % key)
def add_component(self, module_name, module_conf, contribution_name,
results_wavelengths, results_lumin, infos):
def add_module(self, module_name, module_conf):
"""Add a new module information to the SED.
Parametres
----------
module_name : string
Name of the module. This name can be suffixed with anything
using a dot.
module_conf : dictionary
Dictionary containing the module parametres.
TODO: Complete the parametre dictionary with the default values from
the module if they are not present.
"""
Add a new module contribution to the SED
self.modules.append((module_name, module_conf))
The module name and parametres are added to the module list of the
SED. If the module adds some information to the SED, it is added to
the info dictionnary.
def add_contribution(self, contribution_name, results_wavelengths,
results_lumin):
"""
Add a new luminosity contribution to the SED.
The luminosity contribution of the module is added to the contribution
table doing an interpolation between the current wavelength grid and
......@@ -149,12 +162,6 @@ class SED(object):
Parametres
----------
module_name : string
Name of the SED creation module used.
module_conf : dictionary
The dictionnary containing the module parametres.
contribution_name : string
Name of the contribution added. This name is used to retrieve the
luminosity contribution and allows one module to add more than
......@@ -166,16 +173,8 @@ class SED(object):
results_lumin : array of floats
The vector of the Lλ luminosities (in W/nm) of the module results.
infos : dictionary
The dictionary of the informations added by the module to the SED.
If a key of the dictionary is yet present in __info, its value
will be overwritten.
"""
self.modules.append((module_name, module_conf))
self.contribution_names.append(contribution_name)
for key, value in infos.items():
self.add_info(key, value)
# If the SED luminosity table is empty, then there is nothing to
# compute.
......
......@@ -8,6 +8,7 @@ Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
"""
import numpy as np
from . import common
from ...extern.lsst import Sed as lsst
......@@ -21,6 +22,10 @@ class Module(common.SEDCreationModule):
This module is based on the code from the Large Synoptic Survey Telescope.
http://dev.lsstcorp.org/trac/
The parametres added and available for the statistical analysis are:
ccmdust_A_v, ccmdust_ebv and ccmdust_R_v.
"""
parametre_list = {
......@@ -100,20 +105,28 @@ class Module(common.SEDCreationModule):
# We only want to add the extinction as a SED component. We compute
# the difference because extended_flambda and flambda have the same
# wavelength grit (see addCCMDust definition).
# wavelength grid (see addCCMDust definition).
extinction = extended_l_lambda - l_lambda
# If the extinction was applied to a specific contribution flux, we
# suffix the ccmdust contribution name with its own.
ccmdust_contrib_name = 'ccmdust'
if parametres['contribution_name']:
ccmdust_contrib_name += '_' + parametres['contribution_name']
# Integrate the value of the extinction (-1 is because the extinction
# spectrum is negative.
extinction_value = -1 * np.trapz(extinction, wavelen)
# Base name for adding information to the SED.
name = self.name or 'ccmdust'
sed.add_module(name, parametres)
# Add the parametres values to the SED information.
sed.add_info(name + '_A_v', parametres['A_v'])
sed.add_info(name + '_ebv', parametres['ebv'])
sed.add_info(name + '_R_v', parametres['R_v'])
# Add the extinction value to the SED information
sed.add_info(name + '_extinction', extinction_value)
sed.add_component(
'ccmdust',
parametres,
ccmdust_contrib_name,
sed.add_contribution(
name,
wavelen,
extinction,
{}
extinction
)
......@@ -25,11 +25,17 @@ class SEDCreationModule(object):
# instructions for the configuration.
comments = ""
def __init__(self, **kwargs):
def __init__(self, name=None, **kwargs):
"""Instantiate a SED creation module
A name can be given to the module. This can be useful when a same
module is used several times with different parametres in the SED
creation process.
The module parametres values can be passed as keyworded paramatres.
"""
self.name = name
# parametres is a dictionnary containing the actual values for each
# module parametre.
self.parametres = kwargs
......@@ -101,24 +107,30 @@ class SEDCreationModule(object):
self._process(sed, parametres)
def get_module(module_name):
"""Return the main class of the module provided
def get_module(name):
"""Get a SED creation module from its name
Parametres
----------
module_name : string
The name of the module we want to get the class.
The name of the module we want to get the class. This name can be
prefixed by anything using a dot, then the part before the dot is
used to determine the module to load (e.g. 'dh2002_ir.1' will return
the 'dh2002_ir' module).
Returns
-------
module_class : class
a pcigale.sed.modules.Module instance
"""
# Determine the real module name by removing the dotted prefix.
module_name = name.split('.')[0]
try:
# TODO Find a better way to do dynamic import
import_string = 'from . import ' + module_name + ' as module'
exec import_string
return module.Module()
return module.Module(name=name)
except ImportError:
print('Module ' + module_name + ' does not exists!')
raise
......@@ -8,21 +8,21 @@ Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
"""
import numpy as np
from . import common
from pcigale.data import Database
class Module(common.SEDCreationModule):
"""
Compute the infra-red re-emission of absorbed energy using Dale and Helou
(2002) templates.
Module computing the infra-red re-emission corresponding to an amount of
attenuation using the Dale and Helou (2002) templates.
This module integrates the energy absorbed by the dust and normalises the
Dale and Helou (2002) template corresponding to the given α to this amount
of energy. This module can only be used on a SED that has gone through a
dust extinction module that has produced at least one extinction
(negative) contribution.
Given an amount of attenuation (e.g. resulting from the action of a dust
extinction module) this module normalises the Dale and Helou (2002)
template corresponding to a given α to this amount of energy and add it
to the SED.
Information added to the SED: NAME_alpha.
"""
......@@ -33,11 +33,11 @@ class Module(common.SEDCreationModule):
"Alpha slope.",
None
),
'extinction_contrib_names': (
'extinction_value_names': (
'array of strings',
None,
"List of the extinction contributions to process. This module "
"will add a new re-emission contribution for each one.",
"List of extinction value names (in the SED's info dictionary). "
"A new re-emission contribution will be added for each one.",
None
)
}
......@@ -52,7 +52,7 @@ class Module(common.SEDCreationModule):
"""
alpha = parametres["alpha"]
extinction_contrib_names = parametres["extinction_contrib_names"]
extinction_value_names = parametres["extinction_value_names"]
# Get the template set out of the database
database = Database()
......@@ -61,19 +61,18 @@ class Module(common.SEDCreationModule):
ir_template = dh2002.get_template(alpha)
sed.add_info('ir_alpha_slope', alpha)
# Base name for adding information to the SED.
name = self.name or 'dh2002_ir'
for contrib_name in extinction_contrib_names:
# The extinction contribution is negative
extinction = -1. * np.trapz(
sed.get_lumin_contribution(contrib_name),
sed.wavelength_grid
)
sed.add_component(
'dh2002_ir',
parametres,
'dh2002_ir_' + contrib_name,
sed.add_module(name, parametres)
sed.add_info(name + '_alpha', alpha)
# We multiply the extinction * ir_template by -1 because the value
# of the former is positive in the SED's info and its effects are
# negative.
for extinction in extinction_value_names:
sed.add_contribution(
name + '_' + extinction,
dh2002.wavelength_grid,
extinction * ir_template,
{}
-1 * sed.info[extinction] * ir_template
)
......@@ -86,11 +86,16 @@ class Module(common.SEDCreationModule):
igm_effect = red_l_lambda - init_l_lambda
sed.add_component(
'igmattenuation',
parametres,
'igmattenuation',
# Base name for adding information to the SED.
name = self.name or 'igmattenuation'
sed.add_module(name, parametres)
sed.add_info(name + '_redshift', parametres['redshift'])
sed.add_info(name + '_rtau', parametres['rtau'])
sed.add_contribution(
name,
new_wavelen,
igm_effect,
{}
igm_effect
)
......@@ -13,10 +13,11 @@ from . import common
class Module(common.SEDCreationModule):
"""Read a spectrum from a file and add it to the SED.
"""Module reading a spectrum from a file and adding it to the SED.
Note that this module uses the atpy module, which is not automatically
installed when one installs pcigale.
"""
parametre_list = {
......@@ -53,11 +54,14 @@ class Module(common.SEDCreationModule):
"""
filename = parametres['filename']
table = atpy.Table(filename)
sed.add_component(
'loadfile',
parametres,
'loadfile_' + filename,
# Base name for adding information to the SED.
name = self.name or 'loadfile'
sed.add_module(name, parametres)
sed.add_contribution(
name + '_' + filename,
table[parametres['lambda_column']],
table[parametres['l_lambda_column']],
{}
table[parametres['l_lambda_column']]
)
......@@ -13,12 +13,35 @@ from . import common
from pcigale.data import Database
# Time lapse used to compute the average star formation rate. We use a
# constant to keep it easily changeable for advanced user while limiting the
# number of parametres. The value is in Gyr.
AV_LAPSE = 0.1
class Module(common.SEDCreationModule):
"""Population synthesis based on Maraston (2005)
"""Module computing the Star Formation History contribution bases on the
Maraston (2005) models.
Implements the population syntesis based on the SSP described in Maraston,
2005, MNRAS, 362, 799.
Information added to the SED:
- imf, metallicity, galaxy_age
- sfr: star formation rate normalised to 1 solar mass formed at the
age of the galaxy.
- average_sfr: SFR averaged on the last 0.1 Gyr of the galaxy history
- mass_total, mass_alive, mass_white_dwarf,mass_neutrino,
mass_black_hole, mass_turn_off : stellar masses in solar mass.
- old_young_separation_age: age (in Gyr) separating the young and the
old star populations (if 0, there is only one population)
- mass_total_old, mass_alive_old, mass_white_dwarf_old,
mass_neutrino_old, mass_black_hole_old, mass_turn_off_old: old
star population masses.
- mass_total_young, mass_alive_young, mass_white_dwarf_young,
mass_neutrino_young, mass_black_hole_young, mass_turn_off_young:
young star population masses.
"""
parametre_list = {
......@@ -35,10 +58,10 @@ class Module(common.SEDCreationModule):
"normalised to solar values.",
None
),
'sfr': (
'sfh': (
'numpy array of floats',
"Solar mass per year",
"sfr is the star formation rate vector (in solar mass per year) "
"Star formation history (in solar mass per year) as an array "
"based on the age grid used for the Maraston 2005 SSP, i.e. "
"np.arange(1e-3,13.701,1e-3).",
None
......@@ -73,7 +96,7 @@ class Module(common.SEDCreationModule):
imf = self.parametres['imf']
metallicity = self.parametres['metallicity']
sfr = np.copy(self.parametres['sfr'])
sfh = np.copy(self.parametres['sfh'])
oldest_age = self.parametres['oldest_age']
separation_age = self.parametres['separation_age']
......@@ -81,101 +104,85 @@ class Module(common.SEDCreationModule):
database = Database()
ssp = database.get_ssp_m2005(imf, metallicity)
database.session.close_all()
# We check that the sfr array length corresponds to the one of the
# We check that the sfh array length corresponds to the one of the
# SSP age grid.
if len(ssp.age_grid) != len(sfr):
raise ValueError("The SFR array must be based on the same age "
if len(ssp.time_grid) != len(sfh):
raise ValueError("The SFH array must be based on the same age "
"grid as the Maraston 2005 SSPs.")
# We limit the star formation history to the age of the oldest stars
# and make the distinction between the old a the young populations.
# Index of the age in sfr[0] nearest to the age of the oldest stars.
idx_max = np.abs(ssp.age_grid - oldest_age).argmin()
# Index of the age in sfr[0] separating the old and the young stars.
idx_sep = np.abs(ssp.age_grid - oldest_age + separation_age).argmin()
# and set SFR=0 for every age after (we need to keep the same age
# grid for the computations.)
sfh[ssp.time_grid > oldest_age] = 0
# The age of the galaxy corresponds to the age of its oldest stars
# and is taken on the age grid.
galaxy_age = ssp.age_grid[idx_max]
# The age of the galaxy is taken from the age grid.
galaxy_age = np.max(ssp.time_grid[ssp.time_grid <= oldest_age])
# We normalise the SFH to have 1 solar mass formed at the age of the
# galaxy.
sfr = sfr / np.trapz(sfr[0:idx_max] * 1e9,
ssp.age_grid[0:idx_max])
sfh = sfh / np.trapz(sfh * 1e9, ssp.time_grid)
# First, we process the old population
old_sfr = np.copy(sfr)
old_sfr[idx_sep:] = 0
old_masses, old_spectrum = ssp.convolve(old_sfr, galaxy_age)
# First, we process the old population (age greater than the
# separation age).
old_sfh = np.copy(sfh)
old_sfh[(galaxy_age - ssp.time_grid) < separation_age] = 0
old_masses, old_spectrum = ssp.convolve(old_sfh, galaxy_age)
# If there is a separation between young and old galaxies, we process
# the young population.
# the young population else we set everything to 0 for the young
# population (always having a young population can be help full e.g.
# to probe various separation ages including 0).
if separation_age:
sed.add_info('old_young_separation_age', separation_age)
young_sfr = np.copy(sfr)
young_sfr[:idx_sep] = 0
young_sfr[idx_max:] = 0
young_masses, young_spectrum = ssp.convolve(young_sfr, galaxy_age)
# The way we enter the information in the SED object is slightly
# different, depending on if there is the age separation or not.
if not separation_age:
sed.add_component(
'm2005_sfh',
parametres,
'm2005_sfh',
ssp.wavelength_grid,
old_spectrum,
{
'stellar_mass_total': old_masses[0],
'stellar_mass_alive': old_masses[1],
'stellar_mass_white_dwarf': old_masses[2],
'stellar_mass_neutrino': old_masses[3],
'stellar_mass_black_hole': old_masses[4],
'stellar_mass_turn_off': old_masses[5]
}
)
young_sfh = np.copy(sfh)
young_sfh[(galaxy_age - ssp.time_grid) >= separation_age] = 0
young_masses, young_spectrum = ssp.convolve(young_sfh, galaxy_age)
else:
sed.add_info('stellar_mass_total',
old_masses[0] + young_masses[0])
sed.add_info('stellar_mass_alive',
old_masses[1] + young_masses[1])
sed.add_info('stellar_mass_white_dwarf',
old_masses[2] + young_masses[2])
sed.add_info('stellar_mass_neutrino',
old_masses[3] + young_masses[3])
sed.add_info('stellar_mass_black_hole',
old_masses[4] + young_masses[4])
sed.add_info('stellar_mass_turn_off',
old_masses[5] + young_masses[5])
sed.add_component(
'm2005_sfh',
parametres,
'm2005_sfh_old',
ssp.wavelength_grid,
old_spectrum,
{
'stellar_mass_total_old': old_masses[0],
'stellar_mass_alive_old': old_masses[1],
'stellar_mass_white_dwarf_old': old_masses[2],
'stellar_mass_neutrino_old': old_masses[3],
'stellar_mass_black_hole_old': old_masses[4],
'stellar_mass_turn_off_old': old_masses[5]
}
)
sed.add_component(
'm2005_sfh',
parametres,
'm2005_sfh_young',
ssp.wavelength_grid,
young_spectrum,
{
'stellar_mass_total_young': young_masses[0],
'stellar_mass_alive_young': young_masses[1],
'stellar_mass_white_dwarf_young': young_masses[2],
'stellar_mass_neutrino_young': young_masses[3],
'stellar_mass_black_hole_young': young_masses[4],
'stellar_mass_turn_off_young': young_masses[5]
}
)
young_masses = np.zeros(6)
young_spectrum = np.zeros(len(ssp.wavelength_grid))
# SFR of the galaxy
sfr = sfh[ssp.time_grid == galaxy_age][0]
# Average SFR on the last AV_LAPSE Gyr of its history
average_sfr = np.mean(sfh[(ssp.time_grid <= oldest_age) &
(ssp.time_grid >= (oldest_age - AV_LAPSE))])
# Base name for adding information to the SED.
name = self.name or 'm2005_sfh'
self.add_module(name, parametres)
sed.add_info('imf', imf)
sed.add_info('metallicity', metallicity)
sed.add_info('old_young_separation_age', separation_age)
sed.add_info('sfr', sfr)
sed.add_info('average_sfr', average_sfr)
sed.add_info('mass_total_old', old_masses[0])
sed.add_info('mass_alive_old', old_masses[1])
sed.add_info('mass_white_dwarf_old', old_masses[2])
sed.add_info('mass_neutrino_old', old_masses[3])
sed.add_info('mass_black_hole_old', old_masses[4])
sed.add_info('mass_turn_off_old', old_masses[5])
sed.add_info('mass_total_young', young_masses[0])
sed.add_info('mass_alive_young', young_masses[1])
sed.add_info('mass_white_dwarf_young', young_masses[2])
sed.add_info('mass_neutrino_young', young_masses[3])
sed.add_info('mass_black_hole_young', young_masses[4])
sed.add_info('mass_turn_off_young', young_masses[5])
sed.add_info('mass_total', old_masses[0] + young_masses[0])
sed.add_info('mass_alive', old_masses[1] + young_masses[1])
sed.add_info('mass_white_dwarf', old_masses[2] + young_masses[2])
sed.add_info('mass_neutrino', old_masses[3] + young_masses[3])
sed.add_info('mass_black_hole', old_masses[4] + young_masses[4])
sed.add_info('mass_turn_off', old_masses[5] + young_masses[5])
sed.add_contribution(name + '_old',
ssp.wavelength_grid,
old_spectrum)
sed.add_contribution(name + '_young',
ssp.wavelength_grid,
young_spectrum)
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