Arc barycentre

[andLaing]

Fixes a bug in makelors that meant barycentres of LOR positions could have radii outside the xenon.

[jacg]

I've been holding off merging the units of measure work into master, but I've decided to go for it. So I rebased your work on top of the new master which contains the units, and adapted your code so that it works with the units. Then I refactored a few more things.

Please check that I haven't messed up anything or misinterpreted the meaning of what you wrote.

[andLaing]

Looks good, much more elegant too.

[andLaing]

All tests pass with a few compilation warnings in my local

[jacg]

I don't trust the parts relating to cc in the implementation: they use the x-projection of points to represent their angular position. So a non-linear transformation is used before taking a mean and then mapping back to the original space. This will not give us the mean in the original space.

I started writing some proptests to check this function more thoroughly, but it's taking quite some effort to write a decent strategy for generating sets of vertices. Once I manage that, I expect that counterexample that proptest presents us after shrinking, will look too much like the typical float-equality errors. So it's probably best to come up with a hand-picked counterexample that shows a blatant discrepancy.

[andLaing]

I was doubting on that too. I've pushed a fix using x and y to get the mean angle then projecting to the mean radius.

[jacg]

• Should we take this opportunity to think about the desired value of t in the barycentre?
• Should we test this more thoroughly, perhaps using proptest, or is it not worth it, given that this reconstruction algorithm is only a toy one?
• Tangentially ... sooner or later we should port the real LOR reconstruction algorithm from julia to here. Can you give me some idea of how hard that will be?

[andLaing]

Should we take this opportunity to think about the desired value of t in the barycentre?

I mean, it's maybe better as a simple mean or a mean of the earliest hits/first n ns. I don't think it's that difficult to implement, just need to decide which is appropriate for this particular toy example.

Should we test this more thoroughly, perhaps using proptest, or is it not worth it, given that this reconstruction algorithm is only a toy one?

I think it's not worth much hassle at the moment. It's just a way to get LORs for testing that're a bit more realistic, right?

Tangentially ... sooner or later we should port the real LOR reconstruction algorithm from julia to here. Can you give me some idea of how hard that will be?

Depends at what level. If you want just the part that takes the selected events and outputs in the LOR format I think it's really easy to do. If we want the whole selection train, which is by no means 100% defined, it's a bit more complicated but we're not doing anything really complicated. The most delicate part would probably be moving to an event loop and then we could implement something simple that does the basics. We have a lot of extra parameters at the moment which are useful for understanding and defining algorithms and should be checked again once we're sure on the simulation and want to do a proper implementation of the event mixing but they aren't necessary in a first rust implementation since we'd likely continue using Julia for initial studies.

[jacg]

I don't think it's that difficult to implement, just need to decide which is appropriate for this particular toy example.

Oh, absolutely.

If you can think of something that is clearly more appropriate than what we have at present, let's put it in before merging, lest we forget. I'm sure you have far better intuition or insight here, than I do.

Should we test this more thoroughly, perhaps using proptest, or is it not worth it, given that this reconstruction algorithm is only a toy one?

I think it's not worth much hassle at the moment. It's just a way to get LORs for testing that're a bit more realistic, right?

Indeed. But the diminished-radius bug wasted your time.

I'd be tempted to add at least one test that does checks that the weights (dEs) are used vaguely sensibly.

Tangentially ... sooner or later we should port the real LOR reconstruction algorithm from julia to here. Can you give me some idea of how hard that will be?

Depends at what level. If you want just the part that takes the selected events and outputs in the LOR format I think it's really easy to do. If we want the whole selection train, which is by no means 100% defined, it's a bit more complicated but we're not doing anything really complicated. The most delicate part would probably be moving to an event loop and then we could implement something simple that does the basics. We have a lot of extra parameters at the moment which are useful for understanding and defining algorithms and should be checked again once we're sure on the simulation and want to do a proper implementation of the event mixing but they aren't necessary in a first rust implementation since we'd likely continue using Julia for initial studies.

Julia is almost certainly better for experimenting (at least until someone on the team puts some effort into identifying useful crates and processes in Rust, and even then ...), while Rust is almost certainly better for executing the code, and having it remain stable and reliable over time, and having it fit in well to the rest of our ecosystem.

I think we should port stuff which is unlikely to change much, over to Rust, keep the stuff which is still highly experimental in Julia, and gradually port more stuff as our understanding crystallizes.

[andLaing]

If you can think of something that is clearly more appropriate than what we have at present, let's put it in before merging, lest we forget. I'm sure you have far better intuition or insight here, than I do.

Ok, I'll implement something a bit more like what we'd do with reco data.

I'd be tempted to add at least one test that does checks that the weights (dEs) are used vaguely sensibly.

Seems reasonable, I'll try and dedicate a bit of time to it.

I think we should port stuff which is unlikely to change much, over to Rust, keep the stuff which is still highly experimental in Julia, and gradually port more stuff as our understanding crystallizes.

Right, since this isn't something that's really urgent, maybe I can start developing the functions relevant for this gradually. Means I can advance a bit more on RUST understanding (hopefully) and get things moving while we're busy with image reco.

[jacg]

I think we should port stuff which is unlikely to change much, over to Rust, keep the stuff which is still highly experimental in Julia, and gradually port more stuff as our understanding crystallizes.

Right, since this isn't something that's really urgent, maybe I can start developing the functions relevant for this gradually. Means I can advance a bit more on RUST understanding (hopefully) and get things moving while we're busy with image reco.

Not urgent at all, compared to other things. Also, I'd like to be involved (or do it all), as that would ensure that I understand what is going on in that area, which is mostly a complete mystery to me at the moment.

[andLaing]

I noticed a bug in the separation of the two positions in the code, here: https://github.com/jacg/petalo/blob/d36a1ae8d15a5ee43f80b18d5c0e93b835741a9a/src/bin/makelor.rs#L183

We just select using v.volume_id == 0 and separate using v.track_id == 1, there are non-negligible events where volume_id == 0 for vertices associated to other tracks. I think this could be patched in a couple of ways, neither of which is ideal: allowing v.track_id == 1 || v.parent_id == 1 then filtering the track_id == 2 part or only considering vertices directly from the two main tracks v.track_id == 1 then filtering all but track_id == 2 from the b.

[jacg]

there are non-negligible events where volume_id == 0 for vertices associated to other tracks

Hmm, how do they get there? They all appear to be phot or Rayl, not a single compt. And tend to have energies just over 29 keV ... that's some well-known process that generates X-rays in the LXe, isn't it?

Why isn't partitioning on v.track_id == 1 || v.parent_id == 1 sufficient? Because there is a small number with parent_id > 2. Most of these parents are 3 or 4, some more in single-figures and a few in the hundreds and thousands. At this point it becomes a real pain to work out whether they descended from gamma 1 or 2, perhaps even impossible with the incomplete data we store in the MC output files. But these are toy algorithms used for sanity checks, so it's probably not worth the effort.

So we should probably just ignore those that we can't trivially link to the original pair of gammas. How about

• filter on parent_id <= 2
• partition on track_id == 1 || parent_id == 1

?

[andLaing]

Hmm, how do they get there? They all appear to be phot or Rayl, not a single compt. And tend to have energies just over 29 keV ... that's some well-known process that generates X-rays in the LXe, isn't it?

Yeah K-alpha and K-beta lines are there abouts.

So we should probably just ignore those that we can't trivially link to the original pair of gammas. How about

filter on parent_id <= 2
partition on track_id == 1 || parent_id == 1

Reasonable compromise.

[jacg]

LGTM.

Shall I merge, or is there anything else you want to do here?

[andLaing]

think that's all we wanted for now.