11588621403

Committed 30 Oct 2024 07:14AM UTC coverage: 66.568% (-0.005%) from 66.573%

Build # 11588621403

Build Type

Pull #1081

gihtub

Committed by

iprafols

Commit Message

linted new code

Pull Request Pull Request #1081: Raise error if both Z and Z_DLA columns in DLA catalog [bump minor]

Run Details

1416 of 2659 branches covered (53.25%)

Branch coverage included in aggregate %.

1 of 2 new or added lines in 1 file covered. (50.0%)

5 existing lines in 1 file now uncovered.

6680 of 9503 relevant lines covered (70.29%)

2.11 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

80.51

/py/picca/delta_extraction/masks/dla_mask.py

"""This module defines the classes DlaMask and Dla used in the
masking of DLAs"""
import logging

from astropy.table import Table
import fitsio
import numpy as np
from scipy.constants import (
    speed_of_light as SPEED_LIGHT,
    e as ELEMENTARY_CHARGE,
    epsilon_0 as EPSILON_0,
    m_p as PROTON_MASS,
    m_e as ELECTRON_MASS,
    k as BOLTZMAN_CONSTANT_K,
)
from scipy.special import voigt_profile

from picca.delta_extraction.astronomical_objects.forest import Forest
from picca.delta_extraction.errors import MaskError
from picca.delta_extraction.mask import Mask, accepted_options, defaults
from picca.delta_extraction.utils import (
    ABSORBER_IGM, update_accepted_options, update_default_options)

accepted_options = update_accepted_options(accepted_options, [
    "dla mask limit", "los_id name", "mask file", "filename"
])

defaults = update_default_options(
    defaults, {
        "dla mask limit": 0.8,
        "los_id name": "THING_ID",
    })

LAMBDA_LYA = float(ABSORBER_IGM["LYA"]) ## Lya wavelength [A]
LAMBDA_LYB = float(ABSORBER_IGM["LYB"]) ## Lyb wavelength [A]
OSCILLATOR_STRENGTH_LYA = 0.41641
OSCILLATOR_STRENGTH_LYB = 0.079142
GAMMA_LYA = 6.2648e08  # s^-1 damping constant
GAMMA_LYB = 1.6725e8  # s^-1 damping constant

def dla_profile(lambda_, z_abs, nhi):
    """Compute DLA profile

    Arguments
    ---------
    lambda_: array of floats
    Wavelength (in Angs)

    z_abs: float
    Redshift of the absorption

    nhi: float
    DLA column density in log10(cm^-2)
    """
    transmission = np.exp(
        -compute_tau(lambda_, z_abs, nhi, LAMBDA_LYA, OSCILLATOR_STRENGTH_LYA, GAMMA_LYA)
        -compute_tau(lambda_, z_abs, nhi, LAMBDA_LYB, OSCILLATOR_STRENGTH_LYB, GAMMA_LYB)
    )
    return transmission

# constants to compute the optical depth of the DLA absoprtion
GAS_TEMP = 5 * 1e4  # K

# precomputed factors to save time
# compared to equation 36 of Garnett et al 2017 there is a sqrt(2) missing
# this is because otherwise we need to divide this quantity by sqrt(2)
# when calling scipy.special.voigt
GAUSSIAN_BROADENING_B = np.sqrt(BOLTZMAN_CONSTANT_K * GAS_TEMP / PROTON_MASS)
# the 1e-10 appears as the wavelengths are given in Angstroms
LORENTZIAN_BROADENING_GAMMA_PREFACTOR = 1e-10 / (4 * np.pi)
TAU_PREFACTOR = (
    ELEMENTARY_CHARGE**2 * 1e-10 / ELECTRON_MASS / SPEED_LIGHT / 4 / EPSILON_0)

def compute_tau(lambda_, z_abs, log_nhi, lambda_t, oscillator_strength_f, gamma):
    r"""Compute the optical depth for DLA absorption.

    Tau is computed using equations 34 to 36 of Garnett et al. 2017. We add
    a factor 4pi\epsion_0 in the denominator of their equation 34 so that
    dimensions match. The equations we use are:

    \tau(\lambda, z_{abs}, N_{HI}) = N_{HI} \frac {e^{2} f\lambda_{t} }
        {4 \epsion_0 m_{e} c } \phi{\nu, b, \gamma}

    where e is the elementary charge, \lambda_{t} is the transition wavelength
    and f is the oscillator strength of the transition. The line profile
    function \phi is a Voigt profile, where \nu is ther elative velocity

    \nu = c ( \frac{ \lambda } { \lambda_{t}* (1+z_{DLA}) }  ) ,

    b / \sqrt{2} is the standard deviation of the Gaussian (Maxwellian)
    broadening contribution:

    b = \sqrt{ \frac{ 2kT }{ m_{p} } }

    and \gamma is the width of the Lorenztian broadening contribution:

    \gamma = \frac { \Gamma \lambda_{t} } { 4\pi }

    where \Gamma is a damping constant

    Arguments
    ---------
    lambda_: array of floats
    Wavelength (in Angs)

    z_abs: float
    Redshift of the absorption

    log_nhi: float
    DLA column density in log10(cm^-2)

    lambda_t: float
    Transition wavelength, in Angstroms, e.g. for Lya 1215.67

    oscillator_strength_f: float
    Oscillator strength, e.g. f = 0.41611 for Lya

    gamma: float
    Damping constant (in s^-1)

    Return
    ------
    tau: array of float
    The optical depth.
    """
    # compute broadenings for the voight profile
    relative_velocity_nu = SPEED_LIGHT * (lambda_ / (1 + z_abs) / lambda_t - 1)
    lorentzian_broadening_gamma = (
        LORENTZIAN_BROADENING_GAMMA_PREFACTOR * gamma * lambda_t)

    # convert column density to m^2
    nhi = 10**log_nhi * 1e4

    # the 1e-10 converts the wavelength from Angstroms to meters
    tau = TAU_PREFACTOR * nhi * oscillator_strength_f * lambda_t * voigt_profile(
        relative_velocity_nu, GAUSSIAN_BROADENING_B, lorentzian_broadening_gamma)

    return tau

class DlaMask(Mask):
    """Class to mask DLAs

    Methods
    -------
    __init__
    apply_mask

    Attributes
    ----------
    (see Mask in py/picca/delta_extraction/mask.py)

    dla_mask_limit: float
    Lower limit on the DLA transmission. Transmissions below this number are
    masked

    logger: logging.Logger
    Logger object

    mask: astropy.Table
    Table containing specific intervals of wavelength to be masked for DLAs
    """
    def __init__(self, config):
        """Initializes class instance.

        Arguments
        ---------
        config: configparser.SectionProxy
        Parsed options to initialize class

        Raise
        -----
        MaskError if there are missing variables
        MaskError if input file does not have extension DLACAT
        MaskError if input file does not have fields THING_ID, Z, NHI in extension
        DLACAT
        MaskError upon OsError when reading the mask file
        """
        self.logger = logging.getLogger(__name__)

        super().__init__(config)

        # first load the dla catalogue
        filename = config.get("filename")
        if filename is None:
            raise MaskError("Missing argument 'filename' required by DlaMask")

        los_id_name = config.get("los_id name")
        if los_id_name is None:
            raise MaskError(
                "Missing argument 'los_id name' required by DlaMask")

        self.logger.progress(f"Reading DLA catalog from: {filename}")

        accepted_zcolnames = ["Z_DLA", "Z"]
        z_colname = accepted_zcolnames[0]
        try:
            with fitsio.FITS(filename) as hdul:
                hdul_colnames = set(hdul["DLACAT"].get_colnames())
                z_colname = hdul_colnames.intersection(accepted_zcolnames)
                if not z_colname:
                    raise ValueError(f"Z colname has to be one of {', '.join(accepted_zcolnames)}")
                if len(z_colname)>1 :
                    raise ValueError(
                        f"Not clear which column should be used for the DLA redshift among {z_colname}. "
                        "Please remove or rename one of the columns from the DLA fits file.")
                z_colname = z_colname.pop()
                columns_list = [los_id_name, z_colname, "NHI"]
                cat = {col: hdul["DLACAT"][col][:] for col in columns_list}
        except OSError as error:
            raise MaskError(
                f"Error loading DlaMask. File {filename} does "
                "not have extension 'DLACAT'"
            ) from error
        except ValueError as error:
            aux = "', '".join(columns_list)
            raise MaskError(
                f"Error loading DlaMask. File {filename} does "
                f"not have fields '{aux}' in HDU 'DLACAT'"
            ) from error

        # group DLAs on the same line of sight together
        self.los_ids = {}
        for los_id in np.unique(cat[los_id_name]):
            w = los_id == cat[los_id_name]
            self.los_ids[los_id] = list(zip(cat[z_colname][w], cat['NHI'][w]))
        num_dlas = np.sum([len(los_id) for los_id in self.los_ids.values()])

        self.logger.progress(f'In catalog: {num_dlas} DLAs')
        self.logger.progress(f'In catalog: {len(self.los_ids)} forests have a DLA\n')

        # setup transmission limit
        # transmissions below this number are masked
        self.dla_mask_limit = config.getfloat("dla mask limit")
        if self.dla_mask_limit is None:
            raise MaskError("Missing argument 'dla mask limit' "
                            "required by DlaMask")

        # load mask
        mask_file = config.get("mask file")
        if mask_file is not None:
            try:
                self.mask = Table.read(mask_file,
                                       names=('type', 'wave_min', 'wave_max',
                                              'frame'),
                                       format='ascii')
                self.mask = self.mask['frame'] == 'RF_DLA'
            except (OSError, ValueError) as error:
                raise MaskError(
                    f"ERROR: Error while reading mask_file file {mask_file}"
                ) from error
        else:
            self.mask = Table(names=('type', 'wave_min', 'wave_max', 'frame'))

    def apply_mask(self, forest):
        """Apply the mask. The mask is done by removing the affected
        pixels from the arrays in Forest.mask_fields

        Arguments
        ---------
        forest: Forest
        A Forest instance to which the correction is applied

        Raise
        -----
        MaskError if Forest.wave_solution is not 'log'
        """
        lambda_ = 10**forest.log_lambda

        # load DLAs
        if self.los_ids.get(forest.los_id) is not None:
            dla_transmission = np.ones(len(lambda_))
            for (z_abs, nhi) in self.los_ids.get(forest.los_id):
                dla_transmission *= dla_profile(lambda_, z_abs,
                                                nhi)

            # find out which pixels to mask
            w = dla_transmission > self.dla_mask_limit
            if len(self.mask) > 0:
                for mask_range in self.mask:
                    for (z_abs, nhi) in self.los_ids.get(forest.los_id):
                        w &= ((lambda_ / (1. + z_abs) < mask_range['wave_min'])
                              | (lambda_ /
                                 (1. + z_abs) > mask_range['wave_max']))

            # do the actual masking
            forest.transmission_correction *= dla_transmission
            for param in Forest.mask_fields:
                self._masker(forest, param, w)

1	"""This module defines the classes DlaMask and Dla used in the
2	masking of DLAs"""
3	import logging	3✔
4
5	from astropy.table import Table	3✔
6	import fitsio	3✔
7	import numpy as np	3✔
8	from scipy.constants import (	3✔
9	speed_of_light as SPEED_LIGHT,
10	e as ELEMENTARY_CHARGE,
11	epsilon_0 as EPSILON_0,
12	m_p as PROTON_MASS,
13	m_e as ELECTRON_MASS,
14	k as BOLTZMAN_CONSTANT_K,
15	)
16	from scipy.special import voigt_profile	3✔
17
18	from picca.delta_extraction.astronomical_objects.forest import Forest	3✔
19	from picca.delta_extraction.errors import MaskError	3✔
20	from picca.delta_extraction.mask import Mask, accepted_options, defaults	3✔
21	from picca.delta_extraction.utils import (	3✔
22	ABSORBER_IGM, update_accepted_options, update_default_options)
23
24	accepted_options = update_accepted_options(accepted_options, [	3✔
25	"dla mask limit", "los_id name", "mask file", "filename"
26	])
27
28	defaults = update_default_options(	3✔
29	defaults, {
30	"dla mask limit": 0.8,
31	"los_id name": "THING_ID",
32	})
33
34	LAMBDA_LYA = float(ABSORBER_IGM["LYA"]) ## Lya wavelength [A]	3✔
35	LAMBDA_LYB = float(ABSORBER_IGM["LYB"]) ## Lyb wavelength [A]	3✔
36	OSCILLATOR_STRENGTH_LYA = 0.41641	3✔
37	OSCILLATOR_STRENGTH_LYB = 0.079142	3✔
38	GAMMA_LYA = 6.2648e08 # s^-1 damping constant	3✔
39	GAMMA_LYB = 1.6725e8 # s^-1 damping constant	3✔
40
41	def dla_profile(lambda_, z_abs, nhi):	3✔
42	"""Compute DLA profile
43
44	Arguments
45	---------
46	lambda_: array of floats
47	Wavelength (in Angs)
48
49	z_abs: float
50	Redshift of the absorption
51
52	nhi: float
53	DLA column density in log10(cm^-2)
54	"""
55	transmission = np.exp(	3✔
56	-compute_tau(lambda_, z_abs, nhi, LAMBDA_LYA, OSCILLATOR_STRENGTH_LYA, GAMMA_LYA)
57	-compute_tau(lambda_, z_abs, nhi, LAMBDA_LYB, OSCILLATOR_STRENGTH_LYB, GAMMA_LYB)
58	)
59	return transmission	3✔
60
61	# constants to compute the optical depth of the DLA absoprtion
62	GAS_TEMP = 5 * 1e4 # K	3✔
63
64	# precomputed factors to save time
65	# compared to equation 36 of Garnett et al 2017 there is a sqrt(2) missing
66	# this is because otherwise we need to divide this quantity by sqrt(2)
67	# when calling scipy.special.voigt
68	GAUSSIAN_BROADENING_B = np.sqrt(BOLTZMAN_CONSTANT_K * GAS_TEMP / PROTON_MASS)	3✔
69	# the 1e-10 appears as the wavelengths are given in Angstroms
70	LORENTZIAN_BROADENING_GAMMA_PREFACTOR = 1e-10 / (4 * np.pi)	3✔
71	TAU_PREFACTOR = (	3✔
72	ELEMENTARY_CHARGE*2 1e-10 / ELECTRON_MASS / SPEED_LIGHT / 4 / EPSILON_0)
73
74	def compute_tau(lambda_, z_abs, log_nhi, lambda_t, oscillator_strength_f, gamma):	3✔
75	r"""Compute the optical depth for DLA absorption.
76
77	Tau is computed using equations 34 to 36 of Garnett et al. 2017. We add
78	a factor 4pi\epsion_0 in the denominator of their equation 34 so that
79	dimensions match. The equations we use are:
80
81	\tau(\lambda, z_{abs}, N_{HI}) = N_{HI} \frac {e^{2} f\lambda_{t} }
82	{4 \epsion_0 m_{e} c } \phi{\nu, b, \gamma}
83
84	where e is the elementary charge, \lambda_{t} is the transition wavelength
85	and f is the oscillator strength of the transition. The line profile
86	function \phi is a Voigt profile, where \nu is ther elative velocity
87
88	\nu = c ( \frac{ \lambda } { \lambda_{t}* (1+z_{DLA}) } ) ,
89
90	b / \sqrt{2} is the standard deviation of the Gaussian (Maxwellian)
91	broadening contribution:
92
93	b = \sqrt{ \frac{ 2kT }{ m_{p} } }
94
95	and \gamma is the width of the Lorenztian broadening contribution:
96
97	\gamma = \frac { \Gamma \lambda_{t} } { 4\pi }
98
99	where \Gamma is a damping constant
100
101	Arguments
102	---------
103	lambda_: array of floats
104	Wavelength (in Angs)
105
106	z_abs: float
107	Redshift of the absorption
108
109	log_nhi: float
110	DLA column density in log10(cm^-2)
111
112	lambda_t: float
113	Transition wavelength, in Angstroms, e.g. for Lya 1215.67
114
115	oscillator_strength_f: float
116	Oscillator strength, e.g. f = 0.41611 for Lya
117
118	gamma: float
119	Damping constant (in s^-1)
120
121	Return
122	------
123	tau: array of float
124	The optical depth.
125	"""
126	# compute broadenings for the voight profile
127	relative_velocity_nu = SPEED_LIGHT * (lambda_ / (1 + z_abs) / lambda_t - 1)	3✔
128	lorentzian_broadening_gamma = (	3✔
129	LORENTZIAN_BROADENING_GAMMA_PREFACTOR * gamma * lambda_t)
130
131	# convert column density to m^2
132	nhi = 10*log_nhi 1e4	3✔
133
134	# the 1e-10 converts the wavelength from Angstroms to meters
135	tau = TAU_PREFACTOR * nhi * oscillator_strength_f * lambda_t * voigt_profile(	3✔
136	relative_velocity_nu, GAUSSIAN_BROADENING_B, lorentzian_broadening_gamma)
137
138	return tau	3✔
139
140	class DlaMask(Mask):	3✔
141	"""Class to mask DLAs
142
143	Methods
144	-------
145	__init__
146	apply_mask
147
148	Attributes
149	----------
150	(see Mask in py/picca/delta_extraction/mask.py)
151
152	dla_mask_limit: float
153	Lower limit on the DLA transmission. Transmissions below this number are
154	masked
155
156	logger: logging.Logger
157	Logger object
158
159	mask: astropy.Table
160	Table containing specific intervals of wavelength to be masked for DLAs
161	"""
162	def __init__(self, config):	3✔
163	"""Initializes class instance.
164
165	Arguments
166	---------
167	config: configparser.SectionProxy
168	Parsed options to initialize class
169
170	Raise
171	-----
172	MaskError if there are missing variables
173	MaskError if input file does not have extension DLACAT
174	MaskError if input file does not have fields THING_ID, Z, NHI in extension
175	DLACAT
176	MaskError upon OsError when reading the mask file
177	"""
178	self.logger = logging.getLogger(__name__)	3✔
179
180	super().__init__(config)	3✔
181
182	# first load the dla catalogue
183	filename = config.get("filename")	3✔
184	if filename is None:	3✔
185	raise MaskError("Missing argument 'filename' required by DlaMask")	3✔
186
187	los_id_name = config.get("los_id name")	3✔
188	if los_id_name is None:	3✔
189	raise MaskError(	3✔
190	"Missing argument 'los_id name' required by DlaMask")
191
192	self.logger.progress(f"Reading DLA catalog from: {filename}")	3✔
193
194	accepted_zcolnames = ["Z_DLA", "Z"]	3✔
195	z_colname = accepted_zcolnames[0]	3✔
196	try:	3✔
197	with fitsio.FITS(filename) as hdul:	3✔
198	hdul_colnames = set(hdul["DLACAT"].get_colnames())	3✔
199	z_colname = hdul_colnames.intersection(accepted_zcolnames)	3✔
200	if not z_colname:	3!
201	raise ValueError(f"Z colname has to be one of {', '.join(accepted_zcolnames)}")	×
202	if len(z_colname)>1 :	3!
NEW 203	raise ValueError(	×
204	f"Not clear which column should be used for the DLA redshift among {z_colname}. "
205	"Please remove or rename one of the columns from the DLA fits file.")
206	z_colname = z_colname.pop()	3✔
207	columns_list = [los_id_name, z_colname, "NHI"]	3✔
208	cat = {col: hdul["DLACAT"][col][:] for col in columns_list}	3✔
UNCOV 209	except OSError as error:	×
210	raise MaskError(	×
211	f"Error loading DlaMask. File {filename} does "
212	"not have extension 'DLACAT'"
213	) from error
UNCOV 214	except ValueError as error:	×
215	aux = "', '".join(columns_list)	×
216	raise MaskError(	×
217	f"Error loading DlaMask. File {filename} does "
218	f"not have fields '{aux}' in HDU 'DLACAT'"
219	) from error
220
221	# group DLAs on the same line of sight together
222	self.los_ids = {}	3✔
223	for los_id in np.unique(cat[los_id_name]):	3✔
224	w = los_id == cat[los_id_name]	3✔
225	self.los_ids[los_id] = list(zip(cat[z_colname][w], cat['NHI'][w]))	3✔
226	num_dlas = np.sum([len(los_id) for los_id in self.los_ids.values()])	3✔
227
228	self.logger.progress(f'In catalog: {num_dlas} DLAs')	3✔
229	self.logger.progress(f'In catalog: {len(self.los_ids)} forests have a DLA\n')	3✔
230
231	# setup transmission limit
232	# transmissions below this number are masked
233	self.dla_mask_limit = config.getfloat("dla mask limit")	3✔
234	if self.dla_mask_limit is None:	3✔
235	raise MaskError("Missing argument 'dla mask limit' "	3✔
236	"required by DlaMask")
237
238	# load mask
239	mask_file = config.get("mask file")	3✔
240	if mask_file is not None:	3!
UNCOV 241	try:	×
242	self.mask = Table.read(mask_file,	×
243	names=('type', 'wave_min', 'wave_max',
244	'frame'),
245	format='ascii')
UNCOV 246	self.mask = self.mask['frame'] == 'RF_DLA'	×
247	except (OSError, ValueError) as error:	×
248	raise MaskError(	×
249	f"ERROR: Error while reading mask_file file {mask_file}"
250	) from error
251	else:
252	self.mask = Table(names=('type', 'wave_min', 'wave_max', 'frame'))	3✔
253
254	def apply_mask(self, forest):	3✔
255	"""Apply the mask. The mask is done by removing the affected
256	pixels from the arrays in Forest.mask_fields
257
258	Arguments
259	---------
260	forest: Forest
261	A Forest instance to which the correction is applied
262
263	Raise
264	-----
265	MaskError if Forest.wave_solution is not 'log'
266	"""
267	lambda_ = 10**forest.log_lambda	3✔
268
269	# load DLAs
270	if self.los_ids.get(forest.los_id) is not None:	3✔
271	dla_transmission = np.ones(len(lambda_))	3✔
272	for (z_abs, nhi) in self.los_ids.get(forest.los_id):	3✔
273	dla_transmission *= dla_profile(lambda_, z_abs,	3✔
274	nhi)
275
276	# find out which pixels to mask
277	w = dla_transmission > self.dla_mask_limit	3✔
278	if len(self.mask) > 0:	3!
UNCOV 279	for mask_range in self.mask:	×
280	for (z_abs, nhi) in self.los_ids.get(forest.los_id):	×
281	w &= ((lambda_ / (1. + z_abs) < mask_range['wave_min'])	×
282	\| (lambda_ /
283	(1. + z_abs) > mask_range['wave_max']))
284
285	# do the actual masking
286	forest.transmission_correction *= dla_transmission	3✔
287	for param in Forest.mask_fields:	3✔
288	self._masker(forest, param, w)	3✔

igmhub / picca / 11588621403

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous