Hello everyone,
I want to implement in Geant4 a new process that describes the annihilation of a positron with an electron to an invisible final state, through a s-channel diagram:
e+ e- -> A’ -> CHI CHI
where CHI is the (invisible) final state particle and A’ is the intermediate (virtual) particle. I know the cross section for this process, including the finite A’ width.
If this process is implemented as a G4VDiscreteProcess, the method GetMeanFreePath would be used to return the mean free path. This is computed from the total cross section, using the positron energy at the end of the previous step.
I am unsure about how the energy loss across the step due to ionization and Bremmstrahlung can be handled here: due to resonant process, the annihilation cross section depends significantly on the positron energy, and the positron looses energy along the step due to the aforementioned continuous processes.
Dear VI,
I apologize to reply after a long time to the post.
If I understood correctly, according to the G4 physics manual (here), sec. 7.4, the method that G4 is using to solve this issue is based on the idea that the energy-dependent cross section along the step is replaced by an upper limit for the cross section 𝜎𝑚. This is used to determine the mean free path of the particle, before the step takes place.
Then, if the process is selected by G4 as the one delimiting the step length, in the PostStepDoIt() method, the sampling of the final state is performed only with the probability 𝑝 = 𝜎(𝐸𝑓 )/𝜎𝑚, alternatively no interaction happen and tracking of the particle is continued. (cit. from reference manual).
What is not clear to me is why the probability is computed with the cross section calculated at the end-of-step energy, and not via the average cross section along the step, that is the one proportional to the interaction probability, as shown by the Equation reported in sec. 7.4 of the physics manual:
where here Ei and Ef are, respectively, the energy at the beginning and at the end of the step. Since the PostStepDoIt method is called after the step is determined, the integral can be computed, I assume.
To confirm that actually G4 is using the energy at the end of the step, I looked at the source code of the G4VEmProcess class, PostStepDoIt method (here), and I saw that really the energy at the end of the step is used.
May I ask you for a clarification? I am sure there’s something that I am missing here!
Dear colleagues,
I’d like to provide an update on this.
Model
Consider a simple model where there is a pseudo-positron particle moving forward in a material and loosing a constant, deterministic amount of energy per unit length −α. The only other process available for this particle is the annihilation process, with a simplified cross section formula:
The number density of electrons per unit volume is k.
Let’s compute the expected flux of pseudo-positrons at depth z in the material N(z), starting at z=z0 with energy E=E0. This model can be solved exactly, starting from the differential equation:
The energy loss is deterministic, so we have E(z)=E0−αz. Thus:
This differential equation can be solved exactly, giving:
This gives a number of interactions per unit length, −dN/dz, equal to:
I implemented a simple toy simulation, with the algorithm for energy-dependent cross section implemented. I used both the cross section computed at the end of the step, and the average cross section along the step. The result is reported below; the analytical estimate (green) agrees with the calculation using the final energy.
we will have CPU penalty but will not have an exact result if will try to average x-sections, because dynamics of tracking and discrete processes will mess up any such attempt.
In the recent Geant4 11.1, both for EM and hadronic cross sections we introduced G4EmTableType.hh with enumerator assuming 5 different types of cross sections. For hadronic physics we add G4HadXSType with also 5 types. If your cross section fit to any of these types then the tracking will be nearly exact, not dependent on geometry and other circumstances. “Nearly” means that for cross section increased when energy is decreasing it may be needed to add step limitation to guarantee unbiased tracking.
Of course, if you have a narrow peak in the cross section or your cross section does not fit to any type, then you will need to take care on step limits and provide tracking with small steps.
Hi VI,
thanks. May I ask you if you have a reference where the integral approach for the energy dependency of the cross section in MC is described in details?
I tried to look at: V.N.Ivanchenko et al. In Proc. of Int. Conf. MC91: Detector and event simulation in high energy physics, 1991, HEP INDEX 30 (1992) No. 3237, 79–85. Amsterdam, 1992 but no formulas where discussed here.
The relevant part in this reference was:
When particle coming through matter is simulated, all absorption coefficients are calculated, then the total absorption coefficient is found which defines the range of the particle before the interaction point. Depending on cut parameters for the geometrical volume, this way may be subdivided into several steps to obtain necessary accuracy for ionization losses and multiple scattering. For charged particles before generation of secondary particles the cross section of the interaction may be recalculated to take into account the energy dependence of the cross section. Technically it is done by using the flag of empty interaction.
