From 4bc0c5a5b3e2e0b080b6b36bd805df68251c1a8e Mon Sep 17 00:00:00 2001
From: Rolf Verberg <rv43@cornell.edu>
Date: Wed, 2 Sep 2026 11:51:52 -0400
Subject: [PATCH 1/6] style: some trivial style changes

---
 CHAP/edd/processor.py | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/CHAP/edd/processor.py b/CHAP/edd/processor.py
index 70efc0c6..add8e696 100755
--- a/CHAP/edd/processor.py
+++ b/CHAP/edd/processor.py
@@ -399,7 +399,8 @@ def _get_sum_axes_data(self, nxdata, detector_id, sum_axes=True):
                 unique_points.append(point)
                 sum_indices.append([i])
         if not unique_points:
-            return NXdata(NXfield(np.empty((0, data.shape[-1])), 'detector_data'))
+            return NXdata(
+                NXfield(np.empty((0, data.shape[-1])), 'detector_data'))
         unique_points = np.asarray(unique_points).T
         mean_data = np.empty((unique_points.shape[1], data.shape[-1]))
         for i in range(unique_points.shape[1]):
@@ -1326,7 +1327,8 @@ def process(self, data):
         # Get and subtract the detector baselines
         self._subtract_baselines()
 
-        # Calibrate detector channel energies based on fluorescence peaks
+        # Calibrate detector channel energies based on fluorescence
+        # peaks
         self._calibrate()
 
         # Combine the calibration and detector configuration
@@ -2889,7 +2891,8 @@ def process(self, data):
                 else:
                     self.logger.warning('Skip adding points')
             if self.standalone and self.update:
-                # Return nxprocess structure and values separately for writer
+                # Return nxprocess structure and values separately for
+                # writer
                 values = self._get_values(results)
                 ret = [nxsetup, values]
             else:
@@ -3609,9 +3612,6 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             material_names_used = material_names[use_peaks]
             peak_locations_used = peak_locations[use_peaks]
 
-            # Get the selected HKLs, lattice spacings, material names
-            # and peak locations used in the fit
-
             # Get the masked energies and spectra for this detector
             energies_masked = energies[mask]
             intensities_masked = nxdata.nxsignal.nxdata[:,mask]
@@ -3621,10 +3621,6 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             # Get the masked energies and spectra after accounting
             # for unused peaks
             if not all(use_peaks):
-                # FIX the fitting works, but adding pointwise breaks
-                # peaks are skipped. Must remap the intensities_masked
-                # and the best fit and residuals to the same shape as
-                # the initial intensities_masked shape
                 mask_ranges_used, _, _ = \
                     select_mask_and_hkls(
                         energies_masked, mean_data_masked, hkls,
@@ -3702,7 +3698,8 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                 normal_strains.append(strain)
                 det_angles.append(detector.attrs['eta'])
 
-            # Insert the peaks omitted from the fit due to find_peak_cutoff
+            # Insert the peaks omitted from the fit due to
+            # find_peak_cutoff
             amplitudes_vary = np.asarray(fit_results['amplitudes_vary'])
             if num_points > 1:
                 amplitudes_vary = np.moveaxis(amplitudes_vary, -1, 0)
@@ -3712,7 +3709,8 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             amplitudes_vary = np.insert(
                 amplitudes_vary, insert_peak_indices, [False], axis=-1)
 
-            # Remap the results if peaks where masked
+            # Remap the results if peaks where omitted from the fit due
+            # to find_peak_cutoff
             best_fits_masked = np.asarray(fit_results['best_fits'])
             residuals_masked = np.asarray(fit_results['residuals'])
             if all(use_peaks):

From 24808b81e5ab4482e8ceaf459a46313f52bf0c3f Mon Sep 17 00:00:00 2001
From: Rolf Verberg <rv43@cornell.edu>
Date: Wed, 14 Oct 2026 10:54:44 -0400
Subject: [PATCH 2/6] fix: fixes a merge error in CHAP/edd/processor.py

---
 CHAP/edd/processor.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/CHAP/edd/processor.py b/CHAP/edd/processor.py
index 8f2f3064..2f84db1b 100755
--- a/CHAP/edd/processor.py
+++ b/CHAP/edd/processor.py
@@ -3713,9 +3713,6 @@ def _get_use_peaks(energies, mean_data, peak_locations):
 
             # Insert the peaks omitted from the fit due to
             # find_peak_cutoff
-            amplitudes_vary = np.asarray(fit_results['amplitudes_vary'])
-            if num_points > 1:
-                amplitudes_vary = np.moveaxis(amplitudes_vary, -1, 0)
             insert_peak_indices = [
                 vv-ii for ii, vv in enumerate(
                     i for i, v in enumerate(use_peaks) if not v)]

From 304b15d198cb2ce35ab3ea076205f05c15ece889 Mon Sep 17 00:00:00 2001
From: Rolf Verberg <rv43@cornell.edu>
Date: Wed, 14 Oct 2026 16:08:52 -0400
Subject: [PATCH 3/6] add: adds multiple MultpiPeak models for multi component
 fitting

---
 CHAP/edd/processor.py | 119 ++++++++++++++++++++++++++----------------
 CHAP/edd/utils.py     |  53 ++++++++++++++-----
 CHAP/utils/fit.py     |  62 +++++++++++-----------
 CHAP/utils/models.py  |   3 ++
 4 files changed, 149 insertions(+), 88 deletions(-)

diff --git a/CHAP/edd/processor.py b/CHAP/edd/processor.py
index 2f84db1b..274332a1 100755
--- a/CHAP/edd/processor.py
+++ b/CHAP/edd/processor.py
@@ -3319,25 +3319,29 @@ def _get_nxprocess(self, nxentry, calibration_config):
             # Setup the NXdata object under NXdetector for the
             # detector wide peak fitting results
             shape = det_nxdata.intensity.shape
-            det_nxdata.best_fit = NXfield(
-                shape=shape, dtype=np.float64,
-                maxshape=(None, *shape[1:]), chunks=(1, *shape[1:]))
-            det_nxdata.included_peaks = NXfield(
-                shape=[shape[0], len(hkls_fit)], dtype=bool,
-                maxshape=(None, *shape[1:]), chunks=(1, *shape[1:]))
-            det_nxdata.included_peaks.attrs['hkls'] = dumps(
-                peak_fit_info.get('hkls', ''))
-            det_nxdata.included_peaks.attrs['use_peaks'] = peak_fit_info.get(
-                'use_peaks')
-            det_nxdata.residual = NXfield(
-                shape=shape, dtype=np.float64,
-                maxshape=(None, *shape[1:]), chunks=(1, *shape[1:]))
-            det_nxdata.redchi = NXfield(
-                shape=[shape[0]], dtype=np.float64,
-                maxshape=(None,), chunks=(1,))
-            det_nxdata.success = NXfield(
-                shape=[shape[0]], dtype='bool',
-                maxshape=(None,), chunks=(1,))
+            for fit_type in {'uniform', 'unconstrained'}:
+                nxdetector[f'data_{fit_type}'] = NXdata()
+                ddet_nxdata = nxdetector[f'data_{fit_type}']
+                ddet_nxdata.best_fit = NXfield(
+                    shape=shape, dtype=np.float64,
+                    maxshape=(None, *shape[1:]), chunks=(1, *shape[1:]))
+                ddet_nxdata.included_peaks = NXfield(
+                    shape=[shape[0], len(hkls_fit)], dtype=bool,
+                    maxshape=(None, *shape[1:]), chunks=(1, *shape[1:]))
+                ddet_nxdata.included_peaks.attrs['hkls'] = dumps(
+                    peak_fit_info.get('hkls', ''))
+                ddet_nxdata.included_peaks.attrs['use_peaks'] = \
+                    peak_fit_info.get('use_peaks')
+                ddet_nxdata.residual = NXfield(
+                    shape=shape, dtype=np.float64,
+                    maxshape=(None, *shape[1:]), chunks=(1, *shape[1:]))
+                ddet_nxdata.redchi = NXfield(
+                    shape=[shape[0]], dtype=np.float64,
+                    maxshape=(None,), chunks=(1,))
+                ddet_nxdata.success = NXfield(
+                    shape=[shape[0]], dtype='bool',
+                    maxshape=(None,), chunks=(1,))
+
 
             # Setup the NXcollection objects for the per material peak
             # fitting results
@@ -3676,10 +3680,10 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                 material_names_used)
 
             # Perform a component-wise uniform fit
-            fit_results_uniform = {}
+            component_hkl_indices = {}
             if num_peak > 1:
                 for material_name in set(material_names):
-                    component_mask_ranges, component_hkl_indices, _ = \
+                    component_mask_ranges, hkl_indices, _ = \
                         select_mask_and_hkls(
                             energies_masked, mean_data_masked, hkls_used,
                             ds_used, detector.tth_calibrated,
@@ -3689,27 +3693,24 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                             detector_id=detector.get_id(),
                             label='Sum of spectra in the map ({material_name})',
                             interactive=False)
-                    component_mask = np.zeros(
-                        (energies_masked.size)).astype(bool)
-                    indices = np.arange(energies_masked.size, dtype=np.int32)
-                    for min_, max_ in component_mask_ranges:
-                        component_mask = np.logical_or(
-                            component_mask,
-                            np.logical_and(indices >= min_, indices < max_))
-                    self.logger.info(f'Fitting detector {detector.get_id()} '
-                                     f'for the {material_name} component...')
-                    fit_results_uniform[material_name] = get_spectra_fits(
-                        np.squeeze(intensities_masked), energies_masked,
-                        peak_locations_used[component_hkl_indices], detector,
-                        fit_type='uniform', mask=~component_mask,
-                        max_nfev=self.config.max_nfev,
-                        num_proc=self.config.num_proc, **self.run_config)
-                    self.logger.info('... done')
+                    component_hkl_indices[material_name] = hkl_indices
+                self.logger.info(
+                    f'Fitting detector {detector.get_id()} (uniform)...')
+                fit_results_uniform = get_spectra_fits(
+                    np.squeeze(intensities_masked), energies_masked,
+                    peak_locations_used, detector, fit_type='uniform',
+                    component_hkl_indices=component_hkl_indices,
+                    max_nfev=self.config.max_nfev,
+                    num_proc=self.config.num_proc, **self.run_config)
+                self.logger.info('... done')
+                amplitudes_vary_uniform = np.asarray([
+                    fit_results['best_vary'][f'peak{i+1}_amplitude']
+                    for i in range(num_peak)])
 
             # FIX for multiple materials
-            if len(component_strains) == 1 and 'eta' in detector.attrs:
-                normal_strains.append(strains)
-                det_angles.append(detector.attrs['eta'])
+#            if len(component_strains) == 1 and 'eta' in detector.attrs:
+#                normal_strains.append(strains)
+#                det_angles.append(detector.attrs['eta'])
 
             # Insert the peaks omitted from the fit due to
             # find_peak_cutoff
@@ -3719,6 +3720,10 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             amplitudes_vary = np.insert(
                 np.moveaxis(amplitudes_vary, -1, 0), insert_peak_indices,
                 [False], axis=-1)
+            if num_peak > 1:
+                amplitudes_vary_uniform = np.insert(
+                    np.moveaxis(amplitudes_vary_uniform, -1, 0),
+                    insert_peak_indices, [False], axis=-1)
 
             # Remap the results if peaks where masked
             best_fits_masked = fit_results['best_fits']
@@ -3731,6 +3736,17 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                 residuals = np.zeros(map_shape)
                 best_fits[:,mask_used] = best_fits_masked
                 residuals[:,mask_used] = residuals_masked
+            if num_peak > 1:
+                best_fits_uniform_masked = fit_results_uniform['best_fits']
+                residuals_uniform_masked = fit_results_uniform['residuals']
+                if all(use_peaks):
+                    best_fits_uniform = best_fits_uniform_masked
+                    residuals_uniform = residuals_uniform_masked
+                else:
+                    best_fits_uniform = np.zeros(map_shape)
+                    residuals_uniform = np.zeros(map_shape)
+                    best_fits_uniform[:,mask_used] = best_fits_uniform_masked
+                    residuals_uniform[:,mask_used] = residuals_uniform_masked
 
             # Store results as full arrays (leading axis = num_points)
             det_id = detector.get_id()
@@ -3738,6 +3754,9 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             results.update({
                 f'{path}/intensity': nxdata.nxsignal.nxdata[:,mask],
                 f'{path}/norm': intensity_norms,
+            })
+            path = f'{det_id}/data_unconstrained'
+            results.update({
                 f'{path}/best_fit': best_fits,
                 f'{path}/included_peaks': amplitudes_vary,
                 f'{path}/residual': residuals,
@@ -3745,18 +3764,26 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                 f'{path}/success': fit_results['success'],
             })
             if num_peak > 1:
-                for material_name, result_uniform \
-                        in fit_results_uniform.items():
-                    path = f'{det_id}/{material_name}/uniform'
-                    results.update({
-                        f'{path}/strain': result_uniform['strain'],
-                    })
+                path = f'{det_id}/data_uniform'
+                results.update({
+                    f'{path}/best_fit': best_fits_uniform,
+                    f'{path}/included_peaks': amplitudes_vary_uniform,
+                    f'{path}/residual': residuals_uniform,
+                    f'{path}/redchi': fit_results_uniform['redchis'],
+                    f'{path}/success': fit_results_uniform['success'],
+                })
             for material_name, component_strain in component_strains.items():
                 path = f'{det_id}/{material_name}/unconstrained'
                 results.update({
                     f'{path}/strain': component_strain['mean'],
                     f'{path}/strain_stdev': component_strain['stdev'],
                 })
+                if num_peak > 1:
+                    path = f'{det_id}/{material_name}/uniform'
+                    results.update({
+                        f'{path}/strain':
+                            fit_results_uniform['strain'][material_name],
+                    })
             for j, (hkl, material_name) in enumerate(zip(
                     hkls[use_peaks], material_names[use_peaks])):
                 hkl_name = '_'.join(str(hkl)[1:-1].split(' '))
diff --git a/CHAP/edd/utils.py b/CHAP/edd/utils.py
index 6996b331..14d0b438 100755
--- a/CHAP/edd/utils.py
+++ b/CHAP/edd/utils.py
@@ -1269,7 +1269,7 @@ def get_rolling_sum_spectra(
 
 def get_spectra_fits(
         spectra, energies, peak_locations, detector, mask=None,
-        fit_type='unconstrained', **kwargs):
+        fit_type='unconstrained', component_hkl_indices=None, **kwargs):
     """Return the fit results for the map of spectra provided.
 
     :param spectra: Intensity spectra to fit.
@@ -1280,11 +1280,14 @@ def get_spectra_fits(
     :type peak_locations: list[float]
     :param detector: MCA detector element configuration.
     :type detector: MCAElementStrainAnalysisConfig
-    :ivar mask: Mask to apply to the data during fitting (boolean
+    :param mask: Mask to apply to the data during fitting (boolean
         array, `True` meaning mask in the fit).
-    :param mask: numpy.ndarray, optional
+    :type mask: numpy.ndarray, optional
     :param fit_type: Type of fit, defaults to `'unconstrained'`.
     :type fit_type: Literal['uniform', 'unconstrained'], optional
+    :param component_hkl_indices: The component names and its peak
+        indices in peak_locations for individual components.
+    :type component_hkl_indices: dict[str, list[int]], optional
     :returns: The best fit parameters and their errors and vary
         properties, uniform strains, best fits, residuals between the
         best fits and the input spectra, reduced chi, and fit success
@@ -1298,6 +1301,9 @@ def get_spectra_fits(
     from CHAP.pipeline import PipelineData
     from CHAP.utils.fit import FitProcessor
 
+    if component_hkl_indices is None:
+        component_hkl_indices = {}
+
     # Construct the fit model
     models = []
     if detector.background is not None:
@@ -1313,16 +1319,28 @@ def get_spectra_fits(
         for peak in backgroundpeaks:
             peak.prefix = f'bkgd_{peak.prefix}'
         models += backgroundpeaks
-    models.append(
-        {'model_type': 'multipeak', 'centers': list(peak_locations),
-         'fit_type': fit_type, 'peak_models': detector.peak_models,
-         'centers_range': detector.centers_range,
-         'centers_range_fraction': detector.centers_range_fraction,
-         'fwhm_min': detector.fwhm_min, 'fwhm_max': detector.fwhm_max})
+    num_component = len(component_hkl_indices)
+    if num_component > 1:
+        for material_name, hkl_indices in component_hkl_indices.items():
+            models.append(
+                {'model_type': 'multipeak',
+                 'centers': list(peak_locations[hkl_indices]),
+                 'fit_type': fit_type, 'peak_models': detector.peak_models,
+                 'centers_range': detector.centers_range,
+                 'centers_range_fraction': detector.centers_range_fraction,
+                 'fwhm_min': detector.fwhm_min, 'fwhm_max': detector.fwhm_max,
+                 'prefix': f'{material_name}_'})
+    else:
+        models.append(
+            {'model_type': 'multipeak',
+             'centers': list(peak_locations),
+             'fit_type': fit_type, 'peak_models': detector.peak_models,
+             'centers_range': detector.centers_range,
+             'centers_range_fraction': detector.centers_range_fraction,
+             'fwhm_min': detector.fwhm_min, 'fwhm_max': detector.fwhm_max})
     config = {
         'abs_height_cutoff': detector.abs_height_cutoff,
 #        'code': 'lmfit',
-        'force_fitmap': True,
         'models': models,
 #        'plot': True,
 #        'print_report': True,
@@ -1342,6 +1360,18 @@ def get_spectra_fits(
     if mask is not None:
         data.append(PipelineData(name='mask', data=mask))
     fit = FitProcessor.run(data=data, config=config, **kwargs)
+    fit_strain = None
+    if fit_type == 'uniform' and len(peak_locations) > 1:
+        if num_component > 1:
+            fit_strain = {}
+            for name in component_hkl_indices:
+                fit_strain[name] = 1 - fit.best_values[
+                    f'{name}_scale_factor']
+        elif num_component:
+            fit_strain = {str(next(iter(component_hkl_indices))):
+                1 - fit.best_values['scale_factor']}
+        else:
+            fit_strain = 1 - fit.best_values['scale_factor']
     return {
         'best_errors': fit.best_errors,
         'best_fits': fit.best_fit,
@@ -1350,6 +1380,5 @@ def get_spectra_fits(
         'residuals': fit.residual,
         'redchis': fit.redchi,
         'success': fit.success,
-        'strain': 1 - fit.best_values['scale_factor']
-            if fit_type == 'uniform' and len(peak_locations) > 1 else None,
+        'strain': fit_strain,
     }
diff --git a/CHAP/utils/fit.py b/CHAP/utils/fit.py
index 4d16d293..59215a63 100755
--- a/CHAP/utils/fit.py
+++ b/CHAP/utils/fit.py
@@ -176,21 +176,20 @@ def process(self, data):
             return data
 
         # Expand multipeak model if present
-        found_multipeak = False
-        multipeak_info = None
-        for i, model in enumerate(deepcopy(self.config.models)):
+        multipeak_models = []
+        for i, model in reversed(list(enumerate(
+                deepcopy(self.config.models)))):
             if isinstance(model, MultipeakModel):
-                if found_multipeak:
-                    raise ValueError(
-                        f'Invalid parameter models ({self.config.models}) '
-                        '(multiple instances of multipeak not allowed)')
-                parameters, models = self.create_multipeak_model(model)
-                if parameters:
-                    self.config.parameters += parameters
-                self.config.models += models
-                self.config.models.pop(i)
-                found_multipeak = True
-                multipeak_info = model.model_dump()
+                multipeak_models.append(self.config.models.pop(i))
+        num_mpm = len(multipeak_models)
+        for i, model in enumerate(multipeak_models):
+            if num_mpm > 1 and not model.prefix:
+                model.prefix = f'mpm{i+1}_'
+            parameters, models = self.create_multipeak_model(model)
+            if parameters:
+                self.config.parameters += parameters
+            self.config.models += models
+        multipeak_info = [m.model_dump() for m in multipeak_models]
 
         # Instantiate the Fit object and fit the data
         y = np.squeeze(data[1])[None,:] \
@@ -228,20 +227,22 @@ def _uniform_model(
 
             if not model.centers_range_fraction:
                 parameters = [FitParameter(
-                    name='scale_factor', value=1.0, min=FLOAT_MIN)]
+                    name=f'{model.prefix}scale_factor', value=1.0,
+                    min=FLOAT_MIN)]
             else:
                 parameters = [FitParameter(
-                    name='scale_factor', value=1.0,
+                    name=f'{model.prefix}scale_factor', value=1.0,
                     min=1.0-model.centers_range_fraction,
                     max=1.0+model.centers_range_fraction)]
             peak_models = []
             for i, cen in enumerate(model.centers):
                 peak_models.append(peak_model_class(
                     model_type=model.peak_models,
-                    prefix=f'peak{i+1}_',
+                    prefix=f'{model.prefix}peak{i+1}_',
                     parameters=[
                          {'name': 'amplitude', 'min': FLOAT_MIN},
-                         {'name': 'center', 'expr': f'scale_factor*{cen}'},
+                         {'name': 'center',
+                          'expr': f'{model.prefix}scale_factor*{cen}'},
                          {'name': 'sigma', 'min': sig_min, 'max': sig_max}]))
             return parameters, peak_models
 
@@ -253,7 +254,7 @@ def _unconstrained_model(
                 if not (model.centers_range and model.centers_range_fraction):
                     peak_models.append(peak_model_class(
                         model_type=model.peak_models,
-                        prefix=f'peak{i+1}_',
+                        prefix=f'{model.prefix}peak{i+1}_',
                         parameters=[
                              {'name': 'amplitude', 'min': FLOAT_MIN},
                              {'name': 'center', 'value': cen},
@@ -264,7 +265,7 @@ def _unconstrained_model(
                         model.centers_range, cen*model.centers_range_fraction)
                     peak_models.append(peak_model_class(
                         model_type=model.peak_models,
-                        prefix=f'peak{i+1}_',
+                        prefix=f'{model.prefix}peak{i+1}_',
                         parameters=[
                              {'name': 'amplitude', 'min': FLOAT_MIN},
                              {'name': 'center', 'value': cen,
@@ -1387,7 +1388,7 @@ def fit(self, config=None, **kwargs):
         if config is None:
             num_proc = kwargs.pop('num_proc', num_proc_max)
             self._abs_height_cutoff = kwargs.pop('abs_height_cutoff')
-            self._multipeak_info = kwargs.pop('multipeak_info', None)
+            self._multipeak_info = kwargs.pop('multipeak_info', [])
             self._plot = kwargs.pop('plot', False)
             self._plot_init = kwargs.pop('plot_init', False)
             self._print_report = kwargs.pop('print_report', False)
@@ -1744,6 +1745,9 @@ def _create_prefixes(self, models):
                         models[i].prefix = f'{name}{n}_'
 
     def _fit(self, n, current_best_values, return_result=False, **kwargs):
+        # Third party modules
+        from asteval import Interpreter
+
         # Do not attempt a fit if the normalized data is zero or
         # entirely below the cutoff
         y_min = self._y_norm[n].min()
@@ -1831,15 +1835,12 @@ def _fit(self, n, current_best_values, return_result=False, **kwargs):
 
         parameters_save = deepcopy(self._parameters)
         parameters_bounds_save = deepcopy(self._parameter_bounds)
-        if self._multipeak_info is not None:
-            # Third party modules
-            from asteval import Interpreter
-
-            centers = self._multipeak_info.get('centers')
-            centers_range = self._multipeak_info.get('centers_range')
+        for multipeak_info in self._multipeak_info:
+            centers = multipeak_info.get('centers')
+            centers_range = multipeak_info.get('centers_range')
             centers_range_fraction = \
-                self._multipeak_info.get('centers_range_fraction')
-            model_type = self._multipeak_info.get('peak_models')
+                multipeak_info.get('centers_range_fraction')
+            model_type = multipeak_info.get('peak_models')
             min_height = None if self._rel_height_cutoff is None \
                 else y_max*self._rel_height_cutoff
             use_peaks, _, peak_heights, peak_widths = \
@@ -1847,10 +1848,11 @@ def _fit(self, n, current_best_values, return_result=False, **kwargs):
                     self._x, self._y_norm[n], centers, centers_range,
                     centers_range_fraction, min_height=min_height, min_width=5)
 
+            prefix = multipeak_info.get('prefix')
             ast = Interpreter()
             for i, (use_peak, height, width) in enumerate(zip(
                     use_peaks, peak_heights, peak_widths)):
-                name = f'peak{i+1}_amplitude'
+                name = f'{prefix}peak{i+1}_amplitude'
                 ast(f'fwhm = {width}')
                 ast(f'height = {height}')
                 sigma = ast(fwhm_factor[model_type])
diff --git a/CHAP/utils/models.py b/CHAP/utils/models.py
index 8ff0a857..52a85859 100755
--- a/CHAP/utils/models.py
+++ b/CHAP/utils/models.py
@@ -1078,6 +1078,8 @@ class MultipeakModel(CHAPBaseModel):
     :ivar peak_models: Type of peaks, defaults to `'gaussian'`.
     :vartype peak_models: Literal['gaussian', 'lorentzian', 'pvoigt'],
         optional.
+    :ivar prefix: Model prefix, defaults to `''`.
+    :vartype prefix: str, optional
     """
 
     model_type: Literal['multipeak']
@@ -1089,6 +1091,7 @@ class MultipeakModel(CHAPBaseModel):
     fwhm_min: Optional[confloat(allow_inf_nan=False)] = None
     fwhm_max: Optional[confloat(allow_inf_nan=False)] = None
     peak_models: Literal['gaussian', 'lorentzian', 'pvoigt'] = 'gaussian'
+    prefix: Optional[str] = ''
 
 
 class FitConfig(CHAPBaseModel):

From cac20411365c7bbbc1774c4988b0570e6d3f2a79 Mon Sep 17 00:00:00 2001
From: Rolf Verberg <rv43@cornell.edu>
Date: Fri, 16 Oct 2026 10:20:58 -0400
Subject: [PATCH 4/6] feat: adds data_unconstrained, data_uniform and uniform
 fit parameters

---
 CHAP/edd/processor.py | 471 +++++++++++++++++++++++-------------------
 CHAP/edd/utils.py     |  15 +-
 2 files changed, 271 insertions(+), 215 deletions(-)

diff --git a/CHAP/edd/processor.py b/CHAP/edd/processor.py
index 274332a1..ff24a056 100755
--- a/CHAP/edd/processor.py
+++ b/CHAP/edd/processor.py
@@ -2952,78 +2952,82 @@ def _add_fit_nxcollection(
 
         # Peak-by-peak results
         # FIX linkdims doesn't seem to be working
-        for hkl, material_name in zip(hkls, material_names):
-            hkl_name = '_'.join(str(hkl)[1:-1].split(' '))
-            nxcollection = nxdetector[material_name]['unconstrained']
-            nxcollection[hkl_name] = NXparameters()
-            nxparameter = nxcollection[hkl_name]
-            # Create initial centers field
-            nxparameter.center_initial_guess = NXdata()
-            self._linkdims(
-                nxparameter.center_initial_guess, det_nxdata,
-                skip_field_dims=['energy'])
-            # Report HKL peak amplitudes
-            nxparameter.amplitudes = NXdata()
-            self._linkdims(
-                nxparameter.amplitudes, det_nxdata, skip_field_dims=['energy'])
-            nxparameter.amplitudes.values = NXfield(
-                shape=[shape[0]], dtype=np.float64, attrs={'units': 'counts'},
-                maxshape=(None,), chunks=(1,))
-            nxparameter.amplitudes.errors = NXfield(
-                shape=[shape[0]], dtype=np.float64, maxshape=(None,),
-                chunks=(1,))
-            nxparameter.amplitudes.attrs['signal'] = 'values'
-            # Report HKL peak centers
-            nxparameter.centers = NXdata()
-            self._linkdims(
-                nxparameter.centers, det_nxdata, skip_field_dims=['energy'])
-            nxparameter.centers.values = NXfield(
-                shape=[shape[0]], dtype=np.float64, attrs={'units': 'keV'},
-                maxshape=(None,), chunks=(1,))
-            nxparameter.centers.errors = NXfield(
-                shape=[shape[0]], dtype=np.float64, maxshape=(None,),
-                chunks=(1,))
-            nxparameter.centers.attrs['signal'] = 'values'
-            # Report HKL peak FWHMs
-            nxparameter.sigmas = NXdata()
-            self._linkdims(
-                nxparameter.sigmas, det_nxdata, skip_field_dims=['energy'])
-            nxparameter.sigmas.values = NXfield(
-                shape=[shape[0]], dtype=np.float64, attrs={'units': 'keV'},
-                maxshape=(None,), chunks=(1,))
-            nxparameter.sigmas.errors = NXfield(
-                shape=[shape[0]], dtype=np.float64, maxshape=(None,),
-                chunks=(1,))
-            nxparameter.sigmas.attrs['signal'] = 'values'
-            if peak_fit_info.get('peak_models') == 'pvoigt':
-                # Report HKL peak fractions
-                nxparameter.fractions = NXdata()
+        for fit_type in {'uniform', 'unconstrained'}:
+            for hkl, material_name in zip(hkls, material_names):
+                hkl_name = '_'.join(str(hkl)[1:-1].split(' '))
+                nxcollection = nxdetector[material_name][fit_type]
+                nxcollection[hkl_name] = NXparameters()
+                nxparameter = nxcollection[hkl_name]
+                # Create initial centers field
+                nxparameter.center_initial_guess = NXdata()
+                self._linkdims(
+                    nxparameter.center_initial_guess, det_nxdata,
+                    skip_field_dims=['energy'])
+                # Report HKL peak amplitudes
+                nxparameter.amplitudes = NXdata()
                 self._linkdims(
-                    nxparameter.fractions, det_nxdata,
+                    nxparameter.amplitudes, det_nxdata,
                     skip_field_dims=['energy'])
-                nxparameter.fractions.values = NXfield(
+                nxparameter.amplitudes.values = NXfield(
+                    shape=[shape[0]], dtype=np.float64,
+                    attrs={'units': 'counts'}, maxshape=(None,), chunks=(1,))
+                nxparameter.amplitudes.errors = NXfield(
                     shape=[shape[0]], dtype=np.float64, maxshape=(None,),
                     chunks=(1,))
-                nxparameter.fractions.errors = NXfield(
+                nxparameter.amplitudes.attrs['signal'] = 'values'
+                # Report HKL peak centers
+                nxparameter.centers = NXdata()
+                self._linkdims(
+                    nxparameter.centers, det_nxdata,
+                    skip_field_dims=['energy'])
+                nxparameter.centers.values = NXfield(
+                    shape=[shape[0]], dtype=np.float64, attrs={'units': 'keV'},
+                    maxshape=(None,), chunks=(1,))
+                nxparameter.centers.errors = NXfield(
                     shape=[shape[0]], dtype=np.float64, maxshape=(None,),
                     chunks=(1,))
-                nxparameter.fractions.attrs['signal'] = 'values'
-            # Report HKL peak strains
-            nxparameter.strains = NXdata()
-            self._linkdims(
-                nxparameter.strains, det_nxdata,
-                skip_field_dims=['energy'])
-            values = np.full(shape=[shape[0]], fill_value=np.nan)
-            nxparameter.strains.values = NXfield(
-                value=values, shape=[shape[0]], dtype=np.float64,
-                maxshape=(None,), chunks=(1,))
-            nxparameter.strains.errors = NXfield(
-                value=values, shape=[shape[0]], dtype=np.float64,
-                maxshape=(None,), chunks=(1,))
-            nxparameter.strains.residuals = NXfield(
-                value=values, shape=[shape[0]], dtype=np.float64,
-                maxshape=(None,), chunks=(1,))
-            nxparameter.strains.attrs['signal'] = 'values'
+                nxparameter.centers.attrs['signal'] = 'values'
+                # Report HKL peak FWHMs
+                nxparameter.sigmas = NXdata()
+                self._linkdims(
+                    nxparameter.sigmas, det_nxdata, skip_field_dims=['energy'])
+                nxparameter.sigmas.values = NXfield(
+                    shape=[shape[0]], dtype=np.float64, attrs={'units': 'keV'},
+                    maxshape=(None,), chunks=(1,))
+                nxparameter.sigmas.errors = NXfield(
+                    shape=[shape[0]], dtype=np.float64, maxshape=(None,),
+                    chunks=(1,))
+                nxparameter.sigmas.attrs['signal'] = 'values'
+                if peak_fit_info.get('peak_models') == 'pvoigt':
+                    # Report HKL peak fractions
+                    nxparameter.fractions = NXdata()
+                    self._linkdims(
+                        nxparameter.fractions, det_nxdata,
+                        skip_field_dims=['energy'])
+                    nxparameter.fractions.values = NXfield(
+                        shape=[shape[0]], dtype=np.float64, maxshape=(None,),
+                        chunks=(1,))
+                    nxparameter.fractions.errors = NXfield(
+                        shape=[shape[0]], dtype=np.float64, maxshape=(None,),
+                        chunks=(1,))
+                    nxparameter.fractions.attrs['signal'] = 'values'
+                # Report HKL peak strains
+                if fit_type == 'unconstrained':
+                    nxparameter.strains = NXdata()
+                    self._linkdims(
+                        nxparameter.strains, det_nxdata,
+                        skip_field_dims=['energy'])
+                    values = np.full(shape=[shape[0]], fill_value=np.nan)
+                    nxparameter.strains.values = NXfield(
+                        value=values, shape=[shape[0]], dtype=np.float64,
+                        maxshape=(None,), chunks=(1,))
+                    nxparameter.strains.errors = NXfield(
+                        value=values, shape=[shape[0]], dtype=np.float64,
+                        maxshape=(None,), chunks=(1,))
+                    nxparameter.strains.residuals = NXfield(
+                        value=values, shape=[shape[0]], dtype=np.float64,
+                        maxshape=(None,), chunks=(1,))
+                    nxparameter.strains.attrs['signal'] = 'values'
 
     def _populate_peak_fit_info(self):
         """Populate _peak_fit_info with all configured HKLs for each
@@ -3047,7 +3051,8 @@ def _populate_peak_fit_info(self):
             ds_fit = np.asarray([ds[i] for i in detector.hkl_indices])
             material_names_fit = [
                 material_names[i] for i in detector.hkl_indices]
-            peak_locations = get_peak_locations(ds_fit, detector.tth_calibrated)
+            peak_locations = get_peak_locations(
+                ds_fit, detector.tth_calibrated)
             self._peak_fit_info.append({
                 'hkls': [''.join(map(str, hkl)) for hkl in hkls_fit],
                 'materials': material_names_fit,
@@ -3173,7 +3178,7 @@ def _get_values(self, results):
         """
         return [{'data': v, 'path': k} for k, v in results.items()]
 
-    def _get_nxprocess(self, nxentry, calibration_config):
+    def _get_nxprocess(self, nxentry, calibration_config, name):
         """Return a standalone NXprocess for the strain analysis
         results & metadata.
 
@@ -3183,6 +3188,8 @@ def _get_nxprocess(self, nxentry, calibration_config):
         :param calibration_config: 2&theta calibration configuration.
         :type calibration_config:
             CHAP.edd.models.MCATthCalibrationConfig
+        :param name: NXprocess name.
+        :type name: str
         :return: Strain analysis results & associated metadata.
         :rtype: nexusformat.nexus.NXprocess
         """
@@ -3211,7 +3218,7 @@ def _get_nxprocess(self, nxentry, calibration_config):
                 'StrainAnalysis Configuration, or re-run the pipeline with '
                 'the --interactive flag.')
 
-        nxprocess = NXprocess()
+        nxprocess = NXprocess(name=name)
         nxprocess.calibration_config = \
             calibration_config.model_dump_json()
         nxprocess.strain_analysis_config = \
@@ -3349,24 +3356,42 @@ def _get_nxprocess(self, nxentry, calibration_config):
                 nxdetector[name] = NXcollection()
                 nxcollection = nxdetector[name]
                 nxcollection['uniform'] = NXcollection()
-                nxcollection['uniform'].strain = NXfield(
-                    dtype=np.float64,
-                    shape=(num_points,),
-                    attrs={'long_name': 'Strain from peak fit (\u03B5)'},
-                    maxshape=(None,), chunks=(1,))
+                nxcollection['uniform'].strains = NXdata(
+                    NXfield(
+                        dtype=np.float64,
+                        name='values',
+                        shape=(num_points,),
+                        attrs={'long_name':
+                                   'Strain values from peak fit (\u03B5)'},
+                        maxshape=(None,), chunks=(1,)),
+                    NXfield(
+                        dtype=np.float64,
+                        name='errors',
+                        shape=(num_points,),
+                        attrs={'long_name':
+                                   'Errors in strain values from peak fit '
+                                   '(\u03B5)'},
+                        maxshape=(None,), chunks=(1,)))
+                self._linkdims(
+                    nxcollection['uniform'].strains, det_nxdata,
+                    skip_field_dims=['energy'])
                 nxcollection['unconstrained'] = NXcollection()
-                nxcollection['unconstrained'].strain = NXfield(
-                    dtype=np.float64,
-                    shape=(num_points,),
-                    attrs={'long_name': 'Strain from peak fit (\u03B5)'},
-                    maxshape=(None,), chunks=(1,))
-                nxcollection['unconstrained'].strain_stdev = NXfield(
-                    dtype=np.float64,
-                    shape=(num_points,),
-                    attrs={'long_name':
-                               'Standard deviation in strain from peak fit '
-                               '(\u03B5)'},
-                    maxshape=(None,), chunks=(1,))
+                nxcollection['unconstrained'].strains = NXdata(
+                    NXfield(
+                        dtype=np.float64,
+                        name='means',
+                        shape=(num_points,),
+                        attrs={'long_name':
+                                   'Mean strains from peak fit (\u03B5)'},
+                        maxshape=(None,), chunks=(1,)),
+                    NXfield(
+                        dtype=np.float64,
+                        name='stdevs',
+                        shape=(num_points,),
+                        attrs={'long_name':
+                                   'Standard deviations in strain from peak '
+                                   'fit (\u03B5)'},
+                        maxshape=(None,), chunks=(1,)))
 
             # Add the peak fit results
             self._add_fit_nxcollection(
@@ -3397,7 +3422,8 @@ def _get_nxroot(self, nxentry, calibration_config):
 
         nxroot = NXroot()
         nxroot[nxentry.nxname] = nxentry
-        nxprocess = self._get_nxprocess(nxentry, calibration_config)
+        nxprocess = self._get_nxprocess(
+            nxentry, calibration_config, f'{nxentry.nxname}_strainanalysis')
         nxroot[f'{nxentry.nxname}_strainanalysis'] = nxprocess
         return nxroot
 
@@ -3564,6 +3590,113 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                 return np.asarray([]).astype(bool)
             return use_peaks
 
+        def _store_global_fit_results(
+                results, fit_results, insert_peak_indices, mask_used,
+                map_shape, num_peak_used, det_id, component_peak_indices=None,
+                fit_type='unconstrained'):
+            # Insert peaks omitted from the fit due to find_peak_cutoff
+            if component_peak_indices is None:
+                included_peaks = np.asarray([
+                    fit_results['best_vary'][f'peak{i+1}_amplitude']
+                    for i in range(num_peak_used)])
+            else:
+                included_peaks = np.asarray([
+                    fit_results['best_vary'][f'{name}_peak{i+1}_amplitude']
+                    for name, i in component_peak_indices])
+            included_peaks = np.insert(
+                np.moveaxis(included_peaks, -1, 0), insert_peak_indices,
+                [False], axis=-1)
+            # Remap the results if peaks where masked
+            best_fits_masked = fit_results['best_fits']
+            residuals_masked = fit_results['residuals']
+            if mask_used is None:
+                best_fits = best_fits_masked
+                residuals = residuals_masked
+            else:
+                best_fits = np.zeros(map_shape)
+                residuals = np.zeros(map_shape)
+                best_fits[:,mask_used] = best_fits_masked
+                residuals[:,mask_used] = residuals_masked
+            # Store results as full arrays (leading axis = num_points)
+            path = f'{det_id}/data_{fit_type}'
+            results.update({
+                f'{path}/best_fit': best_fits,
+                f'{path}/included_peaks': included_peaks,
+                f'{path}/residual': residuals,
+                f'{path}/redchi': fit_results['redchis'],
+                f'{path}/success': fit_results['success'],
+            })
+            return best_fits_masked
+
+        def _store_peak_fit_results(
+            results, fit_results, hkls, material_names, use_peaks, det_id,
+            component_strains=None, component_peak_indices=None,
+            fit_type='unconstrained'):
+            if fit_type == 'uniform':
+                for name, strain in fit_results['strain'].items():
+                    path = f'{det_id}/{name}/uniform/strains'
+                    results.update({
+                        f'{path}/values': strain['value'],
+                        f'{path}/errors': strain['error'],
+                    })
+            else:
+                for name, strain in component_strains.items():
+                    path = f'{det_id}/{name}/unconstrained/strains'
+                    results.update({
+                        f'{path}/means': strain['mean'],
+                        f'{path}/stdevs': strain['stdev'],
+                    })
+            for j, (hkl, material_name) in enumerate(zip(
+                    hkls[use_peaks], material_names[use_peaks])):
+                hkl_name = '_'.join(str(hkl)[1:-1].split(' '))
+                path = f'{det_id}/{material_name}/{fit_type}/{hkl_name}'
+                if component_peak_indices is None:
+                    prefix = f'peak{j+1}_'
+                else:
+                    name, i = component_peak_indices[j]
+                    assert name == material_name
+                    prefix = f'{name}_peak{i+1}_'
+                results.update({
+                    f'{path}/amplitudes/values':
+                        fit_results['best_values'][f'{prefix}amplitude'],
+                    f'{path}/amplitudes/errors':
+                        fit_results['best_errors'][f'{prefix}amplitude'],
+                    f'{path}/centers/values': centers[j],
+                    f'{path}/centers/errors':
+                        fit_results['best_errors'][f'{prefix}center'],
+                    f'{path}/sigmas/values':
+                        fit_results['best_values'][f'{prefix}sigma'],
+                    f'{path}/sigmas/errors':
+                        fit_results['best_errors'][f'{prefix}sigma'],
+                })
+                if detector.peak_models == 'pvoigt':
+                    results.update({
+                        f'{path}/fractions/values':
+                            fit_results['best_values'][f'{prefix}fraction'],
+                        f'{path}/fractions/errors':
+                            fit_results['best_errors'][f'{prefix}fraction'],
+                    })
+                # Strain values: NaN where unconstrained center is zero
+#                has_center = centers[j].astype(bool)
+#                centers_errors_j = \
+#                    fit_results['best_errors'][f'{prefix}center']
+#                strain_errors_raw = np.asarray([
+#                    (e / c if e is not None else np.nan)
+#                    for e, c in zip(centers_errors_j, centers[j])
+#                ])
+#                results.update({
+#                    f'{path}/strains/values':
+#                        np.where(has_center, strains[j], np.nan),
+#                    f'{path}/strains/errors':
+#                        np.where(has_center, strain_errors_raw, np.nan),
+#                    f'{path}/strains/residuals':
+#                        np.where(
+#                            has_center,
+#                            (component_strains[material_name]['mean'] -
+#                             strains[j]),
+#                            np.nan),
+#                })
+
         # Get and subtract the detector baselines
         self._subtract_baselines()
 
@@ -3628,6 +3761,7 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             ds_used = ds[use_peaks]
             material_names_used = material_names[use_peaks]
             peak_locations_used = peak_locations[use_peaks]
+            num_peak_used = peak_locations_used.size
 
             # Get the masked energies and spectra for this detector
             energies_masked = energies[mask]
@@ -3637,7 +3771,8 @@ def _get_use_peaks(energies, mean_data, peak_locations):
 
             # Get the masked energies and spectra after accounting
             # for unused peaks
-            if not all(use_peaks):
+            mask_used = None
+            if num_peak_used != use_peaks.size:
                 mask_ranges_used, _, _ = \
                     select_mask_and_hkls(
                         energies_masked, mean_data_masked, hkls,
@@ -3659,7 +3794,7 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             intensity_norms = intensities_masked.max(
                 axis=intensities_masked.ndim-1)
 
-            # Perform the fit
+            # Perform the unconstrained fit
             self.logger.info(f'Fitting detector {detector.get_id()} ...')
             fit_results = get_spectra_fits(
                 np.squeeze(intensities_masked), energies_masked,
@@ -3667,21 +3802,10 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                 num_proc=self.config.num_proc, **self.run_config)
             self.logger.info('... done')
 
-            # Get the strains
-            num_peak = peak_locations_used.size
-            centers = np.asarray([
-                fit_results['best_values'][f'peak{i+1}_center']
-                for i in range(num_peak)])
-            amplitudes_vary = np.asarray([
-                fit_results['best_vary'][f'peak{i+1}_amplitude']
-                for i in range(num_peak)])
-            strains, component_strains = _get_component_strains(
-                centers, peak_locations_used[:,None], amplitudes_vary,
-                material_names_used)
-
             # Perform a component-wise uniform fit
             component_hkl_indices = {}
-            if num_peak > 1:
+            component_peak_indices = None
+            if num_peak_used > 1:
                 for material_name in set(material_names):
                     component_mask_ranges, hkl_indices, _ = \
                         select_mask_and_hkls(
@@ -3703,9 +3827,23 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                     max_nfev=self.config.max_nfev,
                     num_proc=self.config.num_proc, **self.run_config)
                 self.logger.info('... done')
-                amplitudes_vary_uniform = np.asarray([
-                    fit_results['best_vary'][f'peak{i+1}_amplitude']
-                    for i in range(num_peak)])
+                if len(component_hkl_indices) > 1:
+                    component_peak_indices = num_peak_used*[None]
+                    for name, hkl_indices in component_hkl_indices.items():
+                        for i, index in enumerate(hkl_indices):
+                            component_peak_indices[index] = (name, i)
+                    assert None not in component_peak_indices
+
+            # Get the strains
+            centers = np.asarray([
+                fit_results['best_values'][f'peak{i+1}_center']
+                for i in range(num_peak_used)])
+            amplitudes_vary = np.asarray([
+                fit_results['best_vary'][f'peak{i+1}_amplitude']
+                for i in range(num_peak_used)])
+            strains, component_strains = _get_component_strains(
+                centers, peak_locations_used[:,None], amplitudes_vary,
+                material_names_used)
 
             # FIX for multiple materials
 #            if len(component_strains) == 1 and 'eta' in detector.attrs:
@@ -3717,36 +3855,6 @@ def _get_use_peaks(energies, mean_data, peak_locations):
             insert_peak_indices = [
                 vv-ii for ii, vv in enumerate(
                     i for i, v in enumerate(use_peaks) if not v)]
-            amplitudes_vary = np.insert(
-                np.moveaxis(amplitudes_vary, -1, 0), insert_peak_indices,
-                [False], axis=-1)
-            if num_peak > 1:
-                amplitudes_vary_uniform = np.insert(
-                    np.moveaxis(amplitudes_vary_uniform, -1, 0),
-                    insert_peak_indices, [False], axis=-1)
-
-            # Remap the results if peaks where masked
-            best_fits_masked = fit_results['best_fits']
-            residuals_masked = fit_results['residuals']
-            if all(use_peaks):
-                best_fits = best_fits_masked
-                residuals = residuals_masked
-            else:
-                best_fits = np.zeros(map_shape)
-                residuals = np.zeros(map_shape)
-                best_fits[:,mask_used] = best_fits_masked
-                residuals[:,mask_used] = residuals_masked
-            if num_peak > 1:
-                best_fits_uniform_masked = fit_results_uniform['best_fits']
-                residuals_uniform_masked = fit_results_uniform['residuals']
-                if all(use_peaks):
-                    best_fits_uniform = best_fits_uniform_masked
-                    residuals_uniform = residuals_uniform_masked
-                else:
-                    best_fits_uniform = np.zeros(map_shape)
-                    residuals_uniform = np.zeros(map_shape)
-                    best_fits_uniform[:,mask_used] = best_fits_uniform_masked
-                    residuals_uniform[:,mask_used] = residuals_uniform_masked
 
             # Store results as full arrays (leading axis = num_points)
             det_id = detector.get_id()
@@ -3755,80 +3863,24 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                 f'{path}/intensity': nxdata.nxsignal.nxdata[:,mask],
                 f'{path}/norm': intensity_norms,
             })
-            path = f'{det_id}/data_unconstrained'
-            results.update({
-                f'{path}/best_fit': best_fits,
-                f'{path}/included_peaks': amplitudes_vary,
-                f'{path}/residual': residuals,
-                f'{path}/redchi': fit_results['redchis'],
-                f'{path}/success': fit_results['success'],
-            })
-            if num_peak > 1:
-                path = f'{det_id}/data_uniform'
-                results.update({
-                    f'{path}/best_fit': best_fits_uniform,
-                    f'{path}/included_peaks': amplitudes_vary_uniform,
-                    f'{path}/residual': residuals_uniform,
-                    f'{path}/redchi': fit_results_uniform['redchis'],
-                    f'{path}/success': fit_results_uniform['success'],
-                })
-            for material_name, component_strain in component_strains.items():
-                path = f'{det_id}/{material_name}/unconstrained'
-                results.update({
-                    f'{path}/strain': component_strain['mean'],
-                    f'{path}/strain_stdev': component_strain['stdev'],
-                })
-                if num_peak > 1:
-                    path = f'{det_id}/{material_name}/uniform'
-                    results.update({
-                        f'{path}/strain':
-                            fit_results_uniform['strain'][material_name],
-                    })
-            for j, (hkl, material_name) in enumerate(zip(
-                    hkls[use_peaks], material_names[use_peaks])):
-                hkl_name = '_'.join(str(hkl)[1:-1].split(' '))
-                path = f'{det_id}/{material_name}/unconstrained/{hkl_name}'
-                prefix = f'peak{j+1}_'
-                results.update({
-                    f'{path}/amplitudes/values':
-                        fit_results['best_values'][f'{prefix}amplitude'],
-                    f'{path}/amplitudes/errors':
-                        fit_results['best_errors'][f'{prefix}amplitude'],
-                    f'{path}/centers/values': centers[j],
-                    f'{path}/centers/errors':
-                        fit_results['best_errors'][f'{prefix}center'],
-                    f'{path}/sigmas/values':
-                        fit_results['best_values'][f'{prefix}sigma'],
-                    f'{path}/sigmas/errors':
-                        fit_results['best_errors'][f'{prefix}sigma'],
-                })
-                if detector.peak_models == 'pvoigt':
-                    results.update({
-                        f'{path}/fractions/values':
-                            fit_results['best_values'][f'{prefix}fraction'],
-                        f'{path}/fractions/errors':
-                            fit_results['best_errors'][f'{prefix}fraction'],
-                    })
-                # Strain values: NaN where unconstrained center is zero
-                has_center = centers[j].astype(bool)
-                centers_errors_j = \
-                    fit_results['best_errors'][f'{prefix}center']
-                strain_errors_raw = np.asarray([
-                    (e / c if e is not None else np.nan)
-                    for e, c in zip(centers_errors_j, centers[j])
-                ])
-                results.update({
-                    f'{path}/strains/values':
-                        np.where(has_center, strains[j], np.nan),
-                    f'{path}/strains/errors':
-                        np.where(has_center, strain_errors_raw, np.nan),
-                    f'{path}/strains/residuals':
-                        np.where(
-                            has_center,
-                            (component_strains[material_name]['mean'] -
-                             strains[j]),
-                            np.nan),
-                })
+            best_fits_masked = _store_global_fit_results(
+                results, fit_results, insert_peak_indices, mask_used,
+                map_shape, num_peak_used, det_id)
+            if num_peak_used > 1:
+                _store_global_fit_results(
+                    results, fit_results_uniform, insert_peak_indices,
+                    mask_used, map_shape, num_peak_used, det_id,
+                    component_peak_indices=component_peak_indices,
+                    fit_type='uniform')
+            _store_peak_fit_results(
+                results, fit_results, hkls, material_names, use_peaks,
+                det_id, component_strains=component_strains)
+            if num_peak_used > 1:
+                _store_peak_fit_results(
+                    results, fit_results_uniform, hkls, material_names,
+                    use_peaks, det_id,
+                    component_peak_indices=component_peak_indices,
+                    fit_type='uniform')
 
             if self.json_results:
                 raise RuntimeError('self.json_results needs testing')
@@ -3849,6 +3901,7 @@ def _get_use_peaks(energies, mean_data, peak_locations):
                                 f'{det_id}/{material_name}/unconstrained/'
                                 f'{hkl_name}/{fitparam}/{fitquant}'
                             )
+
                             results.update({key: placeholder})
 
             # Create an animation of the fit points
diff --git a/CHAP/edd/utils.py b/CHAP/edd/utils.py
index 14d0b438..9dd9b4c7 100755
--- a/CHAP/edd/utils.py
+++ b/CHAP/edd/utils.py
@@ -1365,13 +1365,16 @@ def get_spectra_fits(
         if num_component > 1:
             fit_strain = {}
             for name in component_hkl_indices:
-                fit_strain[name] = 1 - fit.best_values[
-                    f'{name}_scale_factor']
-        elif num_component:
-            fit_strain = {str(next(iter(component_hkl_indices))):
-                1 - fit.best_values['scale_factor']}
+                fit_strain[name] = {
+                    'value': 1 - fit.best_values[f'{name}_scale_factor'],
+                    'error': 1 - fit.best_errors[f'{name}_scale_factor'],
+                }
         else:
-            fit_strain = 1 - fit.best_values['scale_factor']
+            fit_strain = {
+                str(next(iter(component_hkl_indices))):
+                    {'value': 1 - fit.best_values['scale_factor'],
+                     'error': 1 - fit.best_errors['scale_factor']}
+            }
     return {
         'best_errors': fit.best_errors,
         'best_fits': fit.best_fit,

From 9d4ec66a0d9abcf922208bd3a66c84afd5ebb5b2 Mon Sep 17 00:00:00 2001
From: Rolf Verberg <rv43@cornell.edu>
Date: Fri, 16 Oct 2026 15:00:30 -0400
Subject: [PATCH 5/6] fix: fixes bug in the uniform strain errors

---
 CHAP/edd/utils.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/CHAP/edd/utils.py b/CHAP/edd/utils.py
index 9dd9b4c7..b342bb4a 100755
--- a/CHAP/edd/utils.py
+++ b/CHAP/edd/utils.py
@@ -1367,13 +1367,13 @@ def get_spectra_fits(
             for name in component_hkl_indices:
                 fit_strain[name] = {
                     'value': 1 - fit.best_values[f'{name}_scale_factor'],
-                    'error': 1 - fit.best_errors[f'{name}_scale_factor'],
+                    'error': fit.best_errors[f'{name}_scale_factor'],
                 }
         else:
             fit_strain = {
                 str(next(iter(component_hkl_indices))):
                     {'value': 1 - fit.best_values['scale_factor'],
-                     'error': 1 - fit.best_errors['scale_factor']}
+                     'error': fit.best_errors['scale_factor']}
             }
     return {
         'best_errors': fit.best_errors,

From f92e5b3e89c41619edf2c3e6300fe2b96d3a2c1f Mon Sep 17 00:00:00 2001
From: Rolf Verberg <rv43@cornell.edu>
Date: Mon, 19 Oct 2026 12:59:05 -0400
Subject: [PATCH 6/6] fix: removes obsolete Rosette Strain computation from
 Strain analysis

This will be covered later in a dedicated processor.
---
 CHAP/edd/models.py    |  4 ---
 CHAP/edd/processor.py | 83 -------------------------------------------
 2 files changed, 87 deletions(-)

diff --git a/CHAP/edd/models.py b/CHAP/edd/models.py
index d2356458..26eca20c 100755
--- a/CHAP/edd/models.py
+++ b/CHAP/edd/models.py
@@ -859,9 +859,6 @@ class StrainAnalysisConfig(MCACalibrationConfig):
         configuration
         :class:`~CHAP.edd.models.MCADetectorStrainAnalysis`).
     :vartype abs_height_cutoff: int, optional
-    :ivar calc_rosette_strains: Compute and add the rosette strain
-         fields defaults to `False`.
-    :vartype calc_rosette_strains: bool, optional
     :ivar find_peak_cutoff: Use scipy.signal.find_peaks to exclude
         peaks for all spectra for a given detector and user specified
         mask. A particular HKL peak is removed from the set of HKLs,
@@ -892,7 +889,6 @@ class StrainAnalysisConfig(MCACalibrationConfig):
     #:vartype oversampling: FIX
 
     abs_height_cutoff: Optional[conint(gt=0)] = None
-    calc_rosette_strains: Optional[bool] = False
     find_peak_cutoff: Optional[confloat(ge=0, allow_inf_nan=False)] = 0.0
     max_nfev: Optional[conint(gt=0)] = None
     num_proc: Optional[conint(gt=0)] = max(1, os.cpu_count()//4)
diff --git a/CHAP/edd/processor.py b/CHAP/edd/processor.py
index ff24a056..39c44086 100755
--- a/CHAP/edd/processor.py
+++ b/CHAP/edd/processor.py
@@ -3224,32 +3224,6 @@ def _get_nxprocess(self, nxentry, calibration_config, name):
         nxprocess.strain_analysis_config = \
             self.config.model_dump_json()
 
-        # Add the Rosette strains fields
-        if self.config.calc_rosette_strains:
-            num_points = self._nxdata_detectors[0].nxsignal.shape[0]
-            nxprocess.strain_rosette = NXparameters()
-            nxprocess.strain_rosette.e_xx = NXdata(
-                NXfield(
-                    dtype=np.float64, name='values', shape=(num_points,),
-                    maxshape=(None,), chunks=(1,)),
-                NXfield(
-                    dtype=np.float64, name='errors', shape=(num_points,),
-                    maxshape=(None,), chunks=(1,)))
-            nxprocess.strain_rosette.e_yy = NXdata(
-                NXfield(
-                    dtype=np.float64, shape=(num_points,), maxshape=(None,),
-                    chunks=(1,), name='values'),
-                NXfield(
-                    dtype=np.float64, shape=(num_points,), maxshape=(None,),
-                    chunks=(1,), name='errors'))
-            nxprocess.strain_rosette.e_xy = NXdata(
-                NXfield(
-                    dtype=np.float64, shape=(num_points,), maxshape=(None,),
-                    chunks=(1,), name='values'),
-                NXfield(
-                    dtype=np.float64, shape=(num_points,), maxshape=(None,),
-                    chunks=(1,), name='errors'))
-
         if not self._peak_fit_info:
             # FIX this is a temporary fix to be able to run update
             # after setup.
@@ -3845,11 +3819,6 @@ def _store_peak_fit_results(
                 centers, peak_locations_used[:,None], amplitudes_vary,
                 material_names_used)
 
-            # FIX for multiple materials
-#            if len(component_strains) == 1 and 'eta' in detector.attrs:
-#                normal_strains.append(strains)
-#                det_angles.append(detector.attrs['eta'])
-
             # Insert the peaks omitted from the fit due to
             # find_peak_cutoff
             insert_peak_indices = [
@@ -3911,60 +3880,8 @@ def _store_peak_fit_results(
                     nxdata, energies_masked, intensities_masked,
                     intensity_norms, best_fits_masked, detector.get_id())
 
-        # Calculate and add the Rosette strains
-        if self.config.calc_rosette_strains:
-            raise RuntimeError('Rosette strains needs testing')
-            popt, perr = self._fit_strain_rosette(normal_strains, det_angles)
-            if self.json_results:
-                results.update({
-                    'strain_rosette/e_xx/values': np.asarray([popt[0]]),
-                    'strain_rosette/e_xx/errors': np.asarray([perr[0]]),
-                    'strain_rosette/e_yy/values': np.asarray([popt[1]]),
-                    'strain_rosette/e_yy/errors': np.asarray([perr[1]]),
-                    'strain_rosette/e_xy/values': np.asarray([popt[2]]),
-                    'strain_rosette/e_xy/errors': np.asarray([perr[2]]),
-                })
-
         return results
 
-    def _fit_strain_rosette(self, normal_strains, det_angles):
-        # Third party modules
-        from scipy.optimize import curve_fit
-
-        # FIX make an individual processor? Right now not updated
-        def strain_rosette_calc(angle, e_xx, e_yy, e_xy):
-            # Calculate the xx strain at an angle, given the xx, yy,
-            # and xy strains.
-            # Accept scalar or array `angle` (radians).
-            raise ValueError('_fit_strain_rosette needs updating/checking')
-            a = np.asarray(angle)
-            c = np.cos(a)
-            s = np.sin(a)
-            e1 = e_xx * c**2 + e_yy * s**2 + 2.0 * e_xy * c * s
-            return e1
-
-        normal_strains = np.squeeze(normal_strains)
-        assert normal_strains.ndim == 1
-        det_angles = np.asarray(det_angles)
-        assert normal_strains.size == det_angles.size
-        det_angles_rad = np.radians(det_angles)
-        # Use normal strain near 0 degrees as initial guess for e_xx
-        e_xx_guess = normal_strains[np.argmin(np.abs(det_angles - 0))]
-        # Use normal strain near 90 degrees as initial guess for e_yy
-        e_yy_guess = normal_strains[np.argmin(np.abs(det_angles - 90))]
-        e_xy_guess = 0.0
-        bounds_guess = 100*np.max(np.abs(normal_strains))
-        bounds = (
-            [-bounds_guess, -bounds_guess, -bounds_guess],
-            [bounds_guess, bounds_guess, bounds_guess])
-
-        popt, pcov = curve_fit(
-            strain_rosette_calc, det_angles_rad, normal_strains,
-            p0=(e_xx_guess, e_yy_guess, e_xy_guess),
-            bounds=bounds, **{'max_nfev': 10000})
-        perr = np.sqrt(np.diag(pcov))
-        return popt, perr
-
 
 if __name__ == '__main__':
     # Local modules