diff --git a/.github/workflows/python-ci.yml b/.github/workflows/python-ci.yml
index 90967cd9..1d33313f 100644
--- a/.github/workflows/python-ci.yml
+++ b/.github/workflows/python-ci.yml
@@ -61,7 +61,7 @@ jobs:
         name: unit-tests-job
         flags: unittests
         files: ./coverage-unit.xml
-        fail_ci_if_error: true
+        fail_ci_if_error: false
         verbose: true
         token: ${{ secrets.CODECOV_TOKEN }}
         slug: EasyScience/EasyReflectometryLib
diff --git a/CHANGELOG.md b/CHANGELOG.md
index bb6a1617..b786ccd6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,12 @@
+# Version 1.4.0 (unreleased)
+
+Add Mighell-based handling of zero-variance points in fitting (issue #256).
+Zero-variance data points are no longer forcibly discarded; instead, a hybrid
+objective applies a Mighell substitution for zero-variance points while using
+standard weighted least squares for the rest. The previous masking behavior is
+available via `objective='legacy_mask'`. New `objective` parameter on
+`MultiFitter`, `fit()`, and `fit_single_data_set_1d()`.
+
 # Version 1.3.3 (17 June 2025)
 
 Added Chi^2 and fit status to fitting results.
diff --git a/MIGHELL_IMPLEMENTATION.md b/MIGHELL_IMPLEMENTATION.md
new file mode 100644
index 00000000..26a61a9d
--- /dev/null
+++ b/MIGHELL_IMPLEMENTATION.md
@@ -0,0 +1,113 @@
+# Mighell Implementation Notes
+
+## Summary
+
+The `mighell` objective implemented in `easyreflectometry.fitting` is mathematically consistent with the
+`chi^2_gamma` statistic described in `mighell_fitting.pdf`.
+
+However, the paper derives this statistic for Poisson-distributed count data, whereas this library applies it to
+reflectometry values that are typically already normalized and are not raw counts. That distinction is important for
+interpreting the fit behavior.
+
+## Formula From The PDF
+
+Section 3 of `mighell_fitting.pdf` defines the statistic
+
+$$
+\chi^2_\gamma = \sum_i \frac{[n_i + \min(n_i, 1) - m_i]^2}{n_i + 1}
+$$
+
+where:
+
+- `n_i` are Poisson-distributed observed counts
+- `m_i` are model values
+
+The paper proposes this in preference to the modified Neyman statistic for Poisson data.
+
+## Implemented Form
+
+The code in `src/easyreflectometry/fitting.py` implements the same objective in weighted least-squares form.
+
+For `objective='mighell'`, it constructs:
+
+$$
+y_{\mathrm{eff},i} = y_i + \min(y_i, 1)
+$$
+
+and
+
+$$
+\sigma_i = \sqrt{y_i + 1}
+$$
+
+with weights
+
+$$
+w_i = \frac{1}{\sigma_i}
+$$
+
+so the minimized least-squares objective is
+
+$$
+\sum_i \left(\frac{y_{\mathrm{eff},i} - f_i}{\sigma_i}\right)^2
+=
+\sum_i \frac{[y_i + \min(y_i,1) - f_i]^2}{y_i + 1}
+$$
+
+which is exactly the same algebraic form as the paper's `chi^2_gamma` statistic.
+
+
+This has been implemented:
+- shifting the target by `min(y, 1)`
+- using `sqrt(y + 1)` as the effective uncertainty
+- minimizing the resulting weighted residual sum of squares
+
+## Important Scope Limitation
+
+The PDF is explicitly about Poisson-distributed count data.
+
+In `reflectometry-lib`, `mighell` is being applied to reflectometry intensities / reflectivities, which are generally:
+
+- normalized values rather than raw counts
+- potentially processed quantities rather than direct Poisson observations
+- not guaranteed to satisfy the assumptions used in the paper's derivation
+
+
+## The Full Mighell Fit Can Look Worse
+
+The poor visual agreement of the full `mighell` fit against the measured reflectivity is expected under this
+implementation (see results of running `zero_variances_example.py or
+notebooks\zero_variance_fitting.ipynb)
+
+That is because `mighell` here is not only a reweighting of the residuals. It also changes the fitted target:
+
+$$
+y \rightarrow y + \min(y,1)
+$$
+
+For any point with `0 < y < 1`, this doubles the target value.
+
+At low `Q`, where reflectivity values can still be relatively large, the transformed target can be substantially higher
+than the measured reflectivity. The optimizer is therefore solving a different problem than “fit the plotted
+reflectivity values as closely as possible”.
+
+The effect:
+
+- `legacy_mask` and `hybrid` often nearly overlap
+- full `mighell` can visibly deviate from the experimental curve
+- the `mighell` objective value can be small while the classical chi-square is poor
+
+## Hybrid mode
+
+The `hybrid` mode is not taken directly from the PDF.
+
+It is a project-specific adaptation that:
+
+- uses standard weighted least squares where variances are valid
+- applies the Mighell substitution only where variance is zero or invalid
+
+This makes `hybrid` a better (?) reflectometry-oriented compromise rather than a direct reproduction of the paper.
+
+
+I would propose using `hybrid` as default, but give users an option of both masking the zero variance points and applying proper Mighell algorithm.
+
diff --git a/notebooks/zero_variance_fitting.ipynb b/notebooks/zero_variance_fitting.ipynb
new file mode 100644
index 00000000..ea3ede96
--- /dev/null
+++ b/notebooks/zero_variance_fitting.ipynb
@@ -0,0 +1,724 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "1a8a52ac",
+   "metadata": {},
+   "source": [
+    "# Zero-Variance Handling in Reflectometry Fitting\n",
+    "\n",
+    "This notebook demonstrates how `easyreflectometry` handles data points with **zero variance**\n",
+    "during fitting. We compare three objective strategies:\n",
+    "\n",
+    "| Objective | Behaviour |\n",
+    "|---|---|\n",
+    "| `legacy_mask` | Drops zero-variance points entirely (old default) |\n",
+    "| `hybrid` | Keeps all points; applies Mighell substitution only to zero-variance entries (**new default**) |\n",
+    "| `mighell` | Applies the Mighell (1999) transform to every point |\n",
+    "\n",
+    "The experimental data comes from `tests/_static/ref_zero_var.txt`, which contains 192 data points,\n",
+    "6 of which have zero error (= zero variance)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5c6429ac",
+   "metadata": {},
+   "source": [
+    "## 1. Import Required Libraries"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4983acb8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "import warnings\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "from easyreflectometry.calculators import CalculatorFactory\n",
+    "from easyreflectometry.data.measurement import load\n",
+    "from easyreflectometry.fitting import MultiFitter\n",
+    "from easyreflectometry.model import Model\n",
+    "from easyreflectometry.model import PercentageFwhm\n",
+    "from easyreflectometry.sample import Layer\n",
+    "from easyreflectometry.sample import Material\n",
+    "from easyreflectometry.sample import Multilayer\n",
+    "from easyreflectometry.sample import Sample\n",
+    "\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "96e174d4",
+   "metadata": {},
+   "source": [
+    "## 2. Load Experimental Data\n",
+    "\n",
+    "The file `ref_zero_var.txt` is a comma-separated file with columns: `q (Å⁻¹)`, `R`, `error`.\n",
+    "Several points near the high-Q end have `error = 0.0`, meaning zero variance."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "47fb3264",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "DATA_PATH = os.path.join('..', 'tests', '_static', 'ref_zero_var.txt')\n",
+    "\n",
+    "# Load raw data for inspection\n",
+    "raw = np.loadtxt(DATA_PATH, delimiter=',', comments='#')\n",
+    "q_raw, r_raw, err_raw = raw[:, 0], raw[:, 1], raw[:, 2]\n",
+    "\n",
+    "# Load through easyreflectometry (produces a scipp DataGroup)\n",
+    "data = load(DATA_PATH)\n",
+    "data_key = next(iter(data['data'].keys()))\n",
+    "coord_key = next(iter(data['coords'].keys()))\n",
+    "result_model_key = f'{data_key}_model'\n",
+    "\n",
+    "print(f'Total data points : {len(q_raw)}')\n",
+    "print(f'Q range           : [{q_raw.min():.4e}, {q_raw.max():.4e}] Å⁻¹')\n",
+    "print(f'R range           : [{r_raw.min():.4e}, {r_raw.max():.4e}]')\n",
+    "print(f'Data key          : {data_key}')\n",
+    "print(f'Coord key         : {coord_key}')\n",
+    "print(f'Model key         : {result_model_key}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c4495939",
+   "metadata": {},
+   "source": [
+    "## 3. Inspect Data and Identify Zero-Variance Points"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1be2ca78",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "zero_mask = err_raw == 0.0\n",
+    "zero_indices = np.where(zero_mask)[0]\n",
+    "print(f'Number of zero-variance points: {zero_mask.sum()}')\n",
+    "print(f'Indices: {zero_indices}')\n",
+    "print('Q-values with zero variance:')\n",
+    "for idx in zero_indices:\n",
+    "    print(f'  [{idx:3d}]  Q = {q_raw[idx]:.4e} Å⁻¹,  R = {r_raw[idx]:.4e}')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "373cdc62",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "\n",
+    "# Plot all data\n",
+    "valid = ~zero_mask\n",
+    "ax.errorbar(q_raw[valid], r_raw[valid], yerr=err_raw[valid],\n",
+    "            fmt='o', ms=3, color='C0', alpha=0.7, label='Valid points')\n",
+    "ax.plot(q_raw[zero_mask], r_raw[zero_mask],\n",
+    "        'rx', ms=10, mew=2, label=f'Zero-variance points ({zero_mask.sum()})')\n",
+    "\n",
+    "ax.set_yscale('log')\n",
+    "ax.set_xlabel('Q (Å⁻¹)')\n",
+    "ax.set_ylabel('Reflectivity')\n",
+    "ax.set_title('Raw experimental data — zero-variance points highlighted')\n",
+    "ax.legend()\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "428759d3",
+   "metadata": {},
+   "source": [
+    "## 4. Define Materials\n",
+    "\n",
+    "A simple film-on-substrate structure: **Si** substrate / **SiO₂** native oxide / **Film** / **D₂O** superphase."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "006fb3cf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "si = Material(2.07, 0, 'Si')\n",
+    "sio2 = Material(3.47, 0, 'SiO2')\n",
+    "film = Material(2.0, 0, 'Film')\n",
+    "d2o = Material(6.36, 0, 'D2O')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6031b9ad",
+   "metadata": {},
+   "source": [
+    "## 5. Define Layers and Sample Structure"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2643141f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "si_layer = Layer(si, 0, 0, 'Si layer')\n",
+    "sio2_layer = Layer(sio2, 30, 3, 'SiO2 layer')\n",
+    "film_layer = Layer(film, 250, 3, 'Film layer')\n",
+    "superphase = Layer(d2o, 0, 3, 'D2O superphase')\n",
+    "\n",
+    "sample = Sample(\n",
+    "    Multilayer(si_layer),\n",
+    "    Multilayer(sio2_layer),\n",
+    "    Multilayer(film_layer),\n",
+    "    Multilayer(superphase),\n",
+    "    name='Film Structure',\n",
+    ")\n",
+    "print(sample)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "51bcedeb",
+   "metadata": {},
+   "source": [
+    "## 6. Create the Model"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "22e6fbd9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "resolution_function = PercentageFwhm(0.02)\n",
+    "model = Model(sample, 1, 1e-6, resolution_function, 'Film Model')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7eaf6993",
+   "metadata": {},
+   "source": [
+    "## 7. Set Calculator Backend and Compute Initial Curve"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "175afdcc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "interface = CalculatorFactory()\n",
+    "model.interface = interface\n",
+    "\n",
+    "# Compute initial model reflectivity\n",
+    "r_init = interface.fit_func(q_raw, model.unique_name)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fff782bb",
+   "metadata": {},
+   "source": [
+    "## 8. Plot Initial Model vs Experimental Data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a7656b8b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "\n",
+    "ax.errorbar(q_raw[valid], r_raw[valid], yerr=err_raw[valid],\n",
+    "            fmt='o', ms=3, color='C0', alpha=0.6, label='Experiment (valid)')\n",
+    "ax.plot(q_raw[zero_mask], r_raw[zero_mask],\n",
+    "        'rx', ms=10, mew=2, label='Experiment (zero variance)')\n",
+    "ax.plot(q_raw, r_init, '-', color='C1', lw=1.5, label='Initial model')\n",
+    "\n",
+    "ax.set_yscale('log')\n",
+    "ax.set_xlabel('Q (Å⁻¹)')\n",
+    "ax.set_ylabel('Reflectivity')\n",
+    "ax.set_title('Experimental data vs initial model (before fitting)')\n",
+    "ax.legend()\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6c916d41",
+   "metadata": {},
+   "source": [
+    "## 9. Set Fitting Parameters and Constraints"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a85ab7e8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sio2_layer.thickness.fixed = False\n",
+    "sio2_layer.thickness.bounds = (15, 50)\n",
+    "\n",
+    "film_layer.thickness.fixed = False\n",
+    "film_layer.thickness.bounds = (200, 300)\n",
+    "\n",
+    "film.sld.fixed = False\n",
+    "film.sld.bounds = (0.1, 3)\n",
+    "\n",
+    "model.background.fixed = False\n",
+    "model.background.bounds = (1e-7, 1e-5)\n",
+    "\n",
+    "model.scale.fixed = False\n",
+    "model.scale.bounds = (0.5, 1.5)\n",
+    "\n",
+    "print('Free parameters:')\n",
+    "for p in model.get_fit_parameters():\n",
+    "    print(f'  {p.name:20s} = {float(p.value):.4g}  bounds={p.bounds}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b83af903",
+   "metadata": {},
+   "source": [
+    "## 10. Zero-Variance Handling — Three Objective Modes\n",
+    "\n",
+    "We now fit the same data with each objective mode and compare the results.\n",
+    "\n",
+    "### How each mode handles zero-variance points\n",
+    "\n",
+    "- **`legacy_mask`** — zero-variance points are dropped before fitting. The fitter sees fewer\n",
+    "  data points, and those Q-values have no influence on the result.\n",
+    "- **`hybrid`** (default) — valid points use standard weighted least-squares ($w = 1/\\sigma$).\n",
+    "  Zero-variance points get the **Mighell (1999)** substitution:\n",
+    "  $y_\\text{eff} = y + \\min(y, 1)$, $\\sigma = \\sqrt{\\max(y + 1,\\, \\varepsilon)}$.\n",
+    "  This keeps all data in the fit while giving zero-variance points reasonable weight.\n",
+    "- **`mighell`** — applies the Mighell transform to *every* point, not just zero-variance ones.\n",
+    "  This changes both the weighting and the fitted target for the entire dataset.\n",
+    "\n",
+    "Below we report two metric families:\n",
+    "\n",
+    "- **objective chi²** — the quantity actually minimized by the selected objective\n",
+    "- **classical chi²** — a standard variance-weighted chi² computed only on the original\n",
+    "  positive-variance points, for comparison"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8182f0a9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def build_fresh_model():\n",
+    "    \"\"\"Build a fresh model+interface so each fit starts from the same initial state.\"\"\"\n",
+    "    _si = Material(2.07, 0, 'Si')\n",
+    "    _sio2 = Material(3.47, 0, 'SiO2')\n",
+    "    _film = Material(2.0, 0, 'Film')\n",
+    "    _d2o = Material(6.36, 0, 'D2O')\n",
+    "\n",
+    "    _si_layer = Layer(_si, 0, 0, 'Si layer')\n",
+    "    _sio2_layer = Layer(_sio2, 30, 3, 'SiO2 layer')\n",
+    "    _film_layer = Layer(_film, 250, 3, 'Film layer')\n",
+    "    _superphase = Layer(_d2o, 0, 3, 'D2O superphase')\n",
+    "\n",
+    "    _sample = Sample(\n",
+    "        Multilayer(_si_layer),\n",
+    "        Multilayer(_sio2_layer),\n",
+    "        Multilayer(_film_layer),\n",
+    "        Multilayer(_superphase),\n",
+    "        name='Film Structure',\n",
+    "    )\n",
+    "    _resolution = PercentageFwhm(0.02)\n",
+    "    _model = Model(_sample, 1, 1e-6, _resolution, 'Film Model')\n",
+    "\n",
+    "    _sio2_layer.thickness.fixed = False\n",
+    "    _sio2_layer.thickness.bounds = (15, 50)\n",
+    "    _film_layer.thickness.fixed = False\n",
+    "    _film_layer.thickness.bounds = (200, 300)\n",
+    "    _film.sld.fixed = False\n",
+    "    _film.sld.bounds = (0.1, 3)\n",
+    "    _model.background.fixed = False\n",
+    "    _model.background.bounds = (1e-7, 1e-5)\n",
+    "    _model.scale.fixed = False\n",
+    "    _model.scale.bounds = (0.5, 1.5)\n",
+    "\n",
+    "    _model.interface = CalculatorFactory()\n",
+    "    return _model"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fd6f4833",
+   "metadata": {},
+   "source": [
+    "### 11. Fit with `legacy_mask`\n",
+    "\n",
+    "Zero-variance points are simply dropped."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e856e00f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_mask = build_fresh_model()\n",
+    "fitter_mask = MultiFitter(model_mask, objective='legacy_mask')\n",
+    "\n",
+    "with warnings.catch_warnings(record=True) as w_mask:\n",
+    "    warnings.simplefilter('always')\n",
+    "    result_mask = fitter_mask.fit(data)\n",
+    "\n",
+    "for w in w_mask:\n",
+    "    print(f'[WARNING] {w.message}')\n",
+    "print(f'Success                : {result_mask[\"success\"]}')\n",
+    "print(f'Objective reduced chi² : {result_mask[\"objective_reduced_chi\"]:.6f}')\n",
+    "print(f'Objective total chi²   : {fitter_mask.objective_chi2:.6f}')\n",
+    "print(f'Classical reduced chi² : {result_mask[\"classical_reduced_chi\"]:.6f}')\n",
+    "print(f'Classical total chi²   : {fitter_mask.classical_chi2:.6f}')\n",
+    "\n",
+    "r_fit_mask = result_mask[result_model_key].values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "79755052",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "result_mask\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "24cb6fe0",
+   "metadata": {},
+   "source": [
+    "### 12. Fit with `hybrid` (new default)\n",
+    "\n",
+    "All points kept; Mighell substitution applied only to zero-variance entries."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "20c034ba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_hybrid = build_fresh_model()\n",
+    "fitter_hybrid = MultiFitter(model_hybrid, objective='hybrid')\n",
+    "\n",
+    "with warnings.catch_warnings(record=True) as w_hybrid:\n",
+    "    warnings.simplefilter('always')\n",
+    "    result_hybrid = fitter_hybrid.fit(data)\n",
+    "\n",
+    "for w in w_hybrid:\n",
+    "    print(f'[WARNING] {w.message}')\n",
+    "print(f'Success                : {result_hybrid[\"success\"]}')\n",
+    "print(f'Objective reduced chi² : {result_hybrid[\"objective_reduced_chi\"]:.6f}')\n",
+    "print(f'Objective total chi²   : {fitter_hybrid.objective_chi2:.6f}')\n",
+    "print(f'Classical reduced chi² : {result_hybrid[\"classical_reduced_chi\"]:.6f}')\n",
+    "print(f'Classical total chi²   : {fitter_hybrid.classical_chi2:.6f}')\n",
+    "\n",
+    "r_fit_hybrid = result_hybrid[result_model_key].values"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "79258c5f",
+   "metadata": {},
+   "source": [
+    "### 13. Fit with `mighell`\n",
+    "\n",
+    "Mighell transform applied to *all* points — the chi² landscape changes entirely."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "90b08fdd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_mighell = build_fresh_model()\n",
+    "fitter_mighell = MultiFitter(model_mighell, objective='mighell')\n",
+    "\n",
+    "with warnings.catch_warnings(record=True) as w_mighell:\n",
+    "    warnings.simplefilter('always')\n",
+    "    result_mighell = fitter_mighell.fit(data)\n",
+    "\n",
+    "for w in w_mighell:\n",
+    "    print(f'[WARNING] {w.message}')\n",
+    "print(f'Success                : {result_mighell[\"success\"]}')\n",
+    "print(f'Objective reduced chi² : {result_mighell[\"objective_reduced_chi\"]:.6f}')\n",
+    "print(f'Objective total chi²   : {fitter_mighell.objective_chi2:.6f}')\n",
+    "print(f'Classical reduced chi² : {result_mighell[\"classical_reduced_chi\"]:.6f}')\n",
+    "print(f'Classical total chi²   : {fitter_mighell.classical_chi2:.6f}')\n",
+    "\n",
+    "r_fit_mighell = result_mighell[result_model_key].values"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6136693d",
+   "metadata": {},
+   "source": [
+    "## 14. Compare Fit Results\n",
+    "\n",
+    "### Parameter comparison table"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "77c23990",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "models = {\n",
+    "    'legacy_mask': model_mask,\n",
+    "    'hybrid': model_hybrid,\n",
+    "    'mighell': model_mighell,\n",
+    "}\n",
+    "results = {\n",
+    "    'legacy_mask': result_mask,\n",
+    "    'hybrid': result_hybrid,\n",
+    "    'mighell': result_mighell,\n",
+    "}\n",
+    "\n",
+    "# Build descriptive labels to disambiguate duplicates (e.g. two \"thickness\" params)\n",
+    "ref_params = model_mask.get_fit_parameters()\n",
+    "seen = {}\n",
+    "labels = []\n",
+    "for p in ref_params:\n",
+    "    n = p.name\n",
+    "    seen[n] = seen.get(n, 0) + 1\n",
+    "    labels.append(f'{n} ({seen[n]})')\n",
+    "# Simplify if name appears only once\n",
+    "name_counts = {}\n",
+    "for p in ref_params:\n",
+    "    name_counts[p.name] = name_counts.get(p.name, 0) + 1\n",
+    "labels_clean = []\n",
+    "seen2 = {}\n",
+    "for p in ref_params:\n",
+    "    n = p.name\n",
+    "    seen2[n] = seen2.get(n, 0) + 1\n",
+    "    if name_counts[n] > 1:\n",
+    "        labels_clean.append(f'{n}_{seen2[n]}')\n",
+    "    else:\n",
+    "        labels_clean.append(n)\n",
+    "\n",
+    "obj_keys = ['legacy_mask', 'hybrid', 'mighell']\n",
+    "\n",
+    "header = f'{\"Parameter\":<24s} {\"legacy_mask\":>14s} {\"hybrid\":>14s} {\"mighell\":>14s}'\n",
+    "print(header)\n",
+    "print('-' * len(header))\n",
+    "\n",
+    "for i, label in enumerate(labels_clean):\n",
+    "    vals = []\n",
+    "    for obj in obj_keys:\n",
+    "        params = models[obj].get_fit_parameters()\n",
+    "        vals.append(float(params[i].value))\n",
+    "    print(f'{label:<24s} {vals[0]:>14.4g} {vals[1]:>14.4g} {vals[2]:>14.4g}')\n",
+    "\n",
+    "print('-' * len(header))\n",
+    "print(f'{\"objective red. chi²\":<24s} {result_mask[\"objective_reduced_chi\"]:>14.4f} '\n",
+    "      f'{result_hybrid[\"objective_reduced_chi\"]:>14.4f} {result_mighell[\"objective_reduced_chi\"]:>14.6f}')\n",
+    "print(f'{\"classical red. chi²\":<24s} {result_mask[\"classical_reduced_chi\"]:>14.4f} '\n",
+    "      f'{result_hybrid[\"classical_reduced_chi\"]:>14.4f} {result_mighell[\"classical_reduced_chi\"]:>14.4f}')\n",
+    "print(f'{\"success\":<24s} {str(result_mask[\"success\"]):>14s} '\n",
+    "      f'{str(result_hybrid[\"success\"]):>14s} {str(result_mighell[\"success\"]):>14s}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "906b9cdb",
+   "metadata": {},
+   "source": [
+    "### Fitted reflectivity curves — all three objectives"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a4f8f3a7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "\n",
+    "# Experimental data\n",
+    "ax.errorbar(\n",
+    "    q_raw[valid],\n",
+    "    r_raw[valid],\n",
+    "    yerr=err_raw[valid],\n",
+    "    fmt='o',\n",
+    "    ms=2,\n",
+    "    color='grey',\n",
+    "    alpha=0.5,\n",
+    "    label='Experiment',\n",
+    ")\n",
+    "ax.plot(q_raw[zero_mask], r_raw[zero_mask], 'rx', ms=10, mew=2, label='Zero-variance points')\n",
+    "\n",
+    "# Fitted curves\n",
+    "ax.plot(\n",
+    "    q_raw,\n",
+    "    r_fit_mask,\n",
+    "    '-',\n",
+    "    lw=1.5,\n",
+    "    label=(\n",
+    "        f'legacy_mask (obj={result_mask[\"objective_reduced_chi\"]:.3g}, '\n",
+    "        f'class={result_mask[\"classical_reduced_chi\"]:.3g})'\n",
+    "    ),\n",
+    ")\n",
+    "ax.plot(\n",
+    "    q_raw,\n",
+    "    r_fit_hybrid,\n",
+    "    '--',\n",
+    "    lw=1.5,\n",
+    "    label=(\n",
+    "        f'hybrid (obj={result_hybrid[\"objective_reduced_chi\"]:.3g}, '\n",
+    "        f'class={result_hybrid[\"classical_reduced_chi\"]:.3g})'\n",
+    "    ),\n",
+    ")\n",
+    "ax.plot(\n",
+    "    q_raw,\n",
+    "    r_fit_mighell,\n",
+    "    ':',\n",
+    "    lw=1.5,\n",
+    "    label=(\n",
+    "        f'mighell (obj={result_mighell[\"objective_reduced_chi\"]:.3g}, '\n",
+    "        f'class={result_mighell[\"classical_reduced_chi\"]:.3g})'\n",
+    "    ),\n",
+    ")\n",
+    "\n",
+    "ax.set_yscale('log')\n",
+    "ax.set_xlabel('Q (Å⁻¹)')\n",
+    "ax.set_ylabel('Reflectivity')\n",
+    "ax.set_title('Fitted reflectivity — comparison of objective modes')\n",
+    "ax.legend(fontsize=9)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "52821f2d",
+   "metadata": {},
+   "source": [
+    "## 15. Examine Residuals\n",
+    "\n",
+    "We compute normalised residuals $(R_\\text{exp} - R_\\text{model}) / \\sigma$ for each objective.\n",
+    "Zero-variance points are shown separately — for `legacy_mask` they were excluded from\n",
+    "the fit, while `hybrid` and `mighell` included them with transformed weights."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0292ef11",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, axes = plt.subplots(3, 1, figsize=(10, 10), sharex=True)\n",
+    "\n",
+    "fits = {\n",
+    "    'legacy_mask': r_fit_mask,\n",
+    "    'hybrid': r_fit_hybrid,\n",
+    "    'mighell': r_fit_mighell,\n",
+    "}\n",
+    "\n",
+    "for ax, (obj_name, r_fit) in zip(axes, fits.items()):\n",
+    "    # Normalised residuals for valid points\n",
+    "    residuals_valid = (r_raw[valid] - r_fit[valid]) / err_raw[valid]\n",
+    "    ax.plot(q_raw[valid], residuals_valid, 'o', ms=2, alpha=0.6, color='C0', label='Valid points')\n",
+    "\n",
+    "    # Un-normalised residuals at zero-variance points (no sigma to normalise by)\n",
+    "    residuals_zero = r_raw[zero_mask] - r_fit[zero_mask]\n",
+    "    ax.plot(q_raw[zero_mask], np.zeros_like(residuals_zero), 'rx', ms=10, mew=2,\n",
+    "            label='Zero-var Q positions')\n",
+    "\n",
+    "    ax.axhline(0, color='k', lw=0.5)\n",
+    "    ax.set_ylabel('Normalised residual')\n",
+    "    ax.set_title(f'{obj_name}')\n",
+    "    ax.legend(fontsize=8, loc='upper right')\n",
+    "\n",
+    "axes[-1].set_xlabel('Q (Å⁻¹)')\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a2f39948",
+   "metadata": {},
+   "source": [
+    "## Summary\n",
+    "\n",
+    "| Objective | Zero-var handling | Objective chi² interpretation | Classical chi² comparison |\n",
+    "|---|---|---|---|\n",
+    "| `legacy_mask` | Dropped from fit | Same as classical chi² on retained points | Directly comparable |\n",
+    "| `hybrid` | Mighell substitution for zero-var only | Slightly modified objective | Usually close to classical when zero-var fraction is small |\n",
+    "| `mighell` | Mighell transform for **all** points | **Not** a classical chi²; target and weights both change | Can look visibly worse against plotted reflectivity even when the objective is minimized well |\n",
+    "\n",
+    "The `hybrid` mode (new default) is recommended: it keeps all data points in the fit while\n",
+    "preserving a classical chi² comparison on the original positive-variance points.\n",
+    "\n",
+    "The full `mighell` mode matches the Mighell paper's objective mathematically, but that paper is\n",
+    "derived for Poisson count data. Applied to normalized reflectivity, it should be interpreted as a\n",
+    "Poisson-style objective rather than a visually comparable reduced chi² fit."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "era",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/src/easyreflectometry/fitting.py b/src/easyreflectometry/fitting.py
index 86df41e3..ffcfec03 100644
--- a/src/easyreflectometry/fitting.py
+++ b/src/easyreflectometry/fitting.py
@@ -11,14 +11,132 @@
 from easyreflectometry.data import DataSet1D
 from easyreflectometry.model import Model
 
+_VALID_OBJECTIVES = ('legacy_mask', 'mighell', 'hybrid', 'auto')
+_EPS = 1e-30
+
+
+def _validate_objective(objective: str) -> str:
+    """Validate and resolve the objective string.
+
+    :param objective: The objective mode string.
+    :type objective: str
+    :return: Resolved objective string ('auto' becomes 'hybrid').
+    :rtype: str
+    :raises ValueError: If the objective is not one of the valid options.
+    """
+    if objective not in _VALID_OBJECTIVES:
+        raise ValueError(f'Unknown objective {objective!r}. Valid options: {_VALID_OBJECTIVES}')
+    if objective == 'auto':
+        return 'hybrid'
+    return objective
+
+
+def _prepare_fit_arrays(
+    x_vals: np.ndarray,
+    y_vals: np.ndarray,
+    variances: np.ndarray,
+    objective: str,
+) -> tuple[np.ndarray, np.ndarray, np.ndarray, dict]:
+    """Prepare x, y_eff, and weights arrays for fitting based on the objective mode.
+
+    For ``legacy_mask``, zero-variance points are removed from all arrays.
+    For ``hybrid``, valid-variance points use standard WLS while zero-variance
+    points use Mighell-transformed y and weights.
+    For ``mighell``, all points use the Mighell transform.
+
+    Note: ``variances`` here means σ² (the scipp convention), not σ.
+
+    :param x_vals: Independent variable values.
+    :type x_vals: np.ndarray
+    :param y_vals: Observed dependent variable values.
+    :type y_vals: np.ndarray
+    :param variances: Variance (σ²) of each observed point.
+    :type variances: np.ndarray
+    :param objective: One of 'legacy_mask', 'hybrid', 'mighell'.
+    :type objective: str
+    :return: Tuple of (x_out, y_eff, weights, stats) where stats is a dict
+             with keys 'valid', 'mighell_substituted', 'masked'.
+    :rtype: tuple[np.ndarray, np.ndarray, np.ndarray, dict]
+    """
+    n = len(y_vals)
+    zero_mask = variances <= 0.0
+    n_zero = int(np.sum(zero_mask))
+    n_valid = n - n_zero
+
+    if objective == 'legacy_mask':
+        valid = ~zero_mask
+        x_out = x_vals[valid]
+        y_eff = y_vals[valid]
+        if n_valid > 0:
+            weights = 1.0 / np.sqrt(variances[valid])
+        else:
+            weights = np.array([])
+        stats = {'valid': n_valid, 'mighell_substituted': 0, 'masked': n_zero, 'transformed_all_points': False}
+        return x_out, y_eff, weights, stats
+
+    # hybrid or mighell
+    y_eff = np.copy(y_vals)
+    sigma = np.empty(n)
+
+    if objective == 'mighell':
+        apply_mighell = np.ones(n, dtype=bool)
+    else:
+        # hybrid: apply Mighell only to zero-variance points
+        apply_mighell = zero_mask
+
+    # Standard WLS for non-Mighell points
+    standard = ~apply_mighell
+    if np.any(standard):
+        sigma[standard] = np.sqrt(variances[standard])
+
+    # Mighell transform for selected points
+    if np.any(apply_mighell):
+        y_m = y_vals[apply_mighell]
+        delta = np.minimum(y_m, 1.0)
+        y_eff[apply_mighell] = y_m + delta
+        sigma[apply_mighell] = np.sqrt(np.maximum(y_m + 1.0, _EPS))
+
+    weights = 1.0 / sigma
+    n_mighell = int(np.sum(apply_mighell))
+    stats = {
+        'valid': n - n_mighell,
+        'mighell_substituted': n_mighell,
+        'masked': 0,
+        'transformed_all_points': bool(objective == 'mighell'),
+    }
+    return x_vals, y_eff, weights, stats
+
+
+def _compute_weighted_chi2(y_obs: np.ndarray, y_calc: np.ndarray, sigma: np.ndarray) -> float:
+    """Return weighted chi-square for finite, strictly positive uncertainties."""
+    valid = np.isfinite(y_obs) & np.isfinite(y_calc) & np.isfinite(sigma) & (sigma > 0.0)
+    if not np.any(valid):
+        return 0.0
+    residual = (y_obs[valid] - y_calc[valid]) / sigma[valid]
+    return float(np.sum(residual**2))
+
+
+def _compute_reduced_chi2(chi2: float, n_points: int, n_params: int) -> float | None:
+    """Return reduced chi-square or None when degrees of freedom are not positive."""
+    dof = int(n_points) - int(n_params)
+    if dof <= 0:
+        return None
+    return float(chi2 / dof)
+
 
 class MultiFitter:
-    def __init__(self, *args: Model):
-        r"""A convinence class for the :py:class:`easyscience.Fitting.Fitting`
+    def __init__(self, *args: Model, objective: str = 'hybrid'):
+        r"""A convenience class for the :py:class:`easyscience.Fitting.Fitting`
         which will populate the :py:class:`sc.DataGroup` appropriately
         after the fitting is performed.
 
-        :param args: Reflectometry model
+        :param args: Reflectometry model(s).
+        :param objective: Zero-variance handling strategy. One of
+            ``'hybrid'`` (default, Mighell for zero-variance, WLS otherwise),
+            ``'mighell'`` (Mighell transform for all points),
+            ``'legacy_mask'`` (drop zero-variance points),
+            ``'auto'`` (alias for ``'hybrid'``).
+        :type objective: str
         """
 
         # This lets the unique_name be passed with the fit_func.
@@ -32,22 +150,35 @@ def wrapped(*args, **kwargs):
         self._models = args
         self.easy_science_multi_fitter = EasyScienceMultiFitter(args, self._fit_func)
         self._fit_results: list[FitResults] | None = None
-
-    def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup:
+        self._classical_fit_metrics: list[dict] | None = None
+        self._objective = _validate_objective(objective)
+
+    def fit(self, data: sc.DataGroup, id: int = 0, objective: str | None = None) -> sc.DataGroup:
+        """Perform the fitting and populate the DataGroups with the result.
+
+        :param data: DataGroup to be fitted to and populated.
+        :type data: sc.DataGroup
+        :param id: Unused parameter kept for backward compatibility.
+        :type id: int
+        :param objective: Per-call override for the zero-variance objective.
+            If ``None``, uses the instance default set at construction.
+        :type objective: str or None
+        :return: A new DataGroup with fitted model curves, SLD profiles, and fit statistics.
+        :rtype: sc.DataGroup
+
+        :note: Under the ``mighell`` objective all points are transformed,
+               so ``reduced_chi`` is not a classical chi-square statistic.
+               Under ``hybrid``, only zero-variance points are transformed;
+               when they are a small fraction of the data the chi-square
+               remains approximately classical.
         """
-        Perform the fitting and populate the DataGroups with the result.
-
-        :param data: DataGroup to be fitted to and populated
-        :param method: Optimisation method
+        obj = _validate_objective(objective) if objective is not None else self._objective
 
-        :note: Points with zero variance in the data will be automatically masked
-               out during fitting. A warning will be issued if any such points
-               are found, indicating the number of points masked per reflectivity.
-        """
         refl_nums = [k[3:] for k in data['coords'].keys() if 'Qz' == k[:2]]
         x = []
         y = []
         dy = []
+        original_arrays = []
 
         # Process each reflectivity dataset
         for i in refl_nums:
@@ -55,34 +186,40 @@ def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup:
             y_vals = data['data'][f'R_{i}'].values
             variances = data['data'][f'R_{i}'].variances
 
-            # Find points with non-zero variance
-            zero_variance_mask = variances == 0.0
-            num_zero_variance = np.sum(zero_variance_mask)
+            x_out, y_eff, weights, stats = _prepare_fit_arrays(x_vals, y_vals, variances, obj)
 
-            if num_zero_variance > 0:
+            if stats['masked'] > 0:
                 warnings.warn(
-                    f'Masked {num_zero_variance} data point(s) in reflectivity {i} due to zero variance during fitting.',
+                    f'Masked {stats["masked"]} data point(s) in reflectivity {i} '
+                    'due to zero variance during fitting.',
+                    UserWarning,
+                )
+            if stats.get('transformed_all_points'):
+                warnings.warn(
+                    f'Applied Mighell transform to all {len(y_vals)} point(s) '
+                    f'in reflectivity {i} during fitting.',
+                    UserWarning,
+                )
+            elif stats['mighell_substituted'] > 0:
+                warnings.warn(
+                    f'Applied Mighell substitution to {stats["mighell_substituted"]} '
+                    f'zero-variance point(s) in reflectivity {i} during fitting.',
                     UserWarning,
                 )
 
-            # Keep only points with non-zero variances
-            valid_mask = ~zero_variance_mask
-            x_vals_masked = x_vals[valid_mask]
-            y_vals_masked = y_vals[valid_mask]
-            variances_masked = variances[valid_mask]
-
-            x.append(x_vals_masked)
-            y.append(y_vals_masked)
-            dy.append(1 / np.sqrt(variances_masked))
+            x.append(x_out)
+            y.append(y_eff)
+            dy.append(weights)
+            original_arrays.append({'x': x_vals, 'y': y_vals, 'variances': variances})
 
         result = self.easy_science_multi_fitter.fit(x, y, weights=dy)
         self._fit_results = result
+        self._classical_fit_metrics = []
         new_data = data.copy()
         for i, _ in enumerate(result):
             id = refl_nums[i]
-            new_data[f'R_{id}_model'] = sc.array(
-                dims=[f'Qz_{id}'], values=self._fit_func[i](data['coords'][f'Qz_{id}'].values)
-            )
+            model_curve = self._fit_func[i](data['coords'][f'Qz_{id}'].values)
+            new_data[f'R_{id}_model'] = sc.array(dims=[f'Qz_{id}'], values=model_curve)
             sld_profile = self.easy_science_multi_fitter._fit_objects[i].interface.sld_profile(self._models[i].unique_name)
             new_data[f'SLD_{id}'] = sc.array(dims=[f'z_{id}'], values=sld_profile[1] * 1e-6, unit=sc.Unit('1/angstrom') ** 2)
             if 'attrs' in new_data:
@@ -90,41 +227,90 @@ def fit(self, data: sc.DataGroup, id: int = 0) -> sc.DataGroup:
             new_data['coords'][f'z_{id}'] = sc.array(
                 dims=[f'z_{id}'], values=sld_profile[0], unit=(1 / new_data['coords'][f'Qz_{id}'].unit).unit
             )
+            original = original_arrays[i]
+            sigma_classical = np.sqrt(np.clip(original['variances'], 0.0, None))
+            n_classical_points = int(np.sum(original['variances'] > 0.0))
+            classical_chi2 = _compute_weighted_chi2(original['y'], model_curve, sigma_classical)
+            classical_reduced_chi = _compute_reduced_chi2(classical_chi2, n_classical_points, result[i].n_pars)
+            objective_chi2 = float(result[i].chi2)
+            objective_reduced_chi = float(result[i].reduced_chi)
+
+            self._classical_fit_metrics.append(
+                {
+                    'classical_chi2': classical_chi2,
+                    'classical_reduced_chi': classical_reduced_chi,
+                    'objective_chi2': objective_chi2,
+                    'objective_reduced_chi': objective_reduced_chi,
+                    'n_classical_points': n_classical_points,
+                }
+            )
+
+            new_data['objective_chi2'] = objective_chi2
+            new_data['objective_reduced_chi'] = objective_reduced_chi
+            new_data['classical_chi2'] = classical_chi2
+            new_data['classical_reduced_chi'] = classical_reduced_chi
             new_data['reduced_chi'] = float(result[i].reduced_chi)
             new_data['success'] = result[i].success
         return new_data
 
-    def fit_single_data_set_1d(self, data: DataSet1D) -> FitResults:
+    def fit_single_data_set_1d(self, data: DataSet1D, objective: str | None = None) -> FitResults:
+        """Perform fitting on a single 1D dataset.
+
+        :param data: The 1D dataset to fit. Note that ``data.ye`` stores
+            variances (σ²), not standard deviations.
+        :type data: DataSet1D
+        :param objective: Per-call override for the zero-variance objective.
+            If ``None``, uses the instance default set at construction.
+        :type objective: str or None
+        :return: Fit results from the minimizer.
+        :rtype: FitResults
         """
-        Perform the fitting and populate the DataGroups with the result.
+        obj = _validate_objective(objective) if objective is not None else self._objective
 
-        :param data: DataGroup to be fitted to and populated
-        :param method: Optimisation method
-        """
         x_vals = np.asarray(data.x)
         y_vals = np.asarray(data.y)
         variances = np.asarray(data.ye)
 
-        zero_variance_mask = variances == 0.0
-        num_zero_variance = int(np.sum(zero_variance_mask))
+        x_out, y_eff, weights, stats = _prepare_fit_arrays(x_vals, y_vals, variances, obj)
 
-        if num_zero_variance > 0:
+        if stats['masked'] > 0:
+            warnings.warn(
+                f'Masked {stats["masked"]} data point(s) in single-dataset fit '
+                'due to zero variance during fitting.',
+                UserWarning,
+            )
+        if stats.get('transformed_all_points'):
+            warnings.warn(
+                f'Applied Mighell transform to all {len(y_vals)} point(s) '
+                'in single-dataset fit during fitting.',
+                UserWarning,
+            )
+        elif stats['mighell_substituted'] > 0:
             warnings.warn(
-                f'Masked {num_zero_variance} data point(s) in single-dataset fit due to zero variance during fitting.',
+                f'Applied Mighell substitution to {stats["mighell_substituted"]} '
+                'zero-variance point(s) in single-dataset fit during fitting.',
                 UserWarning,
             )
 
-        valid_mask = ~zero_variance_mask
-        if not np.any(valid_mask):
+        if obj == 'legacy_mask' and len(x_out) == 0:
             raise ValueError('Cannot fit single dataset: all points have zero variance.')
 
-        x_vals_masked = x_vals[valid_mask]
-        y_vals_masked = y_vals[valid_mask]
-        variances_masked = variances[valid_mask]
-
-        weights = 1.0 / np.sqrt(variances_masked)
-        result = self.easy_science_multi_fitter.fit(x=[x_vals_masked], y=[y_vals_masked], weights=[weights])[0]
+        result = self.easy_science_multi_fitter.fit(x=[x_out], y=[y_eff], weights=[weights])[0]
         self._fit_results = [result]
+        sigma_classical = np.sqrt(np.clip(variances, 0.0, None))
+        model_curve = self._fit_func[0](x_vals)
+        n_classical_points = int(np.sum(variances > 0.0))
+        classical_chi2 = _compute_weighted_chi2(y_vals, model_curve, sigma_classical)
+        classical_reduced_chi = _compute_reduced_chi2(classical_chi2, n_classical_points, result.n_pars)
+        self._classical_fit_metrics = [
+            {
+                'classical_chi2': classical_chi2,
+                'classical_reduced_chi': classical_reduced_chi,
+                'objective_chi2': float(result.chi2),
+                'objective_reduced_chi': float(result.reduced_chi),
+                'n_classical_points': n_classical_points,
+            }
+        ]
         return result
 
     @property
@@ -149,6 +335,33 @@ def reduced_chi(self) -> float | None:
 
         return total_chi2 / total_dof
 
+    @property
+    def classical_chi2(self) -> float | None:
+        """Classical chi-squared using only points with positive variances."""
+        if self._classical_fit_metrics is None:
+            return None
+        return float(sum(metric['classical_chi2'] for metric in self._classical_fit_metrics))
+
+    @property
+    def classical_reduced_chi(self) -> float | None:
+        """Reduced classical chi-squared using only points with positive variances."""
+        if self._classical_fit_metrics is None or self._fit_results is None:
+            return None
+        total_chi2 = self.classical_chi2
+        total_points = sum(metric['n_classical_points'] for metric in self._classical_fit_metrics)
+        n_params = self._fit_results[0].n_pars
+        return _compute_reduced_chi2(total_chi2, total_points, n_params)
+
+    @property
+    def objective_chi2(self) -> float | None:
+        """Objective-space chi-squared returned by the minimizer."""
+        return self.chi2
+
+    @property
+    def objective_reduced_chi(self) -> float | None:
+        """Objective-space reduced chi-squared returned by the minimizer."""
+        return self.reduced_chi
+
     def switch_minimizer(self, minimizer: AvailableMinimizers) -> None:
         """
         Switch the minimizer for the fitting.
diff --git a/tests/test_fitting.py b/tests/test_fitting.py
index 446f10b9..e4145f45 100644
--- a/tests/test_fitting.py
+++ b/tests/test_fitting.py
@@ -5,6 +5,7 @@
 
 import numpy as np
 import pytest
+import scipp as sc
 from easyscience.fitting.minimizers.factory import AvailableMinimizers
 
 import easyreflectometry
@@ -12,6 +13,8 @@
 from easyreflectometry.data import DataSet1D
 from easyreflectometry.data.measurement import load
 from easyreflectometry.fitting import MultiFitter
+from easyreflectometry.fitting import _prepare_fit_arrays
+from easyreflectometry.fitting import _validate_objective
 from easyreflectometry.model import Model
 from easyreflectometry.model import PercentageFwhm
 from easyreflectometry.sample import Layer
@@ -76,7 +79,7 @@ def test_fitting(minimizer):
 
 
 def test_fitting_with_zero_variance():
-    """Test that zero variance points are properly detected and masked during fitting when present in the data."""
+    """Test that zero variance points are handled via Mighell substitution (hybrid default)."""
     import warnings
 
     import numpy as np
@@ -129,26 +132,24 @@ def test_fitting_with_zero_variance():
     model.interface = interface
     fitter = MultiFitter(model)
 
-    # Capture warnings during fitting - check if zero variance points still exist in the data
-    # and are properly handled by the fitting method
+    # Capture warnings during fitting
     with warnings.catch_warnings(record=True) as w:
         warnings.simplefilter('always')
         analysed = fitter.fit(data)
 
-        # Check if any zero variance warnings were issued during fitting
-        fitting_warnings = [str(warning.message) for warning in w if 'zero variance during fitting' in str(warning.message)]
+        # Under hybrid default, zero variance points trigger Mighell substitution warnings
+        mighell_warnings = [str(warning.message) for warning in w if 'Mighell substitution' in str(warning.message)]
+        mask_warnings = [str(warning.message) for warning in w if 'Masked' in str(warning.message)]
+
+        # Hybrid mode should NOT produce mask warnings
+        assert len(mask_warnings) == 0, f'Unexpected mask warnings under hybrid: {mask_warnings}'
 
-        # The fitting method should handle zero variance points gracefully
-        # If there are any zero variance points remaining in the data, they should be masked
-        # and a warning should be issued
-        if len(fitting_warnings) > 0:
-            # Verify the warning message format and that it mentions masking points
-            for warning_msg in fitting_warnings:
-                assert 'Masked' in warning_msg and 'zero variance during fitting' in warning_msg
-                print(f'Info: {warning_msg}')  # Log for debugging
+        # If there are zero-variance points in the loaded data, Mighell warnings should appear
+        if len(mighell_warnings) > 0:
+            for warning_msg in mighell_warnings:
+                assert 'zero-variance point(s)' in warning_msg
 
     # Basic checks that fitting completed
-    # The keys will be based on the filename, not just '0'
     model_keys = [k for k in analysed.keys() if k.endswith('_model')]
     sld_keys = [k for k in analysed.keys() if k.startswith('SLD_')]
     assert len(model_keys) > 0, f'No model keys found in {list(analysed.keys())}'
@@ -157,7 +158,7 @@ def test_fitting_with_zero_variance():
 
 
 def test_fitting_with_manual_zero_variance():
-    """Test the fit method with manually created zero variance points."""
+    """Test the fit method with manually created zero variance points using hybrid (default)."""
     import warnings
 
     import numpy as np
@@ -217,12 +218,12 @@ def test_fitting_with_manual_zero_variance():
         warnings.simplefilter('always')
         analysed = fitter.fit(data)
 
-        # Check that warnings were issued about zero variance points
-        fitting_warnings = [str(warning.message) for warning in w if 'zero variance during fitting' in str(warning.message)]
+        # Under hybrid default, should get Mighell substitution warning, not masking
+        mighell_warnings = [str(warning.message) for warning in w if 'Mighell substitution' in str(warning.message)]
+
+        assert len(mighell_warnings) == 1, f'Expected 1 Mighell warning, got {len(mighell_warnings)}: {mighell_warnings}'
+        assert '7 zero-variance point(s)' in mighell_warnings[0], f'Unexpected warning content: {mighell_warnings[0]}'
 
-        # Should have one warning about the 7 zero variance points (5 + 2)
-        assert len(fitting_warnings) == 1, f'Expected 1 warning, got {len(fitting_warnings)}: {fitting_warnings}'
-        assert 'Masked 7 data point(s)' in fitting_warnings[0], f'Unexpected warning content: {fitting_warnings[0]}'
     # Basic checks that fitting completed despite zero variance points
     assert 'R_0_model' in analysed.keys()
     assert 'SLD_0' in analysed.keys()
@@ -230,9 +231,10 @@ def test_fitting_with_manual_zero_variance():
 
 
 def test_fit_single_data_set_1d_masks_zero_variance_points():
+    """Legacy mask mode: zero-variance points are dropped."""
     model = Model()
     model.interface = CalculatorFactory()
-    fitter = MultiFitter(model)
+    fitter = MultiFitter(model, objective='legacy_mask')
 
     captured = {}
     mock_result = MagicMock()
@@ -255,7 +257,7 @@ def _fake_fit(*, x, y, weights):
         ye=np.array([0.01, 0.0, 0.04]),
     )
 
-    with pytest.warns(UserWarning, match='Masked 1 data point\(s\) in single-dataset fit'):
+    with pytest.warns(UserWarning, match='Masked 1 data point\\(s\\) in single-dataset fit'):
         result = fitter.fit_single_data_set_1d(data)
 
     assert result is mock_result
@@ -286,9 +288,10 @@ def test_reduced_chi_uses_global_dof_across_fit_results():
 
 
 def test_fit_single_data_set_1d_all_zero_variance_raises():
+    """Legacy mask mode raises when all points have zero variance."""
     model = Model()
     model.interface = CalculatorFactory()
-    fitter = MultiFitter(model)
+    fitter = MultiFitter(model, objective='legacy_mask')
 
     data = DataSet1D(
         name='all_zero',
@@ -376,3 +379,407 @@ def _fake_fit(*, x, y, weights):
     assert result is mock_result
     assert np.allclose(captured['x'][0], np.array([0.01, 0.02, 0.03]))
     assert np.allclose(captured['y'][0], np.array([1.0, 0.8, 0.6]))
+
+
+# --- New tests for objective-based zero-variance handling ---
+
+
+def test_objective_validation_rejects_unknown_value():
+    with pytest.raises(ValueError, match='Unknown objective'):
+        _validate_objective('bad_value')
+
+
+def test_objective_validation_resolves_auto():
+    assert _validate_objective('auto') == 'hybrid'
+    assert _validate_objective('hybrid') == 'hybrid'
+    assert _validate_objective('legacy_mask') == 'legacy_mask'
+    assert _validate_objective('mighell') == 'mighell'
+
+
+def test_prepare_fit_arrays_weights_always_positive_and_finite():
+    """Weights must be strictly positive and finite for all inputs and objectives."""
+    test_cases = [
+        # (y_vals, variances, description)
+        (np.array([0.0]), np.array([0.0]), 'y=0, var=0'),
+        (np.array([-0.5]), np.array([0.0]), 'y=-0.5, var=0'),
+        (np.array([-1.0]), np.array([0.0]), 'y=-1, var=0'),
+        (np.array([1e6]), np.array([0.0]), 'y=1e6, var=0'),
+        (np.array([0.5, 0.3, 0.1]), np.array([0.0, 0.0, 0.0]), 'all-zero variances'),
+        (np.array([0.5, 0.3, 0.1]), np.array([0.01, 0.0, 0.04]), 'mixed variances'),
+        (np.array([0.0, -0.5, -1.0, 1e6]), np.array([0.0, 0.0, 0.0, 0.0]), 'edge y values'),
+    ]
+
+    for objective in ('hybrid', 'mighell'):
+        for y_vals, variances, desc in test_cases:
+            x = np.arange(len(y_vals), dtype=float)
+            _, _, weights, _ = _prepare_fit_arrays(x, y_vals, variances, objective)
+            assert len(weights) == len(y_vals), f'Wrong length for {desc}, {objective}'
+            assert np.all(weights > 0), f'Non-positive weight for {desc}, {objective}: {weights}'
+            assert np.all(np.isfinite(weights)), f'Non-finite weight for {desc}, {objective}: {weights}'
+
+
+def test_prepare_fit_arrays_legacy_mask_drops_zero_variance():
+    x = np.array([0.01, 0.02, 0.03])
+    y = np.array([1.0, 0.8, 0.6])
+    var = np.array([0.01, 0.0, 0.04])
+
+    x_out, y_eff, weights, stats = _prepare_fit_arrays(x, y, var, 'legacy_mask')
+
+    assert np.allclose(x_out, [0.01, 0.03])
+    assert np.allclose(y_eff, [1.0, 0.6])
+    assert np.allclose(weights, [1.0 / np.sqrt(0.01), 1.0 / np.sqrt(0.04)])
+    assert stats == {'valid': 2, 'mighell_substituted': 0, 'masked': 1, 'transformed_all_points': False}
+
+
+def test_prepare_fit_arrays_hybrid_transforms_zero_variance():
+    x = np.array([0.01, 0.02, 0.03])
+    y = np.array([1.0, 0.8, 0.6])
+    var = np.array([0.01, 0.0, 0.04])
+
+    x_out, y_eff, weights, stats = _prepare_fit_arrays(x, y, var, 'hybrid')
+
+    # x unchanged
+    assert np.allclose(x_out, x)
+    # Index 0 and 2: standard WLS (unchanged y)
+    assert y_eff[0] == pytest.approx(1.0)
+    assert y_eff[2] == pytest.approx(0.6)
+    assert weights[0] == pytest.approx(1.0 / np.sqrt(0.01))
+    assert weights[2] == pytest.approx(1.0 / np.sqrt(0.04))
+    # Index 1: Mighell transform — y_eff = y + min(y, 1) = 0.8 + 0.8 = 1.6
+    assert y_eff[1] == pytest.approx(0.8 + 0.8)
+    # sigma = sqrt(y + 1) = sqrt(1.8)
+    assert weights[1] == pytest.approx(1.0 / np.sqrt(1.8))
+    assert stats == {'valid': 2, 'mighell_substituted': 1, 'masked': 0, 'transformed_all_points': False}
+
+
+def test_prepare_fit_arrays_mighell_transforms_all():
+    x = np.array([0.01, 0.02])
+    y = np.array([0.5, 0.3])
+    var = np.array([0.01, 0.04])  # All valid, but mighell transforms everything
+
+    x_out, y_eff, weights, stats = _prepare_fit_arrays(x, y, var, 'mighell')
+
+    assert np.allclose(x_out, x)
+    # y_eff = y + min(y, 1) = y + y (since y < 1)
+    assert y_eff[0] == pytest.approx(0.5 + 0.5)
+    assert y_eff[1] == pytest.approx(0.3 + 0.3)
+    # sigma = sqrt(y + 1)
+    assert weights[0] == pytest.approx(1.0 / np.sqrt(1.5))
+    assert weights[1] == pytest.approx(1.0 / np.sqrt(1.3))
+    assert stats == {'valid': 0, 'mighell_substituted': 2, 'masked': 0, 'transformed_all_points': True}
+
+
+def test_fit_single_data_set_1d_hybrid_keeps_zero_variance_points():
+    """Hybrid mode keeps all points (transforms zero-variance ones)."""
+    model = Model()
+    model.interface = CalculatorFactory()
+    fitter = MultiFitter(model)  # default objective='hybrid'
+
+    captured = {}
+    mock_result = MagicMock()
+    mock_result.chi2 = 1.0
+    mock_result.n_pars = 1
+
+    def _fake_fit(*, x, y, weights):
+        captured['x'] = x
+        captured['y'] = y
+        captured['weights'] = weights
+        return [mock_result]
+
+    fitter.easy_science_multi_fitter = MagicMock()
+    fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit)
+
+    data = DataSet1D(
+        name='hybrid_test',
+        x=np.array([0.01, 0.02, 0.03]),
+        y=np.array([1.0, 0.8, 0.6]),
+        ye=np.array([0.01, 0.0, 0.04]),
+    )
+
+    with pytest.warns(UserWarning, match='Mighell substitution'):
+        result = fitter.fit_single_data_set_1d(data)
+
+    assert result is mock_result
+    # All 3 points should be passed through (not masked)
+    assert len(captured['x'][0]) == 3
+    assert len(captured['y'][0]) == 3
+    assert len(captured['weights'][0]) == 3
+
+
+def test_fit_single_data_set_1d_mighell_warning_mentions_all_points():
+    model = Model()
+    model.interface = CalculatorFactory()
+    fitter = MultiFitter(model, objective='mighell')
+
+    mock_result = MagicMock()
+    mock_result.chi2 = 1.0
+    mock_result.reduced_chi = 0.5
+    mock_result.n_pars = 1
+
+    fitter.easy_science_multi_fitter = MagicMock()
+    fitter.easy_science_multi_fitter.fit = MagicMock(return_value=[mock_result])
+    fitter._fit_func = [lambda x: np.zeros_like(x)]
+
+    data = DataSet1D(
+        name='mighell_warning',
+        x=np.array([0.01, 0.02, 0.03]),
+        y=np.array([1.0, 0.8, 0.6]),
+        ye=np.array([0.01, 0.02, 0.04]),
+    )
+
+    with pytest.warns(UserWarning, match=r'Applied Mighell transform to all 3 point\(s\)'):
+        fitter.fit_single_data_set_1d(data)
+
+
+def test_classical_and_objective_chi_are_split_for_fit_results():
+    model = Model()
+    model.interface = CalculatorFactory()
+    fitter = MultiFitter(model, objective='mighell')
+
+    fit_result = MagicMock()
+    fit_result.chi2 = 0.25
+    fit_result.reduced_chi = 0.125
+    fit_result.n_pars = 1
+    fit_result.x = np.array([0.01, 0.02, 0.03])
+
+    fitter.easy_science_multi_fitter = MagicMock()
+    fitter.easy_science_multi_fitter.fit = MagicMock(return_value=[fit_result])
+    fitter.easy_science_multi_fitter._fit_objects = [MagicMock(interface=MagicMock())]
+    fitter.easy_science_multi_fitter._fit_objects[0].interface.sld_profile.return_value = \
+        (np.array([0.0, 1.0]), np.array([1.0, 2.0]))
+
+    fitter._models = [MagicMock(unique_name='model_0',
+                                as_dict=MagicMock(return_value={'name': 'model_0'}))]
+    fitter._fit_func = [lambda x: np.array([0.8, 0.75, 0.7])]
+
+    data = sc.DataGroup(
+        {
+            'coords': {'Qz_0': sc.array(dims=['Qz_0'], values=np.array([0.01, 0.02, 0.03]), unit=sc.Unit('1/angstrom'))},
+            'data': {'R_0': sc.array(dims=['Qz_0'], values=np.array([1.0, 0.9, 0.7]), variances=np.array([0.01, 0.0, 0.04]))},
+            'attrs': {},
+        }
+    )
+
+    analysed = fitter.fit(data)
+
+    expected_classical_chi2 = ((1.0 - 0.8) / 0.1) ** 2 + ((0.7 - 0.7) / 0.2) ** 2
+    expected_classical_reduced = expected_classical_chi2 / (2 - fit_result.n_pars)
+
+    assert analysed['objective_chi2'] == pytest.approx(0.25)
+    assert analysed['objective_reduced_chi'] == pytest.approx(0.125)
+    assert analysed['classical_chi2'] == pytest.approx(expected_classical_chi2)
+    assert analysed['classical_reduced_chi'] == pytest.approx(expected_classical_reduced)
+    assert fitter.objective_chi2 == pytest.approx(0.25)
+    assert fitter.objective_reduced_chi == pytest.approx(0.125)
+    assert fitter.classical_chi2 == pytest.approx(expected_classical_chi2)
+    assert fitter.classical_reduced_chi == pytest.approx(expected_classical_reduced)
+
+
+def test_fit_single_data_set_1d_all_zero_variance_hybrid_does_not_raise():
+    """Hybrid mode handles all-zero-variance data without raising."""
+    model = Model()
+    model.interface = CalculatorFactory()
+    fitter = MultiFitter(model)  # default objective='hybrid'
+
+    captured = {}
+    mock_result = MagicMock()
+    mock_result.chi2 = 1.0
+    mock_result.n_pars = 1
+
+    def _fake_fit(*, x, y, weights):
+        captured['x'] = x
+        captured['y'] = y
+        captured['weights'] = weights
+        return [mock_result]
+
+    fitter.easy_science_multi_fitter = MagicMock()
+    fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit)
+
+    data = DataSet1D(
+        name='all_zero_hybrid',
+        x=np.array([0.01, 0.02, 0.03]),
+        y=np.array([1.0, 0.8, 0.6]),
+        ye=np.array([0.0, 0.0, 0.0]),
+    )
+
+    with pytest.warns(UserWarning, match='Mighell substitution'):
+        result = fitter.fit_single_data_set_1d(data)
+
+    assert result is mock_result
+    assert len(captured['x'][0]) == 3
+
+
+def test_fit_single_data_set_1d_legacy_mask_preserves_old_behavior():
+    """Legacy mask mode drops zero-variance points and warns with old message."""
+    model = Model()
+    model.interface = CalculatorFactory()
+    fitter = MultiFitter(model, objective='legacy_mask')
+
+    captured = {}
+    mock_result = MagicMock()
+    mock_result.chi2 = 1.0
+    mock_result.n_pars = 1
+
+    def _fake_fit(*, x, y, weights):
+        captured['x'] = x
+        captured['y'] = y
+        captured['weights'] = weights
+        return [mock_result]
+
+    fitter.easy_science_multi_fitter = MagicMock()
+    fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit)
+
+    data = DataSet1D(
+        name='legacy_test',
+        x=np.array([0.01, 0.02, 0.03]),
+        y=np.array([1.0, 0.8, 0.6]),
+        ye=np.array([0.01, 0.0, 0.04]),
+    )
+
+    with pytest.warns(UserWarning, match='Masked 1 data point'):
+        result = fitter.fit_single_data_set_1d(data)
+
+    assert result is mock_result
+    assert np.allclose(captured['x'][0], np.array([0.01, 0.03]))
+    assert np.allclose(captured['y'][0], np.array([1.0, 0.6]))
+
+
+def test_fit_multi_dataset_hybrid_uses_transformed_y_and_weights():
+    """Multi-dataset fit with hybrid objective transforms zero-variance points."""
+    import scipp as sc
+
+    qz_values = np.linspace(0.01, 0.3, 10)
+    r_values = np.exp(-qz_values * 50)
+    variances = np.ones_like(r_values) * 0.01
+    variances[3:5] = 0.0  # 2 zero-variance points
+
+    data = sc.DataGroup(
+        {
+            'coords': {'Qz_0': sc.array(dims=['Qz_0'], values=qz_values)},
+            'data': {'R_0': sc.array(dims=['Qz_0'], values=r_values, variances=variances)},
+        }
+    )
+
+    model = Model()
+    model.interface = CalculatorFactory()
+    fitter = MultiFitter(model)
+
+    captured = {}
+
+    def _fake_fit(x, y, weights):
+        captured['x'] = x
+        captured['y'] = y
+        captured['weights'] = weights
+        mock_r = MagicMock()
+        mock_r.reduced_chi = 1.0
+        mock_r.success = True
+        mock_r.chi2 = 1.0
+        mock_r.n_pars = 1
+        mock_r.x = x[0]
+        return [mock_r]
+
+    fitter.easy_science_multi_fitter = MagicMock()
+    fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit)
+    fitter.easy_science_multi_fitter._fit_objects = [MagicMock()]
+    fitter.easy_science_multi_fitter._fit_objects[0].interface.sld_profile.return_value = (
+        np.linspace(0, 100, 5),
+        np.ones(5),
+    )
+
+    import warnings
+
+    with warnings.catch_warnings(record=True) as w:
+        warnings.simplefilter('always')
+        fitter.fit(data)
+
+    # All 10 points should be present (not masked)
+    assert len(captured['x'][0]) == 10
+    assert len(captured['y'][0]) == 10
+    assert len(captured['weights'][0]) == 10
+
+    # Zero-variance points should have Mighell-transformed y
+    for idx in [3, 4]:
+        y_orig = r_values[idx]
+        expected_y_eff = y_orig + min(y_orig, 1.0)
+        assert captured['y'][0][idx] == pytest.approx(expected_y_eff)
+
+    # Check that Mighell warning was emitted
+    mighell_warnings = [str(ww.message) for ww in w if 'Mighell substitution' in str(ww.message)]
+    assert len(mighell_warnings) == 1
+    assert '2 zero-variance point(s)' in mighell_warnings[0]
+
+
+def test_fit_warnings_objective_specific():
+    """Verify that each objective mode produces the correct warning type."""
+    import warnings
+
+    model = Model()
+    model.interface = CalculatorFactory()
+
+    mock_result = MagicMock()
+    mock_result.chi2 = 1.0
+    mock_result.n_pars = 1
+
+    data = DataSet1D(
+        name='warn_test',
+        x=np.array([0.01, 0.02, 0.03]),
+        y=np.array([1.0, 0.8, 0.6]),
+        ye=np.array([0.01, 0.0, 0.04]),
+    )
+
+    for obj, expected_fragment in [
+        ('legacy_mask', 'Masked 1 data point(s)'),
+        ('hybrid', 'Mighell substitution'),
+        ('mighell', 'Applied Mighell transform to all 3 point(s)'),
+    ]:
+        fitter = MultiFitter(model, objective=obj)
+        fitter.easy_science_multi_fitter = MagicMock()
+        fitter.easy_science_multi_fitter.fit = MagicMock(return_value=[mock_result])
+
+        with warnings.catch_warnings(record=True) as w:
+            warnings.simplefilter('always')
+            fitter.fit_single_data_set_1d(data)
+
+        matching = [str(ww.message) for ww in w if expected_fragment in str(ww.message)]
+        assert len(matching) > 0, f'No warning containing {expected_fragment!r} for objective={obj}'
+
+
+def test_multifitter_constructor_rejects_bad_objective():
+    model = Model()
+    model.interface = CalculatorFactory()
+    with pytest.raises(ValueError, match='Unknown objective'):
+        MultiFitter(model, objective='nonsense')
+
+
+def test_fit_per_call_objective_override():
+    """Per-call objective override in fit_single_data_set_1d works."""
+    model = Model()
+    model.interface = CalculatorFactory()
+    fitter = MultiFitter(model, objective='hybrid')  # default
+
+    captured = {}
+    mock_result = MagicMock()
+    mock_result.chi2 = 1.0
+    mock_result.n_pars = 1
+
+    def _fake_fit(*, x, y, weights):
+        captured['x'] = x
+        captured['y'] = y
+        captured['weights'] = weights
+        return [mock_result]
+
+    fitter.easy_science_multi_fitter = MagicMock()
+    fitter.easy_science_multi_fitter.fit = MagicMock(side_effect=_fake_fit)
+
+    data = DataSet1D(
+        name='override_test',
+        x=np.array([0.01, 0.02, 0.03]),
+        y=np.array([1.0, 0.8, 0.6]),
+        ye=np.array([0.01, 0.0, 0.04]),
+    )
+
+    # Override to legacy_mask — should drop the zero-variance point
+    with pytest.warns(UserWarning, match='Masked 1 data point'):
+        fitter.fit_single_data_set_1d(data, objective='legacy_mask')
+
+    assert len(captured['x'][0]) == 2  # one point dropped
diff --git a/zero_variances_example.py b/zero_variances_example.py
new file mode 100644
index 00000000..fb2b42c6
--- /dev/null
+++ b/zero_variances_example.py
@@ -0,0 +1,222 @@
+"""
+Example: Zero-variance handling with the three objective modes.
+
+This script demonstrates how MultiFitter handles data points that have
+zero variance under each of the three objective strategies:
+
+  - legacy_mask : drops zero-variance points (old behaviour)
+  - hybrid      : keeps all points; applies Mighell substitution only to
+                   zero-variance entries (new default)
+  - mighell     : applies the Mighell transform to every point
+
+The example builds a simple film-on-substrate reflectometry model,
+creates synthetic data with a few zero-variance points injected, and
+fits with each mode so you can compare results side-by-side.
+"""
+
+import warnings
+
+import matplotlib.pyplot as plt
+import numpy as np
+import scipp as sc
+
+from easyreflectometry.calculators import CalculatorFactory
+from easyreflectometry.fitting import MultiFitter
+from easyreflectometry.model import Model
+from easyreflectometry.model import PercentageFwhm
+from easyreflectometry.sample import Layer
+from easyreflectometry.sample import Material
+from easyreflectometry.sample import Multilayer
+from easyreflectometry.sample import Sample
+
+
+def build_sample_and_model():
+    """Return a fresh (sample, model) pair with fittable parameters."""
+    si = Material(2.07, 0, 'Si')
+    sio2 = Material(3.47, 0, 'SiO2')
+    film = Material(2.0, 0, 'Film')
+    d2o = Material(6.36, 0, 'D2O')
+
+    si_layer = Layer(si, 0, 0, 'Si layer')
+    sio2_layer = Layer(sio2, 30, 3, 'SiO2 layer')
+    film_layer = Layer(film, 250, 3, 'Film Layer')
+    superphase = Layer(d2o, 0, 3, 'D2O Subphase')
+
+    sample = Sample(
+        Multilayer(si_layer),
+        Multilayer(sio2_layer),
+        Multilayer(film_layer),
+        Multilayer(superphase),
+        name='Film Structure',
+    )
+
+    resolution_function = PercentageFwhm(0.02)
+    model = Model(sample, 1, 1e-6, resolution_function, 'Film Model')
+
+    # Free the parameters we want to fit
+    sio2_layer.thickness.fixed = False
+    sio2_layer.thickness.bounds = (15, 50)
+
+    film_layer.thickness.fixed = False
+    film_layer.thickness.bounds = (200, 300)
+
+    film.sld.fixed = False
+    film.sld.bounds = (0.1, 3)
+
+    model.background.fixed = False
+    model.background.bounds = (1e-7, 1e-5)
+
+    model.scale.fixed = False
+    model.scale.bounds = (0.5, 1.5)
+
+    model.interface = CalculatorFactory()
+
+    return sample, model
+
+
+def make_synthetic_data_with_zero_variances(model):
+    """Generate a DataGroup from the model, then inject zero variances."""
+    qz = np.linspace(0.01, 0.30, 80)
+    r_true = model.interface.fit_func(qz, model.unique_name)
+
+    # Add a little Gaussian noise for realism
+    rng = np.random.default_rng(42)
+    noise_sigma = 0.02 * r_true + 1e-4
+    r_noisy = r_true + rng.normal(0, noise_sigma)
+
+    # Variance = sigma^2
+    variances = noise_sigma**2
+
+    # Inject 6 points with zero variance (simulating detector dead-spots, etc.)
+    zero_indices = [5, 15, 30, 45, 60, 75]
+    variances[zero_indices] = 0.0
+
+    data = sc.DataGroup(
+        {
+            'coords': {
+                'Qz_0': sc.array(dims=['Qz_0'], values=qz, unit=sc.Unit('1/angstrom')),
+            },
+            'data': {
+                'R_0': sc.array(dims=['Qz_0'], values=r_noisy, variances=variances),
+            },
+        }
+    )
+    return data, zero_indices
+
+
+def run_fit(objective, data):
+    """Build a fresh model and fit *data* using the given objective."""
+    _, model = build_sample_and_model()
+    fitter = MultiFitter(model, objective=objective)
+
+    with warnings.catch_warnings(record=True) as caught:
+        warnings.simplefilter('always')
+        result = fitter.fit(data)
+
+    return result, fitter, caught
+
+
+def plot_fit_comparison(data, fit_summaries, zero_indices):
+    """Plot the experiment and all fitted curves on a single matplotlib chart."""
+    qz = data['coords']['Qz_0'].values
+    reflectivity = data['data']['R_0'].values
+    variances = data['data']['R_0'].variances
+    yerr = np.sqrt(np.clip(variances, 0.0, None)) if variances is not None else None
+
+    fig, ax = plt.subplots(figsize=(10, 6))
+    ax.errorbar(
+        qz,
+        reflectivity,
+        yerr=yerr,
+        fmt='o',
+        ms=4,
+        color='black',
+        alpha=0.65,
+        capsize=2,
+        label='Experiment',
+    )
+    ax.scatter(
+        qz[zero_indices],
+        reflectivity[zero_indices],
+        s=80,
+        facecolors='none',
+        edgecolors='crimson',
+        linewidths=1.5,
+        label='Zero-variance points',
+        zorder=4,
+    )
+
+    for summary in fit_summaries:
+        ax.plot(
+            qz,
+            summary['model_curve'],
+            linewidth=2,
+            label=(
+                f"{summary['objective']} fit "
+                f"(obj={summary['objective_reduced_chi']:.3g}, class={summary['classical_reduced_chi']:.3g})"
+            ),
+        )
+
+    ax.set_title('Zero-Variance Objective Comparison')
+    ax.set_xlabel('Qz (1/angstrom)')
+    ax.set_ylabel('Reflectivity')
+    ax.set_yscale('log')
+    ax.grid(True, which='both', alpha=0.25)
+    ax.legend()
+    fig.tight_layout()
+    plt.show()
+
+
+def main():
+    # --- setup ---------------------------------------------------------------
+    _, seed_model = build_sample_and_model()
+    data, zero_idx = make_synthetic_data_with_zero_variances(seed_model)
+
+    total_points = len(data['data']['R_0'].values)
+    print(f'Dataset: {total_points} points, {len(zero_idx)} with zero variance\n')
+
+    # --- fit with each objective --------------------------------------------
+    objectives = ['legacy_mask', 'hybrid', 'mighell']
+    fit_summaries = []
+
+    for obj in objectives:
+        print(f'{"=" * 60}')
+        print(f'Objective: {obj}')
+        print(f'{"=" * 60}')
+
+        result, fitter, caught = run_fit(obj, data)
+
+        # Show warnings emitted during fitting
+        for w in caught:
+            print(f'  [WARNING] {w.message}')
+
+        print(f'  Success      : {result["success"]}')
+        print(f'  Objective reduced chi2 : {result["objective_reduced_chi"]:.6f}')
+        print(f'  Objective total chi2   : {fitter.objective_chi2:.6f}')
+        print(f'  Classical reduced chi2 : {result["classical_reduced_chi"]:.6f}')
+        print(f'  Classical total chi2   : {fitter.classical_chi2:.6f}')
+        print()
+
+        fit_summaries.append(
+            {
+                'objective': obj,
+                'model_curve': result['R_0_model'].values,
+                'objective_reduced_chi': float(result['objective_reduced_chi']),
+                'classical_reduced_chi': float(result['classical_reduced_chi']),
+            }
+        )
+
+    # --- interpretation note -------------------------------------------------
+    print('-' * 60)
+    print('NOTE: Under the mighell objective ALL points are transformed,')
+    print('so objective chi2 is not a classical chi-square statistic.')
+    print('Under hybrid, only zero-variance points are transformed;')
+    print('the classical chi2 is still computed from the original')
+    print('positive-variance points for comparison.')
+    print('-' * 60)
+
+    plot_fit_comparison(data, fit_summaries, zero_idx)
+
+
+if __name__ == '__main__':
+    main()