Photon only toy example

[SimonBiggs]

Issue to track the Photon toy example as discussed in https://github.com/nrc-cnrc/EGSnrc/discussions/658#discussioncomment-258317 by @ftessier. The text from there is copied in below.

---

Accessing the data files is perhaps not the next logical step. In light of the discussion regarding vectorization, it makes more sense to me now to focus on the vectorization logic.

Let's consider only photons for clarity, since for electrons a number of complications aris if one uses multiple-scattering. Say you have a list of photons, in an infinite medium, and 4 interaction (Rayleigh, photoelectric, Compton, pair production). We need to show significant efficiency gains with JIT for the following sequence:

• [ ] generate a random numbers for each photon (already vectorized in prng)
• [ ] sample the exponential distribution (each with its own cross-section value, depending on energy) for the distance to the next interaction (easily vectorized I presume)
• [ ] update photon positions (already vectorized in foregoing examples, could be done at the end)
• [ ] pick the interaction based on a random number and the relative interaction cross-sections (easily vectorized I presume)
• [ ] perhaps sort or mask the photon array if necessary to vectorize the interaction function call (next step)
• [ ] sample the interaction, using bogus function calls for now (can @jit efficiently vectorize (pure) function calls, i.e., map func() to array?)
• [ ] update the photons' direction and energy
• [ ] create an array of all the outgoing particles (except for Rayleigh).

[SimonBiggs]

sample the exponential distribution (each with its own cross-section value, depending on energy) for the distance to the next interaction (easily vectorized I presume)

@ftessier is there a good reference data file within EGSnrc that details how the exponential distribution is to be varied depending on energy?

[ftessier]

It is all in the cross-section tables; you can generate the actual curves using the examin with the xcom data set. But for trying it out you could cook up any dependency: the important aspect is that the exponential depends en energy (and position, if you think about geometries), so that the interpreter cannot cache the results etc.

[SimonBiggs]

Sorry for all the noob questions... where would I be able to find the specification for the layout/data structure of the following files?

• https://github.com/nrc-cnrc/EGSnrc/blob/master/HEN_HOUSE/data/xcom_compton.data
• https://github.com/nrc-cnrc/EGSnrc/blob/master/HEN_HOUSE/data/xcom_pair.data
• https://github.com/nrc-cnrc/EGSnrc/blob/master/HEN_HOUSE/data/xcom_photo.data
• https://github.com/nrc-cnrc/EGSnrc/blob/master/HEN_HOUSE/data/xcom_rayleigh.data

Anyway... time to head off to work. I might get a bit of a chance to play with this tonight. I'm trying to see if I can work it out by reading the examin.mortran file https://github.com/nrc-cnrc/EGSnrc/blob/b735e209236878963c2133f08a77c10c79f12308/HEN_HOUSE/user_codes/examin/examin.mortran#L234-L297 ... but no immediate wins.

[ftessier]

I don't know! 😂 I don't think formats or schemas were written for those, it is essentially hard-coded (!) and creating proper documented formats certainly is near the top of the list in an EGSnrc rewrite. The way I get by is to use the EGSnrc examin application to extract the cross section data. I can probably help with that...

[SimonBiggs]

🙂 glad at least it's not just me. I was beginning to feel a little silly 🙂

Regarding format choices, are you okay with plain old CSV files? That's what I opted to turn the brems data file into:

https://github.com/SimonBiggs/pyegsnrc/tree/main/pyegsnrc/data/bremsstrahlung

[SimonBiggs]

I can probably help with that...

If you're able to get the cross section data into a csv format that'd be massively helpful :) If it helps, here is the notebook I used to convert the brems data file:

https://github.com/SimonBiggs/pyegsnrc/blob/main/notebooks/reformat-nrc-brems-data.ipynb

[darcymason]

I'll need this too, also can help. At a quick glance at a number of data files, there is often a number of lines having varying count of numbers, then the consistent data set. Probably really need to know what that 'header' info is. I was going to try to track down where each is loaded and try to decipher from there.

[mchamberland]

xcom_photo.data was last modified by @mainegra in this commit

@ftessier, perhaps Ernesto can shed some light on those data files?

[ftessier]

Yes indeed, Ernesto knows these formats really well, as he created some of these files not so long ago. I am sending him a note. Typically, all these files have a similar format, with a header specifying the domain grid, and then the function values as a flat array, possibly in more than 1 dimension.

[ftessier]

Also see the EGSnrc manual, section 2.2.5:

The user can also supply their own photon cross section data. This can be accomplished by placing data files named xxxx photo.data, xxxx pair.data, xxxx triplet.data and xxxx rayleigh.data in the $HEN HOUSE/data folder, where xxxx stands for any string different from epdl, xcom or si. These cross sections will then be used if the variable photon xsections is set to xxxx. The format of the ASCII data files is very simple: for each element between 1 and 100 one has to provide the number N of data points in the tabulation followed by N pairs of the logarithm of the photon energy in MeV and the logarithm of the cross section in barn (see also one of the data files provided with EGSnrc).

[ftessier]

From @mainegra:

I would recommend to follow the code in subroutine HATCH and the subroutines called from within HATCH. That's how I check the formats on all those files. I agree that eventually there should be a proper description of all these formats. Maybe the original PEGS4 manual has some info on some of the files? But many have changed over the years.

[SimonBiggs]

Say you have a list of photons, in an infinite medium, and 4 interaction (Rayleigh, photoelectric, Compton, pair production)

Would it be okay to change that to three interactions, triplet, pair and compton? Reasoning being, those are the xcom cross section tables that I was able to reformat. See:

https://github.com/SimonBiggs/pyegsnrc/pull/3#issuecomment-757397960

[ftessier]

Yes, of course, that's fine! We can worry about Rayleigh later... (or not at all 😂 )

[SimonBiggs]

Just a heads up, I usually only get a couple of "project night" evenings a week, the next few project nights I'm going to utilise to address the following maintenance task within PyMedPhys:

https://pymedphys.discourse.group/t/user-guide-or-manuals/75/4?u=simonbiggs#post_4:~:text=I%20can%20pull%20that%20code%20back,and%20tweak%20it%20to%20your%20preferences%3F

In general, I'll prioritise PyMedPhys maintenance tasks before this here. Feel free to run ahead without me, but I shall likely be back to this a few project nights from now :slightly_smiling_face: