Analysis Tools: Fit models to data¶
The tool
module contains the analysis tools in the PDR Toolbox. All tools are derived from ToolBase
.
For examples how to use LineRatioFit, see the notebooks PDRT_Example_Find_n_G0_Single_Pixel.ipynb and PDRT_Example_Make_n_G0_maps.ipynb.
For an example how to use H2ExcitationFit and ExcitationPlot see the notebook PDRT_Example_H2_Excitation.ipynb.
ToolBase¶
The base class of all tools. Tools have a run()
method both of which subclasses must define.
- class pdrtpy.tool.toolbase.ToolBase[source]¶
Bases:
object
Base class object for PDR Toolbox tools. This class implements a simple interface with a run method. Tools will generally do some set up such as reading in observational data before run() can be invoked.
- Attributes:
has_maps
Are the Measurements used map-based?.
has_scalar
Are the Measurements used scalars.
has_vectors
Are the Measurements used a Nx1 vector, e.g.
Methods
run
()Runs the tool.
- property has_maps¶
Are the Measurements used map-based?. (i.e., have 2 spatial axes)
- Returns:
True, if the observational inputs are spatial maps, False otherwise
- Return type:
- property has_scalar¶
Are the Measurements used scalars.
- Returns:
True, if the observational inputs are scalars, False otherwise
- Return type:
- property has_vectors¶
Are the Measurements used a Nx1 vector, e.g. read in from a table with
from_table()
.- Returns:
True, if the observational inputs are a vector, False otherwise
- Return type:
Excitation Diagram Fitting¶
H2ExcitationFit
is a tool for fitting temperature, column density, and ortho-to-para ratio in \(H_2\) excitation diagrams. A two temperature model is assumed, and the fit will find \(T_{hot}, T_{cold}, N_{hot}(H_2), N_{cold}(H_2),\) and optionally OPR. The base class ExcitationFit
can be used to create a tool to fit a different molecule.
- class pdrtpy.tool.h2excitation.ExcitationFit(measurements=None, constantsfile=None)[source]¶
Bases:
ToolBase
Base class for creating excitation fitting tools for various species.
- Parameters:
measurements (list of
Measurement
.) – Input measurements to be fit.- Attributes:
has_maps
Are the Measurements used map-based?.
has_scalar
Are the Measurements used scalars.
has_vectors
Are the Measurements used a Nx1 vector, e.g.
Methods
Add an intensity Measurement to internal dictionary used to compute the excitation diagram.
remove_measurement
(identifier)Delete a measurement from the internal dictionary used to compute column densities.
Safely replace an existing intensity Measurement.
run
()Runs the tool.
- add_measurement(m)[source]¶
Add an intensity Measurement to internal dictionary used to compute the excitation diagram. This method can also be used to safely replace an existing intensity Measurement.
- Parameters:
m – A
Measurement
instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\)
- remove_measurement(identifier)[source]¶
Delete a measurement from the internal dictionary used to compute column densities. Any associated column density will also be removed.
- replace_measurement(m)[source]¶
Safely replace an existing intensity Measurement. Do not change a Measurement in place, use this method. Otherwise, the column densities will be inconsistent.
- Parameters:
m – A
Measurement
instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\)
- class pdrtpy.tool.h2excitation.H2ExcitationFit(measurements=None, constantsfile='RoueffEtAl.tab')[source]¶
Bases:
ExcitationFit
Tool for fitting temperatures, column densities, and ortho-to-para ratio(OPR) from an \(H_2\) excitation diagram. It takes as input a set of \(H_2\) rovibrational line observations with errors represented as
Measurement
.Often, excitation diagrams show evidence of both “hot” and “cold” gas components, where the cold gas dominates the intensity in the low J transitions and the hot gas dominates in the high J transitions. Given data over several transitions, one can fit for \(T_{cold}, T_{hot}, N_{total} = N_{cold}+ N_{hot}\), and optionally OPR. One needs at least 5 points to fit the temperatures and column densities (slope and intercept \(\times 2\)), though one could compute (not fit) them with only 4 points. To additionally fit OPR, one should have 6 points (5 degrees of freedom).
Once the fit is done,
ExcitationPlot
can be used to view the results.- Parameters:
measurements (list of
Measurement
.) – Input \(H_2\) measurements to be fit.- Attributes:
cold_colden
The fitted cold gas total column density
fit_result
The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables.
has_maps
Are the Measurements used map-based?.
has_scalar
Are the Measurements used scalars.
has_vectors
Are the Measurements used a Nx1 vector, e.g.
hot_colden
The fitted hot gas total column density
intensities
The stored intensities.
numcomponents
Number of temperature components in the fit
opr
The ortho-to-para ratio (OPR)
opr_fitted
Was the ortho-to-para ratio fitted?
tcold
The fitted cold gas excitation temperature
temperature
The fitted gas temperatures, returned in a dictionary with keys ‘hot’ and ‘cold’.
thot
The fitted hot gas excitation temperature
total_colden
The fitted total column density
Methods
add_measurement
(m)Add an intensity Measurement to internal dictionary used to compute the excitation diagram.
average_column_density
([position, size, ...])Compute the average column density over a spatial box.
colden
(component)The column density of hot or cold gas component, or total column density.
column_densities
([norm, unit, line])The computed upper state column densities of stored intensities
energies
([line])Upper state energies of stored intensities, in K.
gu
(id, opr)Get the upper state statistical weight $g_u$ for the given transition identifer, and, if the transition is odd-$J$, scale the result by the given ortho-to-para ratio.
intensity
(colden)Given an upper state column density \(N_u\), compute the intensity \(I\).
remove_measurement
(identifier)Delete a measurement from the internal dictionary used to compute column densities.
replace_measurement
(m)Safely replace an existing intensity Measurement.
run
([position, size, fit_opr])Fit the \(log N_u-E\) diagram with two excitation temperatures, a
hot
\(T_{ex}\) and acold
\(T_{ex}\).upper_colden
(intensity, unit)Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.
- average_column_density(position=None, size=None, norm=True, unit=Unit("1 / cm2"), line=True, clip=<Quantity -1.e+40 1 / cm2>)[source]¶
Compute the average column density over a spatial box. The box is created using
astropy.nddata.utils.Cutout2D
.- Parameters:
position (tuple) – The position of the cutout array’s center with respect to the data array. The position can be specified either as a (x, y) tuple of pixel coordinates.
size (int, array_like`) – The size of the cutout array along each axis. If size is a scalar number or a scalar
Quantity
, then a square cutout of size will be created. If size has two elements, they should be in (nx,ny) order [this is the opposite of Cutout2D signature]. Scalar numbers in size are assumed to be in units of pixels. Default value of None means use all pixels (position is ignored)norm (bool) – if True, normalize the column densities by the statistical weight of the upper state, \(g_u\). For ortho-$H_2$ $g_u = OPR times (2J+1)$, for para-$H_2$ $g_u=2J+1$. In LTE, $OPR = 3$.
unit (str or
astropy.units.Unit
) – The units in which to return the column density. Default: \({\rm cm}^{-2}\)line (bool) – if True, the returned dictionary index is the Line name, otherwise it is the upper state \(J\) number.
clip (
astropy.units.Quantity
) – Column density value at which to clip pixels. Pixels with column densities below this value will not be used in the average. Default: a large negative number, which translates to no clipping.
- Returns:
dictionary of column density Measurements, with keys as \(J\) number or Line name
- Return type:
- property cold_colden¶
The fitted cold gas total column density
- Return type:
- colden(component)[source]¶
The column density of hot or cold gas component, or total column density.
- Parameters:
component (str) – ‘hot’, ‘cold’, or ‘total
- Return type:
- column_densities(norm=False, unit=Unit('1 / cm2'), line=True)[source]¶
The computed upper state column densities of stored intensities
- Parameters:
norm (bool) – if True, normalize the column densities by the statistical weight of the upper state, \(g_u\). Default: False
unit (str or
astropy.units.Unit
) – The units in which to return the column density. Default: \({\\rm }cm^{-2}\)line (bool) – if True, the dictionary index is the Line name, otherwise it is the upper state \(J\) number. Default: False
- Returns:
dictionary of column densities indexed by upper state \(J\) number or Line name. Default: False means return indexed by \(J\).
- Return type:
- property fit_result¶
The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables.
- Return type:
- gu(id, opr)[source]¶
Get the upper state statistical weight $g_u$ for the given transition identifer, and, if the transition is odd-$J$, scale the result by the given ortho-to-para ratio. If the transition is even-$J$, the LTE value is returned.
- property hot_colden¶
The fitted hot gas total column density
- Return type:
- property intensities¶
The stored intensities. See
add_measurement()
- Return type:
list of
Measurement
- intensity(colden)[source]¶
Given an upper state column density \(N_u\), compute the intensity \(I\).
\[I = {A \Delta E~N_u \over 4\pi}\]where \(A\) is the Einstein A coefficient and \(\Delta E\) is the energy of the transition.
- Parameters:
colden (
Measurement
) – upper state column density- Returns:
optically thin intensity
- Return type:
- property opr¶
The ortho-to-para ratio (OPR)
- Returns:
The fitted OPR is it was determined in the fit, otherwise the canonical LTE OPR
- Return type:
- property opr_fitted¶
Was the ortho-to-para ratio fitted?
- Returns:
True if OPR was fitted, False if canonical LTE value was used
- Return type:
- run(position=None, size=None, fit_opr=False, **kwargs)[source]¶
Fit the \(log N_u-E\) diagram with two excitation temperatures, a
hot
\(T_{ex}\) and acold
\(T_{ex}\).If
position
andsize
are given, the data will be averaged over a spatial box before fitting. The box is created usingastropy.nddata.utils.Cutout2D
. If position or size is None, the data are averaged over all pixels. If the Measurements are single values, these arguments are ignored.- Parameters:
position (tuple) – The position of the cutout array’s center with respect to the data array. The position can be specified either as a (x, y) tuple of pixel coordinates.
size (int, array_like`) – The size of the cutout array along each axis in pixels. If size is a scalar number or a scalar
Quantity
, then a square cutout of size will be created. If size has two elements, they should be in (nx, ny) order [this is the opposite of Cutout2D signature]. Scalar numbers in size are assumed to be in units of pixels. Default value of None means use all pixels (position is ignored)fit_opr (bool) – Whether to fit the ortho-to-para ratio or not. If True, the OPR will be varied to determine the best value. If False, the OPR is fixed at the canonical LTE value of 3.
- property tcold¶
The fitted cold gas excitation temperature
- Return type:
- property temperature¶
The fitted gas temperatures, returned in a dictionary with keys ‘hot’ and ‘cold’. :rtype: dict
- property thot¶
The fitted hot gas excitation temperature
- Return type:
- property total_colden¶
The fitted total column density
- Return type:
- upper_colden(intensity, unit)[source]¶
Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission. Units of \(I\) need to be equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\).
\[ \begin{align}\begin{aligned}I &= {A \Delta E~N_u \over 4\pi}\\N_u &= 4\pi {I\over A\Delta E}\end{aligned}\end{align} \]where \(A\) is the Einstein A coefficient and \(\Delta E\) is the energy of the transition.
- Parameters:
intensity (
Measurement
) – AMeasurement
instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\)unit (str or
astropy.units.Unit
) – The units in which to return the column density. Default: \({\rm }cm^{-2}\)
- Returns:
a
Measurement
of the column density.- Return type:
LineRatioFit¶
LineRatioFit
is a tool for determining photodissociation region external radiation field and hydrogen nucleus density (commonly given as \(G_0\) and \(n\)) from measured spectral line intensity ratios.
- class pdrtpy.tool.lineratiofit.LineRatioFit(modelset=<pdrtpy.modelset.ModelSet object>, measurements=None)[source]¶
Bases:
ToolBase
LineRatioFit is a tool to fit observations of intensity ratios to a set of PDR models. It takes as input a set of observations with errors represented as
Measurement
andModelSet
for the models to which the data will be fitted. The observations should be spectral line or continuum intensities. They can be spatial maps or single pixel values. They should have the same spatial resolution.The models to be fit are stored as intensity ratios. The input observations will be use to create ratios that correspond to models. From there a minimization fit is done to determine the density and radiation field that best fit the data. At least 3 observations are needed in order to make at least 2 ratios. With fewer ratios, no fitting can be done. More ratios generally means better determined density and radiation field, assuming the data are consistent with each other.
Once the fit is done,
LineRatioPlot
can be used to view the results.- Parameters:
modelset (
ModelSet
) – The set of PDR models to use for fitting.measurements (list or dict of
Measurement
. If dict, the keys should be the Measurement identifiers.) – Input measurements to be fit.
- Attributes:
density
The computed hydrogen nucleus density value(s).
fit_result
The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables.
has_maps
Are the Measurements used map-based?.
has_scalar
Are the Measurements used scalars.
has_vectors
Are the Measurements used a Nx1 vector, e.g.
measurementIDs
The stored measurement IDs, which are strings.
measurements
The stored
measurements
as a dictionary with Measurement IDs as keysmodelset
The underlying
ModelSet
observed_ratios
The list of the observed line ratios that have been input so far.
radiation_field
The computed radiation field value(s).
ratiocount
The number of ratios that match models available in the current
ModelSet
given the current set of measurementstable
Construct the table of input Measurements, and if the fit has been run, the density, radiation field, and \(\chi^2\) values
Methods
Add a Measurement to internal dictionary used to compute ratios.
chisq
([min])The computed chisquare value(s).
read_models
([unit])Given a list of measurement IDs, find and open the FITS files that have matching ratios and populate the _modelratios dictionary.
reduced_chisq
([min])The computed reduced chisquare value(s).
Delete a measurement from the internal dictionary used to compute ratios.
run
(**kwargs)Run the full computation using all the
observations
added.write_chisq
([chi, rchi, overwrite])Write the chisq and reduced-chisq data to a file
- add_measurement(m)[source]¶
Add a Measurement to internal dictionary used to compute ratios. This measurement may be intensity units (erg \({\rm s}^{-1}\) \({\rm cm}^{-2}\)) or integrated intensity (K km/s).
- Parameters:
m (
Measurement
.) – a Measurement instance to be added to this tool
- chisq(min=False)[source]¶
The computed chisquare value(s).
- Parameters:
min (bool) – If True return the minimum reduced \(\chi^2\). In the case of map inputs this will be a spatial map of mininum \(\chi^2\). If False with map inputs the entire \(\chi^2\) hypercube is returned. If True with single pixel inputs, a single value is returned. If False with single pixel inputs, \(\chi^2\) as a function of density and radiation field is returned.
- Return type:
- property density¶
The computed hydrogen nucleus density value(s).
- Return type:
- property fit_result¶
The result of the fitting procedure which includes fit statistics, variable values and uncertainties, and correlations between variables. For each pixel, there will be one instance of
lmfit.minimizer.MinimizerResult
.- Return type:
FitMap
- property measurementIDs¶
The stored measurement IDs, which are strings.
- Return type:
dict_keys
- property measurements¶
The stored
measurements
as a dictionary with Measurement IDs as keys- Return type:
dict of
Measurement
- property observed_ratios¶
The list of the observed line ratios that have been input so far.
- property radiation_field¶
The computed radiation field value(s).
- Return type:
- property ratiocount¶
The number of ratios that match models available in the current
ModelSet
given the current set of measurements- Return type:
- read_models(unit=Unit(dimensionless))[source]¶
Given a list of measurement IDs, find and open the FITS files that have matching ratios and populate the _modelratios dictionary. Uses
pdrtpy.measurement.Measurement
as a storage mechanism.- param m:
list of measurement IDS (string)
- type m:
list
- param unit:
units of the data
- type unit:
string or astropy.Unit
- reduced_chisq(min=False)[source]¶
The computed reduced chisquare value(s).
- Parameters:
min (bool) – If True return the minimum reduced \(\chi_\nu^2\). In the case of map inputs this will be a spatial map of mininum \(\chi_\nu^2\). If False with map inputs the entire \(\chi_\nu^2\) hypercube is returned. If True with single pixel inputs, a single value is returned. If False with single pixel inputs, \(\chi_\nu^2\) as a function of density and radiation field is returned.
- Return type:
- remove_measurement(id)[source]¶
Delete a measurement from the internal dictionary used to compute ratios.
- run(**kwargs)[source]¶
Run the full computation using all the
observations
added. This will check compatibility of input observations (e.g., beam parameters, coordinate types, axes lengths) and raise exceptions if the observations don’t match each other.- param mask:
- Indicate how to mask image observations (Measurements) before computing the density
and radiation field. Possible values are:
- [‘mad’, multiplier] - compute standard deviation using median absolute deviation (astropy.mad_std),
and mask out values between +/- multiplier*mad_std. Example: [‘mad’,1.0]
[‘data’, (low,high)] - mask based on data values, mask out data between low and high
[‘clip’, (low,high)] - mask based on data values, mask out data below low and above high
- [‘error’, (low,high)] - mask based on uncertainty plane, mask out data where the corresponding error pixel value
is below low or above high
None - Don’t mask data. This is the default.
- type mask:
list or None
- param method:
the fitting method to be used. The default is ‘leastsq’, which is Levenberg-Marquardt least squares. For other options see https://lmfit-py.readthedocs.io/en/latest/fitting.html#fit-methods-table
- type method:
str
- param nan_policy:
Specifies action if fit returns NaN values. One of: * ’raise’ : a ValueError is raised [Default] * ’propagate’ : the values returned from userfcn are un-altered * ’omit’ : non-finite values are filtered
- type nan_policy:
str
- raises Exception:
if no models match the input observations, observations are not compatible, or on unrecognized parameters, or NaN encountered.
- property table¶
Construct the table of input Measurements, and if the fit has been run, the density, radiation field, and \(\chi^2\) values
- Return type:
FitMap¶
When fitting either single pixels or spatial maps, the fit results are stored per pixel in an NDData
object that contains ModelResult
objects for H2ExcitationFit
or MinimizerResult
objects for LineRatioFit
. The user can thus examine in detail the fit at any pixel.
- class pdrtpy.tool.fitmap.FitMap(data, *args, **kwargs)[source]¶
Bases:
NDData
- Attributes:
data
~numpy.ndarray-like : The stored dataset.
mask
any type : Mask for the dataset, if any.
- meta
dict-like : Additional meta information about the dataset.
name
The name of this FitMap
psf
Image representation of the PSF for the dataset.
uncertainty
any type : Uncertainty in the dataset, if any.
unit
~astropy.units.Unit : Unit for the dataset, if any.
wcs
any type : A world coordinate system (WCS) for the dataset, if any.
Methods
get_pixel
(world_x, world_y)Return the nearest pixel coordinates to the input world coordinates
- get_pixel(world_x, world_y)[source]¶
Return the nearest pixel coordinates to the input world coordinates
- property name¶
The name of this FitMap :rtype: str