Root data energy information is converted to keV instead of MeV

[robbietuk]

Addresses #356. @ludovicabrusaferri I cannot figure out how you handled this issue in #722.

This is done when properties of the list mode event are queeried to be between upper/low energy window (keV).

Edit: Our work around by setting upper/low energy window (keV) to an MeV value (e.g. 0.65) no longer works. This energy information is needed for scatter estimation code. Using lm_to_projdata to unlist a root file, the upper/low energy window (keV) information in the output sinogram is set from the hroot file, overwriting that in the template sinogram.

[KrisThielemans]

could you put this on release_4 please?

[robbietuk]

We are assuming the root data is in MeV. I cannot find a way to convert the energy information output by GATE to keV. However, other software may write the energy information in keV, which would cause this to fail... Additionally, (at least) my root files do not contain information regarding the energy unit.

One final silly suggestion would be to check the values of the energy information, for both thresholds and actual values, and if less than 1, we assume the energy is in MeV, and greater than 1, the energy information is keV.

STIR should use keV as the energy unit, which includes .hroot and .hs files. To include this PR, we should document that we expect root file energy units to be MeV. I cannot see and easy way to standardise everything to the same unit.

[ludovicabrusaferri]

Hi, yes the root data is in MeV.
Also, the current STIR code sets the threshold from the root header file (in MeV) and not from the scanner template info.
This is all a bit tricky.

With this PR you did, if you try to access the energy info (for instance during the unlisting), that will still be in MeV, as it's been only multiplied during the thresholding but it's not set to keV. This could all be fine in the case we don't need to access the energy info anymore i guess? not sure

I have done this on my branch : https://github.com/ludovicabrusaferri/STIR/blob/5470dcf08d0f1b584e2764795e21a0d3a9f4de5f/src/include/stir/listmode/CListRecordROOT.h#L137-L153.

If this can be of any help, I am also setting the energy info from the scanner template and not from root header:
https://github.com/ludovicabrusaferri/STIR/blob/5470dcf08d0f1b584e2764795e21a0d3a9f4de5f/src/listmode_buildblock/LmToProjData.cxx#L575-L582

Keep me posted! I am happy to help

[robbietuk]

Thank you for the input @ludovicabrusaferri

https://github.com/ludovicabrusaferri/STIR/blob/5470dcf08d0f1b584e2764795e21a0d3a9f4de5f/src/include/stir/listmode/CListRecordROOT.h#L137-L153.

Im sorry but I cant work out why you create energyA and energyB in #722, it is a copy of energy1 and energy2?

With this PR you did, if you try to access the energy info (for instance during the unlisting), that will still be in MeV, as it's been only multiplied during the thresholding but it's not set to keV. This could all be fine in the case we don't need to access the energy info anymore i guess? not sure

As far as I know, the energy information, from list mode, is only ever accessed in these window checks (in the master). Inspired by your code, perhaps we could set the energy1 and energy2 to private variables, ranaming to energy1_MeV/energy2_MeV and instead create

    float get_energy1_keV() const
    { return energy1_MeV * 1000; };
    float get_energy2_keV() const
    { return energy2_MeV * 1000; };

in https://github.com/UCL/STIR/blob/master/src/include/stir/IO/InputStreamFromROOTFile.h. That way, everything is clean. We can then convert your energy window checks to a standardised keV units? What I don't know is the this suggestion's impact on computation, i.e. the multiplcation or the call of a function multiple times for each list mode event.

I am also setting the energy info from the scanner template and not from root header:

I initially thought we could do this, but I think this is dangerous in the long run as lm_to_projdata is not the only list mode utility that uses energy window information, e.g. list_lm_events example.hroot. As the hroot does not have sinogram information, it relys on the energy information window information from the hroot. To get the energy information from two different sources for different utilities seems like it will cause problems.

[KrisThielemans]

I certainly agree that any STIR class should only return numbers in keV. If ROOT stores in MeV, we convert internally before we pass it on (a multiplication with 1000 isn't going to kill us, compared to all the rest!)

For @ludovicabrusaferri 's work, there are 2 things:

• reading energy info from listmode (i.e. correct interpretation of ROOT or whatever files)
• when unlisting, select which events you want. (i.e. use energy window info from the template for that selection).

However, for listmode recon, you either have some kind of filtering of events, or you need to write projectors that are energy aware (which we don't have). At present, the ROOT code allows you to do that filtering, so we cannot remove it until we put something more general in the listmode recon (similarly, Nikos created some things in the TOF branch to specify the TOF window width in the .hroot, such that you can simulate different DACs.)

Another somewhat tricky part is that in contrast to ROOT, in most scanners, listmode wouldn't give you an energy-pair, but an energy window pair index, where the energy window info sits in some header. (i.e. energy info is discretised, such as TOF or indeed ring/crystal info is discretised). So a CListEnergy probably should be able to return indices as well (in addition to the current get_energyA_in_keV(). But that's a discussion for #722, not for here.

Anyway, after all this, I believe that the current PR is good (except that I'd like to see it on release_4 but if it's too much of a hassle, we just list it as a known-bug in release_4.1.htm (on the release_4 branch). Things from #722 seem independent of this.

(I certainly think that breaking #722 in smaller PRs would be a good strategy. Depends if it can be done without having interdependencies between PRs.)

[robbietuk]

I think it is on release_4

[KrisThielemans]

Looks good to me (aside from 1 line in the release_4.1.htm).

@ludovicabrusaferri are you ok with this?

[ludovicabrusaferri]

@KrisThielemans Sounds good to me. I have energyA and energyB because I kept the same strategy used for time in https://github.com/ludovicabrusaferri/STIR/blob/5470dcf08d0f1b584e2764795e21a0d3a9f4de5f/src/include/stir/listmode/CListRecordROOT.h#L86 (https://github.com/ludovicabrusaferri/STIR/blob/5470dcf08d0f1b584e2764795e21a0d3a9f4de5f/src/include/stir/listmode/CListRecordROOT.h#L139). It made sense in my mind to create a copy of the time1/time2 or energy1/energy2 but it's probably unneeded?

[KrisThielemans]

@ludovicabrusaferri not 100% sure really how @NikEfth set-up the ROOT parsing (I knew, but forgot!). It's a special type of handling of listmode data, as other CListRecord classes just get whatever is stored in the record, while for ROOT, @NikEfth copied the data into class members. Both work, not 100% sure what is fastest (likely application dependent). In any case, it'd make sense for you to copy energy info as well then.

I did wonder why your chose A and B, as in other places we use 1 and 2, e.g. here but maybe we're not consistent...

anyway these comments should really sit in #722 I guess 😉