diff --git a/.gitignore b/.gitignore index 318e444d..eb22feb4 100644 --- a/.gitignore +++ b/.gitignore @@ -2,5 +2,6 @@ *.jl.cov *.jl.*.cov *.jl.mem +*.lh5 .ipynb_checkpoints Manifest.toml diff --git a/docs/src/extensions.md b/docs/src/extensions.md index c2e797b9..2fe19470 100644 --- a/docs/src/extensions.md +++ b/docs/src/extensions.md @@ -36,7 +36,7 @@ dsp = read_ldata(l200, :jldsp, :cal, :p03, :r000, det) ``` ## `SolidStateDetectors` extension -LegendDataManagment provides an extension for [SolidStateDetectors](https://github.com/JuliaPhysics/SolidStateDetectors.jl). This makes it possible to create `SolidStateDetector` instances from LEGEND metadata. +LegendDataManagment provides an extension for [SolidStateDetectors](https://github.com/JuliaPhysics/SolidStateDetectors.jl). This makes it possible to create `SolidStateDetector` and `Simulation` instances from LEGEND metadata. Example (requires a `$LEGEND_DATA_CONFIG` environment variable pointing to a legend data-config file): @@ -52,6 +52,21 @@ A detector can also be constructed using the filename of the LEGEND metadata det det = SolidStateDetector(LegendData, "V99000A.json") ``` +In addition, when creating a `Simulation`, all simulation functions in SolidStateDetectors.jl can be applied. As usual, all fields stored in the `Simulation` can be written and read using `LegendHDF5IO`: + +```julia +using LegendDataManagement +using SolidStateDetectors + +sim = Simulation(LegendData, "V99000A.json") +simulate!(sim) # calculate electric field and weighting potentials + +using LegendHDF5IO +ssd_write("V99000A.lh5", sim) +sim_in = ssd_read("V99000A.lh5", Simulation) +``` + + The following code will generate an overview plot of every 5th LEGEND detector (requires the actual LEGEND metadata instead of the metadata in legend-testdata): ```julia diff --git a/ext/LegendDataManagementSolidStateDetectorsExt.jl b/ext/LegendDataManagementSolidStateDetectorsExt.jl index dd1b5306..500a9bc5 100644 --- a/ext/LegendDataManagementSolidStateDetectorsExt.jl +++ b/ext/LegendDataManagementSolidStateDetectorsExt.jl @@ -3,38 +3,33 @@ module LegendDataManagementSolidStateDetectorsExt using SolidStateDetectors -import SolidStateDetectors.ConstructiveSolidGeometry as CSG - using LegendDataManagement using Unitful using PropDicts - const _SSDDefaultNumtype = Float32 + """ - SolidStateDetector[{T<:Real}](data::LegendData, detector::DetectorIdLike - SolidStateDetector[{T<:Real}(::Type{LegendData}, detector_props::AbstractDict) - SolidStateDetector[{T<:Real}(::Type{LegendData}, json_filename::AbstractString) + SolidStateDetector[{T<:AbstractFloat}](data::LegendData, detector::DetectorIdLike) + SolidStateDetector[{T<:AbstractFloat}(::Type{LegendData}, detector_props::AbstractDict) + SolidStateDetector[{T<:AbstractFloat}(::Type{LegendData}, json_filename::AbstractString) LegendDataManagement provides an extension for SolidStateDetectors, a -`SolidStateDetector` can be constructed from LEGEND metadata using the +`SolidStateDetector` can be constructed from LEGEND metadata using the methods above. """ function SolidStateDetectors.SolidStateDetector(data::LegendData, detector::DetectorIdLike) SolidStateDetectors.SolidStateDetector{_SSDDefaultNumtype}(data, detector) end -function SolidStateDetectors.SolidStateDetector{T}(data::LegendData, detector::DetectorIdLike) where {T<:Real} +function SolidStateDetectors.SolidStateDetector{T}(data::LegendData, detector::DetectorIdLike) where {T<:AbstractFloat} detector_props = getproperty(data.metadata.hardware.detectors.germanium.diodes, Symbol(detector)) - SolidStateDetector{T}(LegendData, detector_props) + xtal_props = getproperty(data.metadata.hardware.detectors.germanium.crystals, Symbol(string(detector)[1:end-1])) + SolidStateDetector{T}(LegendData, detector_props, xtal_props) end - -to_SSD_units(::Type{T}, x, unit) where {T} = T(SolidStateDetectors.to_internal_units(x*unit)) - - -function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, filename::String) where {T<:Real} +function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, filename::String) where {T<:AbstractFloat} SolidStateDetector{T}(LegendData, readprops(filename, subst_pathvar = false, subst_env = false, trim_null = false)) end @@ -46,11 +41,58 @@ function SolidStateDetectors.SolidStateDetector(::Type{LegendData}, meta::Abstra SolidStateDetectors.SolidStateDetector{_SSDDefaultNumtype}(LegendData, meta) end -function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::AbstractDict) where {T<:Real} - SolidStateDetectors.SolidStateDetector{T}(LegendData, convert(PropDict, meta)) +function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::AbstractDict) where {T<:AbstractFloat} + SolidStateDetectors.SolidStateDetector{T}(LegendData, convert(PropDict, meta), LegendDataManagement.NoSuchPropsDBEntry("",[])) +end + +function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::PropDict, xtal_meta::Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}) where {T<:AbstractFloat} + config_dict = create_SSD_config_dict_from_LEGEND_metadata(meta, xtal_meta) + return SolidStateDetector{T}(config_dict, SolidStateDetectors.construct_units(config_dict)) +end + +""" + Simulation[{T<:AbstractFloat}](data::LegendData, detector::DetectorIdLike) + Simulation[{T<:AbstractFloat}(::Type{LegendData}, detector_props::AbstractDict) + Simulation[{T<:AbstractFloat}(::Type{LegendData}, json_filename::AbstractString) + +LegendDataManagement provides an extension for SolidStateDetectors, a +`Simulation` can be constructed from LEGEND metadata using the +methods above. +""" +function SolidStateDetectors.Simulation(data::LegendData, detector::DetectorIdLike) + SolidStateDetectors.Simulation{_SSDDefaultNumtype}(data, detector) +end + +function SolidStateDetectors.Simulation{T}(data::LegendData, detector::DetectorIdLike) where {T<:AbstractFloat} + detector_props = getproperty(data.metadata.hardware.detectors.germanium.diodes, Symbol(detector)) + xtal_props = getproperty(data.metadata.hardware.detectors.germanium.crystals, Symbol(string(detector)[1:end-1])) + Simulation{T}(LegendData, detector_props, xtal_props) +end + +function SolidStateDetectors.Simulation{T}(::Type{LegendData}, filename::String) where {T<:AbstractFloat} + Simulation{T}(LegendData, readprops(filename, subst_pathvar = false, subst_env = false, trim_null = false)) +end + +function SolidStateDetectors.Simulation(::Type{LegendData}, filename::String) + Simulation{_SSDDefaultNumtype}(LegendData, filename) +end + +function SolidStateDetectors.Simulation(::Type{LegendData}, meta::AbstractDict) + SolidStateDetectors.Simulation{_SSDDefaultNumtype}(LegendData, meta) end -function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::PropDict) where {T<:Real} +function SolidStateDetectors.Simulation{T}(::Type{LegendData}, meta::AbstractDict) where {T<:AbstractFloat} + SolidStateDetectors.Simulation{T}(LegendData, convert(PropDict, meta), LegendDataManagement.NoSuchPropsDBEntry("", [])) +end + +function SolidStateDetectors.Simulation{T}(::Type{LegendData}, meta::PropDict, xtal_meta::Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}) where {T<:AbstractFloat} + config_dict = create_SSD_config_dict_from_LEGEND_metadata(meta, xtal_meta) + return Simulation{T}(config_dict) +end + + +function create_SSD_config_dict_from_LEGEND_metadata(meta::PropDict, xtal_meta::X; dicttype = Dict{String,Any}) where {X <: Union{PropDict, LegendDataManagement.NoSuchPropsDBEntry}} + # Not all possible configurations are yet implemented! # https://github.com/legend-exp/legend-metadata/blob/main/hardware/detectors/detector-metadata_1.pdf # https://github.com/legend-exp/legend-metadata/blob/main/hardware/detectors/detector-metadata_2.pdf @@ -60,148 +102,224 @@ function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::Pro # https://github.com/legend-exp/legend-metadata/blob/main/hardware/detectors/detector-metadata_6.pdf # https://github.com/legend-exp/legend-metadata/blob/main/hardware/detectors/detector-metadata_7.pdf - gap = to_SSD_units(T, 1, u"mm") + gap = 1.0 dl_thickness_in_mm = :dl_thickness_in_mm in keys(meta.geometry) ? meta.geometry.dl_thickness_in_mm : 0 - li_thickness = to_SSD_units(T, dl_thickness_in_mm, u"mm") + li_thickness = dl_thickness_in_mm - crystal_radius = to_SSD_units(T, meta.geometry.radius_in_mm, u"mm") - crystal_height = to_SSD_units(T, meta.geometry.height_in_mm, u"mm") + crystal_radius = meta.geometry.radius_in_mm + crystal_height = meta.geometry.height_in_mm is_coax = meta.type == "coax" - # main crystal - semiconductor_geometry = CSG.Cone{T}(CSG.ClosedPrimitive; - r = crystal_radius, - hZ = crystal_height / 2, - origin = CartesianPoint{T}(0, 0, crystal_height / 2) + config_dict = dicttype( + "name" => meta.name, + "units" => dicttype( + "length" => "mm", + "potential" => "V", + "angle" => "deg", + "temperature" => "K" + ), + "grid" => dicttype( + "coordinates" => "cylindrical", + "axes" => dicttype( + "r" => dicttype( + "to" => crystal_radius * 1.2, + "boundaries" => "inf" + ), + "phi" => dicttype( + "from" => 0, + "to" => 0, + "boundaries" => "reflecting" + ), + "z" => dicttype( + "from" => -0.2 * crystal_height, + "to" => 1.2 * crystal_height, + "boundaries" => "inf" + ) + ) + ), + "medium" => "vacuum", + "detectors" => [] ) - # borehole - has_borehole = haskey(meta.geometry, :borehole) - if is_coax && !has_borehole - error("Coax detectors should have boreholes") - end - if has_borehole - borehole_depth = to_SSD_units(T, meta.geometry.borehole.depth_in_mm, u"mm") - borehole_radius = to_SSD_units(T, meta.geometry.borehole.radius_in_mm, u"mm") - semiconductor_geometry -= CSG.Cone{T}(CSG.ClosedPrimitive; - r = borehole_radius, - hZ = borehole_depth / 2 + gap, - origin = CartesianPoint{T}(0, 0, is_coax ? borehole_depth/2 - gap : crystal_height - borehole_depth/2 + gap) - ) - end + push!(config_dict["detectors"], dicttype( + "semiconductor" => dicttype( + "material" => "HPGe", + "charge_drift_model" => dicttype( + "include" => joinpath(SolidStateDetectors.get_path_to_example_config_files(), "ADLChargeDriftModel", "drift_velocity_config.yaml"), + ), + # "impurity_density" => dicttype("parameters" => Vector()), + "geometry" => dicttype(), + "temperature" => 78 + ), + "contacts" => [] + )) - # borehole taper - has_borehole_taper = haskey(meta.geometry.taper, :borehole) - if has_borehole_taper - borehole_taper_height = to_SSD_units(T, meta.geometry.taper.borehole.height_in_mm, u"mm") - if haskey(meta.geometry.taper.borehole, :radius_in_mm) - borehole_taper_radius = to_SSD_units(T, meta.geometry.taper.borehole.radius_in_mm, u"mm") - borehole_taper_angle = atan(borehole_taper_radius, borehole_taper_height) - elseif haskey(meta.geometry.taper.borehole, :angle_in_deg) - borehole_taper_angle = to_SSD_units(T, meta.geometry.taper.borehole.angle_in_deg, u"°") - borehole_taper_radius = borehole_taper_height * tan(borehole_taper_angle) - else - error("The borehole taper needs either radius_in_mm or angle_in_deg") - end - has_borehole_taper = borehole_taper_height > 0 && borehole_taper_angle > 0 - if has_borehole_taper && !has_borehole - error("A detector without a borehole cannot have a borehole taper.") + # main crystal + semiconductor_geometry_basis = dicttype("cone" => dicttype( + "r" => crystal_radius, + "h" => crystal_height, + "origin" => [0, 0, crystal_height / 2] + )) + semiconductor_geometry_subtractions = [] + begin + # borehole + has_borehole = haskey(meta.geometry, :borehole) + if is_coax && !has_borehole + error("Coax detectors should have boreholes") end - if has_borehole_taper && is_coax - error("Coax detectors should not have borehole tapers") + if has_borehole + borehole_depth = meta.geometry.borehole.depth_in_mm + borehole_radius = meta.geometry.borehole.radius_in_mm + push!(semiconductor_geometry_subtractions, dicttype("cone" => dicttype( + "r" => borehole_radius, + "h" => borehole_depth + 2*gap, + "origin" => [0, 0, is_coax ? borehole_depth/2 - gap : crystal_height - borehole_depth/2 + gap] + ))) end + + # borehole taper + has_borehole_taper = haskey(meta.geometry.taper, :borehole) if has_borehole_taper - r_center = borehole_radius + borehole_taper_height * tan(borehole_taper_angle) / 2 - hZ = borehole_taper_height/2 - Δr = hZ * tan(borehole_taper_angle) - r_out_bot = r_center - Δr - r_out_top = r_center + Δr * (1 + 2*gap/hZ) - r = ((r_out_bot,), (r_out_top,)) - semiconductor_geometry -= CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ + gap, - origin = CartesianPoint{T}(0, 0, crystal_height - borehole_taper_height/2 + gap) - ) + borehole_taper_height = meta.geometry.taper.borehole.height_in_mm + if haskey(meta.geometry.taper.borehole, :radius_in_mm) + borehole_taper_radius = meta.geometry.taper.borehole.radius_in_mm + borehole_taper_angle = atand(borehole_taper_radius, borehole_taper_height) + elseif haskey(meta.geometry.taper.borehole, :angle_in_deg) + borehole_taper_angle = meta.geometry.taper.borehole.angle_in_deg + borehole_taper_radius = borehole_taper_height * tand(borehole_taper_angle) + else + error("The borehole taper needs either radius_in_mm or angle_in_deg") + end + has_borehole_taper = borehole_taper_height > 0 && borehole_taper_angle > 0 + if has_borehole_taper && !has_borehole + error("A detector without a borehole cannot have a borehole taper.") + end + if has_borehole_taper && is_coax + error("Coax detectors should not have borehole tapers") + end + if has_borehole_taper + r_center = borehole_radius + borehole_taper_radius / 2 + hZ = borehole_taper_height/2 + Δr = hZ * tand(borehole_taper_angle) + r_out_bot = r_center - Δr + r_out_top = r_center + Δr * (1 + 2*gap/hZ) + push!(semiconductor_geometry_subtractions, dicttype("cone" => dicttype( + "r" => dicttype( + "from" => r_out_bot, + "to" => r_out_top + ), + "h" => 2 * hZ, + "origin" => [0, 0, crystal_height - borehole_taper_height / 2 + gap] + ))) + end + end + + # top taper + if haskey(meta.geometry.taper, :top) + top_taper_height = meta.geometry.taper.top.height_in_mm + if haskey(meta.geometry.taper.top, :radius_in_mm) + top_taper_radius = meta.geometry.taper.top.radius_in_mm + top_taper_angle = atand(top_taper_radius, top_taper_height) + elseif haskey(meta.geometry.taper.top, :angle_in_deg) + top_taper_angle = meta.geometry.taper.top.angle_in_deg + top_taper_radius = top_taper_height * tand(top_taper_angle) + else + error("The top taper needs either radius_in_mm or angle_in_deg") + end + has_top_taper = top_taper_height > 0 && top_taper_angle > 0 + if has_top_taper + r_center = crystal_radius - top_taper_radius / 2 + hZ = top_taper_height/2 + 1gap + h = 2 * hZ + Δr = hZ * tand(top_taper_angle) + r_in_bot = r_center + Δr + r_in_top = r_center - Δr + r_out = max(r_in_top, r_in_bot) + gap # ensure that r_out is always bigger as r_in + push!(semiconductor_geometry_subtractions, dicttype("cone" => dicttype( + "r" => dicttype( + "bottom" => dicttype( + "from" => r_in_bot, + "to" => r_out + ), + "top" => dicttype( + "from" => r_in_top, + "to" => r_out + ) + ), + "h" => h, + "origin" => [0, 0, crystal_height - top_taper_height / 2] + ))) + end end - end - # top taper - if haskey(meta.geometry.taper, :top) - top_taper_height = to_SSD_units(T, meta.geometry.taper.top.height_in_mm, u"mm") - if haskey(meta.geometry.taper.top, :radius_in_mm) - top_taper_radius = to_SSD_units(T, meta.geometry.taper.top.radius_in_mm, u"mm") - top_taper_angle = atan(top_taper_radius, top_taper_height) - elseif haskey(meta.geometry.taper.top, :angle_in_deg) - top_taper_angle = to_SSD_units(T, meta.geometry.taper.top.angle_in_deg, u"°") - top_taper_radius = top_taper_height * tan(top_taper_angle) + # bot outer taper + bot_taper_height = meta.geometry.taper.bottom.height_in_mm + if :radius_in_mm in keys(meta.geometry.taper.bottom) + bot_taper_radius = meta.geometry.taper.bottom.radius_in_mm + bot_taper_angle = atand(bot_taper_radius, bot_taper_height) + elseif :angle_in_deg in keys(meta.geometry.taper.bottom) + bot_taper_angle = meta.geometry.taper.bottom.angle_in_deg + bot_taper_radius = bot_taper_height * tand(bot_taper_angle) else - error("The top taper needs either radius_in_mm or angle_in_deg") + error("The bottom outer tape needs either radius_in_mm or angle_in_deg") end - has_top_taper = top_taper_height > 0 && top_taper_angle > 0 - if has_top_taper - r_center = crystal_radius - top_taper_height * tan(top_taper_angle) / 2 - hZ = top_taper_height/2 + 1gap - Δr = hZ * tan(top_taper_angle) - r_in_bot = r_center + Δr - r_in_top = r_center - Δr + has_bot_taper = bot_taper_height > 0 && bot_taper_angle > 0 + if has_bot_taper + r_center = crystal_radius - bot_taper_radius / 2 + hZ = bot_taper_height/2 + 1gap + Δr = hZ * tand(bot_taper_angle) + r_in_bot = r_center - Δr + r_in_top = r_center + Δr r_out = max(r_in_top, r_in_bot) + gap # ensure that r_out is always bigger as r_in - r = ((r_in_bot, r_out),(r_in_top, r_out)) - semiconductor_geometry -= CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, crystal_height - top_taper_height/2) - ) + push!(semiconductor_geometry_subtractions, dicttype("cone" => dicttype( + "r" => dicttype( + "bottom" => dicttype( + "from" => r_in_bot, + "to" => r_out + ), + "top" => dicttype( + "from" => r_in_top, + "to" => r_out + ) + ), + "h" => 2 * hZ, + "origin" => [0, 0, bot_taper_height / 2] + ))) + end + + # groove + has_groove = haskey(meta.geometry, :groove) + if has_groove + groove_inner_radius = meta.geometry.groove.radius_in_mm.inner + groove_outer_radius = meta.geometry.groove.radius_in_mm.outer + groove_depth = meta.geometry.groove.depth_in_mm + has_groove = groove_outer_radius > 0 && groove_depth > 0 && groove_inner_radius > 0 + if has_groove + hZ = groove_depth / 2 + gap + r_in = groove_inner_radius + r_out = groove_outer_radius + push!(semiconductor_geometry_subtractions, dicttype("cone" => dicttype( + "r" => dicttype( + "from" => r_in, + "to" => r_out + ), + "h" => 2 * hZ, + "origin" => [0, 0, groove_depth / 2 - gap] + ))) + end end end - # bot outer taper - bot_taper_height = to_SSD_units(T, meta.geometry.taper.bottom.height_in_mm, u"mm") - if :radius_in_mm in keys(meta.geometry.taper.bottom) - bot_taper_radius = to_SSD_units(T, meta.geometry.taper.bottom.radius_in_mm, u"mm") - bot_taper_angle = atan(bot_taper_radius, bot_taper_height) - elseif :angle_in_deg in keys(meta.geometry.taper.bottom) - bot_taper_angle = to_SSD_units(T, meta.geometry.taper.bottom.angle_in_deg, u"°") - bot_taper_radius = bot_taper_height * tan(bot_taper_angle) + if isempty(semiconductor_geometry_subtractions) + config_dict["detectors"][1]["semiconductor"]["geometry"] = semiconductor_geometry_basis else - error("The bottom outer tape needs either radius_in_mm or angle_in_deg") - end - has_bot_taper = bot_taper_height > 0 && bot_taper_angle > 0 - if has_bot_taper - r_center = crystal_radius - bot_taper_height * tan(bot_taper_angle) / 2 - hZ = bot_taper_height/2 + 1gap - Δr = hZ * tan(bot_taper_angle) - r_in_bot = r_center - Δr - r_in_top = r_center + Δr - r_out = max(r_in_top, r_in_bot) + gap # ensure that r_out is always bigger as r_in - r = ((r_in_bot, r_out),(r_in_top, r_out)) - semiconductor_geometry -= CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, bot_taper_height/2) + config_dict["detectors"][1]["semiconductor"]["geometry"] = dicttype( + "difference" => [semiconductor_geometry_basis, semiconductor_geometry_subtractions...] ) end - # groove - has_groove = haskey(meta.geometry, :groove) - if has_groove - groove_inner_radius = to_SSD_units(T, meta.geometry.groove.radius_in_mm.inner, u"mm") - groove_outer_radius = to_SSD_units(T, meta.geometry.groove.radius_in_mm.outer, u"mm") - groove_depth = to_SSD_units(T, meta.geometry.groove.depth_in_mm, u"mm") - has_groove = groove_outer_radius > 0 && groove_depth > 0 && groove_inner_radius > 0 - if has_groove - hZ = groove_depth / 2 + gap - r_in = groove_inner_radius - r_out = groove_outer_radius - r = ((r_in, r_out), (r_in, r_out)) - semiconductor_geometry -= CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, groove_depth / 2 - gap) - ) - end - end # bulletization @@ -214,33 +332,56 @@ function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::Pro ### P+ CONTACT ### - pp_radius = to_SSD_units(T, meta.geometry.pp_contact.radius_in_mm, u"mm") - pp_depth = to_SSD_units(T, meta.geometry.pp_contact.depth_in_mm, u"mm") - pp_contact_geometry = if is_coax - CSG.Cone{T}(CSG.ClosedPrimitive; - r = ((borehole_radius, borehole_radius), (borehole_radius, borehole_radius)), - hZ = borehole_depth/2, - origin = CartesianPoint{T}(0, 0, borehole_depth / 2) - ) + CSG.Cone{T}(CSG.ClosedPrimitive; - r = borehole_radius, - hZ = 0, - origin = CartesianPoint{T}(0, 0, borehole_depth) - ) + CSG.Cone{T}(CSG.ClosedPrimitive; - r = ((borehole_radius, pp_radius), (borehole_radius, pp_radius)), - hZ = 0 - ) + pp_radius = meta.geometry.pp_contact.radius_in_mm + pp_depth = meta.geometry.pp_contact.depth_in_mm + push!(config_dict["detectors"][1]["contacts"], dicttype( + "material" => "HPGe", + "geometry" => dicttype(), + "id" => 1, + "potential" => 0 + )) + config_dict["detectors"][1]["contacts"][1]["geometry"] = if is_coax + dicttype("union" => [ + dicttype("cone" => dicttype( + "r" => dicttype( + "from" => borehole_radius, + "to" => borehole_radius + ), + "h" => borehole_depth, + "origin" => [0, 0, borehole_depth / 2] + )), + dicttype("cone" => dicttype( + "r" => borehole_radius, + "h" => 0, + "origin" => [0, 0, borehole_depth] + )), + dicttype("cone" => dicttype( + "r" => dicttype( + "from" => borehole_radius, + "to" => pp_radius + ), + "h" => 0 + )) + ]) else - CSG.Cone{T}(CSG.ClosedPrimitive; - r = pp_radius, - hZ = pp_depth / 2, - origin = CartesianPoint{T}(0, 0, pp_depth / 2) - ) + dicttype("cone" => dicttype( + "r" => pp_radius, + "h" => pp_depth, + "origin" => [0, 0, pp_depth / 2] + )) end ### MANTLE CONTACT ### - - mantle_contact_geometry = begin # top plate + + push!(config_dict["detectors"][1]["contacts"], dicttype( + "material" => "HPGe", + "geometry" => dicttype("union" => []), + "id" => 2, + "potential" => meta.characterization.manufacturer.recommended_voltage_in_V + )) + config_dict["detectors"][1]["contacts"][2]["geometry"]["union"] = begin + mantle_contact_parts = [] top_plate = begin r = if !has_borehole || is_coax !has_top_taper ? crystal_radius : crystal_radius - top_taper_radius @@ -249,149 +390,153 @@ function SolidStateDetectors.SolidStateDetector{T}(::Type{LegendData}, meta::Pro r_out = crystal_radius if has_borehole_taper r_in += borehole_taper_radius end if has_top_taper r_out -= top_taper_radius end - ((r_in, r_out), (r_in, r_out)) + dicttype("from" => r_in, "to" => r_out) end - CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = li_thickness / 2, - origin = CartesianPoint{T}(0, 0, crystal_height - li_thickness / 2) - ) + dicttype("cone" => dicttype( + "r" => r, + "h" => li_thickness, + "origin" => [0, 0, crystal_height - li_thickness / 2] + )) end - mc_geometry = top_plate - - # borehole at outer taper + push!(mantle_contact_parts, top_plate) + if has_top_taper - Δr_li_thickness = li_thickness / cos(top_taper_angle) - hZ = top_taper_height/2 + Δr_li_thickness = li_thickness / cosd(top_taper_angle) + h = top_taper_height r_bot = crystal_radius r_top = crystal_radius - top_taper_radius - r = ((r_bot - Δr_li_thickness, r_bot),(r_top - Δr_li_thickness, r_top)) - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, crystal_height - top_taper_height/2) - ) + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => dicttype( + "bottom" => dicttype( + "from" => r_bot - Δr_li_thickness, + "to" => r_bot + ), + "top" => dicttype( + "from" => r_top - Δr_li_thickness, + "to" => r_top + ) + ), + "h" => h, + "origin" => [0, 0, crystal_height - top_taper_height / 2] + ))) end - # contact in borehole if has_borehole_taper - Δr_li_thickness = li_thickness / cos(borehole_taper_angle) - hZ = borehole_taper_height/2 + Δr_li_thickness = li_thickness / cosd(borehole_taper_angle) + h = borehole_taper_height r_bot = borehole_radius r_top = borehole_radius + borehole_taper_radius - r = ((r_bot, r_bot+Δr_li_thickness),(r_top, r_top+Δr_li_thickness)) - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, crystal_height - borehole_taper_height/2) - ) - - hZ = (borehole_depth - borehole_taper_height) / 2 - r = ((borehole_radius, borehole_radius+Δr_li_thickness),(borehole_radius, borehole_radius+Δr_li_thickness)) - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, crystal_height - borehole_taper_height - hZ) - ) + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => dicttype( + "bottom" => dicttype( + "from" => r_bot, + "to" => r_bot + Δr_li_thickness + ), + "top" => dicttype( + "from" => r_top, + "to" => r_top + Δr_li_thickness + ) + ), + "h" => h, + "origin" => [0, 0, crystal_height - borehole_taper_height / 2] + ))) + + h = (borehole_depth - borehole_taper_height) + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => dicttype( + "from" => borehole_radius, + "to" => borehole_radius + Δr_li_thickness + ), + "h" => h, + "origin" => [0, 0, crystal_height - borehole_taper_height - h / 2] + ))) elseif has_borehole && !is_coax # but no borehole taper - hZ = borehole_depth / 2 - r = ((borehole_radius, borehole_radius+li_thickness),(borehole_radius, borehole_radius+li_thickness)) - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, crystal_height - hZ) - ) + h = borehole_depth + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => dicttype( + "from" => borehole_radius, + "to" => borehole_radius + li_thickness + ), + "h" => h, + "origin" => [0, 0, crystal_height - h / 2] + ))) end if has_borehole && !is_coax r = borehole_radius + li_thickness - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = li_thickness / 2, - origin = CartesianPoint{T}(0, 0, crystal_height - borehole_depth - li_thickness / 2) - ) + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => r, + "h" => li_thickness / 2, + "origin" => [0, 0, crystal_height - borehole_depth - li_thickness / 2] + ))) end - # outer surface of mantle contact begin - r = ((crystal_radius-li_thickness, crystal_radius),(crystal_radius-li_thickness, crystal_radius)) - hZ = crystal_height - if has_top_taper hZ -= top_taper_height end - z_origin = hZ/2 + h = crystal_height + if has_top_taper h -= top_taper_height end + z_origin = h/2 if has_bot_taper - hZ -= bot_taper_height + h -= bot_taper_height z_origin += bot_taper_height/2 end - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ / 2, - origin = CartesianPoint{T}(0, 0, z_origin) - ) + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => dicttype( + "from" => crystal_radius - li_thickness, + "to" => crystal_radius + ), + "h" => h, + "origin" => [0, 0, z_origin] + ))) end - # bottom outer taper contact if has_bot_taper - Δr_li_thickness = li_thickness / cos(bot_taper_angle) - hZ = bot_taper_height/2 + Δr_li_thickness = li_thickness / cosd(bot_taper_angle) + h = bot_taper_height r_bot = crystal_radius - bot_taper_radius r_top = crystal_radius - r = ((r_bot - Δr_li_thickness, r_bot),(r_top - Δr_li_thickness, r_top)) - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = hZ, - origin = CartesianPoint{T}(0, 0, hZ) - ) - end + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => dicttype( + "bottom" => dicttype( + "from" => r_bot - Δr_li_thickness, + "to" => r_bot + ), + "top" => dicttype( + "from" => r_top - Δr_li_thickness, + "to" => r_top + ) + ), + "h" => h, + "origin" => [0, 0, h / 2] + ))) + end - # bottom surface of mantle contact (only if it has a groove ?) if has_groove && groove_outer_radius > 0 r_in = groove_outer_radius r_out = crystal_radius if has_bot_taper r_out -= bot_taper_radius end - r = ((r_in, r_out), (r_in, r_out)) - mc_geometry += CSG.Cone{T}(CSG.ClosedPrimitive; - r = r, - hZ = li_thickness / 2, - origin = CartesianPoint{T}(0, 0, li_thickness / 2) - ) + push!(mantle_contact_parts, dicttype("cone" => dicttype( + "r" => dicttype( + "from" => r_in, + "to" => r_out + ), + "h" => li_thickness, + "origin" => [0, 0, li_thickness / 2] + ))) end - - mc_geometry + + mantle_contact_parts end - - # Hardcoded parameter values: In future, should be defined in config file - temperature = T(78) - material = SolidStateDetectors.material_properties[:HPGe] - - # Impurity Model: Information are stored in `meta.production.impcc` - # For now: Constant impurity density: - # n-type: positive impurity density - # p-type: negative impurity density - # Assume p-type - constant_impurity_density = ustrip(uconvert(u"m^-3", T(-1e9) * u"cm^-3")) - impurity_density_model = SolidStateDetectors.CylindricalImpurityDensity{T}( - (0, 0, constant_impurity_density), # offsets - (0, 0, 0) # linear slopes + + config_dict["detectors"][1]["semiconductor"]["impurity_density"] = dicttype( + "name" => "constant", + "value" => "-1e9cm^-3" ) - # Charge Drift Model: - # Use example ADL charge drift model from SSD (Crystal axis <100> is at φ = 0): - adl_charge_drift_config_file = joinpath(dirname(dirname(pathof(SolidStateDetectors))), - "examples/example_config_files/ADLChargeDriftModel/drift_velocity_config.yaml") - charge_drift_model = SolidStateDetectors.ADLChargeDriftModel{T}(adl_charge_drift_config_file); - - semiconductor = SolidStateDetectors.Semiconductor(temperature, material, impurity_density_model, charge_drift_model, semiconductor_geometry) - - operation_voltage = T(meta.characterization.manufacturer.recommended_voltage_in_V) - pp_contact = SolidStateDetectors.Contact( zero(T), material, 1, "Point Contact", pp_contact_geometry ) - mantle_contact = SolidStateDetectors.Contact( operation_voltage, material, 2, "Mantle Contact", mantle_contact_geometry ) - - semiconductor, (pp_contact, mantle_contact) + # evaluate "include" statements - needed for the charge drift model + SolidStateDetectors.scan_and_merge_included_json_files!(config_dict, "") - passives = missing # possible holding structure around the detector - virtual_drift_volumes = missing - SolidStateDetector{T}( meta.name, semiconductor, [pp_contact, mantle_contact], passives, virtual_drift_volumes ) + return config_dict end end # module LegendDataManagementSolidStateDetectorsExt diff --git a/test/test_ext_ssd.jl b/test/test_ext_ssd.jl index 03707048..a4affa05 100644 --- a/test/test_ext_ssd.jl +++ b/test/test_ext_ssd.jl @@ -2,91 +2,37 @@ using LegendDataManagement using Test + using SolidStateDetectors +using SolidStateDetectors: ConstantImpurityDensity +using Unitful +using LegendHDF5IO include("testing_utils.jl") -function SolidStateDetectors.Simulation(l200, detname; medium::Union{String, Symbol} = "vacuum") - - config_dict = Dict() - radius = l200.metadata.hardware.detectors.germanium.diodes[Symbol(detname)].geometry.radius_in_mm / 1000 - height = l200.metadata.hardware.detectors.germanium.diodes[Symbol(detname)].geometry.height_in_mm / 1000 - det = SolidStateDetector(l200, detname); - - function get_material_name(material::NamedTuple) - collect(keys(SolidStateDetectors.material_properties))[findfirst(getindex.(values(SolidStateDetectors.material_properties), :name) .== material.name)] - end - - config_dict["name"] = "$(detname)" - config_dict["units"] = Dict( - "length" => "m", - "angle" => "deg", - "potential" => "V", - "temperature" => "K" - ) - config_dict["grid"] = Dict( - "coordinates" => "cylindrical", - "axes" => Dict( - "r" => Dict( - "to" => "$(radius + 0.01)", - "boundaries" => "inf" - ), - "phi" => Dict( - "from" => 0, - "to" => 0, - "boundaries" => "periodic" - ), - "z" => Dict( - "from" => "-0.01", - "to" => "$(height + 0.01)", - "boundaries" => Dict( - "left" => "inf", - "right" => "inf" - ) - ) - ) - ) - config_dict["medium"] = medium - config_dict["detectors"] = [ - Dict( - "semiconductor" => Dict( - "material" => "$(get_material_name(det.semiconductor.material))", - "temperature" => det.semiconductor.temperature, - "geometry" => SolidStateDetectors.ConstructiveSolidGeometry.Dictionary(det.semiconductor.geometry) - # charge_drift_model ? - # impurity_density ? - ), - "contacts" => [ - Dict( - "material" => "$(get_material_name(c.material))", - "id" => c.id, - "potential" => c.potential, - "geometry" => SolidStateDetectors.ConstructiveSolidGeometry.Dictionary(c.geometry) - ) - for c in det.contacts - ] - ) - ] - T = SolidStateDetectors.get_precision_type(det) - Simulation{T}(config_dict) -end - @testset "test_ext_ssd" begin l200 = LegendData(:l200) for detname in (:V99000A, :B99000A, :C99000A, :P99000A) @testset "$(detname)" begin - det = SolidStateDetector(l200, detname) + det = SolidStateDetector{Float64}(l200, detname) @test det isa SolidStateDetector - sim = Simulation(l200, detname) + sim = Simulation{Float64}(l200, detname) @test sim isa Simulation # Compare active volume from SSD to active volume from LegendDataManagement - sim.detector = SolidStateDetector(sim.detector, contact_id = 1, contact_potential = 0) # Set potential to 0V to avoid long simulation times - sim.detector = SolidStateDetector(sim.detector, contact_id = 2, contact_potential = 0) # (we just need the PointTypes for the active volume) - calculate_electric_potential!(sim, max_tick_distance = 0.1u"mm", refinement_limits = missing, verbose = false) + SolidStateDetectors.apply_initial_state!(sim, ElectricPotential, Grid(sim, max_tick_distance = 0.1u"mm")) active_volume_ssd = SolidStateDetectors.get_active_volume(sim.point_types) active_volume_ldm = LegendDataManagement.get_active_volume(l200.metadata.hardware.detectors.germanium.diodes[Symbol(detname)], 0.0) @test isapprox(active_volume_ssd, active_volume_ldm, rtol = 0.01) + + # The creation via config files allows to save Simulations to files using LegendHDF5IO + lh5name = "$(detname).lh5" + isfile(lh5name) && rm(lh5name) + @test_nowarn ssd_write(lh5name, sim) + @test isfile(lh5name) + @test sim == ssd_read(lh5name, Simulation) + @test_nowarn rm(lh5name) + @test !isfile(lh5name) end end end