[estevesjh]

Problem

SelectionFunction.S_i and first_moment_i integrated λ_tr with one 64-node Gauss–Legendre rule over the whole bracket [max(0, μ−8σ), μ+8σ]. The bin kernel steps at every richness edge, centred where μ(λ_tr) = λ_e, over the kernel width (3 richness units for EmgParams(-1.5, 3.0, 0.3, 0.12)). Meanwhile the bracket of a 10¹⁵ h⁻¹M☉ halo spans ~900. So the rule's error oscillated as the nodes slid past the steps:

• LogNormalMor + that kernel, z = 0.4: S_i(M) of [60, 200) was off by up to 1.3e-2 of its peak (7.4e-4 at 128 nodes, 6.7e-7 at 256), a ripple of period ~0.3 in ln M. No ln M rule can average it out: it kept ClusterCounts at ~1e-3.
• HodMor's density jumps from zero at λ_cen = 1. One rule across a jump converges only as O(1/n): 5e-4 of the [20, 30) peak at 64 nodes, still 1.3e-4 at 256, with the fixed kernel.

Change

The bracket logic is unchanged. The bracket is now cut into panels, each with n_quad GL nodes (new default 20 per panel), at:

1. λ_tr_e ± 5 w_e around every edge's step. The new richness_edge_steps gives centre (λ_e − Δμ)/b_μ and width σ/b_μ, exact for the Y3 fit's b_μ ≠ 1. Each step sits mid-panel and hands its neighbours only a Φ(−5) remnant. Breaking at the step centres alone was measured and rejected: it leaves a half-step at the end of a panel hundreds of units long (2e-5 at 32 nodes per panel).
2. The relation's mean. P(λ_tr|M) is integrated in two halves, with its peak at a panel end.
3. HodMor.support_min, a new optional MOR hook (λ_cen). LogNormalMor needs none: its density vanishes continuously at zero.

Breaks outside the bracket leave empty panels with zero weight, and the integrand is evaluated only on live nodes. That keeps the kernel off λ_tr = 0, where the Y3 σ ∝ λ^{a_σ} vanishes, and fixes a ValueError for HodMor + Y3 below M_min, where the bracket is the single point zero. New public panel_edges shows the layout.

Measured

Against adaptive quadrature (quad_vec with breakpoints): 2 relations × {fixed, Y3} kernel × 3 binnings ([20,30,45,60,200), [20,30,45,1000), [19.5,20.5)) × z ∈ {0.2, 0.4, 0.6}, 201 masses in [10¹³, 10^15.5] h⁻¹M☉.

per-panel order	max error / bin peak	max relative error where S_i > 1e-5	live nodes / cell
16	3.1e-7	4.5e-6	116
20 (default)	3.4e-9	4.5e-8	145
24	2.1e-10	4.2e-9	174

The old single 64-node rule had up to 1.3e-2 / bin peak on the same sweep, and up to 55% relative error in the small-S_i tails.

• Counts: with the (separate, still uncommitted) GL ClusterCounts patch applied, convergence() on number_counts.py's example goes from 9.6e-4 to 1.3e-6 at the default 96-node ln M rule. It is unchanged at 40 nodes per panel, so the residual belongs to the counts rule.
• Cost: ~2.3× the old node count. Per 96-mass call, S_i goes 1.8 → 5.8 ms and first_moment_i 27 → 65–82 ms. The time is spent in the kernel itself, not in the panel bookkeeping.

Tests

13 new tests in tests/test_selection.py:

• exponential convergence in the panel order (three setups);
• adaptive-quadrature checks at the ripple masses and for HodMor.buzzard() + Y3;
• first_moment_i convergence;
• panel layout, the HOD jump as a panel edge, richness_edge_steps;
• the degenerate bracket and n_quad validation.

7 of the accuracy tests fail on the old rule. Full suite: 905 passed, 22 skipped, 33 xfailed (also passes with the GL-counts patch applied: 910 passed).

Not in this PR

• Inner λ_ob rule: richness_bin_first_moment's own fixed 32-node λ_ob rule is unconverged across bins much wider than the kernel: 2.2% in [60, 200) at σ = 3, worse in [45, 1000); Y3 is fine in the DES bins. That bounds first_moment_i / mean_richness there whatever the λ_tr rule does.
• Docs: docs/selection_function.md still shows the old 64-node output block (S_[60,200)(10¹⁵, z=0.3) is now 0.1533, not 0.1618) and says "one Gauss–Legendre rule".
• Changed meaning of n_quad: it now counts nodes per panel. Notes scripts passing n_quad=96 still work but are slower than needed; dropping the argument gives the default.

🤖 Generated with Claude Code