def _fit_helper(self, fit_spectra, reject_outliers, fit_lb, fit_ub,
                    fit_func):
        """
        This is a helper function for fitting different segments of the spectrum
        with Gaussian functions (GLX and GABA).

        Parameters
        ----------
        fit_spectra : ndarray
           The data to fit

        reject_outliers : float or bool
            Z score for outlier rejection. If set to `False`, not outlier
            rejection.

        fit_lb : float
            The lower bound of the part of the ppm scale for which the Gaussian
            is fit.

        fit_ub : float
            The upper bound of the part of the scale fit.

        fit_func: callable
           e.g. `fit_gaussian`

        Returns
        -------
        choose_transients : tuple
            Indices into the original data's transients dimension to select

        non-outlier transients. If reject_outliers is set to `False`, this is
            all the  transients

        model : ndarray
            The model predicition in each transient, based on the fit.

        signal : ndarray
            The original signal in this part of the difference spectrum.

        params : ndarray
            The Gaussian parameters in each transient as fit.

        this_idx : slice object
            A slice into the part of the spectrum that is fit
        """
        # fit_idx should already be set from fitting the creatine params:
        model, signal, params = fit_func(fit_spectra,
                                         self.f_ppm,
                                         lb=fit_lb,
                                         ub=fit_ub)

        # We'll use these indices to reject outliers (or not):
        ii = np.ones(signal.shape[0], dtype=bool)
        # Reject outliers:
        if reject_outliers:
            model, signal, params, ii = self._outlier_rejection(params,
                                                                model,
                                                                signal,
                                                                ii)

        choose_transients = np.where(ii)
        this_idx = ut.make_idx(self.f_ppm, fit_lb, fit_ub)
        return choose_transients, model, signal, params, this_idx