Icrp107 database

[BishopWolf]

This PR addresses the issue #590

Basically replace radar database by the ICRP107 database

[BishopWolf]

@tbaudier This PR in principle does not change the behaviour of the tests, but the results may change a little bit as a more refined database is used.

Regarding the json, I tried to push it minimized (11MB), but the formatting bot expanded it (41MB).

[BishopWolf]

@tbaudier the precommit bot is doing fancy things. I think I have succeeded this time.

[dsarrut]

very good idea, thanks @BishopWolf ! @alexis-pereda can you have a look please to check this is ok for you ?

[BishopWolf]

@dsarrut Everything is working now. I will make a final push with the database schema.

[BishopWolf]

@alexis-pereda The change in the database makes some test yield a slightly different results, however those should be very close to before. The difference in time it's because the database is bigger, so once it's loaded it should work the same for larger simulations. If some of these presets are not respected, then there is a bug somewhere.

A solution could be to separate the icrp107.json into components. An icrp107 folder containing 1252 JSON, one for each radionuclide. This way we don't need to read the entire database. Let me know your thoughts. We can off course optimize it.

[BishopWolf]

@dsarrut @alexis-pereda I just discovered that icrp107 DO have beta+ spectras. It is just located in the same b-spectra field. therefore you don't even need a separate database for this.

Take as example the F18.

[BishopWolf]

Here it is the database separated into components, in binary format to prevent the issues with git transfer
icrp107.zip

you can read it like

isotopePath = Path("./icrp107") / "Lu-177.json"
with open(isotopePath, "rb") as i:
    isot = json.loads(json.load(i))
    print(isot["name"], isot["emissions"])

[BishopWolf]

The comparison of icrp107 vs X. Mougeot for F-18 beta+ spectra:

icrp107 mean energy 249.776 keV yield 0.9673
Mougeot mean energy 250.50 (21) keV

spectra icrp107 (better sampling at lower energies, worse sampling at higher energies)

I think, being Mougeot work more recent it shall be the default and we can rely on icrp107 for the missing isotopes.

isotopePath = Path("./icrp107") / "F-18.json"
with open(isotopePath, "rb") as i:
    isot = json.loads(json.load(i))
    print(isot["emissions"]["beta+"])
    for line in isot["emissions"]["b-spectra"]:
        print(f"| {line[0]*1000:.2f} | {line[1]/1000:.3e} |")

[alexis-pereda]

I looked at it. It is great to have everything and I think this will allow more flexibility the way its done.

However, it breaks the test010_generic_source_energy_spectrum_histogram.py.
When the source energy type is "spectrum_histogram", the "energy_bin_edges" option must be used instead of "energies".
This must contain energy bin edges instead of directly energies, so it will need one more value than the number of weights.

Also, after discussing it with other people, we think we should rename the function get_icrp107_spectrum to simply get_spectrum instead, but with an extra argument (with default value "icrp107") so that it is more generic and we can possibly add other sources for this kind of data.

[BishopWolf]

I have made some tests, separating by isotope seems to make everything faster. If you agree, I will push this change, and then we can work from there.

Basically, we need a function to calculate the energy bins from the energy slots, or maybe drop the last weight (it's always zero) and consider the energies as bins

[BishopWolf]

so it will need one more value than the number of weights.

@alexis-pereda So if I understand correctly, I must remove the last zero weight if the spectrum type is "b-spectra". Please confirm, this is a very trivial modification.

[BishopWolf]

@alexis-pereda Would somethign like this works ? Is it enough? Do you need spectrum_energies in histogram spectrum ?

def set_source_icrp107_energy_spectrum(source, rad):
    if (
        source.particle == "beta-"
        or source.particle == "e-"
        or source.particle == "beta+"
        or source.particle == "e+"
    ):
        rad_spectrum = get_icrp107_spectrum(rad, "b-spectra")
        source.energy.type = "spectrum_histogram"
        source.energy.spectrum_weights = rad_spectrum.weights[:-1]
        source.energy.spectrum_energy_bin_edges= rad_spectrum.energies
    else:
        rad_spectrum = get_icrp107_spectrum(rad, source.particle)
        source.energy.type = "spectrum_discrete"
        source.energy.spectrum_weights = rad_spectrum.weights
        source.energy.spectrum_energies = rad_spectrum.energies

[alexis-pereda]

If your rad_spectrum.energies in the histogram case is correctly providing energy bin edges and not energies, yes.
But I would rather have it named differently as it actually contains different things, even though it is an internal part of the code.

[alexis-pereda]

About the trivial "remove the last 0-weight", I tested it with the test010_generic_source_energy_spectrum_histogram.py test, and it still fails.
I am still trying to figure out why because the test expects to get the same values in the root output as it uses for the source.
(so it does not depend on the old data)

[BishopWolf]

I obtained 4x speed increase in the tests using individual json files compared to a monolithic database. So I have pushed that. It is also easier to debug and maintain.

[BishopWolf]

I am still trying to figure out why because the test expects to get the same values in the root output as it uses for the source.
(so it does not depend on the old data)

The test shall replace this function to use the proper initialization

def add_source_energy_spectrum_histogram(sim, phsp, interpolation: str = None):

    source = sim.add_source("GenericSource", "beam")
    source.attached_to = phsp.name
    source.particle = "e-"
    source.n = 5e5 / sim.number_of_threads
    source.position.type = "point"
    source.direction.type = "iso"

    spectrum = set_source_icrp107_energy_spectrum(source, "Lu177") # After defining the particle
    source.energy.spectrum_histogram_interpolation = interpolation

    return source

[alexis-pereda]

This is not the issue. I did a equivalent change in the test (and to be extra sure, I also tried with this function).
The test runs, it produces a root output, but we get bad results:

(you can replicate and find this output graph in opengate/tests/output/test010/test010_plot_energy_spectrum_histogram_None.png alongside the produced root file)

[BishopWolf]

The problem is this

The old test was made with gammas!!

[alexis-pereda]

I already fixed this and used "b-spectra" on my local tests.
Do you get correct results with the test on your side?

[BishopWolf]

Do you get correct results with the test on your side?

The old test is wrong, it uses gammas with the energy spectra of electrons. Everything else works fine. I will push the latests changes we have discussed.

[alexis-pereda]

Okay, I'll look again with your fix.

(and just to extend a bit on why I am surprised at the tests not passing:
the error in the test in master branch is that it emits gamma particles, but it correctly uses a beta spectrum.
So, while it has to be fixed for users that would use it as a reference, the tests still passes as it only looks at how the energy values are distributed.
Given that, because it only specifically try to match output energies to input energies, I fail to see how it could depend on what exact data is given to the source (as it is the source data that is used in the testing part as "input").
But again, I will see with your update, perhaps I am just missing something obvious here)

[BishopWolf]

I see the test 010 passing in all systems now. It remains the test 073 and this is the only thing different:

WARNING Compare stats
Runs:         4 4 --> OK? True
Events:       10620854 10146490 : +4.68 %  (tol = 1.00 %) --> OK? False
Tracks:       18329508 17500862 : +4.73 %  (tol = 1.00 %) --> OK? False
Steps:        52593098 50195644 : +4.78 %  (tol = 1.00 %) --> OK? False
PPS:          80573.0 589880.0 : -86.3%    speedup = x0.1 --> OK? True
TPS:          139053.0 1017437.0 : -86.3%    speedup = x0.1 --> OK? True
SPS:          398989.0 2918196.0 : -86.3%    speedup = x0.1 --> OK? True
Track e-      44004 42619 : +3.2% --> OK? False
Track gamma   18285504 17458243 : +4.7% --> OK? False
Tracks      : {'e-': 44004, 'gamma': 18285504} --> OK? True
Tracks (ref): {'e-': '42619', 'gamma': '17458243'} --> OK? True
Tracks vs track_types : 18329508 18329508 --> OK? True

In this case it is indeed a database change issue, not much can be done. The difference is not that mch as expected.

I can't explain the slow down as the comparison in both databases yields basically the same data

radar

"Lu177": {
    "energies": [
      0.0716418,
      0.1129498,
      0.1367245,
      0.2083662,
      0.2496742,
      0.3213159
    ],
    "weights": [
      0.001726,
      0.0620,
      0.000470,
      0.1038,
      0.002012,
      0.00216
    ]
  }

icrp107

energy (MeV)	weight
0.071646	0.001540
0.112950	0.064000
0.136725	0.000480
0.208366	0.110000
0.249674	0.002120
0.321316	0.002190

[BishopWolf]

Test 046 fails on these:

Test yield Tc99m: 0.8907654364665489 0.885 False --> OK? False
Test yield Lu177: 0.18033000000000002 0.172168 False --> OK? False
Test yield In111: 1.8474 1.847315 False --> OK? False
Test yield I131: 1.0077017706999998 1.0024600000000004 False --> OK? False

It really needs a tolerance, the difference is just the data provenance.

[BishopWolf]

test 046 is no longer failing

[BishopWolf]

Also, after discussing it with other people, we think we should rename the function get_icrp107_spectrum to simply get_spectrum instead, but with an extra argument (with default value "icrp107") so that it is more generic and we can possibly add other sources for this kind of data.

@alexis-pereda I have pushed this change too.

[BishopWolf]

@alexis-pereda test 73 discovery Lu177 will always fail since the icrp107 library has 4.74% extra gamma yield for Lu177. This IS a database issue, not a test issue. Therefore, either you rerun the test with the new database or we switch back to radar database. I mean set explicitly

set_source_energy_spectrum(source, "Lu177", "radar")

[alexis-pereda]

I think we have to ask @dsarrut for this, as he wrote the test.

@dsarrut is it possible for you to recreate the reference root file for the test test073_test5_discovery_lu177_mt.py
or do you think it is better to keep providing the RADAR source and continue to use it for this test?

[BishopWolf]

@dsarrut @alexis-pereda the tests test073_test2_intevo_mt and test073_test4_intevo_lu177_mt also have this issue. We can leave test2 using radar, but I think it will be good if both test5 are updated to use icrp107 instead. You can load the root files directly from the test simulations.

[BishopWolf]

@alexis-pereda we need a new rebase

[BishopWolf]

I have created an issue in radioactivedacay package to investigate the failing test053

[BishopWolf]

@alexis-pereda @dsarrut the package radioactivedecay has been updated with a fix. We must use version 0.6.1

[BishopWolf]

It remains this error in MacOS with python 3.12, I have created this issue

File "/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/opengate/utility.py", line 9, in <module>
    import pkg_resources
ModuleNotFoundError: No module named 'pkg_resources'

Apparently, the setuptools package is not properly installed in this case. Read this thread

[2023 update] 
The relevant library is removed from Python 3.12 and so a lot more people will hit this answer and not know why. 
The reason is that the latest version of python does not install this by default anymore and you're expected to switch to import_lib