Implemented differential attenuation

Author: cschreib

README.md

@@ -411,6 +411,18 @@ List of available quantities:
 Some of these parameters are defined and used in [Schreiber et al. (2018)](http://adsabs.harvard.edu/cgi-bin/nph-data_query?bibcode=2018arXiv180702523S&db_key=PRE&link_type=ABSTRACT&high=5b4379bf1108433). Look at this paper for more information.
 
 
+## Differential attenuation
+The default dust model of FAST is that of the "uniform screen", in which all stars in a galaxy are placed behind a screen of dust of uniform depth (e.g., Calzetti et al. 2000). This is the simplest possible model, however it is also less physically-motivated. An alternative is the Charlot & Fall (2000) model which separates stars in two categories based on their ages: young stars are expected to still reside within their birth clouds and should thus be more attenuated by dust, while old stars have migrated out of their birth clouds and are only seen through the normal inter-stellar medium dust. This model is more refined, but also has two more free parameters: the extra amount of extinction in birth clouds, and the age at which stars exit their birth clouds.
+
+It is possible to switch to the two-component dust model by enabling the ```DIFFERENTIAL_A_V``` option. Currently, FAST++ allows you to create a grid of birth cloud extinction using the ```A_V_BC_MIN```, ```A_V_BC_MAX```, and ```A_V_BC_STEP``` parameters, as for the usual extinction. However, the age at which stars exit their birth clouds is a fixed parameter defined by ```LOG_BC_AGE_MAX``` (it may be possible to vary this parameter in a future version). When this is enabled, stars younger than ```LOG_BC_AGE_MAX``` will have their Av defined as the *sum* of the values of ```Av``` and ```Av_bc``` as taken from the grid.
+
+For example, for a grid of ```Av = [0, 1, 2]``` and ```Av_bc = [0, 0.5]```, stars in their birth clouds can have effective attenuations of ```[0, 0.5, 1, 1.5, 2, 2.5]```. And if you set ```Av = [0]``` and ```Av_bc = [0, 1, 2]```, then only the young stars will have attenuation.
+
+It is also possible to enable the differential attenuation model but still fix ```Av_bc = 0```. The models will be the same as with the standard dust model, however this allows splitting the predicted ```ldust``` (infrared luminosity) in two components, birth clouds and cirrus. This can be useful when applying infrared luminosity priors (see next section).
+
+Important note: the differential attenuation model can only be enabled when using custom star formation histories (```CUSTOM_SFH```); an error message will be shown if you try to enable it on some other SFH model.
+
+
 ## Using priors on the infrared luminosity
 One of the main degeneracy that arises when fitting UV-to-NIR data is that of dust versus age. When a galaxy has a red SED, unless the signal to noise and the wavelength sampling are high, it is very difficult to say if this is caused by a large amount of dust, or by an older stellar population, or a combination of both. It is for this reason that SFRs obtained from such fits are very uncertain, and that SFRs determined from the far-IR are preferred.
 
@@ -437,6 +449,9 @@ In this example, the fit for source ```2``` will assume Gaussian errors on the l
 
 Which value of ```LOG``` you should use depends on what probability distribution you estimate for the FIR luminosity. I give here a few guidelines to help you decide. If you are limited by the S/N on your flux measurements (i.e., for non-detections or weak detections) then ```0``` is preferable. If you are limited by the uncertainty in the conversion from flux to luminosity (i.e., when Tdust is not known) then ```1``` is preferable. Lastly, if you have high S/N fluxes and a well constrained Tdust, then either ```0``` or ```1``` will be fine.
 
+**Differential attenuation.** When ```DIFFERENTIAL_A_V``` is enabled, it is possible to apply the dust luminosity prior to different components of the galaxy, namely birth clouds or cirrus (or the sum of the two). The component that the prior applies to is then determined by an extra column in the ```.lir``` catalog. This column must be called ```COMP```, and possible values are ```all``` (birth clouds + cirrus), ```bc``` (birth clouds), or ```cirrus``` (cirrus). When available, a prior on the birth cloud luminosity provides a more direct constraint on the current SFR of the galaxy.
+
+
 ## Better treatment of spectra
 The file format for spectra in FAST-IDL is not well defined. The spectra must be given on a wavelength grid where only the central wavelength of each spectral element is specified. FAST-IDL then assumes that this wavelength grid is contiguous (no gap) and uniform (all spectral elements are defined on the same delta_lambda). No check is made to ensure this is true.
 
@@ -467,6 +482,7 @@ In this example the information in both catalogs in the same. But the new syntax
 
 Even when using the old format, the treatment of these spectrum files is also more correct in FAST++. The ```bin``` column correctly combines multiple data points into a single measurement (using inverse variance weighting) rather than simply using the first value of a bin (why this is implemented in this way in FAST-IDL, I do not know). The order of the columns in the file do not matter, while FAST-IDL assumes a fixed format (but does not tell you).
 
+
 ## Velocity broadening
 By default, the template spectra produced by FAST++ assume no velocity dispersion of the stars inside the galaxy. This is a good approximation when only fitting broadband photometry, but it becomes incorrect when fitting high spectral resolution data, such as narrow band photometry or spectra, which can sample the spectral profile of the various absorption lines that compose the spectrum (i.e., mostly the Balmer series).
 

example/cosmos-20115/fast.param

@@ -390,33 +390,48 @@ CUSTOM_SFH_LOOKBACK = 1              # 0 / 1
 #   You can choose more than one by defining an array.
 #
 # o If a grid has already been made for a specific LIBRARY, RESOLUTION,
-#   IMF, tau, age, z, A_v, metallicity, and filter set and/or spectral
-#   elements, the grid will be automatically read from the cache unless
-#   the NO_CACHE option is set to "1". This option will also prevent the
-#   code from writing this cache in the first place.
+#   IMF, tau, age, z, A_v, A_v_bc, metallicity, and filter set and/or
+#   spectral elements, the grid will be automatically read from the cache
+#   unless the NO_CACHE option is set to "1". This option will also
+#   prevent the code from writing this cache in the first place.
+#
+# o DIFFERENTIAL_A_V: By default (0), all stars suffer the same amount
+#   of extinction by dust. If DIFFERENTIAL_A_V is set to "1", young
+#   stars will be place under a thicker screen of dust, as controlled by
+#   A_V_BC_MIN, A_V_BC_MAX, A_V_BC_STEP, and LOG_BC_AGE_MAX.
+#
+# o LOG_BC_AGE_MAX: When DIFFERENTIAL_A_V is set to "1", stars in their
+#   birth cloud will be given extra attenuation. Only stars younger than
+#   LOG_BC_AGE_MAX are considered to be in their birthcloud. The default
+#   of 10 Myr is taken from Charlot&Fall (2001).
 #
 # o NO_MAX_AGE: By default (0), ages that exceed the age of the universe
 #   are not allowed. However, when NO_MAX_AGE is put to "1" you can have
 #   older ages.
 #
 #-----------------------------------------------------------------------
 
-LOG_TAU_MIN    = 8.5            # log [yr]
-LOG_TAU_MAX    = 10.            # log [yr]
-LOG_TAU_STEP   = 0.5            # log [yr], min 0.1
-LOG_AGE_MIN    = 9.21931        # log [yr]
-LOG_AGE_MAX    = 9.21931        # log [yr]
-LOG_AGE_STEP   = 0.05           # log [yr]
-NO_MAX_AGE     = 1              # 0 / 1
-Z_MIN          = 3.71453        # Cannot be 0.
-Z_MAX          = 3.71453
-Z_STEP         = 0.1
-Z_STEP_TYPE    = 0              # 0: Z_STEP, 1: Z_STEP*(1+z)
-A_V_MIN        = 0.0            # [mag]
-A_V_MAX        = 1.0            # [mag]
-A_V_STEP       = 0.1            # [mag]
-METAL          = 0.02           # [0.0096,0.019,0.03]
-NO_CACHE       = 1              # 0 / 1
+LOG_TAU_MIN      = 8.5            # log [yr]
+LOG_TAU_MAX      = 10.            # log [yr]
+LOG_TAU_STEP     = 0.5            # log [yr], min 0.1
+LOG_AGE_MIN      = 9.21931        # log [yr]
+LOG_AGE_MAX      = 9.21931        # log [yr]
+LOG_AGE_STEP     = 0.05           # log [yr]
+NO_MAX_AGE       = 1              # 0 / 1
+Z_MIN            = 3.71453        # Cannot be 0.
+Z_MAX            = 3.71453
+Z_STEP           = 0.1
+Z_STEP_TYPE      = 0              # 0: Z_STEP, 1: Z_STEP*(1+z)
+A_V_MIN          = 0.0            # [mag]
+A_V_MAX          = 1.0            # [mag]
+A_V_STEP         = 0.1            # [mag]
+A_V_BC_MIN       = 0.             # [mag]
+A_V_BC_MAX       = 0.             # [mag]
+A_V_BC_STEP      = 0.1            # [mag]
+DIFFERENTIAL_A_V = 0              # 0 / 1
+LOG_BC_AGE_MAX   = 7.0            # log [yr]
+METAL            = 0.02           # [0.0096,0.019,0.03]
+NO_CACHE         = 1              # 0 / 1
 
 #--- CUSTOM GRID PARAMETERS --------------------------------------------
 

example/hdfn_fs99/fast.param

@@ -393,33 +393,48 @@ CUSTOM_SFH_LOOKBACK = 0              # 0 / 1
 #   You can choose more than one by defining an array.
 #
 # o If a grid has already been made for a specific LIBRARY, RESOLUTION,
-#   IMF, tau, age, z, A_v, metallicity, and filter set and/or spectral
-#   elements, the grid will be automatically read from the cache unless
-#   the NO_CACHE option is set to "1". This option will also prevent the
-#   code from writing this cache in the first place.
+#   IMF, tau, age, z, A_v, A_v_bc, metallicity, and filter set and/or
+#   spectral elements, the grid will be automatically read from the cache
+#   unless the NO_CACHE option is set to "1". This option will also
+#   prevent the code from writing this cache in the first place.
+#
+# o DIFFERENTIAL_A_V: By default (0), all stars suffer the same amount
+#   of extinction by dust. If DIFFERENTIAL_A_V is set to "1", young
+#   stars will be place under a thicker screen of dust, as controlled by
+#   A_V_BC_MIN, A_V_BC_MAX, A_V_BC_STEP, and LOG_BC_AGE_MAX.
+#
+# o LOG_BC_AGE_MAX: When DIFFERENTIAL_A_V is set to "1", stars in their
+#   birth cloud will be given extra attenuation. Only stars younger than
+#   LOG_BC_AGE_MAX are considered to be in their birthcloud. The default
+#   of 10 Myr is taken from Charlot&Fall (2001).
 #
 # o NO_MAX_AGE: By default (0), ages that exceed the age of the universe
 #   are not allowed. However, when NO_MAX_AGE is put to "1" you can have
 #   older ages.
 #
 #-----------------------------------------------------------------------
 
-LOG_TAU_MIN    = 8.5            # log [yr]
-LOG_TAU_MAX    = 10.            # log [yr]
-LOG_TAU_STEP   = 0.5            # log [yr], min 0.1
-LOG_AGE_MIN    = 8.0            # log [yr]
-LOG_AGE_MAX    = 10.0           # log [yr]
-LOG_AGE_STEP   = 0.2            # log [yr]
-NO_MAX_AGE     = 0              # 0 / 1
-Z_MIN          = 0.01           # Cannot be 0.0
-Z_MAX          = 6.00
-Z_STEP         = 0.05
-Z_STEP_TYPE    = 0              # 0: Z_STEP, 1: Z_STEP*(1+z)
-A_V_MIN        = 0.             # [mag]
-A_V_MAX        = 3.             # [mag]
-A_V_STEP       = 0.1            # [mag]
-METAL          = [0.02]         # [0.0096,0.019,0.03]
-NO_CACHE       = 0              # 0 / 1
+LOG_TAU_MIN      = 8.5            # log [yr]
+LOG_TAU_MAX      = 10.            # log [yr]
+LOG_TAU_STEP     = 0.5            # log [yr], min 0.1
+LOG_AGE_MIN      = 8.0            # log [yr]
+LOG_AGE_MAX      = 10.0           # log [yr]
+LOG_AGE_STEP     = 0.2            # log [yr]
+NO_MAX_AGE       = 0              # 0 / 1
+Z_MIN            = 0.01           # Cannot be 0.0
+Z_MAX            = 6.00
+Z_STEP           = 0.05
+Z_STEP_TYPE      = 0              # 0: Z_STEP, 1: Z_STEP*(1+z)
+A_V_MIN          = 0.             # [mag]
+A_V_MAX          = 3.             # [mag]
+A_V_STEP         = 0.1            # [mag]
+A_V_BC_MIN       = 0.             # [mag]
+A_V_BC_MAX       = 0.             # [mag]
+A_V_BC_STEP      = 0.1            # [mag]
+DIFFERENTIAL_A_V = 0              # 0 / 1
+LOG_BC_AGE_MAX   = 7.0            # log [yr]
+METAL            = [0.02]         # [0.0096,0.019,0.03]
+NO_CACHE         = 0              # 0 / 1
 
 
 #--- COSMOLOGY ---------------------------------------------------------

src/fast++-fitter.cpp

@@ -471,7 +471,17 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
             // Model LIR with Gaussian likelihood
             wsp.weight[0] = 1.0/input.lir_err.safe[is];
             wsp.wflux[0] = input.lir.safe[is]*wsp.weight[0];
-            wsp.wmodel[0] = model.props.safe[prop_id::ldust]*wsp.weight[0];
+
+            double lir_model = 0.0;
+            if (input.lir_comp.safe[is] == lir_component::all) {
+                lir_model = model.props.safe[prop_id::ldust];
+            } else if (input.lir_comp.safe[is] == lir_component::bc) {
+                lir_model = model.props.safe[prop_id::ldust_bc];
+            } else if (input.lir_comp.safe[is] == lir_component::cirrus) {
+                lir_model = model.props.safe[prop_id::ldust] - model.props.safe[prop_id::ldust_bc];
+            }
+
+            wsp.wmodel[0] = lir_model*wsp.weight[0];
 
             wfm += wsp.wmodel[0]*wsp.wflux[0];
             wmm += sqr(wsp.wmodel[0]);
@@ -533,7 +543,16 @@ void fitter_t::fit_galaxies(const model_t& model, uint_t i0, uint_t i1) {
 
         // Add LIR as a contribution to chi2 (if given in log units)
         if (!input.lir.empty() && is_finite(input.lir.safe[is]) && input.lir_log.safe[is]) {
-            double log_model = log10(scale*max(model.props.safe[prop_id::ldust], 1e-4));
+            double lir_model = 0.0;
+            if (input.lir_comp.safe[is] == lir_component::all) {
+                lir_model = model.props.safe[prop_id::ldust];
+            } else if (input.lir_comp.safe[is] == lir_component::bc) {
+                lir_model = model.props.safe[prop_id::ldust_bc];
+            } else if (input.lir_comp.safe[is] == lir_component::cirrus) {
+                lir_model = model.props.safe[prop_id::ldust] - model.props.safe[prop_id::ldust_bc];
+            }
+
+            double log_model = log10(scale*max(lir_model, 1e-4));
             tchi2 += sqr((input.lir.safe[is] - log_model)/input.lir_err.safe[is]);
         }