How LUCI Works

LUCI is a general purpose line fitting pipeline designed to unveil the inner workings of how fits of SITELLE IFU datacubes are calculated. In this section, I will explain what we are trying to calculate, how it is calculated, and how you can personalize the fits.

What we calculate

The three primary quantities of interest are the amplitude of the line, the position of the line (often described as the velocity and quoted in km/s), and the broadening of the line (often described as the velocity dispersion).

Velocity

The velocity of a line is calculated using the following equation: $v [km/s] = c [km/s] * \Delta \lambda$

c is the speed of light in kilometers per second. Delta lambda is the shift of the measured line. Although the line position is calculated in units of [cm-1], we translate it into nanometers since $\lambda [nm] = \frac{1e7}{\lambda[cm-1]}$ . $\Delta \lambda = (line\_pos[nm]-line\_ref[nm])/line\_ref[nm]$

where $line\_ref$ is the natural position of the line (for example; $line\_ref=656.28 nm$ for Halpha.)

Velocity Dispersion

The velocity dispersion of a line is calculated using the following equation: \Delta v = \frac{3e5 [km/s] * \sigma}{v [km/s]}

where $\sigma$ is the calculated width of a the fitted Gaussian.

Flux

Similarly, we define the flux for each fitting function as the following:

Flux for a Gaussian Function:

$Flux [erg/s/cm^2/Ang] = \sqrt{2\pi}p_0p_2$

Flux for a Sinc Function:

$Flux [erg/s/cm^2/Ang] = \pi p_0p_2$

Flux for a SincGauss Function:

$Flux [erg/s/cm^2/Ang] = \text{coeff} * p_0\frac{\sqrt{2\pi}p_2}{erf(\frac{p_2}{\sqrt{2}\sigma})}$

Note that $text{coeff}=frac{1.20671}{pi*FWHM_COEFF}$ where the FWHM_COEFF equals $2sqrt{2log{2}}, the $pi$ is there because of the sinc function’s definition, and the 1.20671 is the factor used to go between FWHM and sigma.

Flux calibration

Every expression above is an amplitude times a width, so the flux inherits whatever units the data cube is in. The erg/s/cm^2/Ang written above is only true if the cube itself is flux calibrated.

That is not always the case. ORB’s newer level 3 cubes are stored in counts, with the conversion held in the header rather than applied to the data. ORB documents the level scheme in HDFCube.get_level:

LUCI handles this for you. On reading a cube she checks the flux_calibration header keyword (or BUNIT on the older DR1-style headers) and, if the data are in counts, multiplies the cube by flambda / step_nb / exposure_time before anything else happens. The flux, amplitude and continuum maps are then genuinely in erg/cm^2/s/Ang.

cube = SitelleCube(cube_path=cube_path, output_dir=output_dir,
                   object_name=object_name, redshift=redshift, resolution=resolution)
cube.flux_calibrated   # True once the conversion has been applied

If you would rather fit the raw counts, turn it off:

cube = SitelleCube(..., flux_calibration=False)

LUCI then logs a warning so that the units of the resulting maps are not a surprise. She also warns if a cube says it is uncalibrated but carries no usable flambda, in which case the maps stay in counts and the axis labels are wrong.

Note

ORB is not self-consistent here: HDFCube.to_fits multiplies by flambda alone, leaving out the /dimz/exposure_time, which for a typical SN4 cube is a factor of 17822. LUCI follows the get_level docstring, which is the form that reproduces published surface brightnesses.

How we calculate

Once we have a spectrum, we do two things: we normalize the spectrum by the maximum amplitude and we apply a redshift correction (wavelength = wavelength*(1+redshift)). We do this primarily to constrain the velocity to be between -500 and 500 km/s. This allows our machine learning technique to obtain better initial guess estimates.

Initial Guess

Having a good initial guess is crucial for the success (and speed) of the fitting algorithm. In order to obtain a good initial guess for the fit parameters (line amplitude, line position, and line broadening), we apply a machine learning technique described in Rhea et al. 2020a (disclosure: the author of this code is also the author of this paper). The method uses pre-trained convolutional neural networks to estimate the velocity and broadening of the line in km/s. These are then translated into the line position and broadening. Next, the amplitude is taken as the height of the line corresponding to the shifted line position. We note that the machine learning model has only been trained to obtain velocities between -500 and 500 km/s. Similarly, the model was trained to obtain broadening values between 10 and 200 km/s. You can find more information on this at https://sitelle-signals.github.io/Pamplemousse/index.html <https://sitelle-signals.github.io/Pamplemousse/index.html>. We estimate the amplitude by taking the maximum value of spectrum corresponding to the estimated position plus or minus 2 channels.

The networks are specific to a filter and a resolution, paired with a reference spectrum Reference-Spectrum-R<resolution>-<FILTER>.fits that defines the axis every spectrum gets interpolated onto. If you need one for a filter or resolution that LUCI does not ship, see Adding a New Filter – ML/TrainPredictor.py will build the reference spectrum, generate the synthetic training set, and train the network for you.

How inference runs

Inference runs on ONNX Runtime, not TensorFlow. The networks were trained in Keras, but LUCI ships them converted to ONNX in ML/onnx/ as R<resolution>-PREDICTOR-I-<FILTER>.onnx (and ...-PREDICTOR-I-MDN-<FILTER>.onnx for the mixture-density variants). Fitting therefore does not import TensorFlow at all.

Three things follow from this, all of which matter in practice:

Speed. The predictor and its ONNX session are loaded once per process and reused for every spectrum. Previously the Keras model was re-read from disk for each pixel, which dominated the runtime of a full cube (measured: 5.7 s/pixel before, 3.0 s/pixel after, the remainder being the actual fit).
Install weight. ONNX artifacts are roughly a third the size of the Keras SavedModels, and onnxruntime is a far lighter dependency than TensorFlow.
Fidelity. Every converted model is checked against its Keras original on 200 spectra before being shipped; all 39 agree to float32 precision (worst deviation 6.7e-4 km/s). Any model that fails that gate is reported and not published.

If you retrain a network, convert it with:

uv run tools/convert_models_to_onnx.py --all --validate

That script pins its own legacy-TensorFlow environment via a PEP 723 header, so it builds what it needs on the fly and your project environment never sees TensorFlow.

Internally the predictors sit behind a small interface (LUCI.ml.ParameterPredictor), so the fitting code asks for a velocity/broadening estimate without knowing or caring which backend answers.

If machine learning is not your cup of tea

Pass ML_bool=False and LUCI estimates the initial guesses directly from the data: it locates the brightest peak in the fit window to seed the velocity and uses a default broadening, so the optimiser always starts somewhere sensible. The same fallback kicks in automatically if no trained network exists for your filter/resolution combination – LUCI tells you it is doing so rather than failing.

You can also supply your own priors explicitly with initial_values=[velocity, broadening], which freezes those parameters instead of fitting them.

Fitting Function

The fitting function utilizes scipy.optimize.minimize. Currently, we are using the trust-constr <https://docs.scipy.org/doc/scipy/reference/optimize.minimize-trustconstr.html> optimization algorithm. Before fitting the spectrum, we normalize the spectrum by the maximum amplitude – this makes the fitting process simpler. The fit returns the amplitude of the line (which we then scale to be correct for the un-normalized spectrum), the velocity in km/s, and the velocity dispersion in km/s. If the user chooses, the line velocities and velocity dispersions can be coupled. Additionally, we automatically include the constraint on the NII-doublet flux ratio (setting NII_6583 = 3*NII_6548) using the nii_cons boolean. This can be changed by adding nii_cons=False as an argument to any of the fitting functions.

Available Models

We have implemented three functions: gaussian, sinc, and sincgauss.

We assume a standard form of a Gaussian: $p_0*exp{(-(x-p_1)**2)/(2*p_2**2)}$

We solve for p_0, p_1, and p_2 (x is the wavelength channel and is thus provided). $p_0$ is the amplitude, $p_1$ is the position of the line, and $p_2$ is the broadening.

We adopt the following form $u=(x-p_1)/p_2$ $p_0*sinc(\pi*u)/(\pi*u)$

Note that $p_2$ is FIXED for the sinc function as 1/(2*MPD) (where MPD is the maximum path difference).

$p_0*exp(-b*^2)*((erf(a-i*b)+erf(a+i*b))/(2*erf(a)))$

where $a = p_2/(\sqrt{2}*\sigma)$ $b = (x-p_1)/(\sqrt{2}*\sigma)$

We also have the Dawson integral calculation of the sincgauss function:

$p_0*( (dawson(ai+b)exp(2iab) + dawson(ai-b)exp(-2iab)) / 2*dawson(ai) )$

where sigma is 1/(2*MPD).

Therefore, when using a sincgauss, we have to calculate the MPD. We can adopt the following definition: $MPD = \cos{\theta}\delta_x N$ where $\cos{\theta}$ is the cosine angle defined as $\cos{\theta} = \frac{\lambda_{ref}}{\lambda_{ij}}$ . $\lambda_{ref}$ is the wavelength of the calibration laser and $\lambda_{ij}$ is the measured calibration wavelength of a given pixel (thus $\theta$ is a function of the pixel).

Also note that we divide the sinc width ( $p_2$ ) by $\pi$ based on our definition of the sinc width above.

If you are interested in the broadening, we strongly suggest you use the sincgauss function :)

Transmission

We take into account the transmission of the SITTELLE filters (SN1, SN2, SN3, SN4, C3, and C4). We take the true transmission as the mean of the transmission at different filter angles; the raw data can be found [here](https://www.cfht.hawaii.edu/Instruments/Sitelle/SITELLE_filters.php). The transmission is then applied to the spectrum in the following manner: if the transmission is above 0.5, then we multiply the spectrum by the transmission percentage. Otherwise, we set it to zero. Note that we calculate the noise before applying the transmission.