Muon stop pileup cz restriction

[michaelmackenzie]

This restricts the mu stop pileup cos(theta_z), and accordingly adjusts the rates for the processes so the pileup module will produce sensible distributions when only looking in this cos(theta_z) region. This significantly speeds up Run 1B pileup production for rough estimates that don't need exact precision.

[FNALbuild]

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

• EventGenerator

which require these tests: build.

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

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

About FNALbuild. Code review on Mu2e/Offline.

[FNALbuild]

☀️ The build tests passed at 12c5d86.

Test	Result	Details
test with	✅	Command did not list any other PRs to include
merge	✅	Merged 12c5d86 at e229837
build (prof)	✅	Log file. Build time: 08 min 27 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 6 files
clang-tidy	➡️	12 errors 23 warnings
whitespace check	✅	no whitespace errors found

N.B. These results were obtained from a build of this Pull Request at 12c5d86 after being merged into the base branch at e229837.

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.

[oksuzian]

Yes, here's my review of PR #1787. I found several issues:

---

Issues Found

1. 🐛 Inconsistent rate-scaling approach in DIOGenerator vs. other generators

The DIOGenerator uses a fundamentally different mechanism to account for the cos(θ_z) restriction compared to the other four generators, and it has a subtle bug:

const double r = (_czmax - _czmin)/2.;
if(_fireAll || _randFlat->fire() < r) {

The DIO generator is called once per muon stop (not via Poisson sampling), so it uses a RandFlat draw against r = (czmax - czmin)/2 to probabilistically skip events. However, RandFlat::fire() returns a value in [0, 1). When no cz restriction is applied (default czmin=-1, czmax=1), r = 1.0, and the condition _randFlat->fire() < 1.0 is almost always true — but strictly it can never equal 1.0, so it will fire ~100% of the time. This is functionally correct but semantically fragile. More importantly, every call now always consumes a random number from _randFlat even when there's no restriction, which changes the random number sequence compared to the baseline (before this PR), potentially affecting reproducibility of existing production configurations.

2. 🔤 Typo: "accouunt" in comments (4 occurrences)

In MuCapDeuteronGenerator_tool.cc, MuCapNeutronGenerator_tool.cc, MuCapPhotonGenerator_tool.cc, and MuCapProtonGenerator_tool.cc, the comment reads:

const double rate = _rate * (_czMax - _czMin)/2.; // accouunt for potential cz selection in the produced rates

"accouunt" → "account"

3. ⚠️ Random number sequence change for default configurations

For the Poisson-based generators (MuCapDeuteron, MuCapNeutron, MuCapPhoton, MuCapProton), the RandPoissonQ is initialized at construction with the adjusted rate:

_rate = GlobalConstantsHandle<PhysicsParams>()->getCaptureProtonRate(material);
const double rate = _rate * (_czMax - _czMin)/2.;
_randomUnitSphere = std::make_unique<RandomUnitSphere>(eng, _czMin, _czMax);
_randomPoissonQ = std::make_unique<CLHEP::RandPoissonQ>(eng, rate);

With defaults (czMin=-1, czMax=1), rate = _rate * 1.0, so the value is unchanged. However, passing _czMin, _czMax to the RandomUnitSphere constructor (previously constructed with no cz args) might still change the internal state of the random engine depending on RandomUnitSphere's implementation. This should be verified to not break reproducibility for existing default configs.

4. 🔍 DIOGenerator stores _czMin/_czMax as const but other generators don't

const double _czmin;
const double _czmax;
const bool _fireAll;

In all other generators, _czMin and _czMax are non-const. This is a minor inconsistency but worth aligning for code cleanliness — the other generators should also mark these const since they are set once in the constructor and never modified.

5. 💡 DIOGenerator debug output uses std::cout directly

std::cout << "Cos(theta_z) min " << _czmin << " max " << _czmax << std::endl;
std::cout << "Restricted Spectrum min " << _spectrum.getAbscissa(0) << " max " << _spectrum.getAbscissa(_spectrum.getNbins()-1) << std::endl;
// ... more std::cout lines

Lines 64–71 add 8 lines of std::cout diagnostic output that will print on every job. This looks like debug instrumentation that should either be removed before merge or gated behind a verbosity flag / replaced with the framework's mf::LogInfo messaging.

---

Summary

#	Severity	File	Issue
1	Medium	DIOGenerator_tool.cc	Unconditional RandFlat draw changes random sequence for default (no-restriction) configs
2	Low	4 files	Typo "accouunt" in comments
3	Medium	4 MuCap generators	Passing cz args to RandomUnitSphere ctor may change random sequence even for defaults — needs verification
4	Low	4 MuCap generators	_czMin/_czMax should be const for consistency with DIOGenerator
5	Medium	DIOGenerator_tool.cc	Unconditional std::cout debug output should be removed or use mf::LogInfo

[michaelmackenzie]

@FNALbuild run build test

[FNALbuild]

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

[FNALbuild]

☀️ The build tests passed at 12c5d86.

Test	Result	Details
test with	✅	Command did not list any other PRs to include
merge	✅	Merged 12c5d86 at fa82e62
build (prof)	✅	Log file. Build time: 04 min 06 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 6 files
clang-tidy	➡️	12 errors 23 warnings
whitespace check	✅	no whitespace errors found

N.B. These results were obtained from a build of this Pull Request at 12c5d86 after being merged into the base branch at fa82e62.

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.

[brownd1978]

I'm curious what prompted the addition of fireAll option.

[brownd1978]

What is the use case for fireAll = false? Shouldn't a generator produce an output ever event? Without that, downstream workflows that rely on MCPrimary will fail. You could add a filter in this case, but how would that be different from fireAll = true?

[michaelmackenzie]

@brownd1978 this fireAll is necessary as the module is both used in the muon stop pileup generator and as a standalone DIO generator. In the case of pileup, it would need to produce a DIO only the (czmax - czmin)/2 fraction of the time to keep rates physically consistent. The other pileup tool inputs similarly only produce outputs some fraction of the time, returning an empty collection if they don't fire. In the case where you're specifically simulating this process and not pileup, fireAll should be set (similar to the 1.8 MeV photon line generator) to ensure you always generate your primary of interest.

I can set this to be true for the DIO tail generator config in Production, or I can set it to false in the pileup tool config and have it default to true, whichever is preferred. I can also remove the default and make it explicit in all cases.

[FNALbuild]

📝 The HEAD of main has changed to 68688bb. Tests are now out of date.

[brownd1978]

@michaelmackenzie thanks for the explanation. As it's a formal error for a primary generator to not produce a primary, a way of insuring that in code would be better than a flag that can be misconfigured. How about passing this as a constructor argument, true when the generator is primary, false as part of pileup?

I also think some solution to recording the cz acceptance is required; we should not go back to harvesting numbers from log files to get correct downstream job configurations.