Commit 8ff1cd42 authored by Akke Viitanen's avatar Akke Viitanen
Browse files

Refactor MBH

parent 59269ca5

src/python/catalog_agn.py

+2 −15
Original line number Diff line number Diff line
@@ -628,21 +628,8 @@ class CatalogAGN: 
        return self[key] - a * self["E_BV"]



    def get_occupation_fraction(self):

        """

        Returns the black hole occupation fraction 'f_occ' as a function of

        logMstar. The occupation fraction is the percentage of host galaxies

        that are expected to host SMBH, which could be <1 for dwarf galaxies

        (Miller+15, Burke+25, Zou+).



        The updated AGN duty cycle can be estimated using:



            is_agn_new = occupation_fraction(logMstar) * is_agn_old

        """



        table = Table.read("data/from_zou/table2.csv")

        x = table["logm"]

        y = table["median"]

        return np.interp(self["M"], x, y, left=0.0, right=1.0)

        from mbh import get_occupation_fraction

        return get_occupation_fraction(self["M"])



    def get_lightcurve(

        self,

src/python/mbh.py

+31 −1
Original line number Diff line number Diff line
@@ -4,7 +4,11 @@ 
# Email: akke.viitanen@helsinki.fi

# Date: 2023-04-12 22:41:33



"""Tools for assigning MBH"""

"""

Tools for assigning MBH

"""



from astropy.table import Table

import numpy as np

from scipy.interpolate import CubicSpline

from util import ROOT
@@ -284,6 +288,32 @@ def get_log_mbh(log_mstar, reference="Shankar+16", z=0.0): 
    }[reference]







def get_occupation_fraction(logMstar):



    """

    Returns the black hole occupation fraction 'f_occ' as a function of

    logMstar. The occupation fraction is the percentage of host galaxies

    that are expected to host SMBH, which could be <1 for dwarf galaxies

    (Miller+15, Burke+25, Zou+).



    The updated AGN duty cycle can be estimated using:



        is_agn_new = occupation_fraction(logMstar) * is_agn_old

    """



    table = Table.read("data/from_zou/table2.csv")

    x = table["logm"]

    y = table["median"]



    select = (x.min() <= logMstar) & (logMstar <= x.max())

    if not np.all(select):

        raise ValueError(f"logMstar must be within {x.min()} and {x.max()}")



    return np.interp(logMstar, x, y)







if __name__ == "__main__":



    import matplotlib.pyplot as plt

src/scripts/database_count.py

deleted100644 → 0
+0 −76
Original line number Diff line number Diff line 
#!/usr/bin/env python3

# encoding: utf-8

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2025-07-10 19:28:13



"""

Counts stuff from the database

"""



import argparse

import glob

import math

import os

import random

import re

import subprocess

import sys

import time



import astropy as ap

import astropy.coordinates as c

import astropy.units as u

import fitsio

import matplotlib as mpl

import matplotlib.pyplot as plt

import numpy as np

import scipy as sp



import pandas as pd

import sqlite3



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



    df = pd.read_sql_query(

        """

        SELECT detect_isPrimary, match_candidate, match_truth_type, match_id

        FROM Object

        JOIN MatchesTruth USING (objectId)

        """,

        con

    )

    print(len(df))



    select = df["detect_isPrimary"] == 1

    print(select.sum())



    selectn = select & (df["match_truth_type"] <  0)

    selecti = (df["match_truth_type"] >= 0)

    select0 = select & (df["match_truth_type"] == 0)

    select1 = select & (df["match_truth_type"] == 1)

    select2 = select & (df["match_truth_type"] == 2)

    select3 = select & (df["match_truth_type"] == 3)



    print()

    print(selectn.sum())

    print(selecti.sum())

    print(select0.sum())

    print(select1.sum())

    print(select2.sum())

    print(select3.sum())



    with sqlite3.connect("master.db") as con:

        df = pd.read_sql_query(

            """

            SELECT

                objectId, coord_ra, coord_dec, r_psfFlux

            FROM

                Object

            WHERE

                1 = detect_isPrimary AND

                r_psfFlux > 1000

            """,

            con

        )

    print(len(df))

src/scripts/egg_generate_sed.sh

deleted100755 → 0
+0 −20
Original line number Diff line number Diff line 
#!/bin/sh -x

#

# Generate egg seds on the disk

#

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi



filename_egg=$1

dirname_egg=$(dirname $filename_egg)

filename_egg_sed=$dirname_egg/egg-seds.dat



python3 -c """

import fitsio

ids = fitsio.read(\"$filename_egg\", columns='ID')[0]

for i in ids:

    print(str(i))

""" | while read -r id ; do

    egg-getsed seds=$filename_egg_sed id=$id component=bulge out=$dirname_egg/seds/bulge-$id.fits

    egg-getsed seds=$filename_egg_sed id=$id component=disk  out=$dirname_egg/seds/disk-$id.fits

done

src/scripts/plot_baseline_visit.py

deleted100644 → 0
+0 −102
Original line number Diff line number Diff line 
#!/usr/bin/env python3

# encoding: utf-8

# Author: Akke Viitanen

# Email: akke.viitanen@helsinki.fi

# Date: 2025-07-09 23:55:46



"""

Plot LSST baseline visits

"""





from astropy.table import Table

import numpy as np

import pandas as pd

import sqlite3

import matplotlib as mpl

import matplotlib.pyplot as plt



mpl.style.use("av_presentation")



kwargs_hist = dict(

    bins=np.arange(0, 3660),

    histtype="step",

    cumulative=True

)



def plot_hist(fig, ax, observation_start_mjd, label):

    x = observation_start_mjd - np.min(observation_start_mjd)

    ax.hist(x, label=label, **kwargs_hist)



    select = x < 3 * 365.25

    print(

        "%20s" % label,

        "%6d" % observation_start_mjd.size,

        "%6d" % select.sum()

    )



    return fig, ax





fig = plt.figure()

ax = fig.gca()



# AGILE DR1

table = Table.read("data/baseline_v4.0_10yrs.fits")

observation_ids = np.loadtxt("data/baseline_v4.0_observation_ids_cosmos_3yr.txt")

select = np.isin(table["observationId"], observation_ids)

plot_hist(fig, ax, table["observationStartMJD"][select], label="AGILE DR1")



# DP1 ECDFS

from astropy.coordinates import SkyCoord

import astropy.units as u

table = Table.read("data/dp1/visit.fits")

coord_ecdfs = SkyCoord(53.13, -28.10, unit="deg")

coord_visit = SkyCoord(table["ra"], table["dec"])

select = coord_ecdfs.separation(coord_visit) < 1 * u.deg

plot_hist(fig, ax, table["obsStartMJD"][select], label="DP1 ECDFS")



ax.set_xlim(0, 3.1 * 366)

ax.set_xlabel("MJD")

ax.set_ylabel("Cumulative visits")

ax.legend(loc="lower right")

fig.savefig("fig/baseline_visit_agile_dp1.png")



fig = plt.figure()

ax = fig.gca()



# baseline_v4.0

table = Table.read("data/baseline_v4.0_10yrs.fits")

select = table["target_name"] == "DD:COSMOS"

plot_hist(fig, ax, table["observationStartMJD"][select], label="baseline_v4.0")



# baseline_v4.3.5

with sqlite3.connect("data/baseline_v4.3.5_10yrs.db") as con:

    df = pd.read_sql_query(

        """

        SELECT observationStartMJD

        FROM observations

        WHERE target_name LIKE '%DD:COSMOS%'

        """,

        con

    )

table = Table.from_pandas(df)

plot_hist(fig, ax, table["observationStartMJD"], label="baseline_v4.3.5")



# DDF Ocean6

with sqlite3.connect("data/ddf_ocean_ocean6_v4.3.5_10yrs.db") as con:

    df = pd.read_sql_query(

        """

        SELECT observationStartMJD

        FROM observations

        WHERE target_name='DD:COSMOS'

        """,

        con

    )

table = table.from_pandas(df)

plot_hist(fig, ax, table["observationStartMJD"], label="DDF Ocean6")



ax.set_xlabel("MJD")

ax.set_ylabel("Cumulative visits")

ax.legend(loc="lower right")

fig.savefig("fig/baseline_visit_v4_ocean6.png")