Integration with Unitful.jl

[oschulz]

Geant4.jl currently has it's own implementation of physical units in Geant4.SystemOfUnits. It would be nice to integrate Geant4 with Unitful.jl since it's pretty much the standard for units in the Julia ecosystem.

We could define

using Unitfil

const g4length = u"mm"
const g4area = g4length^2
const g4volume = g4length^3
const g4time = u"ns"
const g4angle = u"rad"
const g4solidangle = u"sr"
const g4frequency = inv(g4time)
const g4charge = u"q"
const g4energy = u"MeV"
const g4mass = g4energy*g4time*g4time/(g4length*g4length) # Simplify?
const g4density = g4mass/g4volume

# Same for: power, force, pressure, current, epot, resistance, capacitance, flux, magnetic_field,
# inductance, temperature, amount, activity, dose, intensity, luminous_flux, illuminance

and maybe also convenience functions

asg4length(x::Unitful.Quantity{<:Real}) = ustrip(g4_length, x)
asg4area(x::Unitful.Quantity{<:Real}) = ustrip(g4area, x)
asg4volume(x::Unitful.Quantity{<:Real}) = ustrip(g4volume, x)
asg4time(x::Unitful.Quantity{<:Real}) = ustrip(g4time, x)
asg4angle(x::Unitful.Quantity{<:Real}) = ustrip(g4angle, x)
asg4solidangle(x::Unitful.Quantity{<:Real}) = ustrip(g4solidangle, x)
asg4frequency(x::Unitful.Quantity{<:Real}) = ustrip(g4frequency, x)
asg4charge(x::Unitful.Quantity{<:Real}) = ustrip(g4charge, x)
asg4energy(x::Unitful.Quantity{<:Real}) = ustrip(g4energy, x)
asg4mass(x::Unitful.Quantity{<:Real}) = ustrip(g4mass, x)
asg4density(x::Unitful.Quantity{<:Real}) = ustrip(g4density, x)
# ...


asg4length(x::Real) = x
asg4area(x::Real) = x
asg4volume(x::Real) = x
asg4time(x::Real) = x
asg4angle(x::Real) = x
asg4solidangle(x::Real) = x
asg4frequency(x::Real) = x
asg4charge(x::Real) = x
asg4energy(x::Real) = x
asg4mass(x::Real) = x
asg4mass(x::Real) = x
# ...

That way, we could use generic code like

const G4RealQuantity{T<:Real} = Union{T,Unitful.Quantity{<:T}}

G4Material(name::String; #=...=#, density::G4RealQuantity, #=...=#) = #=...=# asg4density(density) #=...=#

in "G4JLInterface.jl". Users could then write either

using Unitful
G4Material("quartz", density=2.200*u"g/cm^3", ncomponents=2)

or

using Geant4.SystemOfUnits: g, cm
G4Material("quartz", density=2.200g/cm3, ncomponents=2)

or use unitless numbers directly

using Geant4.SystemOfUnits: g, cm
G4Material("quartz", density=1.3731319963813677e19, ncomponents=2)

Or should we try to replace Geant4.SystemOfUnits by Unitful (aliases) completely?

@giordano, you have a lot of experience with Unitful.jl and derived packages, do you think the design above goes in the right direction?

[giordano]

I haven't thought this through, but it looks a good solution to me.

[oschulz]

@peremato if this looks Ok to you in principle, should I create a PR and we'll iterate on it?

[peremato]

Hi @oschulz. I do not really understand how is going to work. Most of the G4 functions (99.9%) are direct calls to the CxxWrap of the C++ function. The function (constructor) G4Material is an exception to provide a user-friendly form to the various underling calls to C++. For example, to construct a G4Box you need to pass 3 Float64 expressed in mm. If you call it with G4Box("box", 10*u"cm", 10*u"cm", 10*u"cm") you get no method matching G4Box(::String, ::Quantity{Int64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, ::Quantity{Int64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, ::Quantity{Int64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}). Of course you can define G4Box(...,x::G4RealQuantity,...) = G4Box(...,asg4lenght(x),...) as you are suggesting. But, I am afraid that supporting Unitful` will imply to write a new wrapper to all methods of Geant4 (and this is a lot of work). Perhaps it could be automated if we could deduce the unit type for each argument.

[oschulz]

Most of the G4 functions (99.9%) are direct calls to the CxxWrap of the C++ function. [...] Unitful` will imply to write a new wrapper to all methods of Geant4 (and this is a lot of work).

Sorry @peremato, my example above wasn't great - we can obviously only provide the full "luxury" for the few explicitly wrapped (but also frequently used) functions like G4Material(...). Getting this into the auto-generated wrappers would indeed by very challenging.

But in all other functions, users can still write some_g4_functionality(..., asg4length(some_unitful_value), ...). Not quite as short, but very readable and "unit-safe" (esp. if some_unitful_value comes from elsewhere).

As a first step we could maybe add the constants like g4length and convenience functions like asg4length and then "upgrade" select functions like G4Material(...) later on?

[peremato]

I see. Perhaps you can just define a generic conversion of Unitful quantity to what is expected by Geant4. Something like this:

julia> G4unit(x::Unitful.Quantity{<:Real, Unitful.𝐋}) = Float64(ustrip(u"mm", x))
G4unit (generic function with 1 method)
julia> G4unit(x::Unitful.Quantity{<:Real, Unitful.𝐓}) = Float64(ustrip(u"ns", x))
G4unit (generic function with 2 methods)

julia> G4unit(10u"cm")
100.0

julia> G4unit(1u"s")
1.0e9

Calling Box("box", G4unit(1u"cm"), G4unit(1u"cm"), G4unit(1u"cm")) would do the work.

[oschulz]

Perhaps you can just define a generic conversion of Unitful quantity to what is expected by Geant4

I though about that, but then users would actually be allowed to pass a length into a volume argument if they're not careful. The idea behind separate asg4length(x) (and so on) function would be that users would express the dimensions at the "point of use", while the value may come from somewhere else (external configuration, etc.).

[peremato]

The only way the user will not be able to make a mistake or fake a call would be if we could annotate each call of Geant4 with what dimension is expected for each argument. But, this is not easily possible. I see the point of indicating derived dimensions with asg4volume, asg4density, ...
So, please go ahead to provide a PR.

[oschulz]

Thanks, will do!

[oschulz]

Still to-do, sorry ... :-)