nrc-cnrc / EGSnrc

Toolkit for Monte Carlo simulation of ionizing radiation — Trousse d'outils logiciels pour la simulation Monte Carlo du rayonnement ionisant
http://nrc-cnrc.github.io/EGSnrc
GNU Affero General Public License v3.0
243 stars 146 forks source link

Electron beam range too short for energies above 100 MeV #961

Closed camots closed 1 year ago

camots commented 1 year ago

I wanted to simulate a 150 MeV electron pencil beam impinging on a water phantom. Pegs-Data was created up to 500 MeV. There is a big difference compared to Geant4 or Penelope: Looking at the depth-dose curve, the range is about 10cm shorter compared to Geant4 or Penelope.

Thus my question: Is EGSnrc accurate for energies > 100 MeV?

ftessier commented 1 year ago

Interesting question, I have not used EGSnrc above 50 MeV myself! Nominally, EGSnrc is accurate up to 10 GeV. Can you share the ranges given by all three codes, along with the EGSnrc application you are using (DOSXYZnrc, egs_chamber, etc.), and the Monte Carlo parameters your are using?

camots commented 1 year ago

Here is a picture for the depth-dose a 150 MeV pencil beam in water.

The parameters were like this:

#------------------------------------------------------------------------
# Transport parameters
:start MC transport parameter:
   Global ECUT=                    0.600
   Global PCUT=                    0.100
   Global SMAX=                    1e10
   Bound Compton scattering=       On   # Can be Off, On in Regions, 
                                        # Off in regions
   Rayleigh scattering=            Off  # Can be On, On in Regions,
                                        # Off in regions
   Atomic relaxations=             On   # Can be Off, On in Regions,
                                        # Off in regions
   Photoelectron angular sampling= On   # Can be Off, On in Regions,
                                        # Off in regions
   Brems angular sampling=         KM   # Can be Off
   Brems cross sections=           BH   # use Bethe-Heitler (egsnrc default)
   Pair angular sampling=          KM   # Can be Off, Simple or KM
   ESTEPE=                         0.25
   XIMAX=                          0.5
   Skin depth for BCA=             3
   Boundary crossing algorithm=    exact      # Can be PRESTA
   Electron-step algorithm=        PRESTA-II  # Can be PRESTA
   Spin effects=                   On         # Can be off
:stop MC transport parameter:
#------------------------------------------------------------------------

geant4-egs-15MeV.pdf

image

ftessier commented 1 year ago

Wow, quite a difference indeed! What EGSnrc application are you using?

camots commented 1 year ago

I'm using an egs++ (own development) code derived from advanced application.

The same code gives good agreement for energies < 30 MeV. As an example see the depth-dose comparison for a 15 MeV pencil beam (lateraly integrated, as was in the 150 MeV case).

I also verified the pegsless variant - same result. Version 2020 and 2021 - same result.

Best regards!

geant4-egs.pdf

image

ftessier commented 1 year ago

Is it feasible to check with another application, say egs_chamber? I mean is the geometry simple enough, e.g., rectilinear voxels? Although your 15 MeV result kind of confirms that your code is correct. But just to be sure. If you are able to share your geometry input, that would be helpful.

Next check, what are the dose deposition region sizes in the 150 MeV case? There is an internal transport parameter in EGSnrc that has given me some grief before, when the region size became too small compared to the mean free path.

At any rate, it is interesting to investigate this issue, thanks for reporting it and pursuing it!

ftessier commented 1 year ago

Can you share the first graph, but NOT normalized to maximum dose? I would like to know where the discrepancy occurs. I am perplexed by the discrepancy at depth 0 cm.

ftessier commented 1 year ago

According to ESTAR, the CSDA range of 150 MeV electrons in water is 42 cm. You could perform a CSDA calculation with EGSnrc as a baseline, to check that the code is otherwise correct.

kinetic energy CSDA range g / cm $^2$
1.500E+02 4.204E+01

image

camots commented 1 year ago

Thanks for the input!

The voxelsize in the 150 MeV case was 0.5 x 0.5 x 0.5 cm3, 101 x 101 x 400 voxels. The 15 MeV was calculated with a phantom of 0.2 x 0.2 x 0.2 cm3, 200 x 200 x 200 voxels.

Calculate the 15 MeV on the 0.5 x 0.5 x 0.5 cm3, 101 x 101 x 400 for consistency...

Would like to check some the 1

ftessier commented 1 year ago

Note: you can turn Off bound Compton scattering and atomic relaxations if that reduces simulation time.

ftessier commented 1 year ago

The voxel size is fine. It becomes an issue if the voxel size is less than $10^{-8}\lambda$, where $\lambda$ is the particle range (see definition of $EPSEMFP). So not an issue here.

camots commented 1 year ago

Problem solved!!

The problem was caused by a bad density-correction file! Sorry for the trouble!

Thank you and best regards!

geant4-egs-150MeV.pdf geant4-egs-250MeV.pdf

ftessier commented 1 year ago

Great, problem solved! I was quite perplexed because typically transport at higher energies is "easier", given that the relativistic approximations become more accurate. But you never know, always good to check!