Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

change background estimation to gaussian mean+factor*sigma

Browse Source
This commit is contained in:
Thibault Barnouin
2024-03-11 17:25:32 +01:00
parent 38604aa82d
commit 0d8e3a8c47
3 changed files with 14 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
src/FOC_reduction.py
View File

@@ -33,13 +33,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Background estimation
error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51))
subtract_error = 1.20
subtract_error = 1.00
display_bkg = False
display_error = False
# Data binning
rebin = True
pxsize = 0.10
pxsize = 0.05
px_scale = 'arcsec' # pixel, arcsec or full
rebin_operation = 'sum' # sum or average
@@ -50,7 +49,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Smoothing
smoothing_function = 'combine' # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
smoothing_FWHM = 0.2 # If None, no smoothing is done
smoothing_FWHM = 0.10 # If None, no smoothing is done
smoothing_scale = 'arcsec' # pixel or arcsec
# Rotation
@@ -90,7 +89,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
plots_folder = "."
if not path_exists(plots_folder):
system("mkdir -p {0:s} ".format(plots_folder))
infiles = [p[1] for p in prod] # if p[1] not in ['x2rp0202t_c0f.fits', 'x2rp0302t_c0f.fits']]
infiles = [p[1] for p in prod]
data_array, headers = proj_fits.get_obs_data(infiles, data_folder=data_folder, compute_flux=True)
figname = "_".join([target, "FOC"])
@@ -117,13 +116,9 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# Estimate error from data background, estimated from sub-image of desired sub_shape.
background = None
data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_type=error_sub_type, subtract_error=subtract_error, display=display_error, savename="_".join([
data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_type=error_sub_type, subtract_error=subtract_error, display=display_bkg, savename="_".join([
figname, "errors"]), plots_folder=plots_folder, return_background=True)
if display_bkg:
proj_plots.plot_obs(data_array, headers, vmin=data_array[data_array > 0.].min(
)*headers[0]['photflam'], vmax=data_array[data_array > 0.].max()*headers[0]['photflam'], savename="_".join([figname, "bkg"]), plots_folder=plots_folder)
# Align and rescale images with oversampling.
data_array, error_array, headers, data_mask = proj_red.align_data(
data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=False)

16
src/lib/background.py
View File

@@ -157,9 +157,9 @@ def bkg_estimate(img, bins=None, chi2=None, coeff=None):
hist, bin_edges = np.histogram(img[img > 0], bins=bins[-1])
binning = bin_centers(bin_edges)
peak = binning[np.argmax(hist)]
bins_fwhm = binning[hist > hist.max()/2.]
fwhm = bins_fwhm[-1]-bins_fwhm[0]
p0 = [hist.max(), peak, fwhm, 1e-3, 1e-3, 1e-3, 1e-3]
bins_stdev = binning[hist > hist.max()/2.]
stdev = bins_stdev[-1]-bins_stdev[0]
p0 = [hist.max(), peak, stdev, 1e-3, 1e-3, 1e-3, 1e-3]
try:
popt, pcov = curve_fit(gausspol, binning, hist, p0=p0)
except RuntimeError:
@@ -231,7 +231,7 @@ def bkg_fit(data, error, mask, headers, subtract_error=True, display=False, save
weights = 1/chi2**2
weights /= weights.sum()
bkg = np.sum(weights*coeff[:, 1])*subtract_error if subtract_error > 0 else np.sum(weights*coeff[:, 1])
bkg = np.sum(weights*(coeff[:, 1]+np.abs(coeff[:, 2])*subtract_error))
error_bkg[i] *= bkg
@@ -332,12 +332,12 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis
# bkg = np.sqrt(np.sum(image[np.abs(image-hist_max)/hist_max<0.5]**2)/image[np.abs(image-hist_max)/hist_max<0.5].size)
# Fit a gaussian to the log-intensity histogram
bins_fwhm = binning[-1][hist > hist.max()/2.]
fwhm = bins_fwhm[-1]-bins_fwhm[0]
p0 = [hist.max(), binning[-1][np.argmax(hist)], fwhm, 1e-3, 1e-3, 1e-3, 1e-3]
bins_stdev = binning[-1][hist > hist.max()/2.]
stdev = bins_stdev[-1]-bins_stdev[0]
p0 = [hist.max(), binning[-1][np.argmax(hist)], stdev, 1e-3, 1e-3, 1e-3, 1e-3]
popt, pcov = curve_fit(gausspol, binning[-1], hist, p0=p0)
coeff.append(popt)
bkg = popt[1]*subtract_error if subtract_error > 0 else popt[1]
bkg = popt[1]+np.abs(popt[2])*subtract_error
error_bkg[i] *= bkg

2
src/lib/plots.py
View File

@@ -1584,7 +1584,7 @@ class slit(object):
self.angle = angle
self.rect_center = (self.x0, self.y0)-np.dot(rot2D(self.angle), (self.width/2, self.height/2))
self.rect = Rectangle(self.rect_center, self.width, self.height, alpha=0.8, ec='grey', fc='none')
self.rect = Rectangle(self.rect_center, self.width, self.height, angle=self.angle, alpha=0.8, ec='grey', fc='none')
self.ax.add_patch(self.rect)
self.fig.canvas.mpl_connect('button_press_event', self.on_press)