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.


ToolBase

The base class of all tools. Tools have a built-in plotter and 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.

Methods

run()

Runs the tool.

run()[source]

Runs the tool. Each subclass Tool must implement its own run() method.

H2Excitation

Tool for fitting temperatures in \(H_2\) excitation diagrams. A two temperature model is assumed, \(T_{warm}\) and \(T_{cold}\).

class pdrtpy.tool.h2excitation.H2Excitation(measurements=None)[source]

Bases: pdrtpy.tool.toolbase.ToolBase

Tool for fitting temperatures to \(H_2\) Excitation Diagrams

This tool is still under development

Attributes
intensities

The stored intensities.

Methods

add_measurement(m)

Add an intensity Measurement to internal dictionary used to compute the excitation diagram.

average_column_density(norm, x, y, xsize, …)

Compute the average column density over a spatial box.

colden(intensity)

Compute the column density in upper state \(N_u\), given an intensity \(I\) and assuming optically thin emission.

column_densities([norm])

The computed upper state column densities of stored intensities

energies([line])

Upper state energies of stored intensities, in K.

fit_excitation(**kwargs)

Fit the \(log N_u-E\) diagram with two excitation temperatures, a warm \(T_{ex}\) and a cold \(T_{ex}\).

replace_measurement(m)

Safely replace an existing intensity Measurement.

run()

Runs the tool.

excitation_diagram

plot_column_densities

plot_intensities

x_lin

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}}\)

average_column_density(norm, x, y, xsize, ysize, line)[source]

Compute the average column density over a spatial box.

Parameters
  • norm (bool) – if True, normalize the column densities by the statistical weight of the upper state, \(g_u\).

  • x (int) – bottom left corner x

  • y (int) – bottom left corner y

  • xsize (int) – box width, pixels

  • ysize (int) – box height, pixels

  • line (bool) – if True, the returned dictionary index is the Line name, otherwise it is the upper state \(J\) number.

Returns

dictionary of column densities

Return type

dict

colden(intensity)[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

m – A Measurement instance containing intensity in units equivalent to \({\rm erg~cm^{-2}~s^{-1}~sr^{-1}}\)

Returns

a Measurement of the column density.

Return type

Measurement

column_densities(norm=False)[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

Returns

dictionary of column densities indexed by Line name

Return type

dict

energies(line=False)[source]

Upper state energies of stored intensities, in K.

Parameters

line (bool) – if True, the dictionary index is the Line name, otherwise it is the upper state \(J\) number. Default: False

Returns

dictionary indexed by upper state \(J\) number or Line name.

Return type

dict

excitation_diagram(**kwargs)[source]
fit_excitation(**kwargs)[source]

Fit the \(log N_u-E\) diagram with two excitation temperatures, a warm \(T_{ex}\) and a cold \(T_{ex}\). A first pass guess is initially made using data partitioning and two linear fits.

Returns

The fit parameters as in scipy.optimize.curve_fit

Return type

list

property intensities

The stored intensities. See add_measurement()

Return type

list of Measurement

plot_column_densities(**kwargs)[source]
plot_intensities(**kwargs)[source]
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}}\)

x_lin(m1, n1, m2, n2)[source]

LineRatioFit

Tool for determining photodissociation region external radiation field and density (commonly known 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: pdrtpy.tool.toolbase.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 ModelSet 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.

Attributes
density

The computed hydrogen nucleus density value(s).

has_maps

Are the Measurements used map-based or pixel-based?.

measurementIDs

The stored measurement IDs, which are strings.

measurements

The stored measurements as a dictionary with Measurement IDs as keys

modelset

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 measurements

Methods

add_measurement(m)

Add a Measurement to internal dictionary used to compute ratios.

chisq([min])

The computed chisquare value(s).

compute_density_radiation_field()

Compute the best-fit density and radiation field spatial maps by searching for the minimum chi-squared at each spatial pixel.

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).

remove_measurement(id)

Delete a measurement from the internal dictionary used to compute ratios.

run([mask])

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

Measurement

compute_density_radiation_field()[source]

Compute the best-fit density and radiation field spatial maps by searching for the minimum chi-squared at each spatial pixel.

property density

The computed hydrogen nucleus density value(s).

Return type

Measurement

property has_maps

Are the Measurements used map-based or pixel-based?.

Returns

True, if the observational inputs are spatial maps, False if they were single-pixel Measurements

Return type

bool

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 modelset

The underlying ModelSet

property observed_ratios

The list of the observed line ratios that have been input so far.

Return type

list of str

property radiation_field

The computed radiation field value(s).

Return type

Measurement

property ratiocount

The number of ratios that match models available in the current ModelSet given the current set of measurements

Return type

int

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

Measurement

remove_measurement(id)[source]

Delete a measurement from the internal dictionary used to compute ratios.

Parameters

id (str) – the measurement identifier

Raises

KeyError – if id not in existing Measurements

run(mask=['mad', 1.0])[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

[‘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

type mask

list or None

raises Exception

if no models match the input observations, observations are not compatible, or on unrecognized parameters

write_chisq(chi='chisq.fits', rchi='rchisq.fits', overwrite=True)[source]

Write the chisq and reduced-chisq data to a file

Parameters
  • chi (str) – FITS file to write the chisq map to.

  • rchi (str) – FITS file to write the reduced chisq map to.