Puzzling issue: Bragg Curve Changes with Thread Count

Geant4 Version: 11.02 patch 1
Operating System: Ubuntu 22.04 (Linux Kernel v 6.5)
Compiler/Version: GCC 11.4
CMake Version: 3.22.1


Hi all,

I’m having a puzzling problem related to, multithreading, G4PVParameterisation, and parallel worlds. In brief, as I vary the number of simulation threads the proton dose curve changes (which obviously is not physically accurate.).

I have attached figures of this below:


Figure 1: Dose curve for 80 MeV protons in water when simulated using 10 threads.


Figure 2: Dose curve for 80 MeV protons in water when simulated using 2 threads.

I have created a minimal reproducible example of the issue and uploaded it to my GitHub. The issue seems to require a fairly specific set of circumstances to occur. I’m creating parameterised rectangular prism volumes using G4PVParameterisation and placing them in a parallel world. If instead I parameterise my physical world and score using the physical geometry I don’t see the same issue. Additionally, if I replace the G4PVParameterisation with a G4Replica the user code works regardless of whether scoring is done in the physical or parallel world.

I’m sharing this not because I need help with a problem – the issue is fully resolved by using G4Replica. Instead I am sharing this because it might be indicative of some sort of unidentified issue with the Geant4 kernel. I am hopeful that is not the case and instead there is some issue with my user code that I am yet to identify.

Best,

Joseph

Hi,

Thanks for this report. Could you please upload the parallel world construction with G4Replica as well so that we can trace the difference?

BTW, could you please confirm the following lines in your code repository are actually what you intend?

In your RunAction.cc

   analysisManager->CreateH1("Dose", "Dose along phantom profile", _numBins, 0, _numBins);

the last argument is the maximum value of the x-axis. And in your EdepScorer.cc

   analysisManager->FillH1(0, index, edep);

you plot the copy number (not the position) as the x-value. This won’t make the x-axis with unit of mm as shown in your plots.

Hi Joseph,
one thing I notice is that in your parameterisation you introduce a “fudge factor”, multiplying by 0.95 the Z length of the boxes, to “give a little spacing between boxes”… see:
MinimalExample/src/StraightLineParameterisation.cc at main · JDecunha/MinimalExample · GitHub
If this is not matched in your implementation with replicas, then the geometry will be clearly different, as with replicas, in general, there’re NO gaps. Am I missing something?

Hi @makotoasai,

I will update the example to include the implementation with G4Replica as you requested (I will make another post in this thread to let you know once that is complete). Regarding filling the H1, yes, that is as intended, I store the copy number in the H1 and convert the copy number to a physical length when plotting.

@gcosmo, regarding the factor to introduce additional spacing between the boxes. I included that because when working with G4PVParamterisations in the past (typically with curved surfaces) I found that it was necessary to ensure that the faces of the paramaterised surfaces do not touch. I introduced that factor here as part of the troubleshooting. It means that there is a tiny bit of space between the scoring boxes compared to the G4Replica approach, but the boxes are centered around the same positions. I will remove it and re-run, but I believe I previously saw the same results when running without the spacing factor.

Edit: You can disregard what I have written below the dashed line. Removing the gap between parameterised volumes is actually what has resolved the issue and not changing from “GetReplicaNumber” to “GetCopyNumber” in the scorer. For further investigation of the issue, as you requested, I have uploaded a branch which includes an implementation of the user code with G4Replica. You can find that branch here: GitHub - JDecunha/MinimalExample at ReplicaGeometry.

To summarize the current situation. Despite understanding that the issue is related to having a gap between adjacent parameterised volumes. The issue still stands, that for some reason the dose profile returned is substantially different when varying the number of active cores used for the simulation.


@makotoasai I believe I have somewhat figured out the source of the issue. The scorer being used in the user code called called “GetReplicaNumber” rather than “GetCopyNumber”. I didn’t realize there was a meaningful difference between the two functions. I suppose that finding the replica number of a G4PVParameterisation is undefined behavior?

Joseph

@JDecunha, both GetCopyNumber(G4int depth) and GetReplicaNumber(G4int depth) are
valid calls, as calling them from G4TouchableHistory they’re doing exactly the same thing.

I noticed however that you’re using a custom scorer in your code, and in your implemented ProcessHits() there’re no calls to ComputeCurrentSolid(aStep), which I believe is necessary when using parameterisations. This, I guess, might explain the problem you see,