Create new datatype when saving CartesianPoint to HDF5 using LegendHDF5IO

[claudiaalvgar]

This PR allows for CartesianPoint and CylindricalPoint with and without units from #489 to be passed to Event and NBodyChargeCloud ( issue in #183 ). Test code and modes allowed:

# Passing 1 point in Cartesian or Cylindrical
pt = CylindricalPoint(2, deg2rad(20), 1)
e = 10u"keV"
evt1 = Event(pt)
evt2 = Event(pt, e)
println("1-Cyl no units: \n$(evt1), \n$(evt2)\n\n")

pt = CartesianPoint(1, 1.3, 1.2)
evt1 = Event(pt)
evt2 = Event(pt, e)
println("2-Cart no units: \n$(evt1), \n$(evt2)\n\n")

T=Float16
pt = CartesianPoint{T}(1, 1.3, 1.2)
evt1 = Event(pt)
evt2 = Event(pt, e)
println("3-Cart no units, Float16: \n$(evt1), \n$(evt2)\n\n")

pt = CylindricalPoint{T}(1, 1.3, 1.2)
evt1 = Event(pt)
evt2 = Event(pt, e)
println("3.1-Cyl no units, Float16: \n$(evt1), \n$(evt2)\n\n")

T=Float32
pt = CartesianPoint{T}(1, 1.3, 1.2)
evt1 = Event(pt)
evt2 = Event(pt, e)
println("4-Cart no units, Float32: \n$(evt1), \n$(evt2)\n\n")

pt = CylindricalPoint{T}(1, 1.3, 1.2)
evt1 = Event(pt)
evt2 = Event(pt, e)
println("4.1-Cyl no units, Float32: \n$(evt1), \n$(evt2)\n\n")


# Passing 1 point in Cartesian or Cylindrical with units
pt = CartesianPoint(1u"mm", 2u"mm", 3u"mm")
evt1 = Event(pt)
evt2 = Event(pt, e)
println("5-Cart units: \n$(evt1), \n$(evt2)\n\n")

pt = CylindricalPoint(5u"mm", 9u"rad", 2u"mm")
evt1 = Event(pt)
evt2 = Event(pt, e)
println("6-Cyl units: \n$(evt1), \n$(evt2)\n\n")

# Passing 1 vector in Cartesian or Cylindrical
pt1 =  CartesianPoint(1.0, 2.0, 3.0)
pt2 =  CartesianPoint(4.0, 5.0, 6.0)
energy1 = 1.0u"MeV"
energy2 = 2.0u"MeV"
evt1 = Event([pt1,pt2])
evt2 = Event([pt1,pt2], [energy1, energy2])
println("7-Cart no units: \n$(evt1), \n$(evt2)\n\n")

pt1 = CylindricalPoint(5, 9, 2)
pt2 = CylindricalPoint(1, 1.0, 0.0)
evt1 = Event([pt1,pt2])
evt2 = Event([pt1,pt2], [energy1, energy2])
println("8-Cyl no units: \n$(evt1), \n$(evt2)\n\n")

# Passing 1 vector in Cartesian or Cylindrical with units
pt1 =  CartesianPoint(1.0u"mm", 2.0u"mm", 3.0u"mm")
pt2 =  CartesianPoint(4.0u"mm", 5.0u"mm", 6.0u"mm")
evt1 = Event([pt1,pt2])
evt2 = Event([pt1,pt2], [energy1, energy2])
println("9-Cart with units: \n$(evt1), \n$(evt2)\n\n")

pt1 = CylindricalPoint(5u"mm", 9u"rad", 2u"mm")
pt2 = CylindricalPoint(1u"mm", 1.0u"rad", 0.0u"mm")
evt1 = Event([pt1,pt2])
evt2 = Event([pt1,pt2], [energy1, energy2])
println("10-Cyl with units: \n$(evt1), \n$(evt2)\n\n")

# Passing 1 vector of vectors in Cartesian or Cylindrical
pt1 = CartesianPoint(1.0, 2.0, 3.0)
pt2 = CartesianPoint(4.0, 5.0, 6.0)
pt3 = CartesianPoint(7.0, 8.0, 9.0)
pt4 = CartesianPoint(10.0, 11.0, 12.0)

locations = [
    [pt1, pt2],  # first "event"
    [pt3, pt4]   # second "event"
]
energy1 = 1.0u"MeV"
energy2 = 2.0u"MeV"
energy3 = 1.5u"MeV"
energy4 = 2.5u"MeV"

energies = [
    [energy1, energy2],
    [energy3, energy4]
]
evt1 = Event(locations)
evt2 = Event(locations, energies)
println("11-Cart no units: \n$(evt1), \n$(evt2)\n\n")

pt1 = CylindricalPoint(1.0, 2.0, 3.0)
pt2 = CylindricalPoint(4.0, 5.0, 6.0)
pt3 = CylindricalPoint(7.0, 8.0, 9.0)
pt4 = CylindricalPoint(10.0, 11.0, 12.0)

locations = [
    [pt1, pt2],  # first "event"
    [pt3, pt4]   # second "event"
]

energies = [
    [energy1, energy2],
    [energy3, energy4]
]
evt1 = Event(locations)
evt2 = Event(locations, energies)
println("12-Cyl no units: \n$(evt1), \n$(evt2)\n\n")

# Passing 1 vector of vectors in Cartesian or Cylindrical with units
pt1 = CartesianPoint(1.0u"mm", 2.0u"mm", 3.0u"mm")
pt2 = CartesianPoint(4.0u"mm", 5.0u"mm", 6.0u"mm")
pt3 = CartesianPoint(7.0u"mm", 8.0u"mm", 9.0u"mm")
pt4 = CartesianPoint(10.0u"mm", 11.0u"mm", 12.0u"mm")

locations = [
    [pt1, pt2],  # first "event"
    [pt3, pt4]   # second "event"
]

energies = [
    [energy1, energy2],
    [energy3, energy4]
]
evt1 = Event(locations)
evt2 = Event(locations, energies)
println("13-Cart units: \n$(evt1), \n$(evt2)\n\n")

pt1 = CylindricalPoint(1.0u"mm", 2.0u"rad", 3.0u"mm")
pt2 = CylindricalPoint(4.0u"mm", 5.0u"rad", 6.0u"mm")
pt3 = CylindricalPoint(7.0u"mm", 8.0u"rad", 9.0u"mm")
pt4 = CylindricalPoint(10.0u"mm", 11.0u"rad", 12.0u"mm")

locations = [
    [pt1, pt2],  # first "event"
    [pt3, pt4]   # second "event"
]

energies = [
    [energy1, energy2],
    [energy3, energy4]
]
evt1 = Event(locations)
evt2 = Event(locations, energies)
println("14-Cyl units: \n$(evt1), \n$(evt2)\n\n")

#NBodyChargeCloud point no units
center = CartesianPoint(0.2,0.1,0.5)
energy = 1460u"keV"
nbcc = NBodyChargeCloud(center, energy)
evt = Event(nbcc)
println("15-Cart no units: \n$(evt)\n\n")

center = CartesianPoint(0.2,0.1,0.5)
energy = 1460u"keV"
nbcc = NBodyChargeCloud(center, energy, 100)
evt = Event(nbcc)
println("16-Cart no units: \n$(evt)\n\n")

center = CylindricalPoint(0.2,0.1,0.5)
nbcc = NBodyChargeCloud(center, energy)
evt = Event(nbcc)
println("17-Cyl no units: \n$(evt)\n\n")

#NBodyChargeCloud point units
center = CartesianPoint(0.2u"mm",0.3u"mm",5.0u"mm")
nbcc = NBodyChargeCloud(center, energy)
evt = Event(nbcc)
println("18-Cart units: \n$(evt)\n\n")

center = CartesianPoint(0.2u"mm",0.3u"mm",5.0u"mm")
nbcc = NBodyChargeCloud(center, energy, 100)
evt = Event(nbcc)
println("19-Cart units: \n$(evt)\n\n")


center = CylindricalPoint(0.2u"mm",0.0u"rad",2.1u"mm")
nbcc = NBodyChargeCloud(center, energy)
evt = Event(nbcc)
println("20-Cyl units: \n$(evt)\n\n")

[fhagemann]

Thanks so much for taking care of this!!

I would like to have tests to ensure that creating events/charge clouds with units actually works as intended. It looks, for example, that the number of charges in an NBodyChargeCloud are not passed further down.

Can you create tests for each of the functions that you modified/added?

And unrelated, but we should have tests for EVERY function/method that involves CartesianPoint or CylindricalPoint to ensure that there are no bugs after having merged #489 !

[claudiaalvgar]

I implemented the suggestions and also added a fix to run_geant4_simulation so that writedata works with changes in #489 and also test for it in test_geant4.jl

[fhagemann]

With my latest commit be4ec17, we should now be able to write/read CartesianPoint and CylindricalPoint, as well as Vectors containing those points using LegendHDF5IO -- both in the old readdata/writedata way and with the new LH5Array.

Tests failing are unrelated, same issue as before (with the documentation and Geant4 tests expecting pos to have eltype CartesianPoint instead of SVector{3}).

[fhagemann]

I also reverted two of your commits to see if going back to CartesianPoint for the eltype in the Geant4 Table pos column fixes the tests again.

This should be fixed now after reverting the changes made to the Geant4-Table pos column

[fhagemann]

I pushed my requested changes.
The only thing that's still weird to me: we allow for the fields x, y, z of CartesianPoint to have units, whereas for r, φ, z of CylindricalPoint we strip them. This is probably because r/z have different units than φ.

Maybe the most consistent way to proceed would be to also strip the units of CartesianPoint, so that internally it is also unitless, but still allow for all operations with units from outside (e.g. adding a CartesianVector with units, etc.)

[claudiaalvgar]

Yes, in principle, before this PR CylindricalPoint with units wasn't defined. I just added it this way in this PR because sometimes, it could be more convenient to pass the cylindrical coordinate to the Event if your detector is in those coordinates but we can talk about the best way of doing this. Are there any other changes to add?

[oschulz]

Didn't we want to leave units out for now?

[fhagemann]

Didn't we want to leave units out for now?

That’s the question:
I would say, we should for now also strip the units from CartesianPoint, but still allow for the user to create them using units. The coordinates internally should still be floats for now, and let’s add support for units on the internal coordinates themselves in a future PR.

[oschulz]

Sounds good.

[fhagemann]

@claudiaalvgar: could you remove the units from the CartesianPoint coordinates everywhere, such that pt.x, pt.y and pt.z are always unitless numbers in units of meters ?
Then, you should be able to also run the commented new tests for CartesianPoint.

[claudiaalvgar]

Yes, I can add this but I still have the question if the code in Event.jl will stay as it is, if we should delete the changes to account for the units or if there are more changes to implement given this comment:
"

Side comment: I fear a bit that determining T here (and anywhere else) might create some trouble. If we require T of the Event to match T of the Simulation when drifting events, this might not be what we want..

simulate!(evt::Event{T}, sim::Simulation{T},..)

(the T need to match and I fear that this implementation might make it a bit unpredictable, especially because the function does not know the T of the output Event{T} before being called.
"

[fhagemann]

Yes, I can add this but I still have the question if the code in Event.jl will stay as it is,

I don't see why not. Many (all?) methods use to_internal_units which should work perfectly with unitless numbers, as long as those numbers are in meters.

if we should delete the changes to account for the units or

Ideally, all changes should still work, if we use values without units. I would keep the changes, but maybe restrict the type T of AbstractCoordinatePoint{T} to AbstractCoordinatePoint{T <: Real} or AbstractCoordinatePoint{T <: SSDFloat}.
We can then check, if the rest still works or not.

if there are more changes to implement given this comment:

If the event locations are of type T now, they should give the correct T for Event{T}.
We should add tests to ensure this though.

[claudiaalvgar]

My question is: Should this be deleted pt_internal = to_internal_units(location) T = eltype(pt_internal) ? because this is not needed if we don't support units and just do AbstractCoordinatePoint{T <: SSDFloat}

[fhagemann]

We can define @inline to_internal_units(x::AbstractCoordinatePoint{<:Real}) = x if that doesn't exist yet, because then we leave the current code as is, it shouldn't add extra runtime and should result in what we expect.

If we can leave the code as is, it will be easier to support points with eltypes having units in the future, because the majority of the required code modifications should already be there.

But maybe @oschulz has an opinion on this (?)

[claudiaalvgar]

I think I implemented the changes suggested. The modification in MCEventsProcessing.jl was the best I could think of to remove the units and match the output from SolidStateDetectors.simulate_waveforms in SolidStateDetectorsLegendHDF5IOExt.jl since this output is always given without units. Other methods I tried were failing but maybe you have an alternative suggestion for this(?). I also left the CartesianPoint struct as it was since doing this struct CartesianPoint{T <: Real} <: AbstractCartesianPoint{T} was giving the wrong output for CartesianPoint(1, 2, 3) i.e. CartesianPoint{Int64} instead of CartesianPoint{Float64} in CSG_coordinates.jl test