Analysis of Photometric Offsets in Astronomical Data using Python

PhotometricOffsetFigures.ipynb

This Jupyter notebook contains a detailed analysis focused on comparing photometric data from the HUGS and ACS surveys, specifically targeting the NGC 2808 cluster. It involves data cleaning, sampling, and visualization using libraries like pandas, matplotlib, and numpy. Furthermore, it applies fiducial line measurement techniques to discern photometric offsets and undertakes isochrone fitting with Bolometric Corrections to study stellar populations. The notebook concludes with a unique assessment of helium mass fraction effects on turn-off points in color-magnitude diagrams, potentially offering insights into stellar evolution.

snippet.python

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from fidanka.fiducial.fiducial import measure_fiducial_lines
from fidanka.fiducial.utils import bin_color_mag_density
import pickle as pkl
from matplotlib.lines import Line2D
from fidanka.fiducial.fiducial import hull_density
from scipy.interpolate import splrep, BSpline
from scipy.interpolate import interp1d
from fidanka.isochrone.MIST import read_iso
from fidanka.isochrone.isochrone import interp_isochrone_age
from fidanka.bolometric.bctab import BolometricCorrector
import re
from scipy.optimize import curve_fit
from fidanka.isochrone.MIST import read_iso
from fidanka.isochrone.isochrone import interp_isochrone_age

snippet.python

ACSRoot = "../photometry/hlsp_acsggct_hst_acs-wfc_ngc2808_r.rdviq.cal.adj.zpt"
HUGSRoot = "../photometry/HUGS/ngc2808/photometry.pkl"

snippet.python

with open(HUGSRoot, 'rb') as f:
    HUGSdata = pkl.load(f)

snippet.python

ACSdata = pd.read_csv(ACSRoot, sep=r'\s+')
ACSdata = ACSdata[ACSdata['Vvega'] - ACSdata['Ivega'] > 0.5]

snippet.python

n = 20000
ACSSub = ACSdata.sample(n)
ACSmag = ACSSub["Vvega"]
ACScolor = ACSSub["Vvega"] - ACSSub["Ivega"]
 
HUGSSub = HUGSdata[1].sample(n)
HUGSmag = HUGSSub["F606W"]
HUGScolor = HUGSSub["F606W"] - HUGSSub["F814W"]

snippet.python

ACSdensity = hull_density(ACScolor, ACSmag)
HUGSdensity = hull_density(HUGScolor, HUGSmag)

snippet.python

with plt.style.context(pubStyle):
    fig, ax = plt.subplots(1,1,figsize=(10,7))
    ax.scatter(HUGScolor, HUGSmag,s=1,c=HUGSdensity, cmap="Reds")
    ax.scatter(ACScolor, ACSmag, s=1, c=ACSdensity, cmap='Blues', alpha=0.25)
    ax.set_xlim(0,2)
    ax.set_ylim(16, 25)
    ax.invert_yaxis()
 
    ax.set_xlabel("F606W - F814W", fontsize=27)
    ax.set_ylabel("F606W", fontsize=27)
 
 
    legend_elements = [
        Line2D([0], [0], marker='o', color='w', label='HUGS', markerfacecolor='r', markersize=15),
        Line2D([0], [0], marker='o', color='w', label='ACS', markerfacecolor='b', markersize=15),
    ]
    ax.legend(handles = legend_elements, frameon=False, fontsize=27)
    fig.savefig("Figures/photometricOffset.pdf")

snippet.python

ACSFiducial = measure_fiducial_lines(
    ACSSub['Vvega'].values,
    ACSSub['Ivega'].values,
    ACSSub['err'].values,
    ACSSub['err.2'].values,
    binSize='uniformCS',
    targetStat=100,
    mcruns=5,
    convexHullPoints=50,
    nPops=1)

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:58<00:00, 11.61s/it]

snippet.python

HUGSFiducial = measure_fiducial_lines(
    HUGSSub['F606W'].values,
    HUGSSub['F814W'].values,
    HUGSSub['F606W_RMS'].values,
    HUGSSub['F814W_RMS'].values,
    binSize='uniformCS',
    targetStat=100,
    mcruns=5,
    convexHullPoints=50,
    nPops=1)

100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████| 5/5 [00:52<00:00, 10.48s/it]

snippet.python

from fidanka.misc.utils import measusre_perpendicular_distance

snippet.python

with plt.style.context(pubStyle):
    fig, axs = plt.subplots(1,2,figsize=(20,7))
    ax = axs[0]
    ax.scatter(HUGScolor, HUGSmag,s=1,c=HUGSdensity, cmap="Reds")
    ax.scatter(ACScolor, ACSmag, s=1,c=ACSdensity, cmap='Blues', alpha=0.1)
    ax.set_xlim(0,2)
    ax.set_ylim(17, 22)
    ax.invert_yaxis()
 
    ACSSmoothTCK = splrep(ACSFiducial[0].mean[1], ACSFiducial[0].mean[0], s=0.001)
    ax.plot(BSpline(*ACSSmoothTCK)(ACSFiducial[0].mean[1]), ACSFiducial[0].mean[1] , color='g')
    HUGSmoothTCK = splrep(HUGSFiducial[0].mean[1], HUGSFiducial[0].mean[0], s=0.001)
    ax.plot(BSpline(*HUGSmoothTCK)(HUGSFiducial[0].mean[1]), HUGSFiducial[0].mean[1] , color='orange')
    ax.set_xlim(0.6, 1)
 
    ax.set_xlabel("F606W - F814W", fontsize=27)
    ax.set_ylabel("F606W", fontsize=27)
 
 
    legend_elements = [
        Line2D([0], [0], marker='o', color='w', label='HUGS', markerfacecolor='r', markersize=10),
        Line2D([0], [0], marker='o', color='w', label='ACS', markerfacecolor='b', markersize=10),
    ]
    ax.legend(handles = legend_elements, frameon=False, fontsize=27)
 
 
    ax = axs[1]
    ax.set_xlabel("F606W", fontsize=27)
    ax.set_ylabel("HUGS$_{color}$ - ACS$_{color}$", fontsize=27)
 
    domain = np.linspace(19.5, 22, 1000)
    dist = measusre_perpendicular_distance(BSpline(*HUGSmoothTCK), BSpline(*ACSSmoothTCK), domain)
    ax.plot(domain, dist, color='black')
    ax.axhline(np.mean(dist), color='black', linestyle='dashed')
    # ax.axhline(np.median(dist))
    fig.savefig("Figures/photometricOffset.pdf")

snippet.python

plt.plot(np.arange(17, 22, 0.1),measusre_perpendicular_distance(BSpline(*HUGSmoothTCK), BSpline(*ACSSmoothTCK), np.arange(17, 22, 0.1)), 'o')

[<matplotlib.lines.Line2D at 0x7fa884d2c5b0>]

snippet.python

bc = BolometricCorrector("WFC3", 1.13)

snippet.python

rootPath = "../outputs.denseAlpha.fixedLowmass/"

snippet.python

import pathlib

snippet.python

paths = list()
for file in pathlib.Path(rootPath).rglob("alpha-1.901/isochrones.txt"):
    paths.append(file)

snippet.python

isos = list()
testIso = read_iso(paths[0])
header = testIso[list(testIso.keys())[0]].columns
for path in paths:
    isos.append(pd.DataFrame(interp_isochrone_age(read_iso(path), 12.3), columns=header))

snippet.python

isos[0]

<div> <style scoped>

  .dataframe tbody tr th:only-of-type {
      vertical-align: middle;
  }

  .dataframe tbody tr th {
      vertical-align: top;
  }

  .dataframe thead th {
      text-align: right;
  }

</style> <table border=“1” class=“dataframe”

>

<thead>
  <tr style="text-align: right;">
    <th></th>
    <th>EEP</th>
    <th>log10_isochrone_age_yr</th>
    <th>initial_mass</th>
    <th>star_mass</th>
    <th>star_mdot</th>
    <th>he_core_mass</th>
    <th>c_core_mass</th>
    <th>log_L</th>
    <th>log_LH</th>
    <th>log_LHe</th>
    <th>...</th>
    <th>surface_c12</th>
    <th>surface_o16</th>
    <th>log_center_T</th>
    <th>log_center_Rho</th>
    <th>center_gamma</th>
    <th>center_h1</th>
    <th>center_he4</th>
    <th>center_c12</th>
    <th>phase</th>
    <th>age</th>
  </tr>
</thead>
<tbody>
  <tr>
    <th>0</th>
    <td>286.0</td>
    <td>10.089905</td>
    <td>0.235747</td>
    <td>0.235747</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.0</td>
    <td>-1.989956</td>
    <td>-1.989276</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000061</td>
    <td>0.000820</td>
    <td>6.849942</td>
    <td>2.109066</td>
    <td>0.0</td>
    <td>0.688852</td>
    <td>0.309676</td>
    <td>0.000061</td>
    <td>0.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>1</th>
    <td>287.0</td>
    <td>10.089905</td>
    <td>0.260267</td>
    <td>0.260267</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.0</td>
    <td>-1.904474</td>
    <td>-1.903862</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000061</td>
    <td>0.000820</td>
    <td>6.863070</td>
    <td>2.051913</td>
    <td>0.0</td>
    <td>0.688946</td>
    <td>0.309582</td>
    <td>0.000061</td>
    <td>0.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>2</th>
    <td>288.0</td>
    <td>10.089905</td>
    <td>0.262639</td>
    <td>0.262639</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.0</td>
    <td>-1.897316</td>
    <td>-1.896708</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000061</td>
    <td>0.000820</td>
    <td>6.864233</td>
    <td>2.047220</td>
    <td>0.0</td>
    <td>0.688909</td>
    <td>0.309619</td>
    <td>0.000061</td>
    <td>0.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>3</th>
    <td>289.0</td>
    <td>10.089905</td>
    <td>0.265015</td>
    <td>0.265015</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.0</td>
    <td>-1.890155</td>
    <td>-1.889553</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000061</td>
    <td>0.000820</td>
    <td>6.865400</td>
    <td>2.042532</td>
    <td>0.0</td>
    <td>0.688879</td>
    <td>0.309649</td>
    <td>0.000061</td>
    <td>0.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>4</th>
    <td>290.0</td>
    <td>10.089905</td>
    <td>0.267398</td>
    <td>0.267398</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.0</td>
    <td>-1.882988</td>
    <td>-1.882392</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000061</td>
    <td>0.000820</td>
    <td>6.866570</td>
    <td>2.037846</td>
    <td>0.0</td>
    <td>0.688857</td>
    <td>0.309671</td>
    <td>0.000061</td>
    <td>0.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>...</th>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
    <td>...</td>
  </tr>
  <tr>
    <th>187</th>
    <td>601.0</td>
    <td>10.089905</td>
    <td>0.763934</td>
    <td>0.412178</td>
    <td>0.0</td>
    <td>0.332015</td>
    <td>0.0</td>
    <td>2.349763</td>
    <td>2.347601</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000056</td>
    <td>0.000788</td>
    <td>7.712642</td>
    <td>5.599579</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.998341</td>
    <td>0.000004</td>
    <td>2.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>188</th>
    <td>602.0</td>
    <td>10.089905</td>
    <td>0.763948</td>
    <td>0.410825</td>
    <td>0.0</td>
    <td>0.333765</td>
    <td>0.0</td>
    <td>2.364422</td>
    <td>2.362240</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000056</td>
    <td>0.000788</td>
    <td>7.715019</td>
    <td>5.605118</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.998341</td>
    <td>0.000004</td>
    <td>2.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>189</th>
    <td>603.0</td>
    <td>10.089905</td>
    <td>0.763962</td>
    <td>0.409429</td>
    <td>0.0</td>
    <td>0.335503</td>
    <td>0.0</td>
    <td>2.379073</td>
    <td>2.376870</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000056</td>
    <td>0.000788</td>
    <td>7.717450</td>
    <td>5.610761</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.998342</td>
    <td>0.000004</td>
    <td>2.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>190</th>
    <td>604.0</td>
    <td>10.089905</td>
    <td>0.763978</td>
    <td>0.407971</td>
    <td>0.0</td>
    <td>0.337333</td>
    <td>0.0</td>
    <td>2.393712</td>
    <td>2.391482</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000056</td>
    <td>0.000788</td>
    <td>7.719974</td>
    <td>5.616615</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.998342</td>
    <td>0.000004</td>
    <td>2.0</td>
    <td>1.230187e+10</td>
  </tr>
  <tr>
    <th>191</th>
    <td>605.0</td>
    <td>10.089905</td>
    <td>0.764000</td>
    <td>0.406203</td>
    <td>0.0</td>
    <td>0.339454</td>
    <td>0.0</td>
    <td>2.408234</td>
    <td>2.405955</td>
    <td>-10.0</td>
    <td>...</td>
    <td>0.000056</td>
    <td>0.000788</td>
    <td>7.722938</td>
    <td>5.623576</td>
    <td>0.0</td>
    <td>0.000000</td>
    <td>0.998342</td>
    <td>0.000004</td>
    <td>3.0</td>
    <td>1.230187e+10</td>
  </tr>
</tbody>

</table>

<

p>192 rows × 26 columns</p> </div>

snippet.python

mags = list()
for iso in isos:
    mags.append(bc.apparent_mags(
        10**iso['log_Teff'],
        iso['log_g'],
        iso['log_L'],
        filters=("WFC3_UVIS_F606W", "WFC3_UVIS_F814W")
    )
               )

snippet.python

turnOffs = list()
for mag, path in zip(mags, paths):
    color = mag["WFC3_UVIS_F606W"]-mag["WFC3_UVIS_F814W"]
    idStr = re.findall('Pop[A|E]\+0\.\d+', str(path))[0]
    pop = idStr[3]
    Y = float(idStr[5:])
    turnOffColor = min(mag["WFC3_UVIS_F606W"] - mag["WFC3_UVIS_F814W"])
    turnOffMag = mag[color==turnOffColor]["WFC3_UVIS_F606W"]
    if turnOffColor < 0.48:
        turnOffs.append((turnOffColor, turnOffMag.values[0], pop, Y))
 
meanTurnOffColorA = sum([x[0] for x in turnOffs if x[2] == 'A'])/(len(turnOffs)/2)
meanTurnOffMagA = sum([x[1] for x in turnOffs if x[2] == 'A'])/(len(turnOffs)/2)
meanTurnOffColorE = sum([x[0] for x in turnOffs if x[2] == 'E'])/(len(turnOffs)/2)
meanTurnOffMagE = sum([x[1] for x in turnOffs if x[2] == 'E'])/(len(turnOffs)/2)

snippet.python

hevolve = list()
for mag, path, turnOff in zip(mags, paths, turnOffs):
    formater = {'A': 'ko', 'E': 'kx'}
    if turnOff[2] == 'A':
        dist = np.sqrt((turnOff[0]-turnOffs[0][0])**2 + (turnOff[1]-turnOffs[0][1])**2)
        hevolve.append((turnOff[3]-turnOffs[0][3], dist))
hevolve = np.array(hevolve)

snippet.python

with plt.style.context(pubStyle):
    fig, ax = plt.subplots(1,1,figsize=(10,7))
    for h in hevolve:
        ax.plot(h[0], h[1], 'ko')
 
    line = lambda x, m, b:  m*x + b
    fit, covar = curve_fit(line, hevolve[:, 0], hevolve[:, 1])
    X = np.linspace(hevolve[:, 0].min(), hevolve[:, 0].max())
    ax.plot(X, line(X, *fit), 'k--')
    ax.set_xlabel("$\Delta$ Helium Mass Fraction", fontsize=20)
    ax.set_ylabel("Turn off offset from mean [mag]", fontsize=20)
    ax.annotate(f'$O(Y)={fit[0]:0.2f}Y+{fit[1]:0.2f}$', xy=(-0.04, 0.30), color='black', fontsize=20)
 
    fig.savefig("Figures/HeliumMeanOffset.pdf")

snippet.python

with open("BestFit.pkl", 'rb') as f:
    mags = pkl.load(f)

snippet.python

def correct_iso_set(path, age, mu, Av):
    print(os.path.exists(path))
 
    isoAtAge = interp_isochrone_age(iso, age)
 
    isoAMags = bc.apparent_mags(
        10**isoAtAge[:, 10],
        isoAtAge[:, 12],
        isoAtAge[:, 7],
        mu = mu,
        Av = Av,
        filters = ("WFC3_UVIS_F275W", "WFC3_UVIS_F606W", "WFC3_UVIS_F814W")
    )
    return isoAMags

snippet.python

testPath = "/mnt/Astronomy/GraduateSchool/Thesis/GCConsistency/NGC2808/justIso.denseAlpha.fixedLowmass/PopA+0.39/alpha-1.750/isochrones.txt"

snippet.python

isoMags = correct_iso_set(testPath, 17, 1.491e+01, 5.414e-01)

True

snippet.python

fig, ax = plt.subplots(1,1,figsize=(10,10))
ax.scatter(HUGScolor, HUGSmag,s=1,c=HUGSdensity, cmap="Reds")
ax.scatter(ACScolor, ACSmag, s=1,c=ACSdensity, cmap='Blues', alpha=0.1)
ax.set_ylim(17, 22)
ax.invert_yaxis()
ACSSmoothTCK = splrep(ACSFiducial[0].mean[1], ACSFiducial[0].mean[0], s=0.001)
ax.plot(BSpline(*ACSSmoothTCK)(ACSFiducial[0].mean[1]), ACSFiducial[0].mean[1] , color='C0')
HUGSmoothTCK = splrep(HUGSFiducial[0].mean[1], HUGSFiducial[0].mean[0], s=0.001)
ax.plot(BSpline(*HUGSmoothTCK)(HUGSFiducial[0].mean[1]), HUGSFiducial[0].mean[1] , color='C1')
ax.set_xlim(0.6, 1)
 
ax.set_xlabel("F606W - F814W", fontsize=27)
ax.set_ylabel("F606W", fontsize=27)
 
Color = isoMags['WFC3_UVIS_F606W'].values - isoMags["WFC3_UVIS_F814W"].values
Mag = isoMags['WFC3_UVIS_F606W'].values
turnOffColor = min(Color)
turnOffMag = Mag[Color==turnOffColor][0]
print(turnOffMag, turnOffColor)
 
ax.plot(Color, Mag, color='C2', linewidth=10)

19.73960461108486 0.652611121597868





[<matplotlib.lines.Line2D at 0x7fa87ae821c0>]

snippet.python

 
 
 
 
 
<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }
 
    .dataframe tbody tr th {
        vertical-align: top;
    }
 
    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>WFC3_UVIS_F275W</th>
      <th>WFC3_UVIS_F606W</th>
      <th>WFC3_UVIS_F814W</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>33.753726</td>
      <td>26.699337</td>
      <td>24.461482</td>
    </tr>
    <tr>
      <th>1</th>
      <td>33.737848</td>
      <td>26.683081</td>
      <td>24.448594</td>
    </tr>
    <tr>
      <th>2</th>
      <td>33.722102</td>
      <td>26.666965</td>
      <td>24.435789</td>
    </tr>
    <tr>
      <th>3</th>
      <td>33.706476</td>
      <td>26.650976</td>
      <td>24.423060</td>
    </tr>
    <tr>
      <th>4</th>
      <td>33.690944</td>
      <td>26.635087</td>
      <td>24.410393</td>
    </tr>
    <tr>
      <th>...</th>
      <td>...</td>
      <td>...</td>
      <td>...</td>
    </tr>
    <tr>
      <th>313</th>
      <td>23.471173</td>
      <td>14.562014</td>
      <td>13.486965</td>
    </tr>
    <tr>
      <th>314</th>
      <td>23.445802</td>
      <td>14.533347</td>
      <td>13.452625</td>
    </tr>
    <tr>
      <th>315</th>
      <td>23.420594</td>
      <td>14.504854</td>
      <td>13.418372</td>
    </tr>
    <tr>
      <th>316</th>
      <td>23.395397</td>
      <td>14.476370</td>
      <td>13.384128</td>
    </tr>
    <tr>
      <th>317</th>
      <td>23.370260</td>
      <td>14.447949</td>
      <td>13.349919</td>
    </tr>
  </tbody>
</table>
<p>318 rows × 3 columns</p>
</div>