Magnitude and SNR Cuts

from sorcha.modules.PPBrightLimit import PPBrightLimit
from sorcha.modules.PPMagnitudeLimit import PPMagnitudeLimit
from sorcha.modules.PPSNRLimit import PPSNRLimit
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

There are three filters in the survey simulator that perform simple cuts based on SNR or magnitude, the latter to either drop observations too faint to be interesting or too bright and therefore saturating the detector. These are PPSNRLimit, PPMagnitudeLimit, and PPBrightLimit.

This notebook will demonstrate each in turn. First we generate some randomised (and rather unphysical) test data.

rng = np.random.default_rng()
data_size = 1000

filters = ['r', 'g', 'i']

random_magnitudes = 17.5 + (rng.random(data_size) * 6.5)
random_SNRs = 0.5 + (rng.random(data_size) * 7)
random_filters = rng.choice(filters, data_size)

test_data = pd.DataFrame({'PSFMag': random_magnitudes, 'SNR': random_SNRs, 'optFilter': random_filters})

test_data

	PSFMag	SNR	optFilter
0	18.181149	6.134108	r
1	23.134903	2.389628	r
2	22.936785	4.619613	g
3	23.031516	3.607125	g
4	23.069938	2.539725	r
...	...	...	...
995	23.613129	7.040533	r
996	19.478289	5.023488	g
997	20.135023	1.512182	r
998	23.342629	2.482078	r
999	19.186506	6.342341	r

1000 rows × 3 columns

Let’s look at the magnitude filters first. The below shows the original randomised magnitudes.

fig, ax = plt.subplots(1)
ax.hist(random_magnitudes, 50, color='thistle')
ax.set_ylabel('frequency')
ax.set_xlabel('magnitude')
ax.autoscale(axis='x', tight=True)
plt.show()

Testing first the magnitude filter for cutting faint objects, we will cut everything fainter than 21 mag.

test_data_maglimit = PPMagnitudeLimit(test_data, 21.,colname='PSFMag')

maglimit_magnitudes = test_data_maglimit['PSFMag'].values

fig, ax = plt.subplots(1)
ax.hist(maglimit_magnitudes, 50, color='thistle')
ax.set_ylabel('frequency')
ax.set_xlabel('magnitude')
ax.set_xlim((17.5, 24))
ax.axvline(21, linestyle='--', linewidth=1., color='darkslateblue')
plt.show()

Now we will use the brightness limit filter. This allows the user to set a different saturation limit for each filter. As we have data in g, r and i, we will set the saturation limits as 18, 19 and 20 (chosen for clarity, not accuracy!).

test_data_brightlimit = PPBrightLimit(test_data, filters, [18., 19., 20.])

r_magnitudes = test_data_brightlimit[test_data_brightlimit['optFilter'] == 'r']['PSFMag'].values
g_magnitudes = test_data_brightlimit[test_data_brightlimit['optFilter'] == 'g']['PSFMag'].values
i_magnitudes = test_data_brightlimit[test_data_brightlimit['optFilter'] == 'i']['PSFMag'].values

fig, ax = plt.subplots(3, figsize=(5, 12))

ax[0].hist(r_magnitudes, 50, color='lightpink')
ax[0].axvline(18, linestyle='--', linewidth=1., color='crimson')
ax[0].set_ylabel('frequency')

ax[1].hist(g_magnitudes, 50, color='lightgreen')
ax[1].axvline(19, linestyle='--', linewidth=1., color='forestgreen')

ax[2].hist(i_magnitudes, 50, color='paleturquoise')
ax[2].axvline(20, linestyle='--', linewidth=1., color='lightseagreen')

for a in ax:
    a.set_xlabel('magnitude')
    a.set_xlim((17.5, 24))

plt.show()

Alternatively, we can set a single brightness limit of 18.5 across all filters.

test_data_brightlimit_single = PPBrightLimit(test_data, filters, 18.5)

singlelimit_magnitudes = test_data_brightlimit_single['PSFMag'].values

fig, ax = plt.subplots(1)
ax.hist(singlelimit_magnitudes, 50, color='thistle')
ax.set_ylabel('frequency')
ax.set_xlabel('magnitude')
ax.set_xlim((17.5, 24))
ax.axvline(18.5, linestyle='--', linewidth=1., color='darkslateblue')
plt.show()

Now we will look at the SNR limit. First, the original randomised SNRs.

fig, ax = plt.subplots(1)
ax.hist(random_SNRs, 50, color='palegoldenrod')
ax.set_ylabel('frequency')
ax.set_xlabel('SNR')
ax.autoscale(axis='x', tight=True)
plt.show()

Now we will cut everything with an SNR lower than 2 (which is the default in the code).

test_data_SNRlimit = PPSNRLimit(test_data, 2.)

limit_SNRs = test_data_SNRlimit['SNR'].values

fig, ax = plt.subplots(1)
ax.hist(limit_SNRs, 50, color='palegoldenrod')
ax.set_ylabel('frequency')
ax.set_xlabel('SNR')
ax.set_xlim((0.5, 7.5))
ax.axvline(2, linestyle='--', linewidth=1., color='gold')
plt.show()