Add binned stellar mass - BH mass plot

colibre/auto_plotter/black_holes.yml

@@ -67,6 +67,39 @@ stellar_mass_black_hole_mass_50:
     - filename: GalaxyStellarMassBlackHoleMass/Graham2023_S.hdf5
     - filename: GalaxyStellarMassBlackHoleMass/McConnell2013_Fit.hdf5
 
+stellar_mass_black_hole_mass_binned_50:
+  type: "scatter"
+  legend_loc: "upper left"
+  x:
+    quantity: "apertures.mass_star_50_kpc"
+    units: Solar_Mass
+    start: 1e6
+    end: 1e12
+  y:
+    quantity: "black_hole_masses.max"
+    units: Solar_Mass
+    start: 1e3
+    end: 1e11
+  median:
+    plot: true
+    log: true
+    adaptive: true
+    number_of_bins: 30
+    start:
+      value: 1e6
+      units: Solar_Mass
+    end:
+      value: 1e12
+      units: Solar_Mass
+  metadata:
+    title: Stellar Mass-Black Hole Mass relation (50 kpc Stellar Mass)
+    caption: SMBHM relation. The binned observational data are obtained by combining BH mass measurements for galaxies of different morphologies (Graham 2023) and passive vs. active galaxies (Terrazas et al. 2017) with dependencies on stellar mass of morphological types (Moffett et al. 2016, from the GAMA survey) and quenched fractions (Gilbank et al. 2010, SDSS).
+    section: Black Holes
+  observational_data:
+    - filename: GalaxyStellarMassBlackHoleMass/Terrazas2017_binned.hdf5
+    - filename: GalaxyStellarMassBlackHoleMass/Graham2023_binned.hdf5
+    - filename: GalaxyStellarMassBlackHoleMass/McConnell2013_Fit.hdf5
+
 stellar_mass_black_hole_mass_passive_30:
   type: "scatter"
   select_structure_type: 10

colibre/scripts/bh_accretion_evolution.py

@@ -15,7 +15,7 @@
 
 from swiftpipeline.argumentparser import ScriptArgumentParser
 
-from velociraptor.observations import load_observation
+from velociraptor.observations import load_observations
 
 arguments = ScriptArgumentParser(
     description="Creates a BH accretion rate density evolution plot, with added observational data."
@@ -75,15 +75,201 @@
     )
     simulation_labels.append(name)
 
-# Observational data plotting
+path_to_obs_data = f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}"
 
-observational_data = glob.glob(
-    f"{arguments.config.config_directory}/{arguments.config.observational_data_directory}/data/BlackHoleAccretionHistory/*.hdf5"
+observation_lines = []
+observation_labels = []
+
+Aird2015 = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Aird2015.hdf5"]
+)[0]
+observation_lines.append(
+    ax.plot(
+        Aird2015.x.value,
+        Aird2015.y.value,
+        color="brown",
+        zorder=-10000,
+        alpha=0.6,
+        linewidth=1,
+    )[0]
+)
+observation_labels.append(Aird2015.citation)
+
+Annana2020_low = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Annana2020_low.hdf5"]
+)[0]
+observation_lines.append(
+    ax.plot(
+        Annana2020_low.x.value,
+        Annana2020_low.y.value,
+        label=Annana2020_low.citation,
+        zorder=-10000,
+        color="orange",
+        alpha=0.6,
+        linewidth=1,
+    )[0]
+)
+observation_labels.append(Annana2020_low.citation)
+
+Annana2020_high = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Annana2020_high.hdf5"]
+)[0]
+observation_lines.append(
+    ax.plot(
+        Annana2020_high.x.value,
+        Annana2020_high.y.value,
+        label=Annana2020_high.citation,
+        color="olive",
+        zorder=-10000,
+        alpha=0.6,
+        linewidth=1,
+    )[0]
+)
+observation_labels.append(Annana2020_high.citation)
+
+Pouliasis2024_low = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Pouliasis2024_low.hdf5"]
+)[0]
+observation_lines.append(
+    ax.fill_between(
+        Pouliasis2024_low.x.value,
+        Pouliasis2024_low.y.value - Pouliasis2024_low.y_scatter[0].value,
+        Pouliasis2024_low.y.value + Pouliasis2024_low.y_scatter[1].value,
+        color="magenta",
+        label=Pouliasis2024_low.citation,
+        zorder=-10000,
+        alpha=0.2,
+    )
 )
+observation_labels.append(Pouliasis2024_low.citation)
 
-for index, observation in enumerate(observational_data):
-    obs = load_observation(observation)
-    obs.plot_on_axes(ax)
+Pouliasis2024_high = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Pouliasis2024_high.hdf5"]
+)[0]
+observation_lines.append(
+    ax.fill_between(
+        Pouliasis2024_high.x.value,
+        Pouliasis2024_high.y.value - Pouliasis2024_high.y_scatter[0].value,
+        Pouliasis2024_high.y.value + Pouliasis2024_high.y_scatter[1].value,
+        color="teal",
+        label=Pouliasis2024_high.citation,
+        zorder=-10000,
+        alpha=0.2,
+    )
+)
+observation_labels.append(Pouliasis2024_high.citation)
+
+ax.set_xlabel("Redshift $z$")
+ax.set_ylabel(r"BH Accretion Rate Density [M$_\odot$ yr$^{-1}$ Mpc$^{-3}$]")
+
+Delvacchio2014 = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Delvacchio2014.hdf5"]
+)[0]
+observation_lines.append(
+    ax.errorbar(
+        Delvacchio2014.x.value,
+        Delvacchio2014.y.value,
+        xerr=None
+        if Delvacchio2014.x_scatter is None
+        else Delvacchio2014.x_scatter.value,
+        yerr=None
+        if Delvacchio2014.y_scatter is None
+        else Delvacchio2014.y_scatter.value,
+        label=Delvacchio2014.citation,
+        linestyle="none",
+        marker="o",
+        color="maroon",
+        elinewidth=0.75,
+        markeredgecolor="none",
+        markersize=3,
+        zorder=-10,
+    )[0]
+)
+observation_labels.append(Delvacchio2014.citation)
+
+DSilva2023 = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/DSilva2023.hdf5"]
+)[0]
+observation_lines.append(
+    ax.errorbar(
+        DSilva2023.x.value,
+        DSilva2023.y.value,
+        xerr=None if DSilva2023.x_scatter is None else DSilva2023.x_scatter.value,
+        yerr=None if DSilva2023.y_scatter is None else DSilva2023.y_scatter.value,
+        label=DSilva2023.citation,
+        linestyle="none",
+        marker="o",
+        color="darkblue",
+        elinewidth=0.75,
+        markeredgecolor="none",
+        markersize=3,
+        zorder=-10,
+    )[0]
+)
+observation_labels.append(DSilva2023.citation)
+
+Yang2023 = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Yang2023.hdf5"]
+)[0]
+observation_lines.append(
+    ax.errorbar(
+        Yang2023.x.value,
+        Yang2023.y.value,
+        xerr=None if Yang2023.x_scatter is None else Yang2023.x_scatter.value,
+        yerr=None if Yang2023.y_scatter is None else Yang2023.y_scatter.value,
+        label=Yang2023.citation,
+        linestyle="none",
+        marker="o",
+        color="goldenrod",
+        elinewidth=0.75,
+        markeredgecolor="none",
+        markersize=3,
+        zorder=-10,
+    )
+)
+observation_labels.append(Yang2023.citation)
+
+Shen2020_low = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Shen2020_low.hdf5"]
+)[0]
+observation_lines.append(
+    ax.errorbar(
+        Shen2020_low.x.value,
+        Shen2020_low.y.value,
+        xerr=None if Shen2020_low.x_scatter is None else Shen2020_low.x_scatter.value,
+        yerr=None if Shen2020_low.y_scatter is None else Shen2020_low.y_scatter.value,
+        label=Shen2020_low.citation,
+        linestyle="none",
+        marker="o",
+        color="lime",
+        elinewidth=0.5,
+        markeredgecolor="none",
+        markersize=2,
+        zorder=-10,
+    )
+)
+observation_labels.append(Shen2020_low.citation)
+
+Shen2020_high = load_observations(
+    [f"{path_to_obs_data}/data/BlackHoleAccretionHistory/Shen2020_high.hdf5"]
+)[0]
+observation_lines.append(
+    ax.errorbar(
+        Shen2020_high.x.value,
+        Shen2020_high.y.value,
+        xerr=None if Shen2020_high.x_scatter is None else Shen2020_high.x_scatter.value,
+        yerr=None if Shen2020_high.y_scatter is None else Shen2020_high.y_scatter.value,
+        label=Shen2020_high.citation,
+        linestyle="none",
+        marker="o",
+        color="coral",
+        elinewidth=0.5,
+        markeredgecolor="none",
+        markersize=2,
+        zorder=-10,
+    )
+)
+observation_labels.append(Shen2020_high.citation)
 
 ax.set_xlabel("Redshift $z$")
 ax.set_ylabel(r"BH Accretion Rate Density [M$_\odot$ yr$^{-1}$ Mpc$^{-3}$]")
@@ -109,9 +295,17 @@
 ax.tick_params(axis="x", which="minor", bottom=False)
 
 ax.set_xlim(1.02, 0.07)
-ax.set_ylim(6e-7, 6e-4)
+ax.set_ylim(1e-7, 6e-4)
 
-observation_legend = ax.legend(markerfirst=True, loc="lower left")
+# observation_legend = ax.legend(markerfirst=True, loc="lower left")
+observation_legend = ax.legend(
+    observation_lines,
+    observation_labels,
+    markerfirst=True,
+    loc="lower center",
+    fontsize=4,
+    ncol=1,
+)
 
 ax.add_artist(observation_legend)