Updates for ensembling

[sophiemiddleton]

This pull request adds support for Radiative Muon Capture (RMC) 0-nucleon (0N) and 1-nucleon (1N) knockout processes to the normalization and event calculation framework. It introduces new normalization functions, updates the event calculation and template generation scripts to handle these processes, and adds corresponding tests. There are also minor updates to configuration and test files to support these changes.

Normalization and Physics Calculation Enhancements:

• Added new physics constants and functions for RMC 0N and 1N knockout processes, including spectrum integration and normalization functions (rmc_0n_normalization, rmc_1n_normalization, and supporting Plestid spectrum utilities) in normalizations.py [1] [2].
• Extended argument parsing and event calculation in calculateEvents.py to support new RMC 0N and 1N processes, including both internal and external conversion options [1] [2] [3].

Template Generation and Configuration:

• Updated make_template_fcl.py to allow selective calculation of normalization for new RMC 0N/1N processes, with argument parsing and logic for these processes [1] [2].
• Modified FCL configuration to drop the new mu2e::SpectrumConfig product, keeping output files clean.

Testing Improvements:

• Added comprehensive unit tests for RMC 0N and 1N normalization functions, including both internal and external yields, and updated summary printouts in test_normalizations.py [1] [2] [3] [4] [5].
• Adjusted test setup to use more realistic spill times.

These changes collectively enable precise modeling and normalization of RMC 0N and 1N backgrounds, improving the physics reach and flexibility of the event calculation framework.

[FNALbuild]

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

• JobConfig

which require these tests: build.

@Mu2e/fnalbuild-users, @Mu2e/write have access to CI actions on main.

⌛ The following tests have been triggered for 274b68d: build (Build queue - API unavailable)

About FNALbuild. Code review on Mu2e/Offline.

[FNALbuild]

☀️ The build tests passed at 274b68d.

Test	Result	Details
test with	✅	Command did not list any other PRs to include
merge	✅	Merged 274b68d at 5aa34ef
build (prof)	✅	Log file. Build time: 04 min 19 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.
FIXME, TODO	✅	TODO (0) FIXME (0) in 0 files
clang-tidy	✅	0 errors 0 warnings

N.B. These results were obtained from a build of this Pull Request at 274b68d after being merged into the base branch at 5aa34ef.

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]

Overall I think this looks great, thank you for making these updates! My only request/comment is to ensure the RMC normalization is valid for arbitrary values of emin (which may already be true).

[oksuzian]

PR Review Summary — #567

Reviewed at head 274b68d8c205c90201c45c9381bab6b67eb0226b. First pass; @michaelmackenzie's CHANGES_REQUESTED at this same head is carried forward below.

Decision

• 🟡 comment only

Both S1s are in the normalization arithmetic rather than the plumbing, and
neither fires at the committed defaults (Stage1 sets both thresholds to 80 MeV) —
hence comment rather than block. @michaelmackenzie's request is still open, and
his instinct was right; the numbers are below.

Scope understood

• Adds RMC 0-nucleon and 1-nucleon knockout normalizations (rmc_0n_normalization,
  rmc_1n_normalization) plus Plestid spectrum helpers, and wires them through
  calculateEvents.py, make_template_fcl.py and Stages 1/2/3.
• The ensemble's plain RMCInternal/RMCExternal entries are replaced by the
  0N/1N pair, not supplemented — BR_0N_FRAC + BR_1N_FRAC = 0.099 + 0.901 = 1.0,
  so this is the same RMC rate decomposed by knockout mode, and no background is
  dropped. Worth one line in the description, since "adds support" reads as
  additive.
• Also carries an unrelated DBVERSION parsing fix in Stage5_signal.sh and a
  test-livetime change.

Carried forward from @michaelmackenzie's review (state: CHANGES_REQUESTED, no commits since)

• "ensure the RMC normalization is valid for arbitrary values of emin" and the
  two inline comments on normalizations.py:605/:665 — UNADDRESSED, and
  finding 1 below gives the exact failure.
• "Why both RMCN0External and RMCPhaseSpace0NExternal?" — UNADDRESSED in
  code; the author has agreed to unify. See finding 5.
• "This assumes the internal and external photon energies have the same cuts" —
  CLEARED by the author's answer that it is intentional.

Findings

1. 🟠 [S1] For e_min ≤ 75 the spectrum factor is an absolute fraction where a
   relative one is required, under-predicting by 4.4× (0N) and 16× (1N)

   • Evidence: normalizations.py:601-606. The prefactor
     RMC_BR_MUON_CAPTURE * RMC_RATE_GT_57 * RMC_BR_0N_FRAC_GT_57 is already
     BR(0N, E > 57) per stopped muon, so R_spectrum has to be
     R(>e_min)/R(>57). The e_min > 75 branch computes exactly that —
     FRAC_0N_80 / FRAC_0N_57 = 0.03319/0.22887 = 0.14501. The default branch
     instead uses FRAC_0N_57 = 0.22887, whose denominator is the whole
     spectrum, applying a spectrum cut that the prefactor already applied. At
     e_min = 57 the correct factor is 1.0, so the result is low by
     1/0.22887 = 4.37×; for 1N, by 1/0.061620 = 16.2×.

   • Impact: any --rmcn0emin/--rmcn1emin at or below 75 MeV silently yields a
     too-small RMC background. Stage1 hardcodes 80, so today's default path takes
     the correct branch — this is a latent wrong-physics path, not a live one.
     The 75.0 cut is itself discontinuous: 74 MeV gets the E>57 normalization and
     76 MeV gets the E>80 one, with nothing in between.

   • Suggested fix: the machinery to do this exactly is already in the PR, unused.
     plestid_integral is the correctly normalized CDF (I verified its
     antiderivative integrates the plestid_spectrum shape to 1 over [0, kmax]),
     so the whole branch collapses to

     R_spectrum = (plestid_integral(e_min, kmax, kmax, knockout) /
                   plestid_integral(57.0,  kmax, kmax, knockout))

     and it then holds for arbitrary e_min, which is precisely what
     @michaelmackenzie asked for. The two kmax values are recoverable from your
     own constants — see finding 3.

2. 🟠 [S1] Three constants are re-declared locally, and one disagrees with the
   cited shared value

   • Evidence: normalizations.py:62-70 against constants.py, which this same
     module already imports at line 14:
     • RMC_RATE_GT_57 = 1.41e-5 vs RMC_GT_57_PER_CAPTURE = 1.43e-5
       — different values for the same quantity, and the shared one carries a
       source (Phys. Rev. C 59, 2853 (1999)) while the new one carries none.
     • RMC_BR_MUON_CAPTURE = 0.609 duplicates CAPTURES_PER_STOPPED_MUON = 0.609.
     • RMC_INTERNAL_EXTERNAL_RATIO = 0.0069 duplicates INTERNAL_PER_RMC = 0.00690.
   • Impact: rmc_normalization (line 490-495) uses the imported constants while
     rmc_0n_normalization/rmc_1n_normalization use the local copies, so the
     old and new RMC paths now disagree by 1.4% on the same rate. Whichever value
     is right, two homes for it means a future correction reaches only one.
   • Suggested fix: delete all three locals and use the imported names. If
     1.41e-5 is a deliberate update, change it in constants.py with its source
     and let rmc_normalization move with it.

3. 🟡 [S2] plestid_integral and plestid_spectrum are added but never called

   • Evidence: normalizations.py:517-575, 59 lines; the only other occurrences
     of "Plestid" in the file are docstring prose in the two normalization
     functions, which use the hardcoded constants instead.
   • Worth keeping rather than deleting, because the four hardcoded fractions
     are this function: solving plestid_integral(e, kmax, kmax, knockout) for
     the kmax that reproduces each one gives 101.866 MeV from FRAC_0N_57
     and 101.865 MeV from FRAC_0N_80, and 95.449 MeV from both 1N
     values — self-consistent to 0.001 MeV within each mode, and physically
     sensible (1N is lower by roughly the neutron separation energy). Those two
     numbers are the missing input; name them, and finding 1's fix is a two-line
     change. Note neither is the RMC_kmax=90.1 that Stage1 still writes and
     Stage2 still passes as --rmckmax.

4. 🟡 [S2] The four new tests assert only > 0, and CI does not run them

   • Evidence: test_normalizations.py:204-275 — each new test ends in
     assertGreater(..., 0). A 4× or 16× normalization error passes. FNALbuild's
     table for build 658 lists 16 jobs (build, ceSimReco, ceMix, rootOverlaps,
     g4surfaceCheck, …) and no python test target, so nothing runs this file
     automatically either.
   • Suggested fix: pin one expected yield per mode to a few significant figures.
     Given finding 1, a second case at e_min = 57 would have caught it.

5. 🟡 [S2] on_spill_time is scaled by an uncited 88/496, identically for 1BB and 2BB

   • Evidence: test_normalizations.py:28-29 — 4.4e6 * (88/496) = 0.1774, applied
     to both modes. constants.py gives ONEBB_DF = 0.323 and TWOBB_DF = 0.246;
     0.1774 is neither, and the two modes have different duty factors, so one
     shared literal cannot be right for both.
   • Suggested fix: use the named duty-factor constants, or comment where 88/496
     comes from.

6. 🟡 [S2] Stage2_build_sampler.sh: one dataset is validated on disk and then
   listed from tape, and the new file-list check only warns

   • Evidence: line 214 validates every dataset with
     mu2eDatasetFileList "$dataset_name" --disk, but the RPCInternal list at
     line 287 is built with --tape — alone among the six list-building calls,
     the other five pass no flag. Separately, check_file_lists (line 265, used
     at 301 and 363) prints
     ⚠️ WARNING: Some file lists may be empty or incomplete and execution
     continues; an empty list surfaces much later as an IndexError at
     make_template_fcl.py:201 or a ZeroDivisionError at line 164.
   • Suggested fix: drop the stray --tape, and make the file-list check exit 1.
     A check that warns and proceeds reads like a gate but is not one.

7. ⚪ [S3] Three small dead/contradictory items

   • e_threshold = 57.0 is assigned and never read, in both new functions
     (normalizations.py:600, :660).
   • The alias assignments norms["RMCN0External"] = norms["RMCPhaseSpace0NExternal"] = ...
     (make_template_fcl.py:63-73) only ever populate keys nothing reads — the
     consumer loop is for signal in args.prc (line 85) and Stage2 passes only the
     RMCPhaseSpace* names. This is the code side of @michaelmackenzie's
     "why both names?"; dropping the short aliases resolves both.
   • Stage1_initiate_ensemble.sh:24-25 — INCLUDE_RMCN0=1 # ... (default: no).
     The comment contradicts the value; both modes are on by default, so
     --rmcn0/--rmcn1 are opt-out, not opt-in.

Verified, no action needed

• 🟢 The Stage5_signal.sh DBVERSION change is a real fix, not a rewrite. The
  old $(IFS='_'; echo "${PARTS[@]:2}") does not join on IFS — only "${arr[*]}"
  does, "${arr[@]}" expands to separate words that echo joins with spaces. Run
  on MDC2025_best_v1_3 it returns v1 3; the new explicit loop returns v1_3.
• 🟢 make_template_fcl.py:52,55 switching args.tmin → tmin is safe and is
  itself a fix: tmin is defined at line 26 as float(args.tmin), so
  rpc_normalization now receives a float rather than the raw string.
• 🟢 CI is green at this head — FNALbuild build 658, all 16 rows ✅, including
  rootOverlaps, g4surfaceCheck, clang-tidy (0/0) and FIXME/TODO (0/0).

Validation check

• Build/tests run: CI green at 274b68d8 (build 658). The new python tests are
  not part of that suite; I did not execute them. My own checks were arithmetic
  on the constants and a bash reproduction of the DBVERSION parse.
• Config contract check: partial. The Stage1 → Stage2 variable handoff
  (RMC_N0_emin/RMC_N1_emin gating dataset lists, --prc, and sampling
  options) is internally consistent; I did not run the stages end to end.
• Cross-repo consistency: n/a — no Offline or trig-config interface is touched.

Residual risk

• I did not verify that the datasets the new Stage2 branches expect
  (dts.mu2e.RMCPhaseSpace{0,1}N{Internal,External}.MDC2025{au,at}.art) exist with
  enough files, nor that RMCVERSIONINT=au/RMCVERSIONEXT=at are the intended
  versions — that split is asymmetric and unexplained in the diff.
• The physics inputs themselves (0.099/0.901 knockout split, and whether the
  Plestid phase-space shape is the right model here) are outside what I can check;
  I verified only internal consistency and the arithmetic.

Author follow-ups

1. Fix the e_min ≤ 75 branch (finding 1) — using plestid_integral answers
   @michaelmackenzie's arbitrary-e_min request in the same change. The two kmax
   values you need are 101.87 MeV (0N) and 95.45 MeV (1N), recovered from your own
   constants.
2. Reconcile 1.41e-5 against the cited RMC_GT_57_PER_CAPTURE = 1.43e-5 and drop
   the three duplicate constants (finding 2).
3. Pin at least one expected yield in the new tests, and add an e_min = 57 case
   (finding 4).
4. Drop the stray --tape and make the file-list check exit non-zero (finding 6).

[sophiemiddleton]

I have now pushed updates and fixes to your responses

[michaelmackenzie]

I don't think the 80 MeV values are needed here, if emin is an input then in the RMC normalization it should evaluate the fraction above the input emin using the integral function directly

[sophiemiddleton]

I think this is now ready for a merge.

[michaelmackenzie]

This looks good to me. I just had a couple of comments about no longer used code, but these are just comments and not required

[michaelmackenzie]

I don't think these 80 MeV versions are used anymore

[michaelmackenzie]

Are these values used anywhere? Is this function still needed?

[sophiemiddleton]

thanks, I thnk I was using these for validation

[oksuzian]

PR Review Summary — #567

Reviewed at head 3c13c2828abb36da5f001b38079dc3568b55b305. Re-review of 274b68d8.

Decision

• 🔴 request changes

The RMC work is now right, and I verified it numerically rather than taking the
commit message for it — the arbitrary-e_min fix @michaelmackenzie asked for is
in, and it reproduces the production default to better than 0.02%. What blocks is
something else the delta carries: one line was deleted from rpc_normalization,
which multiplies the RPC background by 6.2 and is mentioned nowhere in the PR
description, the commit messages, or any review thread.

I note @michaelmackenzie approved at this head on 2026-08-20 with comments he
explicitly called non-blocking. His review was of the RMC normalization, and I
agree with it; finding 1 below is in a different function that the RMC discussion
never touched.

Scope of this pass

• 274b68d8..3c13c282, four commits ("fixed norms" ×2, "made emin an input arf"
  ×2), across normalizations.py (+140/−89), calculateEvents.py (−16) and a new
  test_rmc_fractions.py (+76).

---

Findings

1. 🔴 [S0] filter_efficiency was deleted from the RPC normalization, raising the RPC yield by a factor of 6.2, in a PR about RMC.

   • Evidence: normalizations.py:535 still computes
     filter_efficiency = float(num_pion_filters) / float(num_pion_stops), but the
     product at :538-544 no longer contains it. On main the same product does
     — normalizations.py:401 there is the filter_efficiency * line this delta
     removes. The variable is now assigned and never read.
   • Impact: with the values the file's own comments record for the table it reads
     (num_pion_stops = 41324703 at :131, num_pion_filters = 6634478 at
     :133), filter_efficiency = 0.1605, so every RPC yield this module produces
     — internal and external, both BB modes — comes out 6.23× larger than
     before. rpc_normalization is live: make_template_fcl.py calls it for the
     ensemble. Nothing in the PR description, the four commit messages, or any
     review thread mentions RPC.
   • Suggested fix: if this is deliberate — say, because the filter efficiency is
     already folded into target_stopped_pions_per_pot and was double-counted —
     then state it in the description, delete the now-dead filter_efficiency
     line, and split it into its own PR so the RPC change is reviewed as an RPC
     change. If it is not deliberate, restore the line. The dead variable left
     behind is what makes me think it is the second.

2. 🟠 [S1] Carried over: three constants are still re-declared locally, and one still disagrees with the cited shared value.
   normalizations.py:62-63 still has RMC_BR_MUON_CAPTURE = 0.609 and
   RMC_RATE_GT_57 = 1.41e-5, against constants.py:9
CAPTURES_PER_STOPPED_MUON = 0.609 and constants.py:11
RMC_GT_57_PER_CAPTURE = 1.43e-5 # Source: Phys. Rev. C 59, 2853 (1999);
   RMC_INTERNAL_EXTERNAL_RATIO = 0.0069 still duplicates
   constants.py:22 INTERNAL_PER_RMC = 0.00690. The old rmc_normalization uses
   the imported names and the new 0N/1N functions use the locals, so the two RMC
   paths still disagree by 1.4% on the same rate, and only one of the two values
   carries a source. Delete the locals, or move 1.41e-5 into constants.py with
   its citation.

3. 🟡 [S2] The plain-RMC branch in calculateEvents.py now computes a yield and prints nothing, and its external twin was deleted.
   calculateEvents.py:41-42 is if(args.prc == "RMC" and int(args.internal) == 1): Yield = rmc_normalization(...) with no print, where every other branch in the
   function prints a <Name>_yield= line for the shell to parse; the
   int(args.internal) == 0 branch is gone entirely. The only remaining caller is
   scripts/old/Stage1_makeinputs.sh:110,112, which would now append nothing to
   ${TAG}.txt instead of failing. Either finish the removal — drop both branches,
   the way the RMCN0*/RMCN1* aliases were dropped — or restore the print.

4. 🟡 [S2] test_rmc_fractions.py is not a test and duplicates the function it checks.
   The new 76-line file re-declares plestid_integral verbatim at :11-33 rather
   than importing it from normalizations.py, so it validates its own copy and
   will silently stop tracking the real one. It contains no assertion — it computes
   four numbers, prints a comparison table and exits 0 regardless — while the
   test_ prefix means unittest and pytest discovery will collect it and record a
   pass. Import plestid_integral and assert the four fractions against the
   literature values, or move it out of test_* and call it a script.

5. 🟡 [S2] The four RMC_SPECTRUM_FRAC_* constants and the function that computes them are now entirely dead — this is @michaelmackenzie's :79 and :222 comments, which he marked non-blocking; I am recording the fuller extent rather than raising the severity.
   The literals at :77-80 are overwritten at import by
   compute_rmc_spectrum_fractions() (called at :260), and nothing then reads
   any of the four: rmc_0n_normalization and rmc_1n_normalization go through
   rmc_spectrum_fraction() instead, and git grep RMC_SPECTRUM_FRAC across
   JobConfig/ensemble finds no other consumer than test_rmc_fractions.py, which
   defines its own. So it is not only the 80 MeV pair — all four, plus the
   compute_rmc_spectrum_fractions function and its global statements, can go.

Prior findings

• 1. 🟠 [was S1] The e_min ≤ 75 branch — FIXED in d61e3be, verified. Both
     functions now call rmc_spectrum_fraction(e_min, RMC_KMAX_*, knockout=*), which
     is plestid_integral(e_min)/plestid_integral(57) — the relative form the
     prefactor requires. I checked the arithmetic rather than the diff: it returns
     exactly 1.000000 at e_min = 57 for both modes, 0.0 above k_max, and there
     is no discontinuity left at 75 MeV. See the verified section for the agreement at
     the production default.
• 2. 🟠 [was S1] Duplicate constants — UNADDRESSED, carried as finding 2.
• 3. 🟡 [was S2] plestid_integral added but never called — FIXED. It is now
     the sole source of the spectrum fractions, and the two k_max values in
     :68-69 (101.8667 and 95.4489 MeV) are the ones I recovered from the hardcoded
     constants last pass, to within the precision I quoted.
• 4. 🟡 [was S2] The four new tests in test_normalizations.py assert only > 0
     — UNADDRESSED; that file is untouched by the delta. Finding 4 above is a
     second instance of the same pattern.
• 5. 🟡 [was S2] on_spill_time scaled by an uncited 88/496 — UNADDRESSED.
• 6. 🟡 [was S2] Stage2_build_sampler.sh: the stray --tape on the RPCInternal
     list, and check_file_lists warning instead of exiting — UNADDRESSED;
     Stage2_build_sampler.sh is untouched by the delta.
• 7. ⚪ [was S3] Of the three small items: the dead e_threshold = 57.0 in both
     new functions is FIXED (removed with the branch). The
     make_template_fcl.py:63-73 aliases and the contradictory
     Stage1_initiate_ensemble.sh:24-25 comment are UNADDRESSED, both files being
     untouched — though the calculateEvents.py side of the alias question is now
     resolved, which was @michaelmackenzie's point.

Verified, no action needed

• 🟢 The new spectrum fraction reproduces the old correct branch at the production
  threshold, so this fix does not move any yield anyone has already computed. At
  e_min = 80, plestid_integral(80)/plestid_integral(57) gives 0.145039 for
  0N against the old 0.03319/0.22887 = 0.145017 — a relative difference of
  1.5e-4 — and 0.021381 for 1N against 0.0013175/0.061620 = 0.021381, to
  1.7e-5.
• 🟢 The two k_max values are consistent with all four literature fractions, not
  just the two the fix uses: computing plestid_integral(E, k_max, k_max, n) / plestid_integral(0, k_max, k_max, n) gives 0.228873 vs 0.22887, 0.061619 vs
  0.061620, 0.033195 vs 0.03319 and 0.0013175 vs 0.0013175 — the largest
  disagreement being 0.016%.
• 🟢 rmc_spectrum_fraction is well behaved outside the range that mattered
  before: e_min above k_max returns 0.0 rather than a negative or a divide
  error, and e_min = 0 returns 4.369 = 1/0.22887, which is the right answer for
  a quantity defined relative to E > 57.

Validation check

• Build/tests run: none complete at this head. mu2e/buildtest has been
  pending at 3c13c282 since 2026-08-20T15:52:51Z and FNALbuild has posted no
  table for it; the green table on this PR is from build 658 at 274b68d8. I did
  not trigger another, since one is already queued — and it would not settle
  anything here, because the Production job list contains no python test target, so
  none of these files run in CI either way. My own checks were the arithmetic
  above and git grep for the dead constants.
• Config contract check: partial, unchanged from the last pass.
• Cross-repo consistency: n/a — no Offline or trig-config interface is touched.

Residual risk

• Whether removing filter_efficiency is physically correct is not something I can
  judge; I can only say it is a 6.2× change that is not described anywhere.
• The physics inputs themselves — the 0.099/0.901 knockout split, and the two
  k_max endpoints as values for Al-27 — remain outside what I can check. I
  verified internal consistency, not the literature.

Author follow-ups

1. Restore filter_efficiency in rpc_normalization, or state the reason for
   removing it and split it out (finding 1). This is the one item I would hold the
   merge for.
2. Reconcile 1.41e-5 against RMC_GT_57_PER_CAPTURE = 1.43e-5 and drop the three
   duplicate constants (finding 2, carried).
3. Finish or revert the plain-RMC branch removal in calculateEvents.py
   (finding 3).
4. Make test_rmc_fractions.py import plestid_integral and assert, or rename it
   (finding 4).
5. Delete the four dead RMC_SPECTRUM_FRAC_* constants and
   compute_rmc_spectrum_fractions (finding 5).

[sophiemiddleton]

Responses:

Author follow-ups

1 Restore filter_efficiency in rpc_normalization, or state the reason for
removing it and split it out (finding 1). This is the one item I would hold the
merge for.

--- > this fixes a bug it needs to be removed

Reconcile 1.41e-5 against RMC_GT_57_PER_CAPTURE = 1.43e-5 and drop the three
duplicate constants (finding 2, carried).

---> can do only one is used 
Finish or revert the plain-RMC branch removal in calculateEvents.py
(finding 3).

---> no reason to do this, it might be useful to look over old data sets 
Make test_rmc_fractions.py import plestid_integral and assert, or rename it
(finding 4).

--> we dont need this file I can remove 

Delete the four dead RMC_SPECTRUM_FRAC_* constants and
compute_rmc_spectrum_fractions (finding 5).

---> I can remove