Neutron generator updates

[michaelmackenzie]

This is so we can simulate the high energy neutron tail component for calo cluster background studies

[FNALbuild]

Hi @michaelmackenzie,
You have proposed changes to files in these packages:

• EventGenerator
• MCDataProducts

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)

About FNALbuild. Code review on Mu2e/Offline.

[oksuzian]

Core change

New helper in the base class centralizes what DIOGenerator did inline:

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;
}

MuCapNeutronGenerator now actually uses the fraction (previously hardcoded 1.), and DIOGenerator was rewritten to call the helper and renamed _energy_fraction → _energyFraction.

Issues / concerns

1. 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 _energyFraction value (and thus normalization) versus current behavior. The comment noted it corrects "the missing prediction near threshold, assuming the rate falls to 0 linearly." Confirm this drop is intentional and acceptable for DIO — if not, the correction should be preserved (e.g. as an optional flag in the helper).

2. pset passed by value but reused after mutation (low, but verify). In both callers fullconfig is retrieved and passed to the helper. The helper takes pset by value and mutates it locally, so callers are safe — good. Just confirm no caller relies on fullconfig being modified afterward.

3. Debug std::cout left in production paths (low). Both generators unconditionally print sampled-fraction lines to stdout on every construction. Consider mf::LogInfo/a verbosity guard instead of raw std::cout for job-log cleanliness.

4. classes_def.xml addition — verify necessity (low). New dictionary entries for std::pair<std::string, RestrictedVar> and its wrapper are added. Confirm these are actually needed for persistence and that a matching entry exists in classes.h (only classes_def.xml is touched here) — mismatches cause genreflex/build failures.

[FNALbuild]

☀️ The build tests passed at fd5c2a7.

Test	Result	Details
test with	✅	Command did not list any other PRs to include
merge	✅	Merged fd5c2a7 at 20895af
build (prof)	✅	Log file. Build time: 08 min 41 sec
ceSimReco	✅	Log file.
g4test_03MT	✅	Log file.
transportOnly	✅	Log file.
POT	✅	Log file.
g4study	✅	Log file.
cosmicSimReco	✅	Log file.
cosmicOffSpill	✅	Log file.
ceSteps	✅	Log file.
ceDigi	✅	Log file.
muDauSteps	✅	Log file.
ceMix	✅	Log file.
rootOverlaps	✅	Log file.
g4surfaceCheck	✅	Log file.
trigger	✅	Log file.
check_cmake	✅	Log file.
FIXME, TODO	✅	TODO (0) FIXME (0) in 3 files
clang-tidy	➡️	4 errors 17 warnings
whitespace check	✅	no whitespace errors found

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.
Build artifacts are deleted after 5 days. If this is not desired, select Keep this build forever on the job page.

[michaelmackenzie]

Responding to the AI comments, after discussion with Dave about the integral calculation:

1. I updated the full spectrum integral to optionally include the low tail linear interpolation estimate, which includes a bug fix where it previously used the restricted PDF information. Before and after this update, I get identical integrals for an unrestricted DIO energy spectrum input of 0.999779 which shows the correction to the full spectrum integral matches. For DIO with E > 95 MeV the integral changes from 3.637e-11 --> 3.636e-11.
2. The intent is to not modify the pset used by the generators, so retrieving by value was intentional.
3. This is consistent with the previous implementation, promising that these fractions can always be cross-checked in the log files if needed.
4. I got warnings about these missing when writing test files, so they appear to be needed.

[brownd1978]

What is the use case for making this correction optional?

[michaelmackenzie]

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.

[brownd1978]

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

[michaelmackenzie]

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.

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 0.5*(<xmin> - <full_var_low>)*<pdf at xmin>/<bin width> so it's a bit more general in the case full_var_low is not 0.

[michaelmackenzie]

I just pushed an update to by default correct the low edge of the full PDF integral if xmin > <provided xmin value>