EData, Data Points, and Target Effects
=======================================

This page describes how experimental data is managed internally and how
target integration, beam convolution, and straggling effects are applied.
This is essential reading for anyone working on improvements to the
convolution routines.

Data Structure Overview
-----------------------

The data hierarchy from top to bottom:

``EData``
   Top-level container. Holds all segments, target effects, and manages the
   overall calculation loop. One ``EData`` object exists per calculation; it is
   cloned for thread-safe parallel fitting.

``ESegment``
   A calculation segment corresponding to a specific entrance/exit particle
   pair, energy range, angle range, and data type. Contains a vector of
   ``EPoint`` objects.

``EPoint``
   A single data point with energy, angle, cross section, and uncertainty.
   When target effects are active, an ``EPoint`` also holds a vector of
   **sub-points** (``integrationPoints_``) used for numerical integration.

``TargetEffect``
   Configuration for experimental effect corrections: Gaussian convolution
   parameters, stopping power equation, density, Q-coefficients, straggling,
   and the number of integration sub-points.

**Key files:**

- ``include/EData.h``, ``src/EData.cpp``
- ``include/ESegment.h``, ``src/ESegment.cpp``
- ``include/EPoint.h``, ``src/EPoint.cpp``
- ``include/TargetEffect.h``, ``src/TargetEffect.cpp``
- ``include/AdaptiveIntegrationGrid.h``, ``src/AdaptiveIntegrationGrid.cpp``

Data Loading and Initialization
-------------------------------

EData::Fill()
^^^^^^^^^^^^^

The ``EData::Fill()`` method (in ``src/EData.cpp``) orchestrates the loading
and initialization of all data:

1. **Read data files** -- for each data segment, parse the four-column data
   file and create ``EPoint`` objects for data points that fall within the
   segment's energy and angle ranges.

2. **Convert to CM frame** -- laboratory energies and angles are converted to
   center-of-mass frame values using the kinematics of the entrance/exit pair.

3. **Map duplicate points** -- if multiple segments reference the same
   physical point (same energy, angle, and reaction), the point is calculated
   only once and other occurrences are mapped to it via ``EnergyMap``.

4. **Create sub-points** -- for each data point that has an associated target
   effect, a grid of sub-points is generated for numerical integration
   (see below).

5. **Initialize sub-points** -- each sub-point's energy-dependent quantities
   (penetrabilities, Coulomb phases, shift functions) are pre-computed via
   ``EPoint::Initialize()``.

Sub-Point Creation
------------------

When a data point has a target effect, the code creates a grid of sub-points
spanning the energy range needed for integration. The sub-points are stored
inside the parent ``EPoint`` in the ``integrationPoints_`` vector.

Energy Range Determination
^^^^^^^^^^^^^^^^^^^^^^^^^^

The sub-point energy range depends on which effects are active:

**Target integration only:**

- From the beam energy (surface) down by the target thickness:

  .. math::

     E_\text{back} = E_\text{CM} - \Delta E_\text{target}

  where :math:`\Delta E_\text{target} = \epsilon(E) \cdot \rho` is the target
  thickness computed from the stopping power :math:`\epsilon(E)` and the
  areal density :math:`\rho`.
- Forward depth is zero (no energy above beam energy).

**Gaussian convolution only:**

- Symmetric range around the beam energy:

  .. math::

     E_\text{range} = E_\text{CM} \pm n_\sigma \cdot \sigma_b

  where :math:`\sigma_b` is the beam energy resolution (Gaussian sigma) and
  :math:`n_\sigma` is the ``convolutionRange`` parameter (typically 5).

**Target integration + convolution:**

- Backward: target thickness plus the convolution tail:

  .. math::

     E_\text{back} = E_\text{CM} - \Delta E_\text{target} - n_\sigma \cdot \sigma_b

- Forward: convolution tail above the beam energy:

  .. math::

     E_\text{forward} = E_\text{CM} + n_\sigma \cdot \sigma_b

**Straggling:**

When straggling is enabled, the backward range is further extended by the
straggling width at the back of the target:

.. math::

   \sigma_\text{straggling} = c_s \sqrt{\Delta E_\text{target}}

where :math:`c_s` is the straggling coefficient (in keV units).

Adaptive Grid Generation
^^^^^^^^^^^^^^^^^^^^^^^^^

Sub-points are placed on an **adaptive energy grid** generated by the
``AdaptiveIntegrationGrid`` class. The grid is denser near resonances and
coarser in smooth regions:

1. **Resonance detection** -- the generator scans all compound nucleus levels
   and identifies resonances that fall within the integration range. For each
   resonance, the total width :math:`\Gamma` is estimated from the reduced
   width amplitudes.

2. **Step size calculation** -- the local step size varies smoothly using a
   Gaussian falloff:

   .. math::

      \Delta E = \Delta E_\text{fine} \cdot e^{-d^2/2\sigma_r^2}
      + \Delta E_\text{base} \cdot (1 - e^{-d^2/2\sigma_r^2})

   where :math:`d` is the distance to the nearest resonance,
   :math:`\Delta E_\text{fine} = \Gamma / N_\text{ppw}` (points per width),
   and :math:`\sigma_r = \Gamma \cdot M / 2` (with :math:`M` being a multiplier).
   This ensures narrow resonances get a fine grid without discontinuities.

3. **Grid construction** -- starting from the highest energy, points are placed
   at adaptively determined intervals down to the lowest energy. The total
   number of points is bounded by ``maxPoints`` (the user-specified integration
   points count).

Each sub-point is an ``EPoint`` object with:

- The same angle as the parent point.
- Energy set to the grid point energy.
- Stopping power pre-computed at that energy (for target integration).

Calculation Flow for Target Effects
------------------------------------

During calculation (in ``AZURECalc``), the processing of a point with target
effects follows this order:

1. **Calculate each sub-point** -- the R-matrix cross section is computed at
   each sub-point energy by calling ``EPoint::Calculate()`` on each sub-point.
   This runs the full FillMatrices/InvertMatrices/CalculateTMatrix/
   CalculateCrossSection pipeline for each sub-point.

2. **Integrate** -- ``EPoint::IntegrateTargetEffect()`` is called on the parent
   point to numerically integrate the sub-point cross sections, yielding the
   effective yield that accounts for target and beam effects.

Integration Methods
-------------------

The integration is implemented in ``EPoint::IntegrateTargetEffect()``
(``src/EPoint.cpp``). All cases use **2-point Gauss-Legendre quadrature** on
each sub-interval, which is exact for cubic polynomials and handles
non-uniform (adaptive) grids correctly.

Cross section values between grid points are obtained by **linear
interpolation** from the two nearest sub-points.

Convolution Only
^^^^^^^^^^^^^^^^

For pure Gaussian convolution (no target integration), the yield at beam
energy :math:`E_0` is:

.. math::

   Y(E_0) = \int \sigma(E') \cdot g(E' - E_0) \, dE'

where :math:`g` is the Gaussian beam profile:

.. math::

   g(E' - E_0) = \frac{1}{\sqrt{2\pi}\,\sigma_b}
   \exp\!\left(-\frac{(E' - E_0)^2}{2\sigma_b^2}\right)

**Implementation:** The code loops over consecutive sub-point pairs
:math:`[E_i, E_{i+1}]`. For each interval, two Gauss-Legendre points are
evaluated. At each Gauss point, the cross section is linearly interpolated
from the bracketing sub-points, multiplied by the Gaussian convolution factor
(``TargetEffect::GetConvolutionFactor()``), and accumulated into the integral.

Target Integration Only
^^^^^^^^^^^^^^^^^^^^^^^^

For target integration without convolution, the yield is:

.. math::

   Y(E_0) = \int_{E_0 - \Delta}^{E_0} \frac{\sigma(E')}{\epsilon(E')} \, dE'

where :math:`\epsilon(E')` is the stopping cross section. The integral runs
from the beam energy at the target surface down to the energy at the back of
the target.

**Implementation:** Similar Gauss-Legendre quadrature over sub-intervals,
but the integrand is :math:`\sigma / \epsilon` (cross section divided by
stopping power). Only intervals within the target range
:math:`[E_\text{back}, E_\text{surface}]` contribute.

Target Integration + Convolution
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The full nested integral is:

.. math::

   Y(E_0) = \frac{1}{\rho} \int_{E_\text{back}}^{E_\text{surface}}
   \left[\int \frac{\sigma(E')}{\epsilon(E')} \cdot g(E' - E_d) \, dE'\right] dE_d

This is a **double integral**: the outer integral runs over target depth
(energy loss), and the inner integral convolves the cross section with the
beam energy distribution at each depth.

**Implementation:** Two nested Gauss-Legendre loops:

1. **Outer loop** (target depth): iterates over sub-point intervals within
   the target range :math:`[E_\text{back}, E_\text{surface}]`. At each Gauss
   point :math:`E_d`, the effective beam sigma is computed (including
   straggling if enabled).

2. **Inner loop** (convolution): for each depth :math:`E_d`, integrates
   :math:`\sigma(E')/\epsilon(E') \cdot g(E' - E_d)` over the range
   :math:`E_d \pm n_\sigma \cdot \sigma_\text{eff}`, where
   :math:`\sigma_\text{eff}` includes both beam resolution and straggling:

   .. math::

      \sigma_\text{eff} = \sqrt{\sigma_b^2 + \sigma_\text{straggling}^2}

   The straggling contribution grows with depth:

   .. math::

      \sigma_\text{straggling} = c_s \sqrt{E_\text{surface} - E_d}

Target Integration + Straggling (no beam convolution)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

When straggling is enabled but Gaussian beam convolution is not, the code
uses the same nested double integral structure, but the convolution kernel
is purely from straggling (no beam sigma component). At the target surface
(:math:`\sigma_\text{straggling} = 0`), the inner integral reduces to a
direct :math:`\sigma/\epsilon` evaluation (no spreading).

Key Helper Functions
--------------------

``TargetEffect::GetConvolutionFactor(energy, centroid)``
   Returns the Gaussian weight :math:`g(E - E_0)` for a fixed sigma.

``TargetEffect::CalculateConvolutionFactor(energy, centroid, config)``
   Returns the Gaussian weight with an energy-dependent sigma (evaluated from
   a user-provided equation).

``TargetEffect::TargetThickness(energy, config)``
   Returns :math:`\epsilon(E) \cdot \rho` -- the product of stopping power
   and density.

``TargetEffect::GetStoppingPowerEq()``
   Returns the ``Equation`` object for the parametrized stopping cross section
   (user-defined functional form with parameters ``a0``, ``a1``, ...).

``AdaptiveIntegrationGrid::GenerateGrid(startEnergy, endEnergy, compound)``
   Generates an adaptive energy grid with finer spacing near resonances.
   Returns a ``std::vector<double>`` of energies from high to low.

Summary: Data Flow Diagram
--------------------------

::

   Data File (lab frame)
       │
       ▼
   EData::Fill()
       ├── Create EPoint objects (convert to CM frame)
       ├── Map duplicate points across segments
       └── For target-effect points:
           ├── Compute target thickness from stopping power
           ├── Determine integration range (target + convolution + straggling)
           ├── Generate adaptive energy grid (AdaptiveIntegrationGrid)
           └── Create sub-point EPoints at grid energies
                   │
                   ▼
   AZURECalc calculation loop (per data point):
       ├── For each sub-point:
       │       └── EPoint::Calculate()
       │           ├── ClearMatrices()
       │           ├── FillMatrices()      ← build A⁻¹ or R matrix
       │           ├── InvertMatrices()    ← GSL LU decomposition
       │           ├── CalculateTMatrix()  ← collision matrix
       │           └── CalculateCrossSection() ← σ(E) at sub-point
       │
       └── EPoint::IntegrateTargetEffect()
           └── Gauss-Legendre quadrature over sub-point grid
               ├── Convolution only: ∫ σ(E') g(E'-E₀) dE'
               ├── Target only: ∫ σ(E')/ε(E') dE'
               └── Both: ∫∫ σ(E')/ε(E') g(E'-E_d) dE' dE_d
                           ▼
                   Final yield stored in parent EPoint