arguments cleanup and plot_absorber() updated

data/datamodel.rst

@@ -41,7 +41,7 @@ Run `qsoabsfind` with the required FITS file. If using a custom constant file, i
 
     qsoabsfind --input <path_to_input_fits_file> [--constant-file <path_to_constant_file>] --output <path_to_output_fits_file>
 
-Replace `<path_to_input_fits_file>` with the path to your input FITS file, `<path_to_constant_file>` with the path to your constant file (if using), and `<path_to_output_fits_file>` with the desired path for your output FITS file. For a quick example run you can run the module on `data/qso_test.fits`. 
+Replace `<path_to_input_fits_file>` with the path to your input FITS file, `<path_to_constant_file>` with the path to your constant file (if using), and `<path_to_output_fits_file>` with the desired path for your output FITS file. For a quick example run you can run the module on `data/qso_test.fits`.
 
 Output FITS File Structure
 --------------------------

qsoabsfind/absfinder.py

@@ -78,7 +78,7 @@ def read_single_spectrum_and_find_absorber(fits_file, spec_index, absorber, **kw
                          and METADATA which must contain keyword Z_QSO.
         spec_index (int): Index of the quasar spectrum to retrieve from the FITS file.
         absorber (str): Name of the absorber to search for (e.g., 'MgII', 'CIV').
-        kwargs (dictionary): search parameters used in convolution_method_absorber_finder_in_QSO_spectra()
+        kwargs (dictionary): search parameters as defined in constants.py
 
     Returns:
         tuple: Contains lists of various parameters related to detected absorbers.
@@ -107,7 +107,7 @@ def read_single_spectrum_and_find_absorber(fits_file, spec_index, absorber, **kw
     lam_obs = spectra.wavelength
 
     # Define the wavelength range for searching the absorber
-    min_wave, max_wave = lam_obs.min() + 200, lam_obs.max() - 200  # avoiding edges
+    min_wave, max_wave = lam_obs.min() + kwargs["lam_edge_sep"], lam_obs.max() - kwargs["lam_edge_sep"]  # avoiding edges
 
     # Retrieve flux and error data, ensuring consistent dtype for Numba compatibility
     residual, error = spectra.flux.astype('float64'), spectra.error.astype('float64')
@@ -124,6 +124,8 @@ def read_single_spectrum_and_find_absorber(fits_file, spec_index, absorber, **kw
     # Verify that the arrays are of equal size
     assert lam_search.size == unmsk_residual.size == unmsk_error.size, "Mismatch in array sizes of lam_search, unmsk_residual, and unmsk_error"
 
+    kwargs.pop("lam_edge_sep") # just remove this keyword as its not used the following function.
+
     (index_spec, pure_z_abs, pure_gauss_fit, pure_gauss_fit_std, pure_ew_first_line_mean, pure_ew_second_line_mean, pure_ew_total_mean, pure_ew_first_line_error, pure_ew_second_line_error, pure_ew_total_error, redshift_err, sn1_all, sn2_all) = convolution_method_absorber_finder_in_QSO_spectra(spec_index, absorber, lam_obs, residual, error, lam_search, unmsk_residual, unmsk_error, **kwargs)
 
     # Print progress for every spectrum processed

qsoabsfind/constants.py

@@ -24,6 +24,7 @@
 ker_width_pixels = [3, 4, 5, 6, 7, 8]  # Pixel size for Gaussian convolution
 pm_pixel = 200  # Pixel size for error calculation to define the threshold for potential absorber features
 mult_resi = 1  # Coefficients to multiply the residual before Gaussian search
+lam_sep = 300 # wavelength shift from edges (in Ang), i.e., lam_min + lam_sep, lam_max - lam_sep, to avoid noisy edges of the spectrum
 
 search_parameters = {
     'MgII': {
@@ -37,6 +38,7 @@
         'use_covariance':False,
         'logwave':True,
         'verbose':False,
+        'lam_edge_sep':lam_sep,
     },
     'CIV': {
         'ker_width_pixels': ker_width_pixels,
@@ -49,6 +51,7 @@
         'use_covariance':False,
         'logwave':True,
         'verbose':False,
+        'lam_edge_sep':lam_sep,
     }
 }
 

qsoabsfind/parallel_convolution.py

@@ -70,7 +70,7 @@ def parse_qso_sequence(qso_sequence):
     # If none of the conditions matched, raise an error
     raise ValueError(f"Invalid QSO sequence format: '{qso_sequence}'. Use 'start-end[:step]' or an integer.")
 
-def parallel_convolution_method_absorber_finder_QSO_spectra(fits_file, spec_indices, absorber='MgII', ker_width_pixels=[3, 4, 5, 6, 7, 8], coeff_sigma=2.5, mult_resi=1, d_pix=0.6, pm_pixel=200, sn_line1=3, sn_line2=2, use_covariance=False, logwave=True, verbose=False, n_jobs=1):
+def parallel_convolution_method_absorber_finder_QSO_spectra(fits_file, spec_indices, absorber, n_jobs, **kwargs):
     """
     Run convolution_method_absorber_finder_in_QSO_spectra in parallel using
     multiprocessing.
@@ -79,23 +79,14 @@ def parallel_convolution_method_absorber_finder_QSO_spectra(fits_file, spec_indi
         fits_file (str): Path to the FITS file containing Normalized QSO spectra.
         spec_indices (list or numpy.array): Indices of quasars in the data matrix.
         absorber (str): Absorber name for searching doublets (MgII, CIV). Default is 'MgII'.
-        ker_width_pixels (list): List of kernel widths in pixels. Default is [3, 4, 5, 6, 7, 8].
-        coeff_sigma (float): Coefficient for sigma to apply threshold in the convolved array. Default is 2.5.
-        mult_resi (float): Factor to shift the residual up or down. Default is 1.
-        d_pix (float): Pixel distance for line separation during Gaussian fitting. Default is 0.6.
-        pm_pixel (int): Pixel parameter for local noise estimation (default 200).
-        sn_line1 (float): Signal-to-noise ratio for thresholding for line1.
-        sn_line2 (float): Signal-to-noise ratio for thresholding for line2.
-        use_covariance (bool): If want to use full covariance of scipy curve_fit for EW error calculation (default is False).
-        logwave (bool): If wavelength on logscale (default True for SDSS).
-        verbose (bool): if True will print lots of output for debugging
         n_jobs (int): Number of parallel jobs to run.
+        kwargs (dictionary): search parameters as defined in constants.py
 
     Returns:
         dict: A dictionary containing combined results from all parallel runs.
     """
 
-    kwargs = {'ker_width_pixels': ker_width_pixels, 'coeff_sigma': coeff_sigma, 'mult_resi': mult_resi, 'd_pix': d_pix, 'pm_pixel': pm_pixel, 'sn_line1': sn_line1, 'sn_line2': sn_line2, 'use_covariance': use_covariance, 'logwave': logwave, 'verbose': verbose}
+    #kwargs = {'ker_width_pixels': ker_width_pixels, 'coeff_sigma': coeff_sigma, 'mult_resi': mult_resi, 'd_pix': d_pix, 'pm_pixel': pm_pixel, 'sn_line1': sn_line1, 'sn_line2': sn_line2, 'use_covariance': use_covariance, 'logwave': logwave, 'verbose': verbose}
 
     params_list = [(fits_file, spec_index, absorber, kwargs) for spec_index in spec_indices]
 
@@ -174,21 +165,21 @@ def main():
     # Add search parameters and package versions to headers
     headers.update({
         'ABSORBER': {"value": args.absorber, "comment": 'Absorber name'},
-        'KERWIDTH': {"value": str(constants.search_parameters[args.absorber]["ker_width_pixels"]), "comment": 'Kernel width in pixels'},
-        'COEFFSIG': {"value": constants.search_parameters[args.absorber]["coeff_sigma"], "comment": 'Coefficient for sigma threshold'},
-        'MULTRE': {"value": constants.search_parameters[args.absorber]["mult_resi"], "comment": 'Multiplicative factor for residuals'},
-        'D_PIX': {"value": constants.search_parameters[args.absorber]["d_pix"], "comment": 'toloerance for line separation (in Ang)'},
-        'PM_PIXEL': {"value": constants.search_parameters[args.absorber]["pm_pixel"], "comment": 'Pixel parameter for local noise estimation'},
-        'SN_LINE1': {"value": constants.search_parameters[args.absorber]["sn_line1"], "comment": 'S/N threshold for first line'},
-        'SN_LINE2': {"value": constants.search_parameters[args.absorber]["sn_line2"], "comment": 'S/N threshold for second line'},
-        'EWCOVAR': {"value": constants.search_parameters[args.absorber]["use_covariance"], "comment": 'Use covariance for EW error calculation'},
-        'LOGWAVE': {"value": constants.search_parameters[args.absorber]["logwave"], "comment": 'Use log wavelength scaling'}
+        'KERWIDTH': {"value": str(constants.search_parameters[args.absorber]["ker_width_pixels"]), "comment": 'Kernel width in pixels (ker_width_pixels)'},
+        'COEFFSIG': {"value": constants.search_parameters[args.absorber]["coeff_sigma"], "comment": 'Coefficient for sigma threshold (coeff_sigma)'},
+        'MULTRE': {"value": constants.search_parameters[args.absorber]["mult_resi"], "comment": 'Multiplicative factor for residuals (mult_resi)'},
+        'D_PIX': {"value": constants.search_parameters[args.absorber]["d_pix"], "comment": 'tolerance for line separation (in Ang) (d_pix)'},
+        'PM_PIXEL': {"value": constants.search_parameters[args.absorber]["pm_pixel"], "comment": 'Pixel parameter for local noise estimation (pm_pixel)'},
+        'SN_LINE1': {"value": constants.search_parameters[args.absorber]["sn_line1"], "comment": 'S/N threshold for first line (sn_line1)'},
+        'SN_LINE2': {"value": constants.search_parameters[args.absorber]["sn_line2"], "comment": 'S/N threshold for second line (sn_line2)'},
+        'EWCOVAR': {"value": constants.search_parameters[args.absorber]["use_covariance"], "comment": 'Use covariance for EW error calculation (use_covariance)'},
+        'LOGWAVE': {"value": constants.search_parameters[args.absorber]["logwave"], "comment": 'Use log wavelength scaling (logwave)'}, 'LAM_ESEP': {"value": constants.search_parameters[args.absorber]["lam_edge_sep"], "comment": 'wavelength edges to avoid noisy regions (edges+/-lam_edge_sep)'}
     })
 
     package_versions = get_package_versions()
     for pkg, ver in package_versions.items():
         headers[pkg.upper()] = {"value": ver, "comment": f'{pkg} version'}
-    headers['QSABFI'] = headers.pop('QSOABSFIND')
+    headers['QSOABFI'] = headers.pop('QSOABSFIND')
     headers['MATPLOT'] = headers.pop('MATPLOTLIB')
 
     # Start timing
@@ -204,15 +195,7 @@ def main():
     # Run the convolution method in parallel
     results = parallel_convolution_method_absorber_finder_QSO_spectra(
         args.input_fits_file, spec_indices, absorber=args.absorber,
-        ker_width_pixels=constants.search_parameters[args.absorber]["ker_width_pixels"],
-        coeff_sigma=constants.search_parameters[args.absorber]["coeff_sigma"],
-        mult_resi=constants.search_parameters[args.absorber]["mult_resi"],
-        d_pix=constants.search_parameters[args.absorber]["d_pix"],
-        pm_pixel=constants.search_parameters[args.absorber]["pm_pixel"],
-        sn_line1=constants.search_parameters[args.absorber]["sn_line1"],
-        sn_line2=constants.search_parameters[args.absorber]["sn_line2"],
-        use_covariance=constants.search_parameters[args.absorber]["use_covariance"],
-        logwave=constants.search_parameters[args.absorber]["logwave"], verbose=constants.search_parameters[args.absorber]["verbose"], n_jobs=args.n_tasks * args.ncpus
+        n_jobs=args.n_tasks * args.ncpus, **constants.search_parameters[args.absorber]
     )
 
     # Save the results to a FITS file

qsoabsfind/utils.py

@@ -278,21 +278,23 @@ def vel_dispersion(c1, c2, sigma1, sigma2, resolution):
 
     return corr_del_v1_sq, corr_del_v2_sq
 
-def plot_absorber(lam, residual, absorber, zabs, xlabel='obs wave (ang)', ylabel='residual', title='QSO', plot_filename=None):
+def plot_absorber(spectra, absorber, zabs, show_error=False, xlabel='obs wave (ang)', ylabel='residual', title='QSO', plot_filename=None):
     """
     Saves a plot of spectra with absorber(s) (full spectrum + zoomed version)
     in the current working directory.
 
     Args:
-        lam (array-like): Observed wavelength.
-        residual (array-like): Residual flux.
+        spectra (object): spectra class, output of QSOSpecRead()
         absorber (str): Type of absorber, e.g., 'MgII', 'CIV'.
         zabs (list, array, or Table): Absorber redshifts, or a Table with 'Z_ABS' and 'GAUSS_FIT' columns.
+        show_error (bool): if error bars should be shown (default False)
         xlabel (str): The label for the x-axis. Default is 'obs wave (ang)'.
         ylabel (str): The label for the y-axis. Default is 'residual'.
         title (str): The super title of the plot. Default is 'QSO'.
         plot_filename (str): If provided, will save the plot to the given filename.
     """
+
+    lam, residual, error = spectra.wavelength, spectra.flux, spectra.error
     # If zabs is a Table or structured array, extract redshifts and fit parameters
     if isinstance(zabs, (Table, np.ndarray)) and ('Z_ABS' in zabs.colnames or 'Z_ABS' in zabs.dtype.names):
         redshifts = zabs['Z_ABS']
@@ -307,16 +309,19 @@ def plot_absorber(lam, residual, absorber, zabs, xlabel='obs wave (ang)', ylabel
     num_absorbers = len(redshifts)
 
     # Create a grid with 2 rows: 1 for the main plot and 1 for zoomed plots
-    fig = plt.figure(figsize=(6 * num_absorbers, 8))
+    fig = plt.figure(figsize=(13.5, 8))
     fig.subplots_adjust(hspace=0.15, wspace=0.15)  # Adjust space between plots
 
     # Super title for the entire figure
     fig.suptitle(title, fontsize=16)
 
     # Create the main plot in the first row
     ax_main = plt.subplot2grid((2, num_absorbers), (0, 0), colspan=num_absorbers)
-    ax_main.plot(lam, residual, ls='-', lw=1.5)
+    ax_main.plot(lam, residual, ls='-', lw=1.5, label='residual')
+    if show_error:
+        ax_main.plot(lam, error, ls='-', lw=1.5, label='error')
     ax_main.set_xlim(3800, 9200)
+    ax_main.legend()
 
     # Determine the absorber line labels
     if absorber == 'MgII':
@@ -345,8 +350,10 @@ def plot_absorber(lam, residual, absorber, zabs, xlabel='obs wave (ang)', ylabel
         ax_zoom = plt.subplot2grid((2, num_absorbers), (1, idx))
         x1, x2 = lines[l1] * shift_z, lines[l2] * shift_z
         mask = (lam > x1 - sep) & (lam < x2 + sep)  # Define zoom range around the lines
-
-        ax_zoom.plot(lam[mask], residual[mask], ls='-', lw=1.5, label='data')
+        if not show_error:
+            ax_zoom.plot(lam[mask], residual[mask], ls='-', lw=1.5, label='data')
+        else:
+            ax_zoom.errorbar(lam[mask], residual[mask], yerr=error[mask], marker='o', color='C0', markersize=6, label='data')
         ax_zoom.axvline(x=x1, color='r', ls='--')
         ax_zoom.axvline(x=x2, color='r', ls='--')
         ax_zoom.set_xlim([x1 - sep, x2 + sep])