Understanding how production cuts are used

[mkelsey]

We are trying to make use of physics-list production cuts to reduce the number of short steps and low-energy secondaries in our simulation (we are interested in net energy deposit, not tracking low-energy ionization electrons or X-rays).

I am seeing unexpected differences in validation jobs we’re using to determine useful cut values. The validation jobs consist of a spherical source of single tracks, aimed inward at a 78 mm R x 25 mm H germanium crystal. I count the total number of energy deposits (steps with GetTotalEnergyDeposit() > 0.) in the crystal, as a function of the primary track energy.

I ran two sets of jobs, one with FTFP_BERT, and one with Shielding_EMZ (EM option4), each set consisting of incident gammas, electrons, or neutrons with a flat energy range from 1 to 300 keV. For each job, we set all the particle cuts to the same value, from 0.01 to 2.00 mm, as follows:

thePL->SetCutValue(cutLength, "gamma");
thePL->SetCutValue(cutLength, "e-");
thePL->SetCutValue(cutLength, "e+");
thePL->SetCutValue(1e-9*mm, "proton");   // See all nuclear recoils
G4EmParameters::GetInstance()->SetApplyCuts(true);

The incident gammas behave exactly as I’d expect, as shown in the two plots below, for Shielding_EMZ and FTFP_BERT, respectively. You can even see the energy thresholds where the cut value kicks in for secondaries (red, green, blue, and yellow, below).

In both cases, the uppermost points in black are for no production cuts at all (SetApplyCuts(false)), and the large error bars reflect that there are tails out to very high numbers of hits (well over 100 in a single event). I don’t fully understand the differences between Shielding_EMZ vs. FTFP_BERT in this case, but I guess they’re due to EM_option4 vs. EMStandard?

Electrons and neutrons both show quite different results, and this is where my main problem lies. Here are the same plots for incident electron primaries:

Notice that all of the different production cut runs overlay one another: there is no reduction in the number of hits, and no threshold effect. There’s no apparent difference between the no-cut job and any of the finite cuts, even for cuts as large as 2 mm. If an incident gamma can show differences of an order of magnitude in the number of hits, why don’t electrons show the same?

The neutron plots look similar to the electron plots, though with an average of just two hits per event, and very long tails up to over 200 hits in a small number of outliers. Again, I don’t understand how the production cuts can have such a strong effect on incident gamms, but not on other primaries.

I’m hoping one of the EM experts can help me to understand either what’s going on, or maybe what I’m doing wrong in configuring these jobs.

[civanch]

Hello Mike,

you need to take into account, that ApplyCuts flag enable/disable secondary thresholds for gamma processes only (photo-effect, Compton,…). For ionisation and bremsstrahlung this flag not used at all, only value of a cut affects number of secondaries.

VI

[mkelsey]

Thanks, Vladimir. This is still confusing me. The documentation and your comment above says,

other processes will use it [production cut] as a default value for their intrinsic limits (e.g., ionisation and bremsstrahlung).

So this means if we set a finite production cut it applies to ionisation and brem unconditionally?

Later in the documentation, it says,

The ParticleDefinition (or ParticleWithCuts) has a Boolean data member: ApplyCut

I see the comment in the documentation at the end of the “Apply Cuts” section,

Note, that for electromagnetic processes for gamma incident a specific ApplyCut option is used which guarantees energy balance and is more efficient because secondary tracks are not produced at all.

This sounds like what you’re talking about above.

Ignore everything below
So, if we want all electron or gamma secondaries restricted by the production cut, should we set the ParticleDefinition::SetApplyCutsFlag() ourselves, for electrons, gammas, etc.? For nuclear recoils, would we set the flag just for GenericIon, or do we need to treat the deuterium, tritium, and He-3 particles individually?

[mkelsey]

I found the G4VUserPhysicsList::SetApplyCuts(flag, "name") function, which is what we need. However, I’ve run into an “order of operations” problem with our physics configuration.

In our framework, we do not have our own physics list class that we use. Instead, we start with one of the reference physics lists (via G4PhysListFactory), and add the additional physics we need by way of RegisterPhysics() calls. This works really well, and allows us to encapsulate different kinds of optional physics in physics-constructors, just the way you’re supposed to.

Unfortunately, that means we don’t have the flexibility to write our own SetCuts() function, and physics constructors don’t have that interface.

If I try to call the SetApplyCuts() function in our setup code (which all happens during PreInit), then I get a segfault because that function queries the G4ParticleTable which is not yet populated.

I don’t mind (and have already set up) carrying the cuts and apply-cuts flag value as local data members for use with our UI commands. But I don’t know how to defer applying those local data members during /run/initialize, without forcing users to put that command into their macros (which we have never required).

[arce]

I have a UI command to do this (before /run/initilize) and I just call
G4EmProcessOptions popt;
popt.SetApplyCuts(GmGenUtils::GetBoolean(newValue));

[mkelsey]

Thanks, Pedro. That’s what I was doing (using G4EmParameters instead of the deprecated G4EmProcessOptions). However, Vladimir explained that this flag only applies to incident gammas (it’s used by Compton, Photoelectric, etc.). It doesn’t apply to incident electrons or other particles.

The physics list SetApplyCuts() set the particle flag that’s used by the Tracking Manager for “all” secondaries, no matter what process created them. My trouble is that the physics list function requires the particle types to have already been created, so it can’t be called before /run/initialize.

[civanch]

Hello Mike,

Unfortunately, description of cuts in the documentation requires some revision. I agree, that there is a space for confusion to users and hope that we will clean-up for the major release of this year.

I would suggest to follow cut definition in extended/electromagnetic examples: TestEm8, TestEm9, where it is shown how to define cuts per region. Cuts for gamma processes are enabled in G4EmStandardPhysics_option1. Cuts may be defined via UI commands or C++ interface.

VI

[mkelsey]

Thanks, Vladimir! I ended up writing our own set of commands, because I couldn’t find one to call the “ApplyCuts” function in the physics list, and I wanted to be sure I could cache the settings during PreInit even if our users switched physics lists via macro command.

Once I figured out that I needed to set both the EM (gamma) SetApplyCuts() as well as the physics list, our event-time performance increased dramatically (order of magnitude or so).

The one setting I haven’t run an extenstive test on is the EM deexcitationIgnoreCuts(false) flag. We do a fair amount of radiogenics, and I think we want to set this to avoid all of the X-rays etc. in the material bulk.

[civanch]

Hi Mike,

deexcitationIgnoreCuts = false by default, it is enabled if a user needs to see all transitions in radioactive decay. In your case it is not needed.

Vladimir