Add photon multigroup cross section generation

[paulromano]

Description

In anticipation of supporting photon random-ray workflows, this PR makes two important updates to multigroup cross section generation: 1) it makes multigroup cross section libraries particle-aware through a particle_type argument, and 2) it adds a PhotonProductionMatrixMGXS class for representing the photon production cross section needed for deterministic photon transport.

For photon libraries, the existing "nu-scatter matrix" quantity represents the expected outgoing photon production matrix. It combines the coherent or incoherent surviving primary photon with secondary photons from atomic relaxation, thick-target bremsstrahlung, and positron annihilation. Photon yield is folded directly into the matrix because a multiplicity defined by dividing production by physical scattering is undefined for reactions that remove the incident photon while producing secondaries.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[GuySten]

How is this feature used with delayed gamma scaling?
What about keff scaling for consistent energy balance?

By using only one PhotonProductionMatrixMGXS we cannot use keff scaling and we need to keep track of if delayed gammas were used in the library.

We might want to separate PhotonProductionMatrixMGXS into PhotonProductionMatrixMGXS, PromptFissionPhotonProductionMatrixMGXS and DelayedFissionPhotonProductionMatrixMGXS to support the two use cases.

Alternatively, if we only care about generating weight windows the consistent energy balance and delayed gamma scaling should not be that important.

[shimwell]

Linking to this issue #3489

[paulromano]

@GuySten Thanks for raising this. This PR is intentionally limited to photon-only transport and generates the photon-to-photon transfer matrix, including secondary photons from normal photon interactions, TTB, and atomic relaxation. It's not intended to handle coupled neutron-photon transport. To make that distinction clearer, I renamed the class from PhotonProductionMatrixXS to PhotonTransferMatrixXS. Delayed-fission scaling and k-eff normalization will need to be handled by a future coupled neutron–photon source treatment.

[jtramm]

This is an exciting feature -- it's been a long term goal to get neutron-photon coupling working with the random ray solver to produce sets of weight windows that work for the coupled MC workflow seamlessly. Great to see progress in this direction, and excited to see the follow up PRs on this line!

I'm still looking through manually, but also did have claude fable take a look and it reported the issues below for your consideration:

Claude Summary

Four major items found:

1. particle_type='neutron' breaks the default neutron library, because
   several MGXS classes index tally filters by position and the branch appends
   a ParticleFilter at the end.
2. The exported photon library is silently wrong in multigroup Monte Carlo
   mode, which is what create_mg_mode builds; a verified fix exists.
3. survival_biasing = True corrupts the production matrix through a
   pre-existing C++ scoring bug, and the Python side hides the NaN.
4. get_slice returns groups in reversed order; pre-existing, but the branch's
   global Tally.get_slice change is what makes matrix slicing run at all, and
   the new unit tests check shapes only.

Plus a list of minor items.

Item 1: particle_type='neutron' breaks positional filter indexing

Symptom. A pin-cell neutron library built with particle_type='neutron'
and otherwise default settings (correction='P0', legendre_order=0) fails
for transport, nu-transport, diffusion-coefficient,
nu-diffusion-coefficient (AttributeError: 'LegendreFilter' object has no attribute 'values') and for consistent scatter matrix,
consistent nu-scatter matrix (Unable to get the bin index ... "P0").
Because the default configuration requires transport, this takes down
create_mg_library, get_condensed_library, build_hdf5_store, get_slice,
get_pandas_dataframe and print_xs for the whole library. The same
particle_type=None library gives identical numbers for every type, so the
particle-aware path is the only difference.

Cause. MGXS.tallies appends ParticleFilter(self.particle_type) after
the class-specific filters (openmc/mgxs/mgxs.py:635). Several classes then
address filters by position:

• TransportXS.rxn_rate_tally and xs_tally: p1_tally.filters[-2]
  (mgxs.py:2813, 2838), expecting the EnergyoutFilter before the
  LegendreFilter.
• DiffusionCoefficient: same pattern at mgxs.py:3027, 3056, 3106.
• ScatterMatrixXS P0-correction paths: filters[-1] = legendre at
  mgxs.py:4290, 4291 and 4394. The simple-formulation case at 4290
  overwrites the ParticleFilter with a one-bin LegendreFilter; the numbers
  come out right only because both filters have one bin, and the derived
  tally's filter list is left corrupted.

Fix. Replace the positional access with find_filter(...) lookups, which
is what the rest of the module does. Inserting the ParticleFilter right
after the domain filter in MGXS.tallies would also unbreak these sites, but
leaves the fragility in place. Either way, add a test that builds a
particle_type='neutron' library in the default configuration and calls
create_mg_library; no current test does.

Item 2: exported photon library is wrong for multigroup Monte Carlo

Symptom. Running the exported library in energy_mode='multi-group'
Monte Carlo on the same infinite-medium problem gives fluxes far below the
continuous-energy answer:

group (high → low energy)	1	2	3	4	5	6	7
rowsum(P) / Σt	1.78	1.54	1.44	1.33	1.19	0.48	0.07
MGMC / CE flux, library as exported	0.88	0.32	0.29	0.25	0.19	0.17	0.18
ν(g) = rowsum(P) / (Σt − Σa)	2.41	2.60	2.92	4.50	11.6	3.72	2.98
MGMC / CE flux, with ν(g) exported	1.00	1.01	0.92	0.98	1.09	1.03	1.00

The residual scatter in the last row is tally noise in the low-flux groups.

Cause. Multigroup Monte Carlo decides survival with absorption / total,
then samples the outgoing group from the stored scatter matrix and multiplies
the weight by the multiplicity, which defaults to 1 when no
multiplicity_matrix is present (src/xsdata.cpp:548). The branch stores the
full production matrix as scatter_matrix, keeps absorption as
photoelectric plus pair production, and writes no multiplicity. Expected
production per collision is therefore (Σt − Σa) / Σt instead of
rowsum(P) / Σt: the outgoing spectrum is right, the magnitude is low by the
factor ν(g). The neutron path handles the analogous "nu-scatter without
multiplicity" case by resetting absorption = total − rowsum at
openmc/mgxs/library.py:1286; the photon branch bypasses that block with
pass at library.py:1274. That adjustment could not work for photons anyway,
because with thick-target bremsstrahlung the production row sum exceeds the
total in most groups (first row of the table), so the adjusted absorption would
be negative. Nothing in the C++ reader checks consistency between total,
absorption and the scatter matrix, so the error is silent, and nothing reads
the new particle_type attribute.

Fix (verified). The PR text says a multiplicity is undefined for photons.
That holds element-wise, but the group-wise value
ν(g) = rowsum(P(g)) / (Σt(g) − Σa(g)) is always defined because
Σt − Σa = Σcoherent + Σincoherent > 0. Exporting it as a row-constant
multiplicity_matrix (via XSdata.set_multiplicity_matrix) makes MGMC
consistent (last table row) and leaves random ray untouched, since
NU_SCATTER returns the stored matrix regardless of multiplicity
(src/scattdata.cpp, get_xs). It also makes MGMC's scatter reaction rate
equal the physical scatter, consistent with absorption. If that is out of
scope for this PR, create_mg_mode and check_library_for_openmc_mgxs should
refuse photon libraries, and the docstrings should say the library is valid
only for solvers that use total plus nu-scatter (random ray).

Item 3: survival_biasing=True silently corrupts the production matrix

Symptom. Same balance problem with settings.survival_biasing = True: the
primary production tally sums to NaN, the reported production matrix loses its
entire primary block (top-group rowsum(P)/Σt drops from 1.78 to 0.56), and no
error is raised. The only trace is a NumPy RuntimeWarning from tally
arithmetic.

Cause. score_general_ce_analog computes the survival-biasing correction
wgt_absorb = wgt_last × neutron_xs(event_nuclide).absorption / neutron_xs(event_nuclide).total for every particle type
(src/tallies/tally_scoring.cpp:1124) and subtracts it from the analog
scatter score (line 1172) and uses it for analog absorption (line 1207).
Photons never receive survival biasing, and in a photon-only run the neutron
cross-section cache is never filled, so the term is 0/0. In a coupled
neutron-photon run the cache holds the parent neutron's values, so by the same
code the score would be finite but wrong; that case is inferred from the code,
not measured. The Python side then hides the damage: hybrid_product maps
NaN to zero (openmc/tallies.py:2338) and _compute_xs does it again
(openmc/mgxs/mgxs.py:1001).

This bug predates the branch, but the branch's headline feature depends on the
analog photon scatter score, and survival biasing is a routine setting in the
fixed-source shielding problems where photon transport is used.

Fix. Guard the correction with settings::survival_biasing && p.type().is_neutron() in score_general_ce_analog. Independently, the MGXS
code should fail loudly if a raw production tally contains NaN instead of
zeroing it.

Item 4: get_slice returns groups in reversed order

Symptom. With real statepoint data, TotalXS.get_slice(groups=[1, 2]).get_xs()
returns [4.81, 0.855] while the full cross section is [0.855, 4.81, 81.9];
PhotonTransferMatrixXS.get_slice(in_groups=[1, 2], out_groups=[2, 3]).get_xs()
returns the expected 2×2 block flipped on both axes.

Cause. Tally.get_slice extracts data in the requested bin order
(openmc/tallies.py:3319) but sorts the filter's bins with np.unique
(tallies.py:3374), so the sliced tally's data and filter disagree whenever
the request is not ascending. MGXS.get_slice and MatrixMGXS.get_slice
always request groups in group order, which is descending energy
(mgxs.py:1611, 2445). Pre-existing for vector MGXS. For matrix MGXS,
develop never got this far: the old isinstance replacement loop in
Tally.get_slice also replaced the EnergyoutFilter (a subclass of
EnergyFilter) with the sliced incoming filter, so ScatterMatrixXS.get_slice
crashed and PhotonTransferMatrixXS.get_slice would too. The branch's one-line
change to type(test_filter) is filter_type (tallies.py:3386) is correct and
consistent with find_filter, and it is required for the new class, but it is
a global behaviour change that the PR does not mention, and it now exposes the
ordering bug.

Related, pre-existing: slicing in_groups without out_groups shrinks
energy_groups to the incoming groups while the data keep all outgoing groups,
so get_xs reshapes wrongly (get_slice(in_groups=[1]) on a 3-group matrix
returns shape (3, 1, 1)).

Fix. Keep the requested order in Tally.get_slice (drop the np.unique,
or sort data and bins together), or make the MGXS slicers request ascending
bins. Strengthen tests/unit_tests/test_photon_mgxs.py:82 to assert values,
not shapes. Note that synthetic-data value checks need one filter object per
tally: MGXS.get_condensed_xs mutates the energy filter in place, and the
tallies built by the filters property share one EnergyFilter, so condensing
un-loaded synthetic tallies gives wrong numbers. Statepoint-loaded tallies have
separate filters and condense correctly.

Minor issues

1. Mixed flux estimators. Library.estimator='collision' gives total
   and absorption a collision-estimated flux while the production matrix
   keeps its tracklength flux (library.py:551 excludes the photon class).
   Only check_library_for_openmc_mgxs rejects it (library.py:1650), and
   only if the user calls create_mg_library. Either reject in
   build_library or let the production's flux tally follow the library
   estimator; only the primary and secondary tallies must be analog.
2. Silently ignored settings. The photon branch of
   check_library_for_openmc_mgxs (library.py:1631) validates num_polar
   and num_azimuthal but not legendre_order, scatter_format or
   histogram_bins. A photon library with legendre_order=3 or
   scatter_format='histogram' builds, passes the check and exports P0
   Legendre data without a message.
3. mgxs_types='all' is resolved at assignment (library.py:215).
   Library(geometry, mgxs_types='all') followed by
   library.particle_type = 'photon' keeps all 607 neutron types and fails
   much later with a long "unsupported" list. Resolve 'all' lazily or
   re-resolve it in the particle_type setter. The setter also resets
   correction to None without saying so in the docstring.
4. Delayed-group MGXS never receive particle_type. build_library
   passes it to MGXS.get_mgxs but not to MDGXS.get_mgxs
   (library.py:538), so a particle-aware library has an inconsistent tally
   set and those tallies cannot merge with the others.
5. PhotonTransferMatrixXS.estimator is a read-only property
   (mgxs.py:5031). Assignment raises AttributeError rather than the
   ValueError every other MGXS raises, build_library has to special-case
   the class, and _valid_estimators = ['analog'] at mgxs.py:4988 is dead.
6. Inaccurate docstring. mgxs.py:4975 says per-nuclide production cannot
   be determined because secondaries come from the bank. Analog scoring
   attributes the events score to p.event_nuclide()
   (src/tallies/tally_scoring.cpp:2338), so per-element production is
   tallyable; the by_nuclide=False restriction is a design choice and should
   be described as one.
7. Fidelity caveat, not a bug. correction=None is forced, so the photon
   library is isotropic P0 with an uncorrected total. Compton scattering is
   strongly forward-peaked at MeV energies, so deep-penetration results will be
   poor until a transport correction is available. The P1 moment of the primary
   scatter score is already tallyable for photons (p.mu() is set for
   coherent and incoherent events), so a transport-corrected total could be
   added later without new C++.