# LARBED Extraction (`pyextal.larbed`)

`pyextal.larbed` is the stage that runs *before* refinement. It turns a raw
EMPAD 4D-STEM acquisition into one 2D rocking curve per reflection — the arrays
that [`LARBEDDiffractionInfo`](ROI.md) and the rest of the pipeline consume.

It was merged in from the standalone
[PyLARBED](https://github.com/rtbusch/PyLARBED) package, so the class and method
names are unchanged; only the import path is new:

```python
from pyextal.larbed import LARBEDAnalysis, LARBEDCalibration
```

---

## Where it sits in the pipeline

```
scan_x256_y256.raw   (EMPAD, 4D: scan_x × scan_y × det_x × det_y)
      │
      ▼
┌──────────────────────────────────────────────┐
│  Stage 0 — Extraction  (pyextal.larbed)      │
│                                              │
│  LARBEDAnalysis                              │
│    → load_data / align_data                  │
│    → assign_gvector + calculate_grid_vectors │
│    → calculate_g_vectors                     │
│    → IntegrateLarbed                         │
│    → to_dinfo(...)  or  save_h5(...)         │
└───────────────┬──────────────────────────────┘
                │  LARBEDDiffractionInfo  /  {sample}.h5  ──▶  /raw
                ▼
        Stage 1 — Pre-refinement   (LARBEDROI, CoarseOptimize, FineOptimize)
                ▼
        Stage 2 — JAX refinement   (jaxEngine)
```

---

## Typical workflow

```python
from pyextal.larbed import LARBEDAnalysis

larbed = LARBEDAnalysis("scan_x256_y256.raw")
larbed.load_data(type=0)
larbed.average_data()
# larbed.align_data()          # optional; only if the direct beam drifts

# 1. Index the reflections. Two basis vectors for a zone-axis pattern …
larbed.assign_gvector((46, 48), (64, 36))          # detector pixel positions
larbed.calculate_grid_vectors(crop=2, order=3, Adjust=True)
larbed.calculate_g_vectors((1, 3, -2), (4, 2, 0))  # matching Miller indices
larbed.plot_grid_vectors(ratio=0.005)              # check the indexing

# 2. Integrate each reflection into a rocking curve
larbed.IntegrateLarbed(ri=3, ro=(6, 8), g=580, m=0.47, Variance=True)
larbed.plot_larbed(ratio=0.5)
```

`assign_gvector` takes **detector pixel positions**; `calculate_g_vectors` takes
the **Miller indices** of those same reflections, in the same order. Supply two
of each for a zone-axis grid, or four to include HOLZ reflections:

```python
larbed.assign_gvector((59.7, 52.3), (52.1, 63.2), (53, 32.2), (62, 62))
larbed.calculate_grid_vectors(order=(8, 8, 0, 0), Adjust=False)
larbed.calculate_g_vectors((-1, 0, 0), (0, 1, 0), (0, 0, 1), (0, 0, 0))
```

`calculate_grid_vectors` dispatches on how many basis vectors it was given: four
selects the HOLZ grid (and `order` must then be a 4-tuple), two selects either
the peak-snapped grid (`Adjust=True`) or the plain lattice grid.

### Aperture and variance

`IntegrateLarbed(ri, ro, g, m, DQE, varB, Variance)` sums each reflection over a
square aperture of half-width `ri`. When `ro=(inner, outer)` is given, the mean
of that annulus is subtracted as background.

With `Variance=True` the counting variance is propagated alongside. The per-pixel
model matches {func}`pyextal.LucyRichardson.sgnoise`: a signal-dependent term
`g*m*I/DQE` plus a constant readout term `g*m*varB`, summed over the aperture.
`varB` defaults to 0.

---

## Handing off to refinement

Two routes, both replacing the older `.npy` round-trip:

```python
# (a) straight into the refinement objects
dinfo = larbed.to_dinfo("si110.dat", thickness=1000.0, gl=37.0)
roi = LARBEDROI(dinfo=dinfo, rotation=214, gx=..., gInclude=[...])

# (b) into the same HDF5 file the later stages append to
larbed.save_h5("sample.h5", prefix="YIG_tilt12", crystal_dat="yig.dat")
```

`save_h5` writes `/raw/dp`, `/raw/variance`, and `/raw/g_vectors` — exactly what
{func}`pyextal.hdf5io.save_raw` defines — so extraction, pre-refinement, and
every JAX refinement run share one file per dataset.

`gl` is the tilt step size, expressed as the number of scan steps spanning the
x-axis reflection. Deriving it is covered in
[Understand how to define/index LARBED dataset](LARBEDindex.md);
{meth}`~pyextal.larbed.LARBEDCalibration.calculate_average_step_length` measures
the pixels-per-step half of it.

{meth}`~pyextal.larbed.LARBEDAnalysis.stack` returns the same three arrays
directly if you would rather handle them yourself, and
{meth}`~pyextal.larbed.LARBEDAnalysis.save_larbed` still writes the loose `.npy`
files.

---

## MTF calibration and deconvolution

```{note}
Legacy. Because `IntegrateLarbed` already integrates over an aperture,
deconvolving against the MTF mainly adds noise, and `deconv_Larbedstack`
discards the propagated variance. It is normally skipped.
```

`LARBEDCalibration` measures the detector point-spread function from a reference
scan taken with no sample and descan off:

```python
calib = LARBEDCalibration("probe_scan.raw")
calib.load_data(type=0)
calib.find_peaks(start=100, end=115)
calib.calculate_average_step_length()
calib.crop_peaks()
calib.interpolate()
calib.calculate_mtf()

larbed.deconv_Larbedstack(mtf_2d=calib.mtf_2d, varB=787, g=580, niter=10)
```

### Two convolution conventions

`calculate_mtf` returns a **Fourier-domain** MTF, and
{func}`pyextal.LucyRichardson.lucy_Richardson` accepts a `mode` argument
selecting how it is applied:

| `mode` | Operator | Use |
| :--- | :--- | :--- |
| `'fourier'` (default) | `ifft2(fft2(img) * mtf)` | Consistent with a Fourier-domain MTF, and what the CBED path has always used. |
| `'real'` | `scipy.signal.convolve(img, mtf, mode='same')` | What the standalone PyLARBED used. Kept so earlier results reproduce exactly. |

The `LARBEDAnalysis` deconvolution methods pass `mode='real'`, so they reproduce
pre-merge PyLARBED output bit for bit. The two operators are genuinely
different — a Fourier-domain array fed to the real-space convolution is not the
same deconvolution — so switching an existing analysis to `'fourier'` will change
its numbers.

---

## Worked example

`example/LARBED/STO_001/STO_001_LARBED.ipynb` runs this whole page against a real
SrTiO₃ [001] dataset and then continues through pre-refinement, using the
indexing parameters validated in the original PyLARBED notebook. It ships with
its crystal file, `STO_001.dat`.

Note the step it takes between extraction and refinement. `STO_001.dat` leaves
both of its beam lists implicit: `sele` picks the matrix beams geometrically at
parse time and `aper` picks the output reflections by g-length, so both move with
the voltage, the tilt, and `xaxis`. The notebook rewrites it as
`STO_001_explicit.dat` with an explicit `hkl` block (all 185 beams) and an
explicit `out` block set to exactly `gInclude`, and refines against that — so the
calculation is reproducible beam for beam.

The output beams are written as a **prefix of the `hkl` list**, so
`hkl[:nout] == hklout` and the engine's `idptr` mapping is the identity. That is
worth preserving: `bloch_parse` resolves each `out` entry to its index in `hkl`,
and any entry it cannot find is *appended* to `hkl` — past the array already
allocated from the declared `nbeams`.

```{warning}
That notebook's `gl` (tilt step size) is a placeholder and its Debye-Waller
factors are literature starting values. Set them for your data before trusting a
refined number.
```

---

## Module layout

| Module | Contents |
| :--- | :--- |
| `pyextal.larbed.analysis` | `LARBEDAnalysis` — the extraction workflow |
| `pyextal.larbed.calibration` | `LARBEDCalibration` — MTF and scan-step calibration |
| `pyextal.larbed.empad` | Raw I/O, virtual apertures, reciprocal-grid construction |
| `pyextal.larbed.peaks` | Least-squares peak models and fitting |

Throughout `pyextal.larbed.empad`, `x` is the **first** (slow) detector axis and
`y` the second — not the screen-space x/y of a plotted image. The grid builders
take a `shape` argument used to reject off-detector nodes; it defaults to
`(124, 124)`, the cropped EMPAD frame.