Neutron generator updates#1878
Conversation
|
Hi @michaelmackenzie,
which require these tests: build. @Mu2e/write, @Mu2e/fnalbuild-users have access to CI actions on main. ⌛ The following tests have been triggered for fd5c2a7: build (Build queue - API unavailable) |
Core changeNew helper in the base class centralizes what double calculateBinnedSpectrumEnergyFraction(fhicl::ParameterSet pset,
std::string var_low = "elow",
std::string var_high = "ehi",
double full_var_low = 0.,
double full_var_high = 0.) const {
BinnedSpectrum spectrum(pset);
pset.erase(var_low);
pset.erase(var_high);
pset.put(var_low, full_var_low);
pset.put(var_high, full_var_high);
BinnedSpectrum fullSpectrum(pset);
...
const double fraction = (fullIntegral > 0.) ? integral / fullIntegral : 0.;
return fraction;
}
Issues / concerns1. Behavior change in DIO: dropped threshold correction (medium). The old DIO code added a near-threshold correction before computing the fraction; the new shared helper omits it: // removed in DIOGenerator_tool.cc
double pmin = _spectrum.getAbscissa(0);
double pdfmin = _spectrum.getPDF(0);
double binsize = _spectrum.getBinWidth();
fullintegral += 0.5*pdfmin*pmin/binsize;This changes the DIO 2. 3. Debug 4. |
|
☀️ The build tests passed at fd5c2a7.
N.B. These results were obtained from a build of this Pull Request at fd5c2a7 after being merged into the base branch at 20895af. For more information, please check the job page here. |
|
Responding to the AI comments, after discussion with Dave about the integral calculation:
|
| } | ||
| // The "full" integral may be missing the lowest tail component in some cases (e.g. DIO) | ||
| // --> apply a linear interpolation correction if requested | ||
| if(correct_full_integral) { |
There was a problem hiding this comment.
What is the use case for making this correction optional?
There was a problem hiding this comment.
Is this correction needed in general? Should the proton/neutron KE spectra also include it? In the case that a different variable was given (perhaps one, like cos(theta), that doesn't go 0 - x) would we need to instead provide the x-value for the absolute edge that isn't included in the BinnedSpectrum PDF?
I only made it optional as I didn't understand if all BinnedSpectrum PDFs are missing a component that can be approximated with a linear interpolation to 0 or not.
|
I suggest taking 1 step at a time:
- fix the bug in the integral calculation of DIO spectrum, where it's use
is unambiguous.
-factorize it out with proper control for the general case (if possible)
simply relying on the user to understand how the spectra are being used and
configure this correctly seems error prone.
…On Sat, Jul 11, 2026 at 10:46 AM Michael MacKenzie ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In EventGenerator/inc/ParticleGeneratorTool.hh
<#1878 (comment)>:
> + pset.put(var_low, full_var_low);
+ pset.put(var_high, full_var_high);
+ BinnedSpectrum fullSpectrum(pset);
+
+ // Calculate the integrals
+ double integral = 0.;
+ double fullIntegral = 0.;
+ for(size_t ibin=0;ibin < spectrum.getNbins();++ibin){
+ integral += spectrum.getPDF(ibin);
+ }
+ for(size_t ibin=0;ibin < fullSpectrum.getNbins();++ibin){
+ fullIntegral += fullSpectrum.getPDF(ibin);
+ }
+ // The "full" integral may be missing the lowest tail component in some cases (e.g. DIO)
+ // --> apply a linear interpolation correction if requested
+ if(correct_full_integral) {
Is this correction needed in general? Should the proton/neutron KE spectra
also include it? In the case that a different variable was given (perhaps
one, like cos(theta), that doesn't go 0 - x) would we need to instead
provide the x-value for the absolute edge that isn't included in the
BinnedSpectrum PDF?
I only made it optional as I didn't understand if all BinnedSpectrum PDFs
are missing a component that can be approximated with a linear
interpolation to 0 or not.
—
Reply to this email directly, view it on GitHub
<#1878?email_source=notifications&email_token=ABAH57YATQ3T7JCETXUQNGT5EJ4PJA5CNFSNUABKM5UWIORPF5TWS5BNNB2WEL2QOVWGYUTFOF2WK43UKJSXM2LFO4XTINRXHA2DMOBUGI42M4TFMFZW63VMORSWC3K7NVSW45DJN5XKKZLWMVXHJLDGN5XXIZLSL5RWY2LDNM#discussion_r3564725389>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ABAH57Y2EYPQYV4TEZJJMKL5EJ4PJAVCNFSNUABFKJSXA33TNF2G64TZHMZDAMRYGA2DMOJSHNEXG43VMU5TIOBVG44DGNJTGQ22C5QC>
.
You are receiving this because you are on a team that was mentioned.Message
ID: ***@***.***>
--
David Brown ***@***.***
Office Phone (510) 486-7261
Lawrence Berkeley National Lab
M/S 50R5008 (50-6026C) Berkeley, CA 94720
|
Is the current implementation following this or do you want the linear interpolation from 0 to bin 0 on by default? I can also update to use |
|
I just pushed an update to by default correct the low edge of the full PDF integral if |
PR Review Summary — "Neutron generator updates" (follow-up at
|
| element | table | _energyFraction |
previously | now |
|---|---|---|---|---|
| C | heeck_finer_binning_2016_szafron-scaled-to-6C.tbl |
2.989e-06 |
all kept (r=1.0) |
~1 per 334,529 stops |
| H | heeck_finer_binning_2016_szafron-scaled-to-1H.tbl |
6.924e-07 |
all kept (r=1.0) |
~1 per 1,444,214 stops |
With PrimaryFilter.MinimumPartMom : 1.0 (IPAMuminusMichel.fcl:8) the job produces effectively nothing. It's also a double-correction in waiting: DIOGenerator_tool.cc:121 already records that same fraction in SpectrumConfig — which is the correct mechanism for a variance-reduced primary, and it is now applied twice.
The root cause is pre-existing and not yours: CompositeMaterialGenerator is a primary generator (it lives in JobConfig/primary/, emits dts.owner.IPAMuminusMichel…) that declares isPrimary=false. It is the odd one out —
| module | passes | semantics |
|---|---|---|
SingleProcessGenerator_module.cc:106 |
true |
primary ✔ |
Pileup_module.cc:119/125, MuStopProductsGun_module.cc:114/119 |
false |
pileup ✔ |
CompositeMaterialGenerator_module.cc:125 |
false |
primary ✘ |
That mislabel was harmless while r ≡ 1.0. This PR makes it fatal, so it has to be resolved here.
Suggested fix: have CompositeMaterialGenerator pass isPrimary=true, and/or state explicitly which generators are meant to get rate scaling. IPAMuminusMichel is the only Production config using this module, so the fix is contained.
2. [S1] Neutron rate silently drops 3.62% with no energy restriction configured
MuCapNeutronGenerator_tool.cc:48 calls the helper with correct_full_integral defaulted to true (ParticleGeneratorTool.hh:48), and :60 folds the result into the Poisson mean.
The live production config — Mu2e/Production JobConfig/pileup/prolog.fcl:294-300 (muonCaptureNeutronGenTool) — is spectrumShape: tabulated with Offline/EventGenerator/data/neutronSpectrum.txt and no elow/ehi. So spectrum and fullSpectrum are built from identical psets and integral == fullIntegral == 298731.0 exactly. The correction at ParticleGeneratorTool.hh:73-78 fires anyway, because xmin = getAbscissa(0) = 0.5 > full_var_low = 0:
corr = 0.5*(22471.0/0.5)*(0.5 - 0) = 11235.5
_energyFraction = 298731/(298731 + 11235.5) = 0.963753 # expected: exactly 1.0
This is live — MuStopPileup.fcl routes the tool through Pileup_module.cc with isPrimary=false. It also poisons the data product: :130 records RestrictedVar(kineticEnergy, 0.9638, xmin=0.25, xmax=100.25), i.e. "96.4% sampled" over what is in fact the full range.
Why the posted validation missed it. The check quoted above ("identical integrals … 0.999779") is DIO/Czarnecki. czarnecki_Al.tbl runs down to 0 → 0, so full.getPDF(0) = weight(0.05) = 0 exactly (CzarneckiSpectrum.cc:38-41 returns _table(0) = 0). With pdfmin = 0 the triangle is identically zero both before and after — which is exactly why those integrals matched. That test structurally cannot see this. The neutron table has pdfmin = 22471 at its low edge and was never exercised.
Why it looks accidental. DIO passes true explicitly at both 81db65d35 and 977788c4b, so flipping the default changed nothing except the neutron:
| commit | default | DIO call site |
|---|---|---|
81db65d35 |
correct_full_integral = false |
(fullconfig, true) — explicit |
977788c4b (head) |
correct_full_integral = true |
(fullconfig, true) — explicit |
The head commit's xmin > full_var_low guard was clearly meant to stop exactly this spurious correction ("correct the low edge of the full PDF integral if xmin > <provided xmin value>"). It doesn't, because of finding 3. Note this PR's own commit message says "…if energy is restricted" — nothing is restricted here.
Suggested fix: revert the default to false (DIO is explicit, so that's a zero-diff change for DIO), and apply finding 3. Fixing finding 3 alone still leaves −1.85%.
3. [S1] The guard and the correction measure from the bin centre, not the bin low edge
ParticleGeneratorTool.hh:73 reads fullSpectrum.getAbscissa(0), but BinnedSpectrum.hh:163 is explicit: "x values for Binned Spectrum is saved as BIN CENTERS not as left edge." The low edge is getXMin() (BinnedSpectrum.hh:143). The code contradicts its own comment on :77, which says "spectrum low edge".
Bin 0's PDF already covers [getXMin(), getXMin()+bw], so extending the triangle to the centre re-adds a half-bin that is already counted. The over-count is exactly pdfmin/4, independent of bin width:
Δ = 0.5*(pdfmin/bw)*[(getXMin()+bw/2) - L] - 0.5*(pdfmin/bw)*[getXMin() - L] = pdfmin/4
For the neutron table that's 22471/4 = 5617.75 — 1.88% of fullIntegral, making the correction exactly 2× too large (base 0.5 instead of 0.25). Denominator too big → fraction too small → rate too low.
It's also why the new guard never fires: for any spectrum whose support already starts at full_var_low, getAbscissa(0) == full_var_low + bw/2 > full_var_low, so the guard always passes.
Suggested fix: const double xmin = fullSpectrum.getXMin();. That makes the guard the exactly-correct predicate — "the spectrum's own lower edge lies above the requested lower bound" ⟺ something is genuinely missing. For the production DIO/Czarnecki config getXMin() == 0.0, so 0 > 0 is false and the correction self-disables, which is right since that spectrum already spans [0, endpoint].
4. [S2] The elow == ehi == 0 "full spectrum" idiom is only valid for tabulated
Concrete evidence for @brownd1978's "proper control for the general case" concern. The docstring (ParticleGeneratorTool.hh:51) promises "if elow == ehi, it does the full spectrum". Against BinnedSpectrum.cc that holds for tabulated (the elo > 0/ehi > 0 guards make 0/0 a genuine no-op) and vacuously for Czarnecki (which never reads ehi from the pset — :33 hardcodes getEndpointEnergy(), so full_var_high is a silent no-op there). Elsewhere it does not:
flat(BinnedSpectrum.cc:37-41):BinWidth = ehi - elow = 0→ throwsBADCONFIG "invalid binWidth = 0".RMC,Bistirlich,ceLeadingLog:initialize<Shape>(0,0,res)iterates zero times → emptyabscissa→BinnedSpectrum.hh:109/:99index an empty vector → undefined behaviour.ejectedProtons: ignoreselow/ehientirely → same phantom-correction class as finding 2.
This matters here because MuCapNeutronGenerator_tool.cc:44 computes _flatSpectrum from spectrumShape == "flat" and reports kFlat at :131 — a flat neutron spectrum is an anticipated configuration — and :48 now calls the helper unconditionally, so such a config would throw at construction. Pre-existing for DIOGenerator; new exposure for MuCapNeutronGenerator. No flat neutron config exists in-tree today, so this is latent rather than observed.
Suggested fix: short-circuit (return 1.) or throw a clear BADCONFIG for shapes that don't support the idiom, and correct the docstring.
5. [S2] DIOGenerator_tool.cc:51's comment is stale — that correction is identically zero
The comment says "correcting for missing 0 - 1 MeV component", but per finding 2 the triangle is exactly 0 for Czarnecki: the full spectrum already includes 0–1 MeV by direct summation, and that region carries ~7e-5 of the spectrum.
Credit where due: the old code took pmin/pdfmin from the restricted spectrum (abscissa 1.05), which was genuinely wrong and double-counted a region fullspect already contained. Switching to fullSpectrum is a real fix. Only the comment needs updating.
6. [S2] The helper doesn't belong on the interface header
ParticleGeneratorTool is a header-only interface (no .cc). A 38-line inline body with loops now pulls BinnedSpectrum.hh into ~14 translation units for the benefit of the 2 that call it, and any BinnedSpectrum.hh edit rebuilds all of them. The function touches no member state, so it isn't interface behaviour at all.
Suggested fix: make it a free function in Mu2eUtilities next to BinnedSpectrum (e.g. BinnedSpectrumUtils.hh/.cc), which also removes the new include from the interface header. Both callers already link Offline::Mu2eUtilities and already include BinnedSpectrum.hh, so the move is zero-cost. Failing that, at least mark it static.
7. [S2] Mu2eXGenerator_tool.cc:59-62 still carries the exact bug DIO just fixed
double pmin = _spectrum.getAbscissa(0); // RESTRICTED spectrum — the old DIO bug
double pdfmin = _spectrum.getPDF(0); // RESTRICTED spectrum
fullintegral += 0.5*pdfmin*pmin/binsize; // unconditionalIt reports _energy_fraction at :126 and didn't receive this fix. Consolidating it onto the new helper is exactly what hoisting the helper was for. (RMCGenerator is the other tool with a genuinely restrictable spectrum — prolog.fcl:373-374 sets elow: 85, ehi: 90.1 — but it needs care: the helper's blind erase/put would desync kMaxUserSet/kMaxUser.)
8. [S3] Minor
DIOGenerator_tool.cc:43-45:integralis now computed and never read — a refactor leftover, along with its "compute normalization" comment. (Checked rather than assumed: this does not break-Werror;integral += …counts as a read for-Wunused-but-set-variable, andfd5c2a789built green with this exact code.) The ctor also now builds threeBinnedSpectrums where the merge-base built two.classes_def.xml:44: understood that the dictionary entries silenced warnings in your test files. Wasart::Wrapper<std::pair<...>>specifically among them?std::pair<std::string,RestrictedVar>is never produced as a data product, and this is the onlyart::Wrapper<std::pair<...>>in the repo — 25+ otherstd::pairselections, none Wrapper'd. Line 43 looks correct and idiomatic either way.ParticleGeneratorTool.hh:47namesfhicl::ParameterSetby value (complete type required) with no#include "fhiclcpp/ParameterSet.h";MuCapNeutronGenerator_tool.cc:49-50addsstd::coutwith no<iostream>. Both compile only via transitive include fromBinnedSpectrum.hh— which finding 6 would remove.std::string var_low/var_highby value →const std::string&. Params are snake_case while the same function's locals are camelCase.full_var_low = 0.is the right floor forkineticEnergy/momentumbut wrong fortotalEnergy(floor = mass). Not reachable today — every current config useskineticEnergy— butMuCapNeutronGenerator's switch can emit"energy"forTOTAL_ENERGY.- Pre-existing but interacting: the neutron fcl omits
BinCenter: true. The table's values (0.25, 0.75, … spacing 0.5) are unambiguously bin centres spanning 0–100, butBinnedSpectrum.cc:95defaults tofalse, shifting support to 0.25–100.25. That shift is what keeps the triangle alive even after thegetXMin()fix — withBinCenter: trueandgetXMin(), the fraction is exactly 1.0. RelatedlyEventGenerator/fcl/prolog.fcl:301still claims "'tabulated' option ignores elow and ehi", contradictingBinnedSpectrum.cc:93-95.
Verified correct — no action needed
- No missing link dependency. Only 2 TUs ODR-use the helper (it's non-virtual, so the other ~12 includers emit no
BinnedSpectrumsymbol), and bothEventGenerator/CMakeLists.txt:537and:581already listOffline::Mu2eUtilities; it also propagates transitively viaLIBRARIES PUBLIC. SCons passes one flat per-directory list includingmu2e_Mu2eUtilities. - Multiplying the rate by the fraction is sound in principle.
physicsParams.Al.capture.neutronRate = 1.2is a total per-capture multiplicity over the table's 0–100 MeV support, and the siblingRMCRatecomment documents that convention ("full rate taken by dividing by closure fraction above 57 MeV"). The mechanism is right; the fraction is wrong (finding 3), and it must not reach primaries (finding 1). Worth annotatingneutronRate = 1.2while here — it carries no comment on what range it spans. - No
_isPrimaryasymmetry between the two tools.MuCapNeutronGenerator_tool.cc:88isn_gen = (_isPrimary) ? 1 : _randomPoissonQ->fire()— the Poisson is never fired when primary, so the fraction on its mean cannot down-scale primary generation. Semantically equivalent to DIO'sif(_isPrimary || …). - The correction is dimensionally consistent — the bare sums are
I/bw, and0.5*(pdfmin/binsize)*(xmin - full_var_low)is the triangle area in the same units. Binning is uniform in all three real tables, so the bare-sum ratio is valid. - The
pseterase/put is safe.spectrumis afhicl::DelegatedParameterin both tools,.get<ParameterSet>()returns a copy, and the helper takespsetby value — the mutation is local, as you said. classes.hneeds no matching edit; the non-conststd::pairspelling is the established Mu2e idiom; no schema-evolution risk.- The restriction mechanism itself works. With
elowset, the neutron fraction tracks the true tabulated fraction (0.41 / 0.22 / 0.078 at 5 / 10 / 20 MeV). The −3.62% is a uniform multiplicative bias on top — the feature this PR is for does work, once the bias is removed.
Validation check
- Build/tests run: not on this head.
mu2e/buildtestreports pending — "This test has not been triggered yet." on977788c4b,81db65d35and27cdacdd0. The only green FNALbuild isfd5c2a789, the first of four commits, and its own description says "Merged fd5c2a7 at 20895af" — i.e. it validated a pure refactor, before either rate change existed. Every behaviour-changing line in this PR is CI-unvalidated. mu2e/muDauStepsis the test that would exercise both blockers (Production/Validation/muDauSteps.fclincludesMuStopPileup.fcl, whose prolog carriesmuonCaptureNeutronGenToolanddioGenTool). It ran only onfd5c2a789.- Config contract check: partial —
elow/ehiare the right keys fortabulated/Czarnecki, but the idiom is unsafe for other shapes (finding 4). - Cross-repo: no
ProductionPR is required — and that's the concern. Both blockers land through unchanged Production configs (pileup/prolog.fcl,primary/IPAMuminusMichel.fcl) with no Production diff to review.mu2e-trig-configis unaffected.
Residual risk
- Any future adopter of the helper with a non-
tabulatedspectrum gets a throw or UB rather than a clear error. SpectrumConfig::fraction_is written but never read anywhere in Offline today, so there's no double-count yet — but the product doesn't record whether the fraction was already applied (it is, only when non-primary), so a future consumer can't tell. Worth stating the contract inSpectrumConfig.hh.
Author follow-ups
IPAMuminusMichel(finding 1) is the urgent one — could you confirm the DIO rate scaling isn't meant to reach primaries, and fixCompositeMaterialGenerator'sisPrimary=false?- Was the
correct_full_integraldefault flip in977788c4bintended to change the neutron rate? The commit message suggests not. - Please use
getXMin()rather thangetAbscissa(0)(finding 3) and revert the default tofalse(finding 2) — together they restore the unrestricted neutron fraction to exactly 1.0. - Cheap confirmation:
mu2e -c EventGenerator/test/Pileup.fclprints[MuCapNeutronGenerator] Sampled spectrum fraction 0.963753where1is expected. Re-running that before/after checks the neutron path directly, which the DIO integral test can't cover. - Could you re-trigger
buildteston the current head, and add intent/scope/validation evidence to the PR description?
|
Addressing the AI comments:
|
|
@FNALbuild run build test |
|
⌛ The following tests have been triggered for 8653622: build (Build queue - API unavailable) |
|
☀️ The build tests passed at 8653622.
N.B. These results were obtained from a build of this Pull Request at 8653622 after being merged into the base branch at dde444e. For more information, please check the job page here. |
4 participants
This is so we can simulate the high energy neutron tail component for calo cluster background studies