Tutorial: Fitting a Cube

Fitting a cube is a little more complicated than fitting a single spectrum, but has been set up to be reasonably flexible and simple to interface with. First, we need to load in the config files

import tomllib

config_file = "config.toml"
local_file = "local.toml"

 with open(config_file, "rb") as f:
     config = tomllib.load(f)

 with open(local_file, "rb") as f:
     local = tomllib.load(f)

Next, a few filenames for various output files. We need the input cube, as well as files for the output parameter maps, fit parameter dictionaries, and fit cube. We also have an error cube and a mask cube that we’ll load in here, but these are not strictly necessary. We also set up some directories for the individual fits to go into at this point.

target = "g316_75"

fit_dict_filename = f"{target}_fit_dict"
mcfine_output_filename = f"{target}_mcfine_output"
fit_cube_filename = f"{target}_fit_cube"
hdu_name = f"{target}.fits"

original_fit_dir = "original"
coherence_forward_dir = "coherence_forward"
coherence_backward_dir = "coherence_backward"

err_hdu_name = hdu_name.replace(".fits", "_noise.fits")
mask_hdu_name = hdu_name.replace(".fits", "_mask.fits")

# Pull out various data we need from HDUs
with fits.open(hdu_name) as hdu, fits.open(err_hdu_name) as err_hdu, fits.open(mask_hdu_name) as mask_hdu:

   data = hdu[0].data
   err = err_hdu[0].data
   mask = mask_hdu[0].data

We’ll set the data up to be read properly into mcfine now. This involves generating error maps and a mask of pixels to fit:

# Add in a calibration uncertainty of 10%
calibration_uncertainty = 0.1
err = np.sqrt(err ** 2 + (calibration_uncertainty * data) ** 2)

# We'll use a pretty unrestrictive mask, which is just any spaxel that has a pixel included in the strict mask
mask = np.nansum(mask, axis=0)
mask[mask < 1] = 0
mask[mask >= 1] = 1

nan_mask = np.where(mask == 0)

We can now throw this all into mcfine! We’ll pass the data here as the path to the .fits file, as that will tell mcfine to load this in using spectral cube (and also extract the velocity axis). Note that when you do this, output maps will be saved as a dictionary of fits files, complete with WCS information.

from mcfine import HyperfineFitter

hf_fitter = HyperfineFitter(data=hdu_name,
                            error=err,
                            mask=mask,
                            config_file=config_file,
                            local_file=local_file,
                            )

We start with a first pass through, fitting all pixels defined by our mask:

print("First-pass fitting")
hf_fitter.multicomponent_fitter(fit_dict_filename=os.path.join(original_fit_dir, fit_dict_filename),
                                mcfine_output_filename=os.path.join(original_fit_dir, mcfine_output_filename),
                                likelihood_filename=os.path.join(original_fit_dir, likelihood_filename),
                                )

This will take a while if you have a lot of fits to do! Go and enjoy your weekend. After this is done, we will perform a coherence pass forwards and backwards. This has the effect of removing potentially bad fits by comparing with neighbours, but typically will only replace 10% or less of the fits

print("Spatial coherence forwards")
hf_fitter.encourage_spatial_coherence(fit_dict_filename=fit_dict_filename,
                                      input_dir=original_fit_dir,
                                      output_dir=coherence_forward_dir,
                                      mcfine_input_filename=os.path.join(original_fit_dir, mcfine_output_filename),
                                      mcfine_output_filename=os.path.join(coherence_forward_dir, mcfine_output_filename),
                                      )
print("Spatial coherence backwards")
hf_fitter.encourage_spatial_coherence(fit_dict_filename=fit_dict_filename,
                                      input_dir=coherence_forward_dir,
                                      output_dir=coherence_backward_dir,
                                      mcfine_input_filename=os.path.join(coherence_forward_dir, mcfine_output_filename),
                                      mcfine_output_filename=os.path.join(coherence_backward_dir, mcfine_output_filename),
                                      reverse_direction=True,
                                      )

Following this, fitting is complete! We will now generate parameter maps and the fit cube

print("Creating maps")
hf_fitter.make_parameter_maps(mcfine_output_filename=os.path.join(coherence_backward_dir, mcfine_output_filename),
                              maps_filename=f"{target}_maps.pkl",
                              )

print("Creating fit cubes")
hf_fitter.create_fit_cube(mcfine_output_filename=os.path.join(coherence_backward_dir, mcfine_output_filename),
                          cube_filename=fit_cube_filename)

This is now everything done! We can now explore the cube in exploring cube fits