Coverage for src/galsbi/ucat_sps/plugins/sample_galaxies_sps.py: 100%
215 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-09 18:05 +0000
1# Copyright (C) 2025 LMU Munich
2# Copyright (C) 2025 ETH Zurich, Institute for Particle Physics and Astrophysics
3# Author: Luca Tortorelli, Silvan Fischbacher
4# created: Apr 2025
6import warnings
8import healpy as hp
9import numpy as np
10import PyCosmo
11from cosmic_toolbox import logger
12from ivy.plugin.base_plugin import BasePlugin
13from ufig import coordinate_util
15from galsbi.ucat import filters_util, galaxy_sampling_util, utils
16from galsbi.ucat.galaxy_population_models.galaxy_position import sample_position_uniform
17from galsbi.ucat.galaxy_population_models.galaxy_shape import (
18 sample_ellipticities_for_galaxy_type,
19)
20from galsbi.ucat.plugins.sample_galaxies_photo import ExtinctionMapEvaluator, in_pos
21from galsbi.ucat_sps.galaxy_population_models.galaxy_agn import (
22 sample_ratio_agn_to_galaxy_bolometric_luminosity,
23)
24from galsbi.ucat_sps.galaxy_population_models.galaxy_dust_attenuation import (
25 sample_dust_attenuation_parameters,
26 stellar_mass_and_ssfr_for_dust_attenuation,
27)
28from galsbi.ucat_sps.galaxy_population_models.galaxy_dust_emission import (
29 sample_Dale2014_dust_emission_parameters,
30)
31from galsbi.ucat_sps.galaxy_population_models.galaxy_gas_ionization import (
32 sample_gas_ionization_Kashino19,
33)
34from galsbi.ucat_sps.galaxy_population_models.galaxy_light_profile import (
35 sample_sersic_for_galaxy_type,
36)
37from galsbi.ucat_sps.galaxy_population_models.galaxy_metallicity_history import (
38 generate_gas_metallicity_history_on_ssp_grid,
39 generate_gas_metallicity_history_snorm_trunc,
40 sample_gas_metallicity,
41)
42from galsbi.ucat_sps.galaxy_population_models.galaxy_size import (
43 sample_r50_from_stellar_mass,
44)
45from galsbi.ucat_sps.config.common_sps import (
46 BURST_SURVIVING_MASS_MODEL,
47 CONDITIONAL_SNORM_BURST_TRUNC,
48 SFH_POPULATION_MODELS,
49 SMOOTH_SURVIVING_MASS_MODEL,
50 TORTORELLI25_SNORM_TRUNC,
51)
52from galsbi.ucat_sps.galaxy_population_models.galaxy_star_formation_history import (
53 evaluate_sfh_on_ssp_grid,
54 load_ssp_time_grid,
55 sample_conditional_sfh_parameters,
56 sample_sfh_snorm_trunc,
57)
58from galsbi.ucat_sps.galaxy_population_models.galaxy_surviving_stellar_mass import (
59 compute_legacy_surviving_stellar_mass,
60 load_surviving_stellar_mass_model,
61 predict_log10_surviving_stellar_mass,
62)
63from galsbi.ucat_sps.galaxy_population_models.galaxy_stellar_mass_function import (
64 initialize_stellar_mass_functions,
65)
66from galsbi.ucat_sps.galaxy_population_models.galaxy_velocity_dispersion import (
67 sample_velocity_dispersion_zahid,
68)
69from galsbi.ucat_sps.sps_sed_generator import MagnitudeGenerator
70from galsbi.ucat_sps.util_sps import (
71 apply_galaxy_observed_magnitude_cut,
72 galaxy_property_seed,
73 get_lookbacktime_from_redshift,
74 stellar_mass_function_seed,
75)
77LOGGER = logger.get_logger(__file__)
78warnings.filterwarnings("once")
81def check_max_mem_error(max_mem_hard_limit=10000):
82 """
83 Check if the catalog does not exceed allowed memory.
85 Prevents job crashes on clusters.
86 """
88 mem_mb_current = utils.memory_usage_psutil()
89 if mem_mb_current > max_mem_hard_limit:
90 raise galaxy_sampling_util.UCatNumGalError(
91 "The sample_galaxies process is taking too much memory: "
92 f"mem_mb_current={mem_mb_current}, max_mem_hard_limit={max_mem_hard_limit}"
93 )
96class Plugin(BasePlugin):
97 """
98 Generate a catalog of galaxies with SPS computed magnitudes in multiple bands.
99 """
101 def __call__(self):
102 par = self.ctx.parameters
104 # Cosmology
105 cosmo = PyCosmo.build()
106 cosmo.set(h=par.h, omega_m=par.omega_m)
108 # Healpix pixelization
109 if par.sampling_mode == "wcs":
110 LOGGER.info("Sampling galaxies based on RA/DEC and pixel scale")
111 # Healpix pixelization
112 w = coordinate_util.wcs_from_parameters(par)
113 self.ctx.pixels = coordinate_util.get_healpix_pixels(
114 par.nside_sampling, w, par.size_x, par.size_y
115 )
116 if len(self.ctx.pixels) < 15:
117 LOGGER.warning(
118 f"Only {len(self.ctx.pixels)} healpy pixels in the footprint,"
119 " consider increasing the nside_sampling"
120 )
121 elif par.sampling_mode == "healpix":
122 LOGGER.info("Sampling galaxies based on healpix pixels")
123 self.ctx.pixels = coordinate_util.get_healpix_pixels_from_map(par)
124 w = None
125 else:
126 raise ValueError(
127 f"Unknown sampling mode: {par.sampling_mode}, must be wcs or healpix"
128 )
129 self.ctx.pixarea = hp.nside2pixarea(par.nside_sampling, degrees=False)
131 # Magnitude calculator
132 all_filters = np.unique(par.filters + [par.lum_fct_filter_band])
134 # backward compatibility - check if full filter names are set
135 if not hasattr(par, "filters_full_names"):
136 par.filters_full_names = filters_util.get_default_full_filter_names(
137 all_filters
138 )
139 warnings.warn(
140 "setting filters to default, this will cause problems if you work"
141 " with filters from different cameras in the same band",
142 stacklevel=1,
143 )
145 # initialize stellar mass functions
146 sm_funcs_lowmass = initialize_stellar_mass_functions(
147 par,
148 cosmo=cosmo,
149 pixarea=self.ctx.pixarea,
150 z_sm_intp=None,
151 mass_key="lowmass", # z_m_intp
152 )
154 sm_funcs_highmass = initialize_stellar_mass_functions(
155 par,
156 cosmo=cosmo,
157 pixarea=self.ctx.pixarea,
158 z_sm_intp=None,
159 mass_key="highmass", # z_m_intp
160 )
162 # setup SPS SED generator
163 sps_sed_generator = MagnitudeGenerator(par)
165 # Select the SFH population model and load the SSP grid once.
166 sfh_population_model = getattr(
167 par, "sfh_population_model", TORTORELLI25_SNORM_TRUNC
168 )
169 if sfh_population_model not in SFH_POPULATION_MODELS:
170 raise ValueError(
171 f"Unknown sfh_population_model={sfh_population_model!r}; expected "
172 f"one of {sorted(SFH_POPULATION_MODELS)}."
173 )
175 ssp_age_yr = None
176 ssp_age_bins_yr = None
177 ssp_age_weights_yr = None
178 smooth_surviving_mass_emulator = None
179 smooth_surviving_mass_limits = None
180 burst_surviving_mass_emulator = None
181 burst_surviving_mass_limits = None
182 if sfh_population_model == TORTORELLI25_SNORM_TRUNC:
183 self.ctx.sfh_time_grid = np.asarray(par.time_grid, dtype=np.float64)
184 else:
185 (
186 ssp_age_yr,
187 ssp_age_bins_yr,
188 ssp_age_weights_yr,
189 ) = load_ssp_time_grid(par.ssp_library_filepath)
190 self.ctx.sfh_time_grid = ssp_age_yr
191 emulator_version = getattr(
192 par, "surviv_stellar_mass_emulator_version", "1.0.0"
193 )
194 (
195 smooth_surviving_mass_emulator,
196 smooth_surviving_mass_limits,
197 ) = load_surviving_stellar_mass_model(
198 par.surviv_stellar_mass_emulator_smooth_directory,
199 SMOOTH_SURVIVING_MASS_MODEL,
200 emulator_version=emulator_version,
201 )
202 if sfh_population_model == CONDITIONAL_SNORM_BURST_TRUNC:
203 (
204 burst_surviving_mass_emulator,
205 burst_surviving_mass_limits,
206 ) = load_surviving_stellar_mass_model(
207 par.surviv_stellar_mass_emulator_burst_directory,
208 BURST_SURVIVING_MASS_MODEL,
209 emulator_version=emulator_version,
210 )
212 # Extinction
213 exctinction_eval = ExtinctionMapEvaluator(par)
215 # Initialize galaxy catalog
216 self.ctx.galaxies = galaxy_sampling_util.Catalog()
217 self.ctx.galaxies.columns = [
218 "id",
219 "z",
220 "log10_stellar_mass",
221 "log10_surviv_stellar_mass",
222 "star_formation_history",
223 "star_formation_rate",
224 "star_formation_rate_10_myr",
225 "mSFR",
226 "mpeak",
227 "mperiod",
228 "mskew",
229 "has_burst",
230 "fburst",
231 "gas_metallicity",
232 "metallicity_history",
233 "gas_ionization",
234 "dust_attenuation",
235 "dust_emission",
236 "int_r50",
237 "r50_phys",
238 "int_e1",
239 "int_e2",
240 "sersic_n",
241 "sigma_vel",
242 "fagn",
243 # "mag",
244 "galaxy_type",
245 "excess_b_v",
246 ]
248 # Columns modified inside loop
249 loop_cols = [
250 "z",
251 "log10_stellar_mass",
252 "log10_surviv_stellar_mass",
253 "star_formation_history",
254 "star_formation_rate",
255 "star_formation_rate_10_myr",
256 "mSFR",
257 "mpeak",
258 "mperiod",
259 "mskew",
260 "has_burst",
261 "fburst",
262 "gas_metallicity",
263 "metallicity_history",
264 "gas_ionization",
265 "dust_attenuation",
266 "dust_emission",
267 "int_r50",
268 "r50_phys",
269 "int_e1",
270 "int_e2",
271 "sersic_n",
272 "sigma_vel",
273 "fagn",
274 # "mag",
275 "galaxy_type",
276 "excess_b_v",
277 ]
278 retain_sfh_histories = getattr(par, "retain_sfh_histories", True)
279 if not retain_sfh_histories:
280 for history_column in ("star_formation_history", "metallicity_history"):
281 self.ctx.galaxies.columns.remove(history_column)
282 loop_cols.remove(history_column)
283 if w is not None:
284 loop_cols += ["x", "y"]
285 self.ctx.galaxies.columns += ["x", "y"]
286 else:
287 loop_cols += ["ra", "dec"]
288 self.ctx.galaxies.columns += ["ra", "dec"]
289 for c in loop_cols:
290 setattr(self.ctx.galaxies, c, [])
292 self.ctx.galaxies.int_magnitude_dict = dict()
293 self.ctx.galaxies.abs_magnitude_dict = dict()
294 for key in par.filters:
295 self.ctx.galaxies.int_magnitude_dict[key] = []
296 self.ctx.galaxies.abs_magnitude_dict[key] = []
297 # loop over stellar mass functions
298 for g in par.galaxy_types:
299 n_gal_type = 0
300 n_gal_type_sampled = 0
301 n_gal_type_max = getattr(par, f"n_gal_max_{g}")
302 max_reached = False
304 for i in LOGGER.progressbar(
305 range(len(self.ctx.pixels)),
306 desc=f"getting {g:<4s} galaxies for healpix pixels",
307 at_level="info",
308 ):
309 # Sample stellar masses and redshifts from double stellar mass function
310 sm_lowmass, z_lowmass = sm_funcs_lowmass[g].sample_z_sm_and_apply_cut(
311 seed_ngal=stellar_mass_function_seed(
312 par,
313 pixel=self.ctx.pixels[i],
314 galaxy_type=g,
315 mass_component="lowmass",
316 ),
317 n_gal_max=n_gal_type_max,
318 )
320 sm_highmass, z_highmass = sm_funcs_highmass[
321 g
322 ].sample_z_sm_and_apply_cut(
323 seed_ngal=stellar_mass_function_seed(
324 par,
325 pixel=self.ctx.pixels[i],
326 galaxy_type=g,
327 mass_component="highmass",
328 ),
329 n_gal_max=n_gal_type_max,
330 )
332 log10_stellar_masses = np.concatenate(
333 (sm_lowmass, sm_highmass)
334 ) # log10(Mass/Msun)
335 z = np.concatenate((z_lowmass, z_highmass))
337 ngal = len(log10_stellar_masses)
338 n_gal_type_sampled += ngal
340 # Check the number of sampled galaxies before generating their
341 # more expensive physical properties.
342 n_gal_type_limit = n_gal_type_max * par.ngal_multiplier
343 if n_gal_type_sampled > n_gal_type_limit:
344 max_reached = True
345 if par.raise_max_num_gal_error:
346 raise galaxy_sampling_util.UCatNumGalError(
347 "exceeded number of "
348 f"{g} galaxies {n_gal_type_sampled}>{n_gal_type_limit}"
349 )
351 # Star formation history
352 magemaxs = par.magemax_z0 - (
353 get_lookbacktime_from_redshift(z, cosmo) / 1.0e9
354 )
355 if sfh_population_model == TORTORELLI25_SNORM_TRUNC:
356 (
357 sfhs_msun_yr,
358 sfrs_msun_yr,
359 mSFR,
360 mpeak,
361 mperiod,
362 mskew,
363 ) = sample_sfh_snorm_trunc(
364 log10_stellar_masses,
365 z,
366 par,
367 cosmo,
368 galaxy_type=g,
369 seed=galaxy_property_seed(
370 par,
371 self.ctx.pixels[i],
372 g,
373 "star_formation_history",
374 ),
375 )
376 star_formation_rates_10_myr = np.full(ngal, np.nan)
377 mburst = np.zeros(ngal, dtype=np.float64)
378 has_burst = np.zeros(ngal, dtype=bool)
379 fburst = np.zeros(ngal, dtype=np.float64)
380 else:
381 include_recent_bursts = (
382 sfh_population_model == CONDITIONAL_SNORM_BURST_TRUNC
383 )
384 (
385 mpeak,
386 mperiod,
387 mskew,
388 has_burst,
389 fburst,
390 ) = sample_conditional_sfh_parameters(
391 log10_stellar_masses,
392 z,
393 magemaxs,
394 par,
395 galaxy_type=g,
396 include_recent_bursts=include_recent_bursts,
397 seed=galaxy_property_seed(
398 par,
399 self.ctx.pixels[i],
400 g,
401 "star_formation_history",
402 ),
403 )
404 (
405 sfhs_msun_yr,
406 sfrs_msun_yr,
407 star_formation_rates_10_myr,
408 mSFR,
409 mburst,
410 ) = evaluate_sfh_on_ssp_grid(
411 log10_stellar_masses,
412 mpeak,
413 mperiod,
414 mskew,
415 magemaxs,
416 fburst,
417 ssp_age_yr,
418 ssp_age_weights_yr,
419 mtrunc=par.mtrunc,
420 chunk_size=getattr(
421 par, "sfh_population_evaluation_chunk_size", 10_000
422 ),
423 retain_histories=retain_sfh_histories,
424 )
426 # gas metallicity history
427 Zgas_final = sample_gas_metallicity(
428 log10_stellar_masses,
429 z,
430 sfrs_msun_yr,
431 par,
432 cosmo,
433 galaxy_type=g,
434 seed=galaxy_property_seed(
435 par, self.ctx.pixels[i], g, "gas_metallicity"
436 ),
437 )
438 if not retain_sfh_histories:
439 Z_gas = None
440 elif sfh_population_model == TORTORELLI25_SNORM_TRUNC:
441 Z_gas = generate_gas_metallicity_history_snorm_trunc(
442 Zgas_final,
443 sfhs_msun_yr,
444 par,
445 Zgas_init=par.Zgas_init,
446 )
447 else:
448 Z_gas = generate_gas_metallicity_history_on_ssp_grid(
449 Zgas_final,
450 sfhs_msun_yr,
451 magemaxs,
452 ssp_age_yr,
453 ssp_age_bins_yr,
454 ssp_age_weights_yr,
455 Zgas_init=par.Zgas_init,
456 chunk_size=getattr(
457 par, "sfh_population_evaluation_chunk_size", 10_000
458 ),
459 )
461 # surviving stellar mass
462 if sfh_population_model == TORTORELLI25_SNORM_TRUNC:
463 log10_surviv_stellar_masses = compute_legacy_surviving_stellar_mass(
464 z,
465 log10_stellar_masses,
466 mSFR,
467 mpeak,
468 mperiod,
469 mskew,
470 Zgas_final,
471 par,
472 )
473 else:
474 log10_surviv_stellar_masses = predict_log10_surviving_stellar_mass(
475 log10_stellar_masses,
476 z,
477 magemaxs,
478 mpeak,
479 mperiod,
480 mskew,
481 Zgas_final,
482 has_burst,
483 fburst,
484 smooth_surviving_mass_emulator,
485 smooth_surviving_mass_limits,
486 burst_surviving_mass_emulator,
487 burst_surviving_mass_limits,
488 )
490 # gas ionization
491 logU = sample_gas_ionization_Kashino19(
492 Zgas_final,
493 sfrs_msun_yr,
494 log10_surviv_stellar_masses,
495 par,
496 seed=galaxy_property_seed(
497 par, self.ctx.pixels[i], g, "gas_ionization"
498 ),
499 )
501 # dust attenuation (internal to the observed galaxy)
502 (
503 dust_stellar_masses,
504 log_specific_star_formation_rate,
505 ) = stellar_mass_and_ssfr_for_dust_attenuation(
506 par,
507 log10_stellar_masses,
508 log10_surviv_stellar_masses,
509 sfrs_msun_yr,
510 )
511 dust_attenuation_params = sample_dust_attenuation_parameters(
512 par,
513 dust_stellar_masses,
514 log_specific_star_formation_rate,
515 z,
516 galaxy_type=g,
517 seed=galaxy_property_seed(
518 par, self.ctx.pixels[i], g, "dust_attenuation"
519 ),
520 )
521 # dust emission
522 dust_emission_params = sample_Dale2014_dust_emission_parameters(
523 par,
524 len(log10_stellar_masses),
525 galaxy_type=g,
526 seed=galaxy_property_seed(
527 par, self.ctx.pixels[i], g, "dust_emission"
528 ),
529 )
531 # intrinsic size z
532 int_r50, int_r50_arcsec, r50_phys = sample_r50_from_stellar_mass(
533 z,
534 log10_surviv_stellar_masses,
535 cosmo,
536 par,
537 galaxy_type=g,
538 seed=galaxy_property_seed(par, self.ctx.pixels[i], g, "size"),
539 )
541 # intrinsic ellipticity
542 np.random.seed(
543 galaxy_property_seed(par, self.ctx.pixels[i], g, "ellipticity")
544 )
545 (
546 int_e1,
547 int_e2,
548 ) = sample_ellipticities_for_galaxy_type(
549 n_gal=ngal, galaxy_type=g, par=par
550 )
552 # sersic index
553 np.random.seed(
554 galaxy_property_seed(par, self.ctx.pixels[i], g, "sersic_index")
555 )
556 sersic_n = sample_sersic_for_galaxy_type(
557 n_gal=ngal,
558 galaxy_type=g,
559 par=par,
560 log10_surviv_stellar_mass=log10_surviv_stellar_masses,
561 )
563 # assign velocity dispersion
564 vel_disp = sample_velocity_dispersion_zahid(
565 par,
566 log10_surviv_stellar_masses,
567 seed=galaxy_property_seed(
568 par, self.ctx.pixels[i], g, "velocity_dispersion"
569 ),
570 )
572 # sample agn luminosity
573 fagn = sample_ratio_agn_to_galaxy_bolometric_luminosity(
574 par,
575 z,
576 log10_stellar_masses,
577 galaxy_type=g,
578 seed=galaxy_property_seed(par, self.ctx.pixels[i], g, "agn"),
579 )
581 # Positions
582 np.random.seed(
583 galaxy_property_seed(par, self.ctx.pixels[i], g, "position")
584 )
585 # x and y for wcs, ra and dec for healpix model
586 x, y = sample_position_uniform(
587 len(z), w, self.ctx.pixels[i], par.nside_sampling
588 )
589 # Make the catalog precision already here to avoid inconsistencies in
590 # the selection
591 x = x.astype(par.catalog_precision)
592 y = y.astype(par.catalog_precision)
594 # apply MW Extinction
595 excess_b_v = exctinction_eval(w, x, y)
596 if len(z) == 1:
597 excess_b_v = np.array([excess_b_v])
599 # select galaxies in image boundaries
600 if w is not None:
601 select = in_pos(x, y, par)
602 else:
603 select = np.ones_like(x, dtype=bool)
604 n_gal = np.count_nonzero(select)
605 n_gal_type += n_gal
607 # store
608 if w is not None:
609 self.ctx.galaxies.x.append(x[select].astype(par.catalog_precision))
610 self.ctx.galaxies.y.append(y[select].astype(par.catalog_precision))
611 else:
612 self.ctx.galaxies.ra.append(x[select].astype(par.catalog_precision))
613 self.ctx.galaxies.dec.append(
614 y[select].astype(par.catalog_precision)
615 )
617 self.ctx.galaxies.z.append(z[select].astype(par.catalog_precision))
619 self.ctx.galaxies.galaxy_type.append(
620 np.ones(n_gal, dtype=np.ushort) * sm_funcs_lowmass[g].galaxy_type
621 )
622 self.ctx.galaxies.excess_b_v.append(
623 excess_b_v[select].astype(par.catalog_precision)
624 )
625 self.ctx.galaxies.log10_stellar_mass.append(
626 log10_stellar_masses[select].astype(par.catalog_precision)
627 )
628 if retain_sfh_histories:
629 self.ctx.galaxies.star_formation_history.append(
630 np.array(sfhs_msun_yr)[select]
631 )
632 self.ctx.galaxies.star_formation_rate.append(sfrs_msun_yr[select])
633 self.ctx.galaxies.star_formation_rate_10_myr.append(
634 star_formation_rates_10_myr[select]
635 )
636 self.ctx.galaxies.mSFR.append(mSFR[select])
637 self.ctx.galaxies.mpeak.append(mpeak[select])
638 self.ctx.galaxies.mperiod.append(mperiod[select])
639 self.ctx.galaxies.mskew.append(mskew[select])
640 self.ctx.galaxies.has_burst.append(has_burst[select])
641 self.ctx.galaxies.fburst.append(fburst[select])
642 self.ctx.galaxies.gas_metallicity.append(
643 np.array(Zgas_final)[select].astype(par.catalog_precision)
644 )
645 if retain_sfh_histories:
646 self.ctx.galaxies.metallicity_history.append(Z_gas[select])
647 self.ctx.galaxies.log10_surviv_stellar_mass.append(
648 log10_surviv_stellar_masses[select]
649 )
650 self.ctx.galaxies.gas_ionization.append(
651 np.array(logU)[select].astype(par.catalog_precision)
652 )
653 self.ctx.galaxies.dust_attenuation.append(
654 np.array(dust_attenuation_params)[select].astype(
655 par.catalog_precision
656 )
657 )
658 self.ctx.galaxies.dust_emission.append(
659 np.array(dust_emission_params)[select].astype(par.catalog_precision)
660 )
661 self.ctx.galaxies.int_r50.append(
662 np.array(int_r50)[select].astype(par.catalog_precision)
663 )
664 self.ctx.galaxies.r50_phys.append(
665 np.array(r50_phys)[select].astype(par.catalog_precision)
666 )
667 self.ctx.galaxies.int_e1.append(
668 np.array(int_e1)[select].astype(par.catalog_precision)
669 )
670 self.ctx.galaxies.int_e2.append(
671 np.array(int_e2)[select].astype(par.catalog_precision)
672 )
673 self.ctx.galaxies.sersic_n.append(
674 np.array(sersic_n)[select].astype(par.catalog_precision)
675 )
676 self.ctx.galaxies.sigma_vel.append(
677 np.array(vel_disp)[select].astype(par.catalog_precision)
678 )
679 self.ctx.galaxies.fagn.append(
680 np.array(fagn)[select].astype(par.catalog_precision)
681 )
683 # check memory footprint
684 check_max_mem_error()
686 if max_reached:
687 break
689 LOGGER.info(
690 f"smfun={g} n_gals={n_gal_type} maximum number of galaxies "
691 f"reached={max_reached} ({n_gal_type_max})"
692 )
694 # Concatenate columns
695 for c in loop_cols:
696 setattr(self.ctx.galaxies, c, np.concatenate(getattr(self.ctx.galaxies, c)))
698 LOGGER.info(
699 f"Keep sources with {par.filters_full_names_prospect_selection} "
700 f"below {par.promage_selection_mag_threshold}"
701 )
703 selected_objects = apply_galaxy_observed_magnitude_cut(par, self.ctx.galaxies)
704 for c in loop_cols:
705 setattr(
706 self.ctx.galaxies,
707 c,
708 getattr(self.ctx.galaxies, c)[selected_objects],
709 )
711 LOGGER.info(f"Computing magnitudes with {par.sed_generator}")
712 output_sed_generator = sps_sed_generator(
713 par,
714 z=self.ctx.galaxies.z,
715 mSFR=self.ctx.galaxies.mSFR,
716 mpeak=self.ctx.galaxies.mpeak,
717 mperiod=self.ctx.galaxies.mperiod,
718 mskew=self.ctx.galaxies.mskew,
719 Zgas_final=self.ctx.galaxies.gas_metallicity,
720 logU=self.ctx.galaxies.gas_ionization,
721 fagn=self.ctx.galaxies.fagn,
722 dust_attenuation_params=self.ctx.galaxies.dust_attenuation,
723 excess_b_v=self.ctx.galaxies.excess_b_v,
724 dust_emission_params=self.ctx.galaxies.dust_emission,
725 vel_disp=self.ctx.galaxies.sigma_vel,
726 ssp_library_filepath=par.ssp_library_filepath,
727 filter_names=par.filters_full_names_prospect,
728 cosmology=cosmo,
729 log10_stellar_mass=self.ctx.galaxies.log10_stellar_mass,
730 has_burst=self.ctx.galaxies.has_burst,
731 fburst=self.ctx.galaxies.fburst,
732 )
734 for j, key in enumerate(par.filters):
735 self.ctx.galaxies.int_magnitude_dict[key] = output_sed_generator[1][
736 :, j
737 ] # maybe save mags before and after mw extinction
738 self.ctx.galaxies.abs_magnitude_dict[key] = output_sed_generator[0][:, j]
740 # Set apparent (lensed) magnitudes, for now equal to intrinsic magnitudes
741 self.ctx.galaxies.magnitude_dict = dict()
742 for band, mag in self.ctx.galaxies.int_magnitude_dict.items():
743 self.ctx.galaxies.magnitude_dict[band] = mag.copy()
745 # Number of galaxies and id
746 self.ctx.numgalaxies = self.ctx.galaxies.z.size
747 self.ctx.galaxies.id = np.arange(self.ctx.numgalaxies)
749 # Backward compatibility
750 self.ctx.galaxies.blue_red = np.ones(len(self.ctx.galaxies.z), dtype=np.ushort)
751 self.ctx.galaxies.blue_red[
752 self.ctx.galaxies.galaxy_type == sm_funcs_lowmass["blue"].galaxy_type
753 ] = 1
754 self.ctx.galaxies.blue_red[
755 self.ctx.galaxies.galaxy_type == sm_funcs_lowmass["red"].galaxy_type
756 ] = 0
758 if par.save_SEDs & (par.sed_generator.lower() == "prospect"):
759 # Store SEDs in the catalog
760 self.ctx.restframe_wavelength_for_SED = output_sed_generator[2]
761 self.ctx.galaxies.sed = output_sed_generator[3]
763 LOGGER.info(
764 f"galaxy counts n_total={self.ctx.numgalaxies} "
765 f"mem_mb_current={utils.memory_usage_psutil():5.1f}"
766 )
768 def __str__(self):
769 return "sample gal photo prospect"