How to very accurately model protons with energy higher than 100 MeV?

Hello

Could you please help me. I need to obtain energy distribution in a layer of liquid water (layer of “G4_WATER” material with 5-20 micrometers thickness) penetrated by the beam of protons with energy up to 400 MeV at nanoscale (on 10-30 nm voxels) as precise as possible.
When the energy of protons is lower than 100 MeV I can use DNA models which (I assume) are the most accurate models available in Geant 4.10.5. But these models cannot be used for energies higher than 100 MeV and that’s where the problem starts. I tried to use combination of DNA and Livermore models by taking physics list from “microdosimetry” example and replacing standard models with models from “G4EmLivermorePhysics” constructor (which seems to be the most accurate after DNA-models-based constructors). Part of PhysicsList file that describes processes for protons:

//G4StepLimiter* fStepLimiter = new G4StepLimiter();
//pmanager->AddDiscreteProcess(fStepLimiter);

G4hMultipleScattering* msc = new G4hMultipleScattering();
msc->AddEmModel(1,new G4WentzelVIModel());
ph->RegisterProcess(msc, particle);

G4CoulombScattering* pcou = new G4CoulombScattering();
ph->RegisterProcess(pcou, particle);

G4hIonisation* hion = new G4hIonisation();
hion->SetEmModel(new G4BraggModel(), 1);
hion->SetEmModel(new G4BetheBlochModel(), 2);
//hion->SetStepFunction(0.01, 1 * nm);
ph->RegisterProcess(hion, particle);

G4hBremsstrahlung* pb = new G4hBremsstrahlung();
ph->RegisterProcess(pb, particle);

G4hPairProduction* pp = new G4hPairProduction();
ph->RegisterProcess(pp, particle);

G4NuclearStopping* pnuc = new G4NuclearStopping();
ph->RegisterProcess(pnuc, particle);

/
//G4DNAExcitation, G4DNAIonisation, G4DNAElastic and G4DNAChargeDecrease described here with models set to G4DummyModel to be activated futher as it is done in “microdosimetry” example
/

I also tried to use G4StepLimiter with the step size of 1 nm and to use SetStepFunction (see commented lines in code). Range cuts for all particles were set to 1 nm. This physics list also uses “G4ElectronCapture” process from “microdosimetry” example that kills all electrons with energies under 7.4 eV.
I attached graphs of resulting dose distributions on voxels for protons with the energies of 99.9 MeV (DNA models are used for primary protons) and 100.1 MeV (non-DNA models are used for primary protons). So each point on graphs represents number of voxels that recieved particular dose during simulation. You can see 3 lines corresponding to 3 sets of parameters used:

  • initial energy of protons = 99.9 MeV, no user limits used (green line)
  • initial energy of protons = 100.1 MeV, no user limits used (red line)
  • initial energy of protons = 100.1 MeV, 1 nm user limits used (blue line)

You can see that using user limits (G4StepLimiter) doesn’t solve the problem: green and blue lines should almost match, but they do not. Number of voxels in high doses area differs several times for these two lines. Such result is not acceptable because this area is especially important for our investigation. SetStepFunction on hIoni process did not affect the results at all.
I would appreciate any ideas how to overcome this discrepancy between models for prorons. The

same problem for electrons doesn’t exist: dose distributions for primary electrons lower (0.99 MeV) and higher (1.01 MeV) than energy limit of applicability of DNA models are very close to each other. So is it possible to somehow get accurate results for protons above 100 MeV?

Hello,

please, try out EM physics Option4. You may check results using simple setup of $G4INSTALL/examples/extended/electromagnetic/TestEm8

VI

Thank you for the comment, but unfortunately “G4EmStandardPhysics_option4” constructor doesn’t solve the problem. Concerning electrons and protons it is actually almost 100% copy of “G4EmLivermorePhysics” constructor with the difference that “G4EmLivermorePhysics” includes nuclear stopping for protons and “G4EmStandardPhysics_option4” does not. The results (graphs) that I posted in previous message were obtained litterally by using combination of “G4EmLivermorePhysics” constructor and DNA models. So the replacement of “G4EmLivermorePhysics” with “G4EmStandardPhysics_option4” doesn’t change anything.
If you suggest not using DNA models at all and use “G4EmStandardPhysics_option4” for all energies of protons and electrons then could you please give more comments about why would it be better for the accuracy of the resulting dose distribution?

Hello,

I may be missunderstand your initial post. Will try to clarify several things.

Firstly, what is the difference between DNA and standard models of hadron ionisation? DNA model is discrete assuming simulation of each ionisation collision . it is technically applicable below 100 MeV. Energy deposition may be computed as total number of ionisations inside the voxel, ionisation may be counted as total number of electrons stopped inside the voxel.

G4hIonisation process is the same in all EM standard PhsyicsLists implements condence history approach. The difference is only in step function, which more tight in Opt4. This process simulates continues ionisation and production of delta-electrons above cut. Continues ionisation at a step is smeared using energy loss fluctuation model. This approach is consistent if the cut for electron energy is low enough allowing sample electron leakage from a voxel correctly.

If the cut or the step limit become too small, there are difficulty connected with limited applicability of fluctuation model for very small steps. In solid or liquid media we estimate safe smalest step as few micrometer. When you ask for nanometer step you may get biased energy depositions.

Another problem - how to score energy deposition in a voxel. More straitforward to me would be making real boundaries - you make cubic voxels of few micron size and compose your geometry from such volumes - in that case, energy deposition per voxel will be accurate and you will have euls for Opt4 similar to DNA at 1 MeV proton beam.

If you indeed need 1 nm voxel, unfortunately, only DNA like approach may work properly.

VI

Thank you very much for the comprehensive explanation.

Unfortunately it’s not possible to use 1 micrometer voxels, since in our model each voxel represents small part of DNA molecule and its size has to be 10-30 nm.

I hoped that there might be some trick with parameters of G4hIonisation process which could allow to obtain results similar to ones obtained with DNA models. As I mentioned in the end of my first post Livermore model for G4eIonisation process for electrons with the energy of 1 MeV gives results that are close to ones obtained with DNA models. So for ions this is not possible, I got it.

Hello,

I would not say , that it i not possible at all but the model of energy loss fluctuations may bias results. How much, I cannot say, it is not factors but some % of shift of mean energy deposit.

VI

Hello!
I am also interested in the question. Is the difference in voxel size only?
For example, with a voxel size Z = 1mm,
physical list (Penelope, Livermore. GEANT-DNA) will show similar results (e.g. dose) ?

Thanks you!