add magnitudes in flux ratio plot

fig/flux_sigma_eta_ratio_20250930.pdf: src/scripts/plots/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 | tee table/table_flux_sigma_eta_ratio_20250930.tex

fig/flux_sigma_eta_ratio_magnitudes_20251006.pdf: src/scripts/plots/plot_flux_sigma_eta_ratio_magnitudes.py

	python3 $< $@ 0 | tee table/table_flux_sigma_eta_ratio_magnitudes_20251006.tex

fig/flux_sigma_eta_ratio_u.pdf: python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 u > table/flux_sigma_eta_u.table

fig/flux_sigma_eta_ratio_g.pdf: python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 g > table/flux_sigma_eta_g.table

fig/flux_sigma_eta_ratio_r.pdf: python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 r > table/flux_sigma_eta_r.table

fig/flux_sigma_eta_ratio_i.pdf: python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 i > table/flux_sigma_eta_i.table

fig/flux_sigma_eta_ratio_z.pdf: python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 z > table/flux_sigma_eta_z.table

fig/flux_sigma_eta_ratio_y.pdf: python/plot_flux_sigma_eta_ratio.py

	python3 $< $@ 0 y > table/flux_sigma_eta_y.table

import os

import sys

import time

import matplotlib.pyplot as plt

import numpy as np

from astropy.table import Table

from lsst_inaf_agile import util

savefig = sys.argv[1]

select_morpho = bool(int(sys.argv[2]))

band = "r"

if len(sys.argv[1:]) > 2:

    band = sys.argv[-1]

# flux_min = 1e-1 * 1000.0

# flux_max = 1e+2 * 1000.0

# NOTE: in uJy

# flux_min = 0.100    # 26.4

# flux_max = 10.00    # 21.4

flux_min = -np.inf

flux_max = +np.inf

# Bins in magnitude space to estimate the statistics in

mag_bins = [

    (20, 22),

    (22, 24),

    (24, 26),

]

if False:

    table = Table.read(

        os.environ.get("PATH_AMUSE", "")

        +

        # "/precario3/viitanen/agile/data/dr1/truth_match_objectTable.csv"

        "/staff2/viitanen/agile/match_detect_isPrimary.csv"

    )

try:

    # Create the match table

    truth_table = Table.read("/staff2/viitanen/agile/data/dr1/catalog.fits")

    object_table = Table.read(

        "/precario3/viitanen/agile/data/dr1/repo_1det/u/viitanen/output/20241211T013547Z/objectTable_tract/9813/objectTable_tract_9813_DC2_u_viitanen_output_20241211T013547Z.csv"

    )

    object_table = object_table[object_table["detect_isPrimary"] == "True"]

    match_table = Table.read(

        "/precario3/viitanen/agile/data/dr1/repo_1det/u/viitanen/output/20241212T094317Z/match_target_truth_summary_objectTable_tract/9813/match_target_truth_summary_objectTable_tract_9813_DC2_u_viitanen_output_20241212T094317Z.fits"

    )

    match_table = match_table[match_table["match_candidate"]]

except FileNotFoundError:

    import sqlite3

    import pandas as pd

    from astropy.table import Table

    with sqlite3.connect("data/catalog/dr1_new_new/db/20250725/master.db") as con:

        select_list = []

        for name, column in [

            (f"lsst-{band}_total_truth", f"lsst-{band}_total"),

            (f"lsst-{band}_point_truth", f"lsst-{band}_point"),

            (f"{band}_psfFlux_object", f"{band}_psfFlux"),

            (f"{band}_calibFlux_object", f"{band}_calibFlux"),

            (f"{band}_cModelFlux_object", f"{band}_cModelFlux"),

            (f"{band}_psfFluxErr_object", f"{band}_psfFluxErr"),

            (f"{band}_calibFluxErr_object", f"{band}_calibFluxErr"),

            (f"{band}_cModelFluxErr_object", f"{band}_cModelFluxErr"),

            (f"{band}_psfFlux_flag_object", f"{band}_psfFlux_flag"),

            (f"{band}_calibFlux_flag_object", f"{band}_calibFlux_flag"),

            (f"{band}_cModel_flag_object", f"{band}_cModel_flag"),

            (f"{band}_extendedness_object", f"{band}_extendedness"),

            ("detect_isPrimary_object", "detect_isPrimary"),

            ("match_truth_type_match", "match_truth_type"),

        ]:

            select_list.append(f'"{column}" AS "{name}"')

        select = ", ".join(select_list[:14])

        t0 = time.time()

        df = pd.read_sql_query(

            f"""

            SELECT

                {select}

            FROM Object AS o

                JOIN MatchesTruth AS m USING (objectId)

                JOIN Truth AS t ON (t.ID = m.match_id)

            WHERE

                1=o.detect_isPrimary AND

                1=m.match_candidate AND

                m.match_truth_type>=0

            """,

            con,

        )

        table = Table.from_pandas(df)

        t1 = time.time()

        print(f"Table loaded in {t1 - t0} seconds", file=sys.stderr)

if False:

    # check that the two arrays are sorted

    assert np.all(truth_table["ID"][:-1] <= truth_table["ID"][1:])

    assert np.all(object_table["objectId"][:-1] <= object_table["objectId"][1:])

    idx1 = np.searchsorted(truth_table["ID"], match_table["match_id"])

    idx2 = np.searchsorted(object_table["objectId"], match_table["objectId"])

    table = Table(

        {

            f"lsst-{band}_total_truth": truth_table[f"lsst-{band}_total"][idx1],

            f"lsst-{band}_point_truth": truth_table[f"lsst-{band}_point"][idx1],

            f"{band}_psfFlux_object": object_table[f"{band}_psfFlux"][idx2],

            f"{band}_calibFlux_object": object_table[f"{band}_calibFlux"][idx2],

            f"{band}_cModelFlux_object": object_table[f"{band}_cModelFlux"][idx2],

            f"{band}_psfFluxErr_object": object_table[f"{band}_psfFluxErr"][idx2],

            f"{band}_calibFluxErr_object": object_table[f"{band}_calibFluxErr"][idx2],

            f"{band}_cModelFluxErr_object": object_table[f"{band}_cModelFluxErr"][idx2],

            f"{band}_psfFlux_flag_object": object_table[f"{band}_psfFlux_flag"][idx2],

            f"{band}_calibFlux_flag_object": object_table[f"{band}_calibFlux_flag"][idx2],

            f"{band}_cModel_flag_object": object_table[f"{band}_cModel_flag"][idx2],

            f"{band}_extendedness_object": object_table[f"{band}_extendedness"][idx2],

            "detect_isPrimary_object": object_table["detect_isPrimary"][idx2],

            "match_truth_type_match": match_table["match_truth_type"],

        }

    )

    table = table[table["detect_isPrimary_object"] == "True"]

    select = np.ones_like(table, dtype=bool)

    for col in table.columns:

        if hasattr(table[col], "mask"):

            select &= ~table[col]

    table = table[select]

def latex_format_number(number):

    """

    Format number in latex style with \\, separating thousands

    """

    ret = f"${number:,}$"

    return ret.replace(",", "\\,")

def do_one(fig, axes, row, col, name_truth, name_flux, is_extended, plot=True):

    # Set the color

    color = f"C{col}"

    uid = None

    if "AGN1" in name_truth:

        uid = 0

    if "AGN2" in name_truth:

        uid = 1

    if "Galaxy" in name_truth:

        uid = 2

    if "Star" in name_truth:

        uid = 3

    # Defines the fluxes, note that we convert nJy from objectTable to uJy

    flux1 = table[f"lsst-{band}_total_truth"]

    flux2 = table[f"{band}_{name_flux}Flux_object"] / 1000.0

    if name_flux == "cModel":

        select_flag = table[f"{band}_{name_flux}_flag_object"] == 0

    else:

        select_flag = table[f"{band}_{name_flux}Flux_flag_object"] == 0

    # Defines the selections

    select_flux = (flux_min <= flux1) * (flux1 < flux_max) * (flux_min <= flux2) * (flux2 < flux_max)

    # Select extended or non-extended objects

    select_extended = (table[f"{band}_extendedness_object"] == 0) | (

        table[f"{band}_extendedness_object"] == 1

    )

    if uid is not None:

        select_uid = table["match_truth_type_match"] == uid

        select = select_flag & select_uid & select_extended

        select_box = select & select_flux

    if select_morpho:

        select_extended &= table[f"{band}_extendedness_object"] == is_extended

        select = select_flag & select_uid & select_extended

        select_box = select_uid & select_flux & select_extended

    # Select pointlike sources

    if "Pointlike" in name_truth:

        select_extended = table[f"{band}_extendedness_object"] == 0

        select = select_flag & select_extended

        select_box = select & select_flux

    if "Extended" in name_truth:

        select_extended = table[f"{band}_extendedness_object"] == 1

        select = select_flag & select_extended

        select_box = select & select_flux

    # AGN with >90%

    if "AGN1 (>90%)" in name_truth:

        flux_point = table[f"lsst-{band}_point_truth"]

        flux_total = table[f"lsst-{band}_total_truth"]

        select_gt_90 = flux_point / flux_total > 0.90

        select &= select_gt_90

        select_box &= select_gt_90

    if "AGN1 (<90%)" in name_truth:

        flux_point = table[f"lsst-{band}_point_truth"]

        flux_total = table[f"lsst-{band}_total_truth"]

        select_gt_90 = flux_point / flux_total < 0.90

        select &= select_gt_90

        select_box &= select_gt_90

    # Estimates the statistics in bins of magnitudes

    ns = []

    sigmas = []

    etas = []

    for mag_min, mag_max in mag_bins:

        mag = util.flux_to_mag(flux1)

        select_box_mag = select_box & (mag_min <= mag) & (mag < mag_max)

        sigma = util.get_sigma_nmad(flux2[select_box_mag], flux1[select_box_mag])

        eta = util.get_fraction_catastrophic_error(flux2[select_box_mag], flux1[select_box_mag])

        ns.append(select_box_mag.sum())

        sigmas.append(sigma)

        etas.append(eta)

    value1 = (name_truth * (row == 0)).replace("%", "\\%")

    value2 = name_flux + "Flux"

    items = [f"{value1:12s}", f"{value2:10s}"]

    for n, sigma, eta in zip(ns, sigmas, etas, strict=False):

        items += [

            f"{latex_format_number(n):10s}",

            f"{sigma:7.3f}",

            f"{eta:7.3f}",

        ]

    print(" & ".join(items), r"\\")

    if not plot:

        return

    # Plot the binned_statistic

    from scipy.stats import binned_statistic

    bins = np.logspace(-2, 3, 16)

    cens = (bins[:-1] * bins[1:]) ** 0.5

    Ntot = binned_statistic(flux1[select], select[select], bins=bins, statistic="sum")[0]

    def get_quantile(q):

        """Return a function that counts the quantile 'q'"""

        def fun(a):

            return np.quantile(a, q)

        return fun

    qs = []

    for q in 0.10, 0.50, 0.90:

        qs.append(

            binned_statistic(

                flux1[select], flux2[select] / flux1[select], bins=bins, statistic=get_quantile(q)

            )[0]

        )

    q10, q50, q90 = qs

    select_gt_5 = Ntot >= 5

    cens_mag = util.flux_to_mag(cens)

    q10_mag = -2.5 * np.log10(q10)

    q50_mag = -2.5 * np.log10(q50)

    q90_mag = -2.5 * np.log10(q90)

    axes[row, col].plot(cens_mag[select_gt_5], q10_mag[select_gt_5], color=color, linestyle="dotted")

    axes[row, col].plot(cens_mag[select_gt_5], q50_mag[select_gt_5], color=color, linestyle="solid")

    axes[row, col].plot(cens_mag[select_gt_5], q90_mag[select_gt_5], color=color, linestyle="dotted")

    # NOTE: magnitude limits are 27.15, 19.65

    # axes[row, col].semilogx()

    # axes[row, col].set_xlim(5e-2, 5e1)

    # axes[row, col].set_xticks([1e-1, 1e0, 1e1])

    axes[row, col].set_xlim(19.6, 27.2)

    axes[row, col].set_xticks([20, 22, 24, 26])

    axes[row, col].set_xlim(axes[row, col].get_xlim()[::-1])

    # NOTE: deltamag limits are -0.44, +0.75

    # axes[row, col].set_ylim(0.50, 1.50)

    # axes[row, col].set_yticks([0.60, 0.80, 1.0, 1.2, 1.4])

    # axes[row, col].set_ylim(-0.50, +0.80)

    axes[row, col].set_ylim(-0.60, +0.60)

    axes[row, col].set_yticks([-0.50, 0.00, +0.50])

    axes[row, col].set_ylim(axes[row, col].get_ylim()[::-1])

    # axes[row, col].axhline(1.0, linestyle="dotted")

    axes[row, col].axhline(0.0, linestyle="dotted")

    if row == 0:

        axes[row, col].set_title(name_truth, fontsize=12)

    # fig.supxlabel("flux truth [uJy]", va="bottom")

    fig.supxlabel(r"$m_\mathrm{truth}$ [ABmag]", va="bottom")

    if col == 0:

        # axes[row, col].set_ylabel(f"{name_flux}Flux / truth")

        mag_str = f"$\\Delta m_{{{name_flux}}}$"

        axes[row, col].set_ylabel(f"{mag_str}")

fig, axes = plt.subplots(

    3,

    5,

    figsize=(8.27, 11.69 / 2),

    sharex=True,

    sharey=True,

    gridspec_kw=dict(hspace=0.03, wspace=0.03),

    dpi=300,

)

# Plot

do_one(fig, axes, 0, 0, "AGN1 (<90%)", "psf", 0, plot=True)

do_one(fig, axes, 1, 0, "AGN1 (<90%)", "calib", 0, plot=True)

do_one(fig, axes, 2, 0, "AGN1 (<90%)", "cModel", 0, plot=True)

do_one(fig, axes, 0, 1, "AGN1 (>90%)", "psf", 0, plot=True)

do_one(fig, axes, 1, 1, "AGN1 (>90%)", "calib", 0, plot=True)

do_one(fig, axes, 2, 1, "AGN1 (>90%)", "cModel", 0, plot=True)

do_one(fig, axes, 0, 2, "AGN2", "psf", 0, plot=True)

do_one(fig, axes, 1, 2, "AGN2", "calib", 0, plot=True)

do_one(fig, axes, 2, 2, "AGN2", "cModel", 0, plot=True)

do_one(fig, axes, 0, 3, "Galaxy", "psf", 0, plot=True)

do_one(fig, axes, 1, 3, "Galaxy", "calib", 0, plot=True)

do_one(fig, axes, 2, 3, "Galaxy", "cModel", 0, plot=True)

do_one(fig, axes, 0, 4, "Star", "psf", 0, plot=True)

do_one(fig, axes, 1, 4, "Star", "calib", 0, plot=True)

do_one(fig, axes, 2, 4, "Star", "cModel", 0, plot=True)

do_one(fig, axes, 0, 4, "Pointlike", "psf", 0, plot=False)

do_one(fig, axes, 1, 4, "Pointlike", "calib", 0, plot=False)

do_one(fig, axes, 2, 4, "Pointlike", "cModel", 0, plot=False)

do_one(fig, axes, 0, 4, "Extended", "psf", 0, plot=False)

do_one(fig, axes, 1, 4, "Extended", "calib", 0, plot=False)

do_one(fig, axes, 2, 4, "Extended", "cModel", 0, plot=False)

fig.savefig(savefig, bbox_inches="tight")