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).

[michaelmackenzie]

Why both RMCN0External and RMCPhaseSpace0NExternal (for example)?

[sophiemiddleton]

I ran my test_normalization with Run-1A and the yield was the same as what you sent me from your calculations. These are just nicknames really, but yes they should match the full name I guess

[michaelmackenzie]

This is set at 75 instead of 80 MeV. And why use hard-coded R_spectrum value instead of doing the integral of the plestid spectrum from emin to kmax? Will this hold for arbitrary emin values?

[michaelmackenzie]

Same comment here as above

[sophiemiddleton]

good spot, I will check this

[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).