Commit fabef1f8 authored by LUSTIG Peter's avatar LUSTIG Peter

update...

parent 4d15fef8
......@@ -249,7 +249,8 @@ class CatalogShifter(object):
transform: callable or None.
If not None, takes coord_in_center and rotation as argument
and transforms them to a :class:`astropy.coordinates.SkyCoord`
and an angle that are used as coord_in_center and rotation.
and an angle that are used as coord_in_center and rotation to
perform transformation.
Returns
-------
catalog : :class:``astropy.table.Table``
......
......@@ -6,21 +6,22 @@ import numpy as np
import sys
from utils import CombineMeasurements
from astropy.io import fits
from astropy.utils.console import ProgressBar
def LoadSingleFile(fname):
_file = open(fname)
flux = u.Quantity(_file[0].header['influx'])
with open(fname) as _file:
flux = u.Quantity(_file[0].header['influx'])
try:
sources = Table.read(_file['DETECTED_SOURCES'])
except KeyError:
print('no sources detected for {}'.format(flux))
sources = None
try:
sources = Table.read(_file['DETECTED_SOURCES'])
except KeyError:
print('no sources detected for {}'.format(flux))
sources = None
fsources = Table.read(_file['FAKE_SOURCES'])
fsources = Table.read(_file['FAKE_SOURCES'])
dthresh = _file[0].header['DTHRESH']
dthresh = _file[0].header['DTHRESH']
return flux, sources, fsources, dthresh
......@@ -116,22 +117,23 @@ def MergeToOneTable(directory, key, outfile, overwrite=False):
header = fits.Header()
print('loading tables')
for i, filename in enumerate(filenames):
print(i)
flux, sources, fsources, dthresh = LoadSingleFile(filename)
header['flux{}'.format(i)] = '{}'.format(flux)
sourceslist.append(sources)
fakesourceslist.append(fsources)
if dthresh > max_dthresh:
max_dthresh = dthresh
with ProgressBar(len(filenames)) as pg:
for i, filename in enumerate(filenames):
flux, sources, fsources, dthresh = LoadSingleFile(filename)
header['flux{}'.format(i)] = '{}'.format(flux)
sourceslist.append(sources)
fakesourceslist.append(fsources)
if dthresh > max_dthresh:
max_dthresh = dthresh
pg.update()
print('done')
header['dthresh'] = max_dthresh
print('combining measurements')
sources, fake_sources = CombineMeasurements(sourceslist, fakesourceslist)
sources['ID'].fill_value = -999999
sources['fake_sources'].fill_value = -999999
sources['group_id'].fill_value = -999999
fake_sources['ID'].fillvalue = -999999
sources['ID'].fill_value = -999999999
sources['fake_sources'].fill_value = -999999999
fake_sources['find_peak'].fill_value = -999999999
fake_sources['ID'].fillvalue = -999999999
print('done')
hdul = HDUList([PrimaryHDU(),
......@@ -149,11 +151,18 @@ directory = 'montecarlo_results/700/'
directory = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'montecarlo_results/photometry/NewRealisation')
# outname = 'NewrealisationPhotometry.fits'
outname = 'NewrealisationPhotometry.fits'
directory = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
'testdir/test/')
directory = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/detection'
directory = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/photometry'
# outname = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/'
# 'PhotometryNewrealisation3.5_MultTables.fits')
directory = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/2mm/realdetection'
outname = '/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/montecarlo_results/Newrealisation/2mm/AllRealdetections.fits'
key = '*'
Combine(directory, key, outname)
# MergeToOneTable(directory, key, outname, overwrite=True)
# Combine(directory, key, outname)
MergeToOneTable(directory, key, outname, overwrite=False)
print('done')
'''
......
This diff is collapsed.
......@@ -25,11 +25,13 @@ from mpl_toolkits.axes_grid1 import make_axes_locatable
import dill as pickle
from matplotlib.ticker import FormatStrFormatter
from collections import OrderedDict
sys.path.append('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness')
from utils import completness_purity_wcs, completness_worker, purity_worker
from utils import find_nearest, DeletePercentiles
from utils import pos_in_mask, Flux1D
from astropy.units import Quantity
# plt.close('all')
class PCEvaluation:
def __init__(self, sources, fake_sources=None, realsources=None,
......@@ -415,7 +417,6 @@ class PCAreaEvaluation:
self.fluxpercentiles = np.zeros((len(sources), 5))
for iflux in range(len(sources)):
print('loop {}'.format(iflux))
_sources = AddXYCoordinatesTo(sources[iflux], wcs)
if fake_sources is not None:
......@@ -812,10 +813,12 @@ if __name__ == '__main__':
sh = data.data.shape
oldwcs = data.wcs
if 1:
if 0:
#completnessareafile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
# 'Completness/allresults_phot_thresh4.fits')
completnessareafile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/allresults_phot_thresh4.fits')
'Completness/PhotometryNewrealisation3.5_MultTables.fits')
pflux, psources, pfsources = UglyLoader(completnessareafile)
exc = 5
......@@ -833,9 +836,8 @@ if __name__ == '__main__':
#bb = PCAreaEvaluation(mask=globmask, sources=sources, wcs=wcs, nsim=nsim,
# realsources=realsources, flux=pflux, fake_sources=None,
# threshold_range=(2.5, 10), threshold_bins=50)
'''
# %%
# % matplotlib tk
plt.figure(figsize=(8, 7))
bb.PlotCounts()
plt.xlim(xmin=2.5, xmax=5)
......@@ -848,11 +850,6 @@ if __name__ == '__main__':
# % matplotlib tk
bb.PlotPercentiles()
plt.show(block=True)
# %%
#bb.FitFluxResolution(10*u.mJy)
#bb.FitFluxResolution(5*u.mJy)
#bb.FitFluxResolution(1*u.mJy)
# plt.show(block=True)
'''
......@@ -872,16 +869,19 @@ if __name__ == '__main__':
aa.PlotDarkcountsThresholds()
plt.show(block=True)
'''
if 0:
if 1:
# %matplotlib tk
fname = ('/home/peter/Dokumente/Uni/Paris/Stage/'
'FirstSteps/Completness/NEWcombined_tables_long.fits')
fname = ('/home/peter/Dokumente/Uni/Paris/Stage/'
'FirstSteps/Completness/NEW_combine_fct_result.fits')
fname = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/PhotometryNewrealisation3.5_MultTables.fits')
FLUX, SOURCE, FSOURCE = UglyLoader(fname)
if 0:
# %matplotlib tk
if 1:
xx = PCEvaluation(SOURCE, FSOURCE, shape=sh, wcs=oldwcs, flux=FLUX,
mapbins=19, threshold_bins=6,
threshold_range=(2.5, 5))
......@@ -894,19 +894,34 @@ if __name__ == '__main__':
ylabel='Purity')
xx.PlotOverview(flux=5*u.mJy)
xx.PlotHitmap(flux=5*u.mJy)
if 0:
plt.close('all')
# %%
%matplotlib tk
# xx.thresholds
fig = plt.figure()
plt.imshow(xx.completness[27, :, :, 2], origin='lower')
50-23
plt.imshow(xx.completness[37, :, :, 2], origin='lower')
if 1:
fluxidxs = [27, 25, 23, 21, 19, 17, 15, 13, 11, 29]
fluxidxs = np.arange(11, 50, 2)
threshidxs = [2]*len(fluxidxs)
_hdul = [fits.PrimaryHDU()]
_hdul.append(fits.ImageHDU(data=xx.completness[-3, :, :, -2],
name='CompletnessMap'))
_hdul[1].header['flux'] = '{}'.format(xx.flux[-3])
_hdul[1].header['dthresh'] = '{}'.format(xx.thresholds[-2])
_hdul[1].header.extend(xx.wcs_3D.sub([1, 2]).to_header())
for ii, (fluxidx, threshidx) in enumerate(zip(fluxidxs, threshidxs)):
_hdul.append(fits.ImageHDU(data=xx.completness[fluxidx, :, :, threshidx],
name='CompletenessBinned{}'.format(ii)))
_hdul[1+ii].header['flux'] = '{}'.format(xx.flux[fluxidx])
_hdul[1+ii].header['dthresh'] = xx.thresholds[threshidx]
_hdul[1+ii].header.extend(xx.wcs_3D.sub([1, 2]).to_header())
_hdul = fits.HDUList(_hdul)
_hdul.writeto('CompletnessBinnedMapExample.fits',
_hdul.writeto('/home/peter/Dokumente/Uni/Paris/Stage/'
'FirstSteps/Completness/'
'CompletnessBinnedMapExampleNew.fits',
overwrite=True)
if 1:
if 0:
yy = PCEvaluation(SOURCE, FSOURCE, shape=sh, wcs=oldwcs, flux=FLUX,
mapbins=9, threshold_bins=6,
threshold_range=(2.5, 5))
......@@ -916,7 +931,7 @@ if __name__ == '__main__':
yy.PlotFixedThreshold(yy.GetPurityBin(4, 4), np.array([3, 5]),
ylabel='Purity')
plt.show(block=True)
# plt.show(block=True)
'''
fname = ('/home/peter/Dokumente/Uni/Paris/Stage/'
......
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)
'''
import numpy as np
import matplotlib.pyplot as plt
import sys
import astropy.units as u
from astropy.table import Table, vstack
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, FluxArea
from completness2d import SaveFigure, FluxArea
from fitmodels import FluxRelativeBoosting, FluxBoosting
from pathlib import Path
# 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/PhotometryNewrealisation3.5.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(3.5, 7, .5)
len(snrmasks)
# %%
'''
class Flux(FluxArea):
pass
'''
totmask = snrmasks[0] | hlsmask
obj = FluxArea(sources, fake_sources, shape=(501, 501), wcs=wcs,
threshold_edges=threshbins, flux_edges=fluxedges,
mask=totmask, flux_centers=fluxcenters)
cp_phot = obj
# %%
if __name__ == '__main__':
# %%
comptbl = Table(names=['maskidx', 'x0', 'x0_err', 'sigma', 'sigma_err'])
phottble = Table(names=['maskidx', 'sigma', 'sigma_err', 'Sd', 'Sd_err'])
for i, snrmask in enumerate(snrmasks):
print(i)
totmask = snrmask | hlsmask
if not np.all(obj.mask == totmask):
obj.ChangeMask(totmask)
(compfct, compfit), _ = obj.FitCompleteness(0)
cx0 = compfct.x0.value
csigma = compfct.sigma.value
cx0_err, csigma_err = np.sqrt(
np.diagonal(compfit.fit_info['param_cov']))
(photfct, photfit), _ = obj.FitFluxboost(mincompleteness=.1)
pSd = photfct.Sd.value
psigma = photfct.sigma.value
pSd_err, psigma_err = np.sqrt(
np.diagonal(photfit.fit_info['param_cov']))
comptbl = vstack([comptbl, Table(data=[[i], [cx0], [cx0_err],
[csigma], [csigma_err]],
names=['maskidx', 'x0', 'x0_err',
'sigma', 'sigma_err'])])
phottble = vstack([phottble, Table(data=[[i], [pSd], [pSd_err],
[psigma], [psigma_err]],
names=['maskidx', 'sigma',
'sigma_err', 'Sd',
'Sd_err'])])
print('Done')
# from pathlib import Path
len(comptbl)
comptbl
here = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
phottble.write(here / '1mm_photfit_results.fits', overwrite=True)
comptbl.write(here / '1mm_compfit_results.fits', overwrite=True)
if 0:
# %%
snrpercentiles = np.linspace(0, 100, 150)
snredges = np.nanpercentile(obj.sources['SNR'], snrpercentiles)
# %matplotlib tk
ax = obj.PlotPercentiles(mean=False)
ax.set_xscale('log')
ax.set_xlabel('Injected Flux [mJy / beam]', fontsize=20)
# SaveFigure(ax.get_figure(), 'fluxboost_median')
# %%
# %matplotlib tk
# obj.FluxPercentiles() # snrbins=snredges)
fitfunc, (ax, axres) = obj.FitFluxboost(weighted=True)
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
# axins = zoomed_inset_axes(ax, 1, loc=2)
fig = ax.get_figure()
axpos = ax.get_position()
axins = fig.add_axes([.5*(axpos.x1+axpos.x0)+0.015, axpos.y0+.06, .35, .15])
axins.set_xscale('log')
axins.axhline(0, color='r')
axins.set_xlim(2, 10)
axins.set_ylim(-.017, .002)
axins.tick_params('both', labelsize=12)
axins.get_xaxis().set_visible(False)
axins.errorbar(obj.percentile_center,
obj.fluxmean-fitfunc(obj.percentile_center),
obj.fluxmean_unc, linestyle='None', marker='o')
# axins.plot([1,2,3,4], 1,2,3,4)
# axins.set_xlim(1, 5)
# axins.set_xlim(0, 5)
from mpl_toolkits.axes_grid1.inset_locator import mark_inset
mark_inset(axres, axins, loc1=1, loc2=3, fc="none", ec="0.5")
# axins.axhline(0)
obj.fluxmean_unc
# SaveFigure(ax.get_figure(), 'fluxboost_fit')
obj.percentile_center
import numpy as np
import matplotlib.pyplot as plt
import sys
import astropy.units as u
from astropy.table import Table, vstack
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
from completness2d import CountsAreaEvaluation, Purity
from pathlib import Path
# 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
# Select mask
emptycountsfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/'
'FalseDetections/180513_123146_NewRealisation_'
'flux0.00mJy_thresh2.5_nsources0_nsim1000_photFalse'
'_detectTrue.fits')
emptycountsfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/'
'FalseDetections/180605_104210_NewRealisation_1mm_'
'flux0.00mJy_thresh2.5_nsources0_nsim1000_'
'photTrue_detectTrue.fits')
emptycountsfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/'
'FalseDetections/180605_164831_NewRealisation_1mm_'
'flux0.00mJy_thresh3.5_nsources0_nsim1000_photTrue'
'_detectTrue.fits')
'''
emptycountsfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/2mm/'
'falsedetections/180531_115409_NewRealisation_2mm_'
'flux0.00mJy_thresh2.5_nsources0_nsim1000_photFalse_'
'detectTrue.fits')
'''
emptysources = Table.read(emptycountsfile)
# %matplotlib tk
realcountsfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/'
'Realmap_dthresh2.5_photTrue.fits')
realcountsfile = ('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/'
'Completness/montecarlo_results/Newrealisation/'
'Realdetections/NewrealisationSources_1mm_threshold3.5.fits')
realsources = Table.read(realcountsfile)
# %%
# %matplotlib tk
threshbins = np.linspace(3.5, 10, 50)
totmask = snrmasks[0] | hlsmask
fc = CountsAreaEvaluation(
emptysources,
shape=shape, wcs=wcs, threshold_edges=threshbins,
mask=totmask, no_sim=1000)
false_counts = fc.Counts()
# fc.PlotCounts()
rc = CountsAreaEvaluation(
realsources,
shape=shape, wcs=wcs, threshold_edges=threshbins,
mask=totmask, no_sim=1)
abs_counts = rc.Counts()
plt.imshow(snrmasks[8] | hlsmask)
# matplotlib tk
plt.imshow(~snrmasks[8]+~snrmasks[7]*2)
np.sum(~(snrmasks[7] | hlsmask))
plt.imshow()
im = np.zeros(snrmasks[0].shape)
plt.close('all')
# fig = plt.figure()
fig, axes = plt.subplots(1, len(snrmasks))
for i, mask in enumerate(snrmasks):
_mask = (~mask).astype(int)
axes[i].imshow(_mask)
_mask *= (i+1)
print(np.amax(_mask))
im += _mask
print(np.sum(~(mask | hlsmask)))
# %matplotlib tk
plt.figure()
plt.imshow(im, origin='lower')
pur = Purity(abs_counts, false_counts, threshbins[0:-1])
fc.Counts()
fc.threshold_edges
# %%
if __name__ == '__main__':
cn_mask = 'Mask'
cn_false = r'$N_{\mathrm{false}}/\mathrm{sim}$'
cn_abs = r'$N_{\mathrm{abs}}$'
cn_pur = 'Purity'
t = Table(names=[cn_mask, cn_false, cn_abs, cn_pur],
dtype=[int, float, int, float])
# t.write('testlatextablenames.fits')
# t.read('testlatextablenames.fits')
for i, snrmask in enumerate(snrmasks):
totmask = snrmask | hlsmask
if not np.all(fc.mask == totmask):
fc.ChangeMask(totmask)
if not np.all(rc.mask == totmask):
rc.ChangeMask(totmask)
_fc = fc.Counts()[0]
_rc = rc.Counts()[0]
tt = Table(names=[cn_mask, cn_false, cn_abs, cn_pur],
data=[[i], [_fc], [_rc], [(_rc-_fc) / _rc]])
t = vstack([t, tt])
here = Path('/home/peter/Dokumente/Uni/Paris/Stage/FirstSteps/Completness/')
t.write(here / '1mm_purity_table.fits', overwrite=True)
if 0:
fc.PlotMask()
plt.close('all')
fig = plt.figure(figsize=(6, 5))
ax = fig.add_axes([.16, .12, .8, .78])
nax = fc.PlotCounts(label='False\nDetections', linewidth=3, ax=ax)
ax = rc.PlotCounts(ax=nax, label='Absolute\nDetections', linestyle='--',
color='r', linewidth=3)
ax.grid()
# %%
pur.Plot()
# SaveFigure(ax.get_figure(), filename='FalseCountsRealCounts')
# plt.close('all')
# %matplotlib tk
if 1:
ax1x0 = .1
ax1y0 = .145
ax1width = .37
ax1height = .75
bothplots = plt.figure(figsize=(10, 4.5))
ax1 = bothplots.add_axes([ax1x0, ax1y0, ax1width, ax1height])
fc.PlotCounts(ax=ax1, label='False\nDetections', linewidth=3)
rc.PlotCounts(ax=ax1, label='Absolute\nDetections', linestyle='--',
color='r', linewidth=3)
ax1.grid()
ax1pos = ax1.get_position()
ax2 = bothplots.add_axes([ax1pos.x1+0.13, ax1pos.y0,
ax1pos.width, ax1pos.height])
pur = Purity(abs_counts, false_counts, threshbins[0:-1])
x = threshbins[0:-1]
y = (abs_counts-false_counts) / abs_counts
completeness35 = np.interp(3.5, x, y)
pur.Plot(ax=ax2, linewidth=3)
ylims = ax2.get_ylim()
xlims = ax2.get_xlim()
ax2.plot([3.5, 3.5], [0, completeness35], color='r', linewidth=3)
ax2.plot([0, 3.5], [completeness35, completeness35], color='r', linewidth=3)
# boxprops = dict(facecolor='white', boxstyle='round, pad')
boxprops = dict(facecolor='white', edgecolor='black', boxstyle='round, pad=.3',
alpha=1)
ax2.set_ylim(.5, ylims[1])
ax2.set_xlim(xmin=xlims[0])
text = r'$\Pi(3.5)={{{:.2f}}}$'.format(completeness35)
ax2.text(.9, .1, text, ha='right', va='bottom', fontsize=20,
transform=ax2.transAxes, bbox=boxprops)
ax2.grid()
# SaveFigure(bothplots, filename='FalseCountsAndPurity')
......@@ -394,21 +394,23 @@ def PlotFixedThreshold(thresholds, bin, completness, allthresholds, flux,
def CombineMeasurements(sourceslist, fakesourceslist):
fake_sources = Table()
sources = Table()
for _fake, _detected in zip(fakesourceslist, sourceslist):
n_fake = len(fake_sources)
n_detected = len(sources)
if _detected is not None:
_detected['ID'] = _detected['ID'] + n_detected
if 'fake_sources' in _detected.keys():
_detected['fake_sources'] = _detected['fake_sources'] + n_fake
sources = vstack([sources, _detected])
if _fake is not None:
_fake['ID'] = _fake['ID'] + n_fake
if 'find_peak' in _fake.keys():
_fake['find_peak'] = _fake['find_peak'] + n_detected
fake_sources = vstack([fake_sources, _fake])
with ProgressBar(len(sourceslist)) as pg:
for _fake, _detected in zip(fakesourceslist, sourceslist):
n_fake = len(fake_sources)
n_detected = len(sources)
if _detected is not None:
_detected['ID'] = _detected['ID'] + n_detected
if 'fake_sources' in _detected.keys():
_detected['fake_sources'] = _detected['fake_sources'] + n_fake
sources = vstack([sources, _detected])
if _fake is not None:
_fake['ID'] = _fake['ID'] + n_fake
if 'find_peak' in _fake.keys():
_fake['find_peak'] = _fake['find_peak'] + n_detected
fake_sources = vstack([fake_sources, _fake])
pg.update()
return sources, fake_sources
......