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:

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.
Convert to CM frame – laboratory energies and angles are converted to center-of-mass frame values using the kinematics of the entrance/exit pair.
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.
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).
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:

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

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

Gaussian convolution only:

Symmetric range around the beam energy:

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

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

Target integration + convolution:

Backward: target thickness plus the convolution tail:

\[E_\text{back} = E_\text{CM} - \Delta E_\text{target} - n_\sigma \cdot \sigma_b\]
Forward: convolution tail above the beam energy:

\[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:

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

where \(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:

Resonance detection – the generator scans all compound nucleus levels and identifies resonances that fall within the integration range. For each resonance, the total width \(\Gamma\) is estimated from the reduced width amplitudes.
Step size calculation – the local step size varies smoothly using a Gaussian falloff:

\[\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 \(d\) is the distance to the nearest resonance, \(\Delta E_\text{fine} = \Gamma / N_\text{ppw}\) (points per width), and \(\sigma_r = \Gamma \cdot M / 2\) (with \(M\) being a multiplier). This ensures narrow resonances get a fine grid without discontinuities.
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:

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.
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 \(E_0\) is:

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

where \(g\) is the Gaussian beam profile:

\[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 \([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:

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

where \(\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 \(\sigma / \epsilon\) (cross section divided by stopping power). Only intervals within the target range \([E_\text{back}, E_\text{surface}]\) contribute.

Target Integration + Convolution

The full nested integral is:

\[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:

Outer loop (target depth): iterates over sub-point intervals within the target range \([E_\text{back}, E_\text{surface}]\). At each Gauss point \(E_d\), the effective beam sigma is computed (including straggling if enabled).
Inner loop (convolution): for each depth \(E_d\), integrates \(\sigma(E')/\epsilon(E') \cdot g(E' - E_d)\) over the range \(E_d \pm n_\sigma \cdot \sigma_\text{eff}\), where \(\sigma_\text{eff}\) includes both beam resolution and straggling:

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

The straggling contribution grows with depth:

\[\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 (\(\sigma_\text{straggling} = 0\)), the inner integral reduces to a direct \(\sigma/\epsilon\) evaluation (no spreading).

Key Helper Functions

TargetEffect::GetConvolutionFactor(energy, centroid): Returns the Gaussian weight \(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 \(\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