script_completeness_2mm.py 2.89 KB
Newer Older
LUSTIG Peter's avatar
LUSTIG Peter committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
import numpy as np
import matplotlib.pyplot as plt
import sys
import astropy.units as u
from astropy.table import Table
from astropy.io import fits
from astropy.wcs import WCS
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness')
from completness2d import CompletenessArea, LogBinEdges, ChiSquaredPerNDF
from completness2d import ModifiedErrorFunction, LevMarLSQFitter, SaveFigure
#%autoreload 2

# Load HLS mask
maskfilename = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
                'Completness/masks.fits')

with fits.open(maskfilename) as hdul:
    hitmask = hdul['hitmask'].data.astype(bool)
    hlsmask = hdul['hls3fwhm'].data.astype(bool)
    wcs = WCS(hdul['hitmask'].header)

# Load SNR masks
anothermaskfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/verify/'
                   'snrmasks_fine_new.fits')
snrmasks = []
with fits.open(anothermaskfile) as hdul:
    i = 0
    while True:
        try:
            hdu = hdul['SNRMASK{}'.format(i)]
            snrmasks.append(hdu.data.astype(bool))
            if i == 0:
                wcs = WCS(hdu.header)
                shape = hdu.data.shape
        except KeyError:
            break
        i += 1

# Load Monte Carlo Data
filename = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
            'Completness/DetectionNewrealisationPixmeans.fits')
#filename = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
#            'montecarlo_results/Newrealisation/2mm/AllRealdetections.fits')
#            DetectionNewrealisationPixmeans.fits

with fits.open(filename) as hdul:
    sources = Table.read(hdul, 'DETECTED_SOURCES')
    fake_sources = Table.read(hdul, 'FAKE_SOURCES')

fluxedges = LogBinEdges(fake_sources['amplitude']) * u.Jy
fluxcenters = u.Quantity(np.squeeze(np.sort(np.unique(
                                fake_sources['amplitude']))))


threshbins = np.arange(2.5, 4.501, .01)
threshbins = np.linspace(2.5, 10, 1001)
threshbins = np.linspace(2.5, 5, 6)
# %%
%matplotlib tk
totmask = snrmasks[0] | hlsmask
obj = CompletenessArea(sources, fake_sources, shape=(501, 501), wcs=wcs,
                       threshold_edges=threshbins, flux_edges=fluxedges,
                       mask=totmask, flux_centers=fluxcenters)
plt.close('all')
fct, (ax0, ax1) = obj.FitCompleteness(idx=2, constant='thresh', weighted=True)
# SaveFigure(ax0.get_figure(), 'completeness_inner_fit')


# %%

fig = plt.figure(figsize=(6, 4))
ax = fig.add_axes([.2, .18, .75, .7])
ax2D = obj.Plot2D(ax=ax)
ax2D.set_ylim([10, 40])
# %%
'''
totmask2 = snrmasks[2] | hlsmask
obj2 = CompletenessArea(sources, fake_sources, shape=(501, 501), wcs=wcs,
                        threshold_edges=threshbins, flux_edges=fluxedges,
                        mask=totmask2, flux_centers=fluxcenters)

# %%
fig = plt.figure(figsize=(6, 4))
ax = fig.add_axes([.2, .18, .75, .7])
ax2D = obj2.Plot2D(ax=ax)
ax2D.set_ylim([10, 40])
obj2.FitCompleteness(idx=2)
'''