convolution function modified

qsoabsfind/absfinder.py

@@ -91,7 +91,7 @@ def read_single_spectrum_and_find_absorber(fits_file, spec_index, absorber, **kw
 
 def convolution_method_absorber_finder_in_QSO_spectra(spec_index, absorber='MgII', lam_obs=None, residual=None, error=None,
 lam_search=None, unmsk_residual=None, unmsk_error=None, 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, resolution=69, logwave=True, verbose=False):
+mult_resi=1, d_pix=0.6, pm_pixel=200, sn_line1=3, sn_line2=2, use_covariance=False, resolution=69, logwave=True, wave_res=0.0001, verbose=False):
     """
     Detect absorbers with doublet properties in SDSS quasar spectra using a
     convolution method. This function identifies potential absorbers based on
@@ -118,6 +118,7 @@ def convolution_method_absorber_finder_in_QSO_spectra(spec_index, absorber='MgII
         use_covariance (bool): if want to use full covariance of scipy curvey_fit for EW error calculation (default is False)
         resolution (float): wavelength resolution of spectrum (in km/s), e.g. SDSS: ~69, DESI: ~70 (also defined in constants)
         logwave (bool): if wavelength on log scale (default True for SDSS)
+        wave_res (float): wavelength pixel size (SDSS: 0.0001 on log scale, DESI: 0.8 on linear scale)
         verbose (bool): if want to print a lot of outputs for debugging (default False)
 
     Returns:
@@ -134,45 +135,50 @@ def convolution_method_absorber_finder_in_QSO_spectra(spec_index, absorber='MgII
             - errors EW total (list): errors on Total Equivalent width of line 1 and line 2 for each absorber
     """
 
-    # Constants
-    if absorber == 'MgII':
-        line1, line2 = lines['MgII_2796'], lines['MgII_2803']
-        f1, f2 = oscillator_parameters['MgII_f1'], oscillator_parameters['MgII_f2']
-    elif absorber == 'CIV':
-        line1, line2 = lines['CIV_1548'], lines['CIV_1550']
-        f1, f2 = oscillator_parameters['CIV_f1'], oscillator_parameters['CIV_f2']
+    # return if there are no wavelength pixels available to search for
+    if lam_search.size==0 or lam_obs.size==0:
+        print(f'INFO: zero wavelength pixels available in search region, spec index = {spec_index}')
+        return ([spec_index], [0], [[0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0]], [0], [0], [0], [0], [0], [0], [0], [0], [0])
+
     else:
-        raise ValueError(f"No support for {absorber}, only supports MgII and CIV")
+        # Constants
+        if absorber == 'MgII':
+            line1, line2 = lines['MgII_2796'], lines['MgII_2803']
+            f1, f2 = oscillator_parameters['MgII_f1'], oscillator_parameters['MgII_f2']
+        elif absorber == 'CIV':
+            line1, line2 = lines['CIV_1548'], lines['CIV_1550']
+            f1, f2 = oscillator_parameters['CIV_f1'], oscillator_parameters['CIV_f2']
+        else:
+            raise ValueError(f"No support for {absorber}, only supports MgII and CIV")
 
-    line_sep = line2 - line1
+        line_sep = line2 - line1
 
-    if not logwave:
-        # average resolution in case wavelength is on linear scale
-        wave_pixel = np.nanmean(lam_search[1:] - lam_search[:-1])
-        resolution  = np.nanmean(wave_pixel/lam_search) * speed_of_light
+        if not logwave:
+            # average resolution in case wavelength is on linear scale
+            wave_pixel = np.nanmean(lam_search[1:] - lam_search[:-1])
+            resolution  = np.nanmean(wave_pixel/lam_search) * speed_of_light
 
-    del_sigma = line1 * resolution / speed_of_light  # in Ang
+        del_sigma = line1 * resolution / speed_of_light  # in Ang
 
-    bd_ct, x_sep = 2.5, 10 # multiple for bound definition (for line centres and widths of line)
+        bd_ct, x_sep = 2.5, 10 # multiple for bound definition (for line centres and widths of line)
 
-    # bounds for gaussian fitting, to avoid very bad candidates
-    bound = ((np.array([2e-2, line1 - bd_ct * d_pix, del_sigma-0.1, 2e-2, line2 - bd_ct * d_pix, del_sigma-0.1])),
-             (np.array([1.11, line1 + bd_ct * d_pix, x_sep * del_sigma+0.1, 1.11, line2 + bd_ct * d_pix, x_sep * del_sigma+0.1])))
+        # bounds for gaussian fitting, to avoid very bad candidates
+        bound = ((np.array([2e-2, line1 - bd_ct * d_pix, del_sigma-0.1, 2e-2, line2 - bd_ct * d_pix, del_sigma-0.1])),
+                 (np.array([1.11, line1 + bd_ct * d_pix, x_sep * del_sigma+0.1, 1.11, line2 + bd_ct * d_pix, x_sep * del_sigma+0.1])))
 
-    # line separation tolerance (fitted line centers should not be outside, centre +/- d_pix)
-    lower_del_lam = line_sep - d_pix
-    upper_del_lam = line_sep + d_pix
+        # line separation tolerance (fitted line centers should not be outside, centre +/- d_pix)
+        lower_del_lam = line_sep - d_pix
+        upper_del_lam = line_sep + d_pix
 
-    # Kernel width computation
-    width_kernel = np.array([ker * resolution * ((f1 * line1 + f2 * line2) / (f1 + f2)) / (speed_of_light * 2.35) for ker in ker_width_pixels])
+        # Kernel width computation
+        width_kernel = np.array([ker * resolution * ((f1 * line1 + f2 * line2) / (f1 + f2)) / (speed_of_light * 2.35) for ker in ker_width_pixels])
 
-    combined_final_our_z = []
+        combined_final_our_z = []
 
-    if len(unmsk_residual) > 0:
         for sig_ker in width_kernel:
             line_centre = (line1 + line2) / 2
 
-            conv_arr = convolution_fun(absorber, mult_resi * unmsk_residual, sig_ker, amp_ratio=0.5, log=logwave, index=spec_index)
+            conv_arr = convolution_fun(absorber, mult_resi * unmsk_residual, sig_ker, amp_ratio=0.5, log=logwave, wave_res=wave_res, index=spec_index)
             sigma_cr = estimate_local_sigma_conv_array(conv_arr, pm_pixel=pm_pixel)
             thr = np.nanmedian(conv_arr) - coeff_sigma * sigma_cr
 
@@ -305,5 +311,3 @@ def convolution_method_absorber_finder_in_QSO_spectra(spec_index, absorber='MgII
                     pure_ew_first_line_error.tolist(), pure_ew_second_line_error.tolist(), pure_ew_total_error.tolist(), redshift_err.tolist(), sn1_all.tolist(), sn2_all.tolist())
         else:
             return ([spec_index], [0], [[0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0]], [0], [0], [0], [0], [0], [0], [0], [0], [0])
-    else:
-        return ([spec_index], [0], [[0, 0, 0, 0, 0, 0]], [[0, 0, 0, 0, 0, 0]], [0], [0], [0], [0], [0], [0], [0], [0], [0])

qsoabsfind/constants.py

@@ -25,7 +25,10 @@
 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
-resolution = 69 # wavelength resolution in km/s for SDSS or DESI
+
+# Some NOTES:
+# SDSS spectra: use resolution 69 (km/s), wave_res = 0.0001, logwave=True
+# DESI spectra: use resolution None, wave_res = 0.8, logwave=False, average resolution will be calculated by code
 
 search_parameters = {
     'MgII': {
@@ -37,10 +40,11 @@
         'sn_line1':3,
         'sn_line2':2,
         'use_covariance':False,
-        'resolution':resolution,
+        'resolution':69,
         'logwave':True,
-        'verbose':False,
+        'wave_res':0.0001,
         'lam_edge_sep':lam_sep,
+        'verbose':False,
     },
     'CIV': {
         'ker_width_pixels': ker_width_pixels,
@@ -51,10 +55,11 @@
         'sn_line1':3,
         'sn_line2':2,
         'use_covariance':False,
-        'resolution':resolution,
+        'resolution':69,
         'logwave':True,
-        'verbose':False,
+        'wave_res':0.0001,
         'lam_edge_sep':lam_sep,
+        'verbose':False,
     }
 }
 

qsoabsfind/utils.py

@@ -96,7 +96,7 @@ def gauss_two_lines_kernel(x, a):
 
     return norm_constant * (-a1 * np.exp(-((x - a[1]) / a[2]) ** 2 / 2) - a2 * np.exp(-((x - a[4]) / a[5]) ** 2 / 2)) * 0.5
 
-def convolution_fun(absorber, residual_arr_after_mask, width, amp_ratio=0.5, log=True, index=None):
+def convolution_fun(absorber, residual_arr_after_mask, width, amp_ratio=0.5, log=True, wave_res=0.0001, index=None):
     """
     Convolves the spectrum with a Gaussian kernel.
 
@@ -106,6 +106,7 @@ def convolution_fun(absorber, residual_arr_after_mask, width, amp_ratio=0.5, log
         width (float): The width of the Gaussian kernel (decide base dupon width of real absorption feature).
         amp_ratio (float): Amplitude ratio for the Gaussian lines (default 0.5).
         log (bool): if log bins should be used for wavelength (dlam = 0.0001, default True)
+        wave_res (float): wavelength pixel size (SDSS: 0.0001 on log scale, DESI: 0.8 on linear scale)
         index (int): QSO index
 
     Returns:
@@ -116,22 +117,22 @@ def convolution_fun(absorber, residual_arr_after_mask, width, amp_ratio=0.5, log
 
     A_main = amplitude_dict[absorber]
     A_secondary = A_main * amp_ratio
-    ct = 5
+    ct = 10
     if absorber == 'MgII':
         ker_parm = np.array([A_main, lines['MgII_2796'], width, A_secondary, lines['MgII_2803'], width])
-        lam_ker_start = lines['MgII_2796']-ct*width # +/- 5sigma , #rest-frame
+        lam_ker_start = lines['MgII_2796']-ct*width # +/- 10sigma , #rest-frame
         lam_ker_end = lines['MgII_2803']+ct*width
     elif absorber == 'CIV':
         ker_parm = np.array([A_main, lines['CIV_1548'], width, A_secondary, lines['CIV_1550'], width])
-        lam_ker_start = lines['CIV_1548']-ct*width # +/- 5sigma , #rest-frame
+        lam_ker_start = lines['CIV_1548']-ct*width # +/- 10sigma , #rest-frame
         lam_ker_end = lines['CIV_1550']+ct*width #
     else:
         raise ValueError(f"Unsupported absorber type for specific Args: {absorber}")
     if log:
-        lam_ker = np.arange(np.log10(lam_ker_start), np.log10(lam_ker_end), 0.0001) #SDSS-like wavelength resolution
+        lam_ker = np.arange(np.log10(lam_ker_start), np.log10(lam_ker_end), wave_res) #SDSS-like wavelength resolution
         lam_ker = 10**lam_ker
     else:
-        lam_ker = np.arange(lam_ker_start, lam_ker_end, 0.8) # DESI-like wavelength resolution
+        lam_ker = np.arange(lam_ker_start, lam_ker_end, wave_res) # DESI-like wavelength resolution
 
     if len(lam_ker)>len(residual_arr_after_mask):
         lam_ker = lam_ker[0: len(residual_arr_after_mask)]
@@ -289,7 +290,7 @@ def combine_fits_files(directory, output_filename='combined.fits'):
 
     # Create HDU list for the combined file
     combined_hdul = fits.HDUList([primary_hdu, combined_absorber_hdu, combined_metadata_hdu])
-    
+
     # Define the full output file path
     output_file_path= output_filename
     if not os.path.isabs(output_file_path):