Improve the optical photon stacking action

[gipert]

At the moment, we send optical photons to the waiting stack and track them at the end only if there is an energy deposition in germanium with the desired energy:

remage/src/RMGGermaniumOutputScheme.cc

Lines 361 to 385 in af65926

	std::optional<G4ClassificationOfNewTrack> RMGGermaniumOutputScheme::StackingActionClassify(
	const G4Track* aTrack,
	int stage
	) {
	// we are only interested in stacking optical photons into stage 1 after stage 0 finished.
	if (stage != 0) return std::nullopt;
	// defer tracking of optical photons.
	if (fDiscardPhotonsIfNoGermaniumEdep &&
	aTrack->GetDefinition() == G4OpticalPhoton::OpticalPhotonDefinition())
	return fWaiting;
	return std::nullopt;
	}
	std::optional<bool> RMGGermaniumOutputScheme::StackingActionNewStage(const int stage) {
	// we are only interested in stacking optical photons into stage 1 after stage 0 finished.
	if (stage != 0) return std::nullopt;
	// if we do not want to discard any photons ourselves, let other output schemes decide (i.e. not
	// force `true` on them).
	if (!fDiscardPhotonsIfNoGermaniumEdep) return std::nullopt;
	const auto event = G4EventManager::GetEventManager()->GetConstCurrentEvent();
	// discard all waiting events, if there was no energy deposition in Germanium.
	return ShouldDiscardEvent(event) ? std::make_optional(false) : std::nullopt;
	}

Still, just stacking itself seems to slow down the simulation quite a bit. We should try to understand if we can be even smarter and avoid stacking photons in some situations. For example, if the energy deposited in argon is above a certain threshold far enough from any other surface.

All of this should of course carefully validated. Any other idea?

[tdixon97]

One thought, we should find what takes time in the stacking/ photon generation? For example if the issue is generating eg. directions for the optical photons maybe we can postpone that to later?

[tdixon97]

Here is a plot to demonstrate the issue.
I used the L-1000 geometry and I ran 226-Ra and 214-Bi simulations using the EdepLow and EdepHigh cuts in Germanium to only run full optical sims for events in a range.

I generated 10^6 primaries and the fraction of full optical simulated events is on the x-axis, while the time is on the y-axis. I also overlay the time for a simulation without optics turned on as the line.

I think this shows the overhead, since the intercept of the points is not the same as a simulation without optics. In addition, you see that the slope of the line is fairly low, indicating that the CPU cost of the optical sims themself is fairly low.

From this plot I would estimate if we can fix the stacking we can have at least a factor of 5 gain in CPU time for direct optical simulations of a narrow range in energy.

[tdixon97]

Update: The situation is dramatically worse for Th, with a factor of around 100x difference. This is also related to the larger energy depositions in LAr

[tdixon97]

Any ideas @ManuelHu ?

[ManuelHu]

I do not really have an idea how to easily fix the stacking. There just seems to be a large performance hit in allocating all the tracks and generally in the stack if it gets large.

For example if the issue is generating eg. directions for the optical photons maybe we can postpone that to later?

unfortunately most of this already happens deep in the respective generating physics processes, so we cannot circumvent the allocation of the track objects (~ particle with momentum, starting point etc already generated) without reimplementing the whole scintillation code...

[tdixon97]

Sure but there is also an option to split the stacking into 3 steps:

• Track primary tracks (i.e betas or gamma from the decay)
• Send secondary tracks in LAr to waiting stack,
• Propagate tracks in Ge, or maybe also LAr close to the Ge (using distance to surface?)
• Only then track the secondary e- tracks in LAr

I think this would work well, the cost is that some tracks in LAr may eventually produce signals in Ge, so we have to be very careful.

I think this would be easy to implement but might take a bit of work to validate.

[EricMEsch]

Sure but there is also an option to split the stacking into 3 steps:

With the current infrastructure that should be possible to add as an optional filter scheme. But we probably would have to revisit the idea of properly staging the tracks. As right now stacks typically only get postponed into the next stage...

It might be an idea to add multiple stages, with each stage going farther away from the detector?

[tdixon97]

Yes this is what im suggesting to have a more sophisticated stacking action mechanism, this could potentially speed up some sims dramatically even without optics

[tdixon97]

@EricMEsch is it easy to have the "distance to the surface" readily available when deciding how to classify tracks this would be very nice

[EricMEsch]

Should be just one call. I am not sure how performant this call is, but the information should be available