From 3dd37fcd2e2302ef6770f0a5becaa1408f85752e Mon Sep 17 00:00:00 2001 From: Shailesh Giri Date: Thu, 17 Oct 2024 15:15:12 +0200 Subject: [PATCH 1/9] Creating extensions --- ext/LegendGeSimExt.jl | 594 +++++++++++++++++++++++++++++++ ext/legend_detector_to_siggen.jl | 263 ++++++++++++++ ext/mjdsiggen_utils.jl | 30 ++ test/Project.toml | 4 +- test/runtests.jl | 6 +- 5 files changed, 894 insertions(+), 3 deletions(-) create mode 100644 ext/LegendGeSimExt.jl create mode 100644 ext/legend_detector_to_siggen.jl create mode 100644 ext/mjdsiggen_utils.jl diff --git a/ext/LegendGeSimExt.jl b/ext/LegendGeSimExt.jl new file mode 100644 index 0000000..f5083bb --- /dev/null +++ b/ext/LegendGeSimExt.jl @@ -0,0 +1,594 @@ +module LegendGeSimExt +#using MJDSigGen: SiggenSimulator +#@info "Error here" +using LegendGeSim +#@info "Error" +using PropDicts +#@info "Error_prop" +using LegendGeSim +#@info "Error_Lgsim" +using LegendGeSim: Environment +#@info "le_env" +#using MJDSigGen: SiggenSimulator +using MJDSigGen +using MJDSigGen: SigGenSetup +#@info "Error_mjd" +#using LegendGeSim: PSSimulator +using LegendGeSim: SiggenSimulator +using ArgCheck +using ArraysOfArrays +using CurveFit +using DelimitedFiles +using Distributions +using DSP +using ElasticArrays +using HDF5 +using IntervalSets +using JSON +using LegendDataManagement +using LegendDataTypes +using LegendHDF5IO +using LegendTextIO +using LinearAlgebra +using Parameters +using Polynomials +using PropDicts +using RadiationDetectorDSP +using RadiationDetectorSignals +using RadiationSpectra +using Random +using Random123 +using RecipesBase +using Requires +using SignalAnalysis +using SolidStateDetectors +using StaticArrays +using Statistics +using StatsBase +using StructArrays +using Tables +using TypedTables +using Unitful + +import LsqFit + +using SolidStateDetectors: ConstructiveSolidGeometry as CSG + +#abstract type PSSimulator end + + +#include("legend_detector_to_siggen.jl") +#@info "Here_then" +#include("mjdsiggen_utils.jl") +#@info "Error_after_files" + + + +#################################### +### FieldGen +#################################### +@info "27" +""" + SiggeSimulator(sim_conf::PropDict) + + LegendGeSimConfig -> SiggenSimulator + +Construct SiggeSimulator instance based on simulation + configuration given in . +""" + +function SiggenSimulator(simulation_settings::PropDict) + # set defaults + inputs = Dict{String, Any}() + for k in (:fieldgen_config, :drift_vel, :impurity_profile, :crystal_metadata_path) + inputs[string(k)] = haskey(simulation_settings, k) ? simulation_settings[k] : "" + end + + # check if provided paths exist + for (param, path) in inputs + if (path != "") && !(isfile(path) || isdir(path)) + throw(ArgumentError("The input for $(param) that you provided does not exist: $(path)")) + end + end + + inputs["offset_in_mm"] = haskey(simulation_settings, :offset_in_mm) ? simulation_settings.offset_in_mm : -1 + + # throw error if both files are provided + if inputs["crystal_metadata_path"] != "" && inputs["impurity_profile"] != "" + throw(ArgumentError("You provided both an .spe/.dat file and path to crystal metadata! Provide one or the other as impurity input for siggen.")) + # warn if none are provided - will use constant impurity density + elseif inputs["crystal_metadata_path"] == "" && inputs["impurity_profile"] == "" + @warn "No impurity inputs given. Simulation with dummy constant impurity density." + # if we get here, one is provided + else + impinput = inputs["impurity_profile"] * inputs["crystal_metadata_path"] + @info "Impurity profile information based on $(impinput)" + end + + # check that offset is provided if .spe/.dat file is + if inputs["impurity_profile"] != "" && inputs["offset_in_mm"] == -1 + throw(ArgumentError("Provide offset_in_mm of this detector corresponding to impurity file $(inputs["impurity_profile"])!")) + end + + + SiggenSimulator( + inputs["fieldgen_config"], + inputs["drift_vel"], + inputs["impurity_profile"], + inputs["offset_in_mm"], + inputs["crystal_metadata_path"], + simulation_settings.cached_name + ) +end + +""" + simulate_detector(det_meta, env, config_name, ps_simulator) + +PropDict, Environment, AbstractString, Siggen -> SigGenSetup + +Construct a fieldgen/siggen configuration file + according to geometry as given in LEGEND metadata + and environment settings given in . +Look up fieldgen generated electric field and weighting potential files + based on given , and if not found, call fieldgen. +""" + +function LegendGeSim.simulate_detector(det_meta::PropDict, env::Environment, simulator::SiggenSimulator; + overwrite::Bool = false) + +# if no cached name given, force simulation from scratch (or else might read past "tmp" file) +if simulator.cached_name == "" + overwrite = true +end + +# returns the name of the resulting siggen config file +# and the name of (already or to be) generated weighting potential file +# ToDo: don't create if already present already at this stage -> check in siggen_config() if name exists and just return name +siggen_config_name, fieldgen_wp_name = siggen_config(det_meta, env, simulator; overwrite) +fieldgen_wp_name = joinpath("cache", fieldgen_wp_name) + +# if the WP file with such a name does not exist yet or want to overwrite it... +if !isfile(fieldgen_wp_name) || overwrite + imp_filename, offset = + if simulator.crystal_metadata_path != "" + # create impurity input on the fly + LegendGeSim.impurity_density_model(det_meta, simulator.crystal_metadata_path, simulator) + else + # user provided .dat file and corresponding offset + # (or if nothing is provided this will be "" and -1) + simulator.impurity_profile, simulator.offset_in_mm + end + #...call fieldgen -> will write the wp file + @info "_|~|_|~|_|~|_ Fieldgen simulation" + fieldgen(siggen_config_name; impurity_profile=imp_filename, offset_mm=offset) + @info "_|~|_|~|_|~|_ Fieldgen simulation complete" +else + #...do nothing, siggen will later read the files based on the generated conifg + @info "Cached simulation found. Reading cached fields from $fieldgen_wp_name" +end +@info "73" +# a SigGenSetup object +SigGenSetup(siggen_config_name) +end + +## ---------- fieldgen + +""" +Convert fit parameters from crystal metadata units (1e9 e/cm^3 VS mm) into Siggen units (1e10 e/cm^3 VS mm) +and create a .dat file (unformatted stream of Float32) to be later used when calling fieldgen(). +""" +function LegendGeSim.impurity_density_model(meta::PropDict, crystal_metadata::PropDict, fit_param::Vector{Float64}, ::SiggenSimulator) + + a,b,c,tau = fit_param + + T = Float32 + dist = T.(crystal_metadata.impurity_measurements.distance_from_seed_end_mm) + # sometimes might not have measurements at seed end? but last measurement can be taken as crystal length + L = dist[end] #mm + + # points in crystal axis from 0 to L with step of 1 mm + # always put 200 points to make all .dat files same length + imp = zeros(Float32, 200) # impurities in 10^10 cm^-3 + + # length is L+1 including both ends + Npoints = Int(ceil(L))+1 # step of 1 mm + zpoints = LinRange(0, L, Npoints) + for i in 1:Npoints + # divide by 10 to convert from 1e9 e/cm^3 to 1e10 e/cm^3 + # invert the order of the values because with siggen z=0 is tail (dirtiest part) + # i.e. z = L in our system (Mirion measurements) + imp[Npoints-i+1] = -(a + b * zpoints[i] + c * exp((zpoints[i] - L)/tau)) / 10. + end + + crystal_name = first(meta.name, length(meta.name)-1) + imp_filename = joinpath("cache", crystal_name*".dat") + open(imp_filename, "w") do io + write(io, imp) + end + @info "Impurity profile for siggen saved in $(imp_filename)" + + # return name of impurity profile file and offset + # offset from seed end + det_z0 = crystal_metadata.slices[Symbol(meta.production.slice)].detector_offset_in_mm + + # siggen offset is from tail i.e. at L offset = 0 + imp_filename, L - det_z0 +end + +@info "123" +""" +Simulation method: siggen +""" + + + +function LegendGeSim.simulate_waveforms(stp_events::Table, detector::SigGenSetup, ::SiggenSimulator) + T = Float32 # This should be somehow defined and be passed properly + println("127") + @info("~.~.~ Siggen") + @info "129" + + nevents = length(stp_events) + wf_array = Array{RDWaveform}(undef, nevents) + for i in 1:nevents + # println("$i/$nevents") + if(i % 100 == 0) println("$i/$nevents") end + signal::Vector{Float32} = simulate_signal(stp_events[i].pos, stp_events[i].edep, detector) + time = range(T(0)u"ns", step = detector.step_time_out*u"ns", length = length(signal)) + + # push!(waveforms, signal) + wf_array[i] = RDWaveform(time, signal) + end + + ## Oliver said this should be the more efficient way of doing it + ## as opposed to an array of separate RDWaveform instances + # Δt = setup.step_time_out*u"ns" + # t_end = length(waveforms[1]) * Δt + # times = fill(0u"ns" : Δt : t_end, length(waveforms)) + # ArrayOfRDWaveforms((times, VectorOfSimilarVectors(waveforms))) + + waveforms = ArrayOfRDWaveforms(wf_array) + # why am I doing this? + waveforms = ArrayOfRDWaveforms((waveforms.time, VectorOfSimilarVectors(waveforms.signal))) + + pss_table = Table( + channel = [1 for idx in 1:length(waveforms)], # lists of ADCs that triggered, 1 for HADES all the time + ievt = stp_events.evtno, + waveform = waveforms + ) + + # using SSD we get mc truth from SSD output + # with siggen, let's just return the input + # maybe mc truth will not be propagated to (sim)raw anyway + pss_truth = Table( + detno = stp_events.detno, + edep = stp_events.edep, + ievt = stp_events.evtno, + pos = stp_events.pos, + thit = stp_events.thit + ) + + pss_table, pss_truth +end + + +#@info "252" + +""" + simulate_wf(pos, edep, siggen_setup) + +AbstractVector, AbstractVector, SigGenSetup -> Vector{Float32} + +Simulate a signal from energy depositions with corresponding positions + based on a given +""" + +function simulate_signal(pos::AbstractVector, edep::AbstractVector, siggen_setup::SigGenSetup) + # this is not so nice, think about it + signal = zeros(Float32, siggen_setup.ntsteps_out) + + # round to siggen crystal grid precision + # a = string(siggen_setup.xtal_grid) # crystal grid precision e.g. 0.1 (mm) + # ndig = length(a[findfirst('.', a)+1:end]) # number of digits after . + + for i in 1:length(pos) + # pos_rounded = round.(ustrip.(pos[i]), digits=ndig) # strip of units and round to siggen precision + # in SSD the output is always in eV -> convert to eV + signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(uconvert(u"eV", edep[i])) + # signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(pos_rounded)) * ustrip(uconvert(u"eV", edep[i])) + + # somehow this doesn't work + # MJDSigGen.get_signal!(signal, siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(uconvert(u"eV", edep[i])) + end + + signal +end + +#Calculating capacitance matrix using Siggen +function LegendGeSim.capacitance_matrix(sim::SigGenSetup) + c = sim.capacitance * u"pF" + [c -c; -c c] +end +#export capacitance_matrix +#@info "288" +#legend_detector_to_siggen +function siggen_config(meta::PropDict, env::Environment, simulator::SiggenSimulator; overwrite::Bool=false) + # ---------------------------------------------------------------------------- + # set up & check if already present in cache + # ---------------------------------------------------------------------------- + + # if no cached name was given by user, make a temporaty name - needed for fieldgen to read from file + # also if no cached name was given means definitely doing from scratch + cached_name = simulator.cached_name + if simulator.cached_name == "" + cached_name = "tmp" + overwrite = true + end + + # append detector name + cached_name = meta.name * "_" * cached_name + + # filenames for fieldgen input/output + # ! no need to add cache/ folder because siggen does it by itself... + fieldgen_names = Dict( + "field_name" => cached_name * "_fieldgen_EV.dat", + "wp_name" => cached_name * "_fieldgen_WP.dat" + ) + + # siggen config filename + sig_conf_name = joinpath("cache", cached_name * "_siggen_config.txt") + + # if field exists, no need to construct config, just read cached simulation -> return + # if config exists, no need to construct again -> return -> maybe just the cached field? not sure + # of course unless asked to overwrite, then construct new config + if ( isfile(joinpath("cache", fieldgen_names["wp_name"])) || isfile(sig_conf_name) ) && !overwrite + return sig_conf_name, fieldgen_names["wp_name"] + end + + @info "Constructing Fieldgen/Siggen configuration file" + + if !ispath("cache") mkpath("cache") end + + + # ---------------------------------------------------------------------------- + # construct config + # ---------------------------------------------------------------------------- + + println("...detector geometry") + + # construct siggen geometry based on LEGEND metadata JSON + environment settings + siggen_geom = meta2siggen(meta, env) + + # collect geometry lines for siggen config file + detector_lines = Vector{String}() + push!(detector_lines, "# detector related parameters") + for (param, value) in siggen_geom + push!(detector_lines, "$param $value") + end + + println("...fieldgen configuration") + + # ---------------------------------------------------------------------------- + # drift velocity, field, wp + # ---------------------------------------------------------------------------- + + # drift velocity input (path to file given by user) + drift_file_path = simulator.drift_vel + # use default if not given by user + if simulator.drift_vel == "" + @warn "No drift velocity input given; using default" + drift_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "drift_vel_tcorr.tab") + end + + # copy to cache + cache_drift_file = cached_name * "_" * basename(drift_file_path) + cp(drift_file_path, joinpath("cache", cache_drift_file), force=true) # later think what to do here + + # add to field fieldgen_names + fieldgen_names["drift_name"] = cache_drift_file + + ## commenting this out: produced a bug + ## (empty file created, next time running empty file read and no field at point error) + # for name in ["field_name", "wp_name"] + # path = joinpath("cache", fieldgen_names[name]) + # if !isfile(path) + # touch(path) + # end + # end + + # add corresponding lines to existing detector geometry lines + push!(detector_lines, "") + push!(detector_lines, "# fieldgen filenames") + for (param, value) in fieldgen_names + push!(detector_lines, "$param $value") + end + + # ---------------------------------------------------------------------------- + # extra settings + # ---------------------------------------------------------------------------- + + # Note: this is similar but not identical to the part above because of the cache folder difference + # -> loop? copy-paste is not nice... + + # read lines from user input for fieldgen + config_file_path = simulator.fieldgen_config + if simulator.fieldgen_config == "" + @warn "No extra fieldgen settings given; using default" + config_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "fieldgen_settings.txt") + end + # copy to cache -> maybe not needed, are appended in main settings? + cache_config_file = joinpath("cache", cached_name * "_" * basename(config_file_path)) + cp(config_file_path, cache_config_file, force=true) # later think what to do here + + fieldgen_lines = readlines(open(cache_config_file)) + fieldgen_lines = replace.(fieldgen_lines, "\t" => " ") + + # unite constructed detector geometry and fieldgen settings + total_lines = vcat(detector_lines, [""], fieldgen_lines) + + # ---------------------------------------------------------------------------- + # final siggen config file + # ---------------------------------------------------------------------------- + + # write a siggen/fieldgen config file + + writedlm(sig_conf_name, total_lines) + println("...fieldgen/siggen config written to $sig_conf_name") + + # return resulting fieldgen/siggen config name + # (currently kinda redundant but that's because we have no idea what's gonna happen later) + sig_conf_name, fieldgen_names["wp_name"] +end +#@info "417" + +""" + meta2siggen(meta, env) + +PropDict, Environment -> Dict + +Construct siggen type config in a form of Dict based on + geometry as given in LEGEND metadata + and environment settings given in . +""" +function meta2siggen(meta::PropDict, env::Environment) + # !!!! this is copy paste from legend_detector_to_SSD !! + # ToDo: rewrite to avoid copy-paste + + # ----------------------------------------------------------- + # DL thickness (quickfix) + + dl_thickness_in_mm = + if env.dl == "vendor" + dl_vendor = meta.characterization.manufacturer.dl_thickness_in_mm + if dl_vendor isa PropDicts.MissingProperty + throw(ArgumentError("No dead layer measurement provided by vendor! Please provide value or skip (default 0).")) + end + dl_vendor + else + env.dl + end + + @info "DL = $(dl_thickness_in_mm)mm" + + # ----------------------------------------------------------- + # operating voltage + + # sometimes recV is null -> in this case, complain and ask to provide (default env.opV is 0 if none given) + operation_voltage = + if env.operating_voltage == 0 + @warn "You did not provide operating voltage in environment settings -> taking recommended voltage from metadata" + recV = meta.characterization.l200_site.recommended_voltage_in_V + if isnothing(recV) + error("Metadata does not provide recommended voltage (null). Please provide operating voltage in settings") + end + recV + else + env.operating_voltage + end + + @info "Simulating at $(operation_voltage)V" + + # ----------------------------------------------------------- + # parameters common for all geometries + + siggen_geometry = Dict( + # -- crystal + "xtal_length" => meta.geometry.height_in_mm, + "xtal_radius" => meta.geometry.radius_in_mm, + # -- p+ contact + "pc_length" => meta.geometry.pp_contact.depth_in_mm, + "pc_radius" => meta.geometry.pp_contact.radius_in_mm, + # -- tapering + + # comment from David Radford: + + # The bottom taper is always at 45 degrees, but the top outer taper is not. + # It's width/angle is defined using either the outer_taper_width parameter or the taper_angle parameter. + # Likewise, the inner taper width/angle is defined using either the inner_taper_width parameter or the taper_angle parameter. + # The inner_taper_length can be anything from zero to the hole_length + + # size of 45-degree taper at bottom of ORTEC-type crystal (for r=z) + # -> for now, ignore bottom taper if not 45 + # -> works for 1) the ones with actual 45 taper, 2) the ones bulletized (saved as 45 deg taper with bullet radius) + # few that have proper non-45 taper - not possible. ToDo: put some warning + "bottom_taper_length" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm : 0, + # current metadata format: always angle + "bottom_taper_width" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm * tan(deg2rad(meta.geometry.taper.bottom.angle_in_deg)) : 0, + # z-length of outside taper for inverted-coax style -> Mariia: you mean top taper here? + "outer_taper_length" => meta.geometry.taper.top.height_in_mm, + "outer_taper_width" => meta.geometry.taper.top.height_in_mm * tan(deg2rad(meta.geometry.taper.top.angle_in_deg)), + # taper angle in degrees, for inner or outer taper -> how? why are they the same? how about borehole? not using this parameter... + # "taper_angle" => meta.geometry.taper.top.outer.angle_in_deg, + # "taper_angle" => meta.geometry.taper.top.inner.angle_in_deg, + # with this setting, it somehow takes taper angle: 45.000038 anyway -> force to zero? + # "taper_angle" => 0, + # depth of full-charge-collection boundary for Li contact + "Li_thickness" => dl_thickness_in_mm, + + # configuration for mjd_fieldgen (calculates electric fields & weighing potentials) + # detector bias for fieldgen, in Volts + "xtal_HV" => operation_voltage, + + # configuration for signal calculation + # crystal temperature in Kelvin + "xtal_temp" => env.crystal_temperature + ) + + # ----------------------------------------------------------- + # non common parameters + + # -- groove (for non PPCs) + if haskey(meta.geometry, :groove) + siggen_geometry["wrap_around_radius"] = meta.geometry.groove.radius_in_mm.outer + siggen_geometry["ditch_depth"] = meta.geometry.groove.depth_in_mm + siggen_geometry["ditch_thickness"] = meta.geometry.groove.radius_in_mm.outer - meta.geometry.groove.radius_in_mm.inner + end + + # -- borehole (for ICPC and Coax) + if haskey(meta.geometry, :borehole) + # length of hole, for inverted-coax style + siggen_geometry["hole_length"] = meta.geometry.borehole.depth_in_mm + # radius of hole, for inverted-coax style + siggen_geometry["hole_radius"] = meta.geometry.borehole.radius_in_mm + # z-length of inside (hole) taper for inverted-coax style + siggen_geometry["inner_taper_length"] = meta.geometry.taper.borehole.height_in_mm + siggen_geometry["inner_taper_width"] = meta.geometry.taper.borehole.height_in_mm * tan(deg2rad(meta.geometry.taper.borehole.angle_in_deg)) + end + + siggen_geometry +end + +#mjdsiggen_utils.jl +#@info "537" + +function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) + E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); + T = eltype(E_pot) + r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid + z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid + ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) + ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) + ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) + axs = (ax1, ax2, ax3) + grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) + data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); + data[:] = E_pot'[:]; + ElectricPotential(data, grid) +end + +function SolidStateDetectors.WeightingPotential(sim::SigGenSetup) + E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); + T = eltype(W_pot) + r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid + z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid + ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) + ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) + ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) + axs = (ax1, ax2, ax3) + grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) + data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); + data[:] = W_pot'[:]; + WeightingPotential(data, grid) +end + +end #module LegendGeSimExt diff --git a/ext/legend_detector_to_siggen.jl b/ext/legend_detector_to_siggen.jl new file mode 100644 index 0000000..7bf9617 --- /dev/null +++ b/ext/legend_detector_to_siggen.jl @@ -0,0 +1,263 @@ + +using PropDicts + +""" + siggen_config(meta, env, siggen_sim, config_name) + +PropDict, Environment, SiggenSimulator, String -> String, String + +Construct Siggen/Fieldgen config based on + geometry as given in LEGEND metadata , + environment settings given in , + and fieldgen settings contained in . + +The resulting siggen config will be cached based on given +The function returns the path to the cached siggen config file, and + the relative path to the fieldgen WP file (to check if already exists later) +""" + +function siggen_config(meta::PropDict, env::Environment, simulator::SiggenSimulator; overwrite::Bool=false) + # ---------------------------------------------------------------------------- + # set up & check if already present in cache + # ---------------------------------------------------------------------------- + + # if no cached name was given by user, make a temporaty name - needed for fieldgen to read from file + # also if no cached name was given means definitely doing from scratch + cached_name = simulator.cached_name + if simulator.cached_name == "" + cached_name = "tmp" + overwrite = true + end + + # append detector name + cached_name = meta.name * "_" * cached_name + + # filenames for fieldgen input/output + # ! no need to add cache/ folder because siggen does it by itself... + fieldgen_names = Dict( + "field_name" => cached_name * "_fieldgen_EV.dat", + "wp_name" => cached_name * "_fieldgen_WP.dat" + ) + + # siggen config filename + sig_conf_name = joinpath("cache", cached_name * "_siggen_config.txt") + + # if field exists, no need to construct config, just read cached simulation -> return + # if config exists, no need to construct again -> return -> maybe just the cached field? not sure + # of course unless asked to overwrite, then construct new config + if ( isfile(joinpath("cache", fieldgen_names["wp_name"])) || isfile(sig_conf_name) ) && !overwrite + return sig_conf_name, fieldgen_names["wp_name"] + end + + @info "Constructing Fieldgen/Siggen configuration file" + + if !ispath("cache") mkpath("cache") end + + + # ---------------------------------------------------------------------------- + # construct config + # ---------------------------------------------------------------------------- + + println("...detector geometry") + + # construct siggen geometry based on LEGEND metadata JSON + environment settings + siggen_geom = meta2siggen(meta, env) + + # collect geometry lines for siggen config file + detector_lines = Vector{String}() + push!(detector_lines, "# detector related parameters") + for (param, value) in siggen_geom + push!(detector_lines, "$param $value") + end + + println("...fieldgen configuration") + + # ---------------------------------------------------------------------------- + # drift velocity, field, wp + # ---------------------------------------------------------------------------- + + # drift velocity input (path to file given by user) + drift_file_path = simulator.drift_vel + # use default if not given by user + if simulator.drift_vel == "" + @warn "No drift velocity input given; using default" + drift_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "drift_vel_tcorr.tab") + end + + # copy to cache + cache_drift_file = cached_name * "_" * basename(drift_file_path) + cp(drift_file_path, joinpath("cache", cache_drift_file), force=true) # later think what to do here + + # add to field fieldgen_names + fieldgen_names["drift_name"] = cache_drift_file + + ## commenting this out: produced a bug + ## (empty file created, next time running empty file read and no field at point error) + # for name in ["field_name", "wp_name"] + # path = joinpath("cache", fieldgen_names[name]) + # if !isfile(path) + # touch(path) + # end + # end + + # add corresponding lines to existing detector geometry lines + push!(detector_lines, "") + push!(detector_lines, "# fieldgen filenames") + for (param, value) in fieldgen_names + push!(detector_lines, "$param $value") + end + + # ---------------------------------------------------------------------------- + # extra settings + # ---------------------------------------------------------------------------- + + # Note: this is similar but not identical to the part above because of the cache folder difference + # -> loop? copy-paste is not nice... + + # read lines from user input for fieldgen + config_file_path = simulator.fieldgen_config + if simulator.fieldgen_config == "" + @warn "No extra fieldgen settings given; using default" + config_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "fieldgen_settings.txt") + end + # copy to cache -> maybe not needed, are appended in main settings? + cache_config_file = joinpath("cache", cached_name * "_" * basename(config_file_path)) + cp(config_file_path, cache_config_file, force=true) # later think what to do here + + fieldgen_lines = readlines(open(cache_config_file)) + fieldgen_lines = replace.(fieldgen_lines, "\t" => " ") + + # unite constructed detector geometry and fieldgen settings + total_lines = vcat(detector_lines, [""], fieldgen_lines) + + # ---------------------------------------------------------------------------- + # final siggen config file + # ---------------------------------------------------------------------------- + + # write a siggen/fieldgen config file + + writedlm(sig_conf_name, total_lines) + println("...fieldgen/siggen config written to $sig_conf_name") + + # return resulting fieldgen/siggen config name + # (currently kinda redundant but that's because we have no idea what's gonna happen later) + sig_conf_name, fieldgen_names["wp_name"] +end + + +""" + meta2siggen(meta, env) + +PropDict, Environment -> Dict + +Construct siggen type config in a form of Dict based on + geometry as given in LEGEND metadata + and environment settings given in . +""" +function meta2siggen(meta::PropDict, env::Environment) + # !!!! this is copy paste from legend_detector_to_SSD !! + # ToDo: rewrite to avoid copy-paste + + # ----------------------------------------------------------- + # DL thickness (quickfix) + + dl_thickness_in_mm = + if env.dl == "vendor" + dl_vendor = meta.characterization.manufacturer.dl_thickness_in_mm + if dl_vendor isa PropDicts.MissingProperty + throw(ArgumentError("No dead layer measurement provided by vendor! Please provide value or skip (default 0).")) + end + dl_vendor + else + env.dl + end + + @info "DL = $(dl_thickness_in_mm)mm" + + # ----------------------------------------------------------- + # operating voltage + + # sometimes recV is null -> in this case, complain and ask to provide (default env.opV is 0 if none given) + operation_voltage = + if env.operating_voltage == 0 + @warn "You did not provide operating voltage in environment settings -> taking recommended voltage from metadata" + recV = meta.characterization.l200_site.recommended_voltage_in_V + if isnothing(recV) + error("Metadata does not provide recommended voltage (null). Please provide operating voltage in settings") + end + recV + else + env.operating_voltage + end + + @info "Simulating at $(operation_voltage)V" + + # ----------------------------------------------------------- + # parameters common for all geometries + + siggen_geometry = Dict( + # -- crystal + "xtal_length" => meta.geometry.height_in_mm, + "xtal_radius" => meta.geometry.radius_in_mm, + # -- p+ contact + "pc_length" => meta.geometry.pp_contact.depth_in_mm, + "pc_radius" => meta.geometry.pp_contact.radius_in_mm, + # -- tapering + + # comment from David Radford: + + # The bottom taper is always at 45 degrees, but the top outer taper is not. + # It's width/angle is defined using either the outer_taper_width parameter or the taper_angle parameter. + # Likewise, the inner taper width/angle is defined using either the inner_taper_width parameter or the taper_angle parameter. + # The inner_taper_length can be anything from zero to the hole_length + + # size of 45-degree taper at bottom of ORTEC-type crystal (for r=z) + # -> for now, ignore bottom taper if not 45 + # -> works for 1) the ones with actual 45 taper, 2) the ones bulletized (saved as 45 deg taper with bullet radius) + # few that have proper non-45 taper - not possible. ToDo: put some warning + "bottom_taper_length" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm : 0, + # current metadata format: always angle + "bottom_taper_width" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm * tan(deg2rad(meta.geometry.taper.bottom.angle_in_deg)) : 0, + # z-length of outside taper for inverted-coax style -> Mariia: you mean top taper here? + "outer_taper_length" => meta.geometry.taper.top.height_in_mm, + "outer_taper_width" => meta.geometry.taper.top.height_in_mm * tan(deg2rad(meta.geometry.taper.top.angle_in_deg)), + # taper angle in degrees, for inner or outer taper -> how? why are they the same? how about borehole? not using this parameter... + # "taper_angle" => meta.geometry.taper.top.outer.angle_in_deg, + # "taper_angle" => meta.geometry.taper.top.inner.angle_in_deg, + # with this setting, it somehow takes taper angle: 45.000038 anyway -> force to zero? + # "taper_angle" => 0, + # depth of full-charge-collection boundary for Li contact + "Li_thickness" => dl_thickness_in_mm, + + # configuration for mjd_fieldgen (calculates electric fields & weighing potentials) + # detector bias for fieldgen, in Volts + "xtal_HV" => operation_voltage, + + # configuration for signal calculation + # crystal temperature in Kelvin + "xtal_temp" => env.crystal_temperature + ) + + # ----------------------------------------------------------- + # non common parameters + + # -- groove (for non PPCs) + if haskey(meta.geometry, :groove) + siggen_geometry["wrap_around_radius"] = meta.geometry.groove.radius_in_mm.outer + siggen_geometry["ditch_depth"] = meta.geometry.groove.depth_in_mm + siggen_geometry["ditch_thickness"] = meta.geometry.groove.radius_in_mm.outer - meta.geometry.groove.radius_in_mm.inner + end + + # -- borehole (for ICPC and Coax) + if haskey(meta.geometry, :borehole) + # length of hole, for inverted-coax style + siggen_geometry["hole_length"] = meta.geometry.borehole.depth_in_mm + # radius of hole, for inverted-coax style + siggen_geometry["hole_radius"] = meta.geometry.borehole.radius_in_mm + # z-length of inside (hole) taper for inverted-coax style + siggen_geometry["inner_taper_length"] = meta.geometry.taper.borehole.height_in_mm + siggen_geometry["inner_taper_width"] = meta.geometry.taper.borehole.height_in_mm * tan(deg2rad(meta.geometry.taper.borehole.angle_in_deg)) + end + + siggen_geometry +end diff --git a/ext/mjdsiggen_utils.jl b/ext/mjdsiggen_utils.jl new file mode 100644 index 0000000..9c1f07a --- /dev/null +++ b/ext/mjdsiggen_utils.jl @@ -0,0 +1,30 @@ +using MJDSigGen: SigGenSetup +function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) + E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); + T = eltype(E_pot) + r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid + z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid + ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) + ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) + ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) + axs = (ax1, ax2, ax3) + grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) + data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); + data[:] = E_pot'[:]; + ElectricPotential(data, grid) +end + +function SolidStateDetectors.WeightingPotential(sim::SigGenSetup) + E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); + T = eltype(W_pot) + r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid + z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid + ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) + ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) + ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) + axs = (ax1, ax2, ax3) + grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) + data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); + data[:] = W_pot'[:]; + WeightingPotential(data, grid) +end \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml index bfaa967..deca42a 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,8 @@ LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0" LegendTestData = "33d6da08-6349-5f7c-b5a4-6ff4d8795aaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +MJDSigGen= "0b925020-e659-485c-98dd-189633cb4be4" [compat] -LegendTestData = "0.2.8" \ No newline at end of file +LegendTestData = "0.2.8" +MJDSigGen = "0.2" \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 08d3107..5ee6203 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,11 +3,12 @@ using Test using LegendGeSim +using MJDSigGen using LegendHDF5IO, HDF5 using LegendTestData using Unitful - +isdir("cache") && rm("cache", force=true, recursive=true) @testset "LegendGeSim: Simulate capacitances of test detectors" begin # note: this part also tests dict settings input as opposed to json file # as well as impurity input @@ -147,5 +148,6 @@ end end # Delete the cached files - rm("cache", force=true, recursive=true) + #rm("cache", force=true, recursive=true) + end From 6229eb1c77fae59cbed44629443f009eb260151a Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 11 Nov 2024 16:53:30 -0500 Subject: [PATCH 2/9] Reorder package imports --- src/LegendGeSim.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/LegendGeSim.jl b/src/LegendGeSim.jl index f42c354..e6e6933 100644 --- a/src/LegendGeSim.jl +++ b/src/LegendGeSim.jl @@ -9,6 +9,7 @@ Template for Julia packages. """ module LegendGeSim +# Registered packages using ArgCheck using ArraysOfArrays using CurveFit @@ -19,12 +20,7 @@ using ElasticArrays using HDF5 using IntervalSets using JSON -using LegendDataManagement -using LegendDataTypes -using LegendHDF5IO -using LegendTextIO using LinearAlgebra -using MJDSigGen using Parameters using Polynomials using PropDicts @@ -45,6 +41,14 @@ using Tables using TypedTables using Unitful +# Unregistered packages +using LegendDataManagement +using LegendDataTypes +using LegendHDF5IO +using LegendTextIO +using MJDSigGen + + import LsqFit using SolidStateDetectors: ConstructiveSolidGeometry as CSG From be830003b488349e024b968f40012a447ea4163a Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Mon, 11 Nov 2024 16:58:43 -0500 Subject: [PATCH 3/9] Adapt LegendGeSim to changes in SolidStateDetectors@0.10.3 --- Project.toml | 2 +- src/legend_detector_to_ssd.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 18fd15e..3b3a759 100644 --- a/Project.toml +++ b/Project.toml @@ -67,7 +67,7 @@ Random123 = "1.2" RecipesBase = "1" Requires = "0.5, 1" SignalAnalysis = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9" -SolidStateDetectors = "0.9.4 - 0.10" +SolidStateDetectors = "0.10.3" StaticArrays = "0.12, 1" Statistics = "1" StatsBase = "0.32, 0.33, 0.34" diff --git a/src/legend_detector_to_ssd.jl b/src/legend_detector_to_ssd.jl index 4af0989..183ca7f 100644 --- a/src/legend_detector_to_ssd.jl +++ b/src/legend_detector_to_ssd.jl @@ -63,7 +63,7 @@ function LEGEND_SolidStateDetector(::Type{T}, meta::PropDict, env::Environment, # "old" semiconductor sc = detector.semiconductor # new semiconductor with our temperature - semiconductor = SolidStateDetectors.Semiconductor(temperature, sc.material, sc.impurity_density_model, sc.charge_drift_model, sc.geometry) + semiconductor = SolidStateDetectors.Semiconductor(temperature, sc.material, sc.impurity_density_model, sc.charge_drift_model, sc.charge_trapping_model, sc.geometry) # new detector with new semiconductor detector = SolidStateDetector{T}(detector.name, semiconductor, detector.contacts, detector.passives, detector.virtual_drift_volumes) From 2e3c70670128a13b483bf125551faefff9af523b Mon Sep 17 00:00:00 2001 From: Shailesh Giri Date: Mon, 18 Nov 2024 11:35:52 -0500 Subject: [PATCH 4/9] ready for review --- Project.toml | 18 ++- ext/LegendGeSimExt.jl | 27 ---- ext/legend_detector_to_siggen.jl | 263 ------------------------------- ext/mjdsiggen_utils.jl | 30 ---- src/LegendGeSim.jl | 4 - src/detector.jl | 54 ------- src/impurity.jl | 46 ------ src/legend_detector_to_siggen.jl | 259 ------------------------------ src/mjdsiggen_utils.jl | 29 ---- src/pss.jl | 195 +++-------------------- src/temp_utils.jl | 7 +- 11 files changed, 39 insertions(+), 893 deletions(-) delete mode 100644 ext/legend_detector_to_siggen.jl delete mode 100644 ext/mjdsiggen_utils.jl delete mode 100644 src/legend_detector_to_siggen.jl delete mode 100644 src/mjdsiggen_utils.jl diff --git a/Project.toml b/Project.toml index 3b3a759..6bb90f8 100644 --- a/Project.toml +++ b/Project.toml @@ -40,12 +40,18 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" +[weakdeps] +MJDSigGen = "0b925020-e659-485c-98dd-189633cb4be4" + +[extensions] +LegendGeSimExt = "MJDSigGen" + [compat] ArgCheck = "1, 2" ArraysOfArrays = "0.5, 0.6" CurveFit = "0.3, 0.5, 0.6" DSP = "0.7" -DelimitedFiles = "1" +DelimitedFiles = " 1" Distributions = "0.23, 0.24, 0.25" ElasticArrays = "1" HDF5 = "0.15, 0.16, 0.17" @@ -57,22 +63,22 @@ LegendHDF5IO = "0.1" LegendTextIO = "0.1" LsqFit = "0.13, 0.14, 0.15" MJDSigGen = "0.2" -Parameters = "0.12" +Parameters = "0.12, 0.13" Polynomials = "2, 3, 4" PropDicts = "0.2" RadiationDetectorDSP = "0.2.1" RadiationDetectorSignals = "0.3" RadiationSpectra = "0.5" -Random123 = "1.2" +Random123 = "1.2, 1.7" # Compatible with 1.7? RecipesBase = "1" Requires = "0.5, 1" SignalAnalysis = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9" -SolidStateDetectors = "0.10.3" +SolidStateDetectors = "0.10.5" StaticArrays = "0.12, 1" Statistics = "1" StatsBase = "0.32, 0.33, 0.34" StructArrays = "0.5, 0.6" Tables = "1" -TypedTables = "1.2" +TypedTables = "1.2, 1.4" #Compatible with 1.4.6? Unitful = "1" -julia = "1.10" +julia = "1.11" diff --git a/ext/LegendGeSimExt.jl b/ext/LegendGeSimExt.jl index f5083bb..f21d147 100644 --- a/ext/LegendGeSimExt.jl +++ b/ext/LegendGeSimExt.jl @@ -1,19 +1,10 @@ module LegendGeSimExt -#using MJDSigGen: SiggenSimulator -#@info "Error here" using LegendGeSim -#@info "Error" using PropDicts -#@info "Error_prop" using LegendGeSim -#@info "Error_Lgsim" using LegendGeSim: Environment -#@info "le_env" -#using MJDSigGen: SiggenSimulator using MJDSigGen using MJDSigGen: SigGenSetup -#@info "Error_mjd" -#using LegendGeSim: PSSimulator using LegendGeSim: SiggenSimulator using ArgCheck using ArraysOfArrays @@ -54,20 +45,10 @@ import LsqFit using SolidStateDetectors: ConstructiveSolidGeometry as CSG -#abstract type PSSimulator end - - -#include("legend_detector_to_siggen.jl") -#@info "Here_then" -#include("mjdsiggen_utils.jl") -#@info "Error_after_files" - - #################################### ### FieldGen #################################### -@info "27" """ SiggeSimulator(sim_conf::PropDict) @@ -166,7 +147,6 @@ else #...do nothing, siggen will later read the files based on the generated conifg @info "Cached simulation found. Reading cached fields from $fieldgen_wp_name" end -@info "73" # a SigGenSetup object SigGenSetup(siggen_config_name) end @@ -215,7 +195,6 @@ function LegendGeSim.impurity_density_model(meta::PropDict, crystal_metadata::Pr imp_filename, L - det_z0 end -@info "123" """ Simulation method: siggen """ @@ -224,9 +203,7 @@ Simulation method: siggen function LegendGeSim.simulate_waveforms(stp_events::Table, detector::SigGenSetup, ::SiggenSimulator) T = Float32 # This should be somehow defined and be passed properly - println("127") @info("~.~.~ Siggen") - @info "129" nevents = length(stp_events) wf_array = Array{RDWaveform}(undef, nevents) @@ -272,7 +249,6 @@ function LegendGeSim.simulate_waveforms(stp_events::Table, detector::SigGenSetup end -#@info "252" """ simulate_wf(pos, edep, siggen_setup) @@ -310,7 +286,6 @@ function LegendGeSim.capacitance_matrix(sim::SigGenSetup) [c -c; -c c] end #export capacitance_matrix -#@info "288" #legend_detector_to_siggen function siggen_config(meta::PropDict, env::Environment, simulator::SiggenSimulator; overwrite::Bool=false) # ---------------------------------------------------------------------------- @@ -439,7 +414,6 @@ function siggen_config(meta::PropDict, env::Environment, simulator::SiggenSimula # (currently kinda redundant but that's because we have no idea what's gonna happen later) sig_conf_name, fieldgen_names["wp_name"] end -#@info "417" """ meta2siggen(meta, env) @@ -559,7 +533,6 @@ function meta2siggen(meta::PropDict, env::Environment) end #mjdsiggen_utils.jl -#@info "537" function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); diff --git a/ext/legend_detector_to_siggen.jl b/ext/legend_detector_to_siggen.jl deleted file mode 100644 index 7bf9617..0000000 --- a/ext/legend_detector_to_siggen.jl +++ /dev/null @@ -1,263 +0,0 @@ - -using PropDicts - -""" - siggen_config(meta, env, siggen_sim, config_name) - -PropDict, Environment, SiggenSimulator, String -> String, String - -Construct Siggen/Fieldgen config based on - geometry as given in LEGEND metadata , - environment settings given in , - and fieldgen settings contained in . - -The resulting siggen config will be cached based on given -The function returns the path to the cached siggen config file, and - the relative path to the fieldgen WP file (to check if already exists later) -""" - -function siggen_config(meta::PropDict, env::Environment, simulator::SiggenSimulator; overwrite::Bool=false) - # ---------------------------------------------------------------------------- - # set up & check if already present in cache - # ---------------------------------------------------------------------------- - - # if no cached name was given by user, make a temporaty name - needed for fieldgen to read from file - # also if no cached name was given means definitely doing from scratch - cached_name = simulator.cached_name - if simulator.cached_name == "" - cached_name = "tmp" - overwrite = true - end - - # append detector name - cached_name = meta.name * "_" * cached_name - - # filenames for fieldgen input/output - # ! no need to add cache/ folder because siggen does it by itself... - fieldgen_names = Dict( - "field_name" => cached_name * "_fieldgen_EV.dat", - "wp_name" => cached_name * "_fieldgen_WP.dat" - ) - - # siggen config filename - sig_conf_name = joinpath("cache", cached_name * "_siggen_config.txt") - - # if field exists, no need to construct config, just read cached simulation -> return - # if config exists, no need to construct again -> return -> maybe just the cached field? not sure - # of course unless asked to overwrite, then construct new config - if ( isfile(joinpath("cache", fieldgen_names["wp_name"])) || isfile(sig_conf_name) ) && !overwrite - return sig_conf_name, fieldgen_names["wp_name"] - end - - @info "Constructing Fieldgen/Siggen configuration file" - - if !ispath("cache") mkpath("cache") end - - - # ---------------------------------------------------------------------------- - # construct config - # ---------------------------------------------------------------------------- - - println("...detector geometry") - - # construct siggen geometry based on LEGEND metadata JSON + environment settings - siggen_geom = meta2siggen(meta, env) - - # collect geometry lines for siggen config file - detector_lines = Vector{String}() - push!(detector_lines, "# detector related parameters") - for (param, value) in siggen_geom - push!(detector_lines, "$param $value") - end - - println("...fieldgen configuration") - - # ---------------------------------------------------------------------------- - # drift velocity, field, wp - # ---------------------------------------------------------------------------- - - # drift velocity input (path to file given by user) - drift_file_path = simulator.drift_vel - # use default if not given by user - if simulator.drift_vel == "" - @warn "No drift velocity input given; using default" - drift_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "drift_vel_tcorr.tab") - end - - # copy to cache - cache_drift_file = cached_name * "_" * basename(drift_file_path) - cp(drift_file_path, joinpath("cache", cache_drift_file), force=true) # later think what to do here - - # add to field fieldgen_names - fieldgen_names["drift_name"] = cache_drift_file - - ## commenting this out: produced a bug - ## (empty file created, next time running empty file read and no field at point error) - # for name in ["field_name", "wp_name"] - # path = joinpath("cache", fieldgen_names[name]) - # if !isfile(path) - # touch(path) - # end - # end - - # add corresponding lines to existing detector geometry lines - push!(detector_lines, "") - push!(detector_lines, "# fieldgen filenames") - for (param, value) in fieldgen_names - push!(detector_lines, "$param $value") - end - - # ---------------------------------------------------------------------------- - # extra settings - # ---------------------------------------------------------------------------- - - # Note: this is similar but not identical to the part above because of the cache folder difference - # -> loop? copy-paste is not nice... - - # read lines from user input for fieldgen - config_file_path = simulator.fieldgen_config - if simulator.fieldgen_config == "" - @warn "No extra fieldgen settings given; using default" - config_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "fieldgen_settings.txt") - end - # copy to cache -> maybe not needed, are appended in main settings? - cache_config_file = joinpath("cache", cached_name * "_" * basename(config_file_path)) - cp(config_file_path, cache_config_file, force=true) # later think what to do here - - fieldgen_lines = readlines(open(cache_config_file)) - fieldgen_lines = replace.(fieldgen_lines, "\t" => " ") - - # unite constructed detector geometry and fieldgen settings - total_lines = vcat(detector_lines, [""], fieldgen_lines) - - # ---------------------------------------------------------------------------- - # final siggen config file - # ---------------------------------------------------------------------------- - - # write a siggen/fieldgen config file - - writedlm(sig_conf_name, total_lines) - println("...fieldgen/siggen config written to $sig_conf_name") - - # return resulting fieldgen/siggen config name - # (currently kinda redundant but that's because we have no idea what's gonna happen later) - sig_conf_name, fieldgen_names["wp_name"] -end - - -""" - meta2siggen(meta, env) - -PropDict, Environment -> Dict - -Construct siggen type config in a form of Dict based on - geometry as given in LEGEND metadata - and environment settings given in . -""" -function meta2siggen(meta::PropDict, env::Environment) - # !!!! this is copy paste from legend_detector_to_SSD !! - # ToDo: rewrite to avoid copy-paste - - # ----------------------------------------------------------- - # DL thickness (quickfix) - - dl_thickness_in_mm = - if env.dl == "vendor" - dl_vendor = meta.characterization.manufacturer.dl_thickness_in_mm - if dl_vendor isa PropDicts.MissingProperty - throw(ArgumentError("No dead layer measurement provided by vendor! Please provide value or skip (default 0).")) - end - dl_vendor - else - env.dl - end - - @info "DL = $(dl_thickness_in_mm)mm" - - # ----------------------------------------------------------- - # operating voltage - - # sometimes recV is null -> in this case, complain and ask to provide (default env.opV is 0 if none given) - operation_voltage = - if env.operating_voltage == 0 - @warn "You did not provide operating voltage in environment settings -> taking recommended voltage from metadata" - recV = meta.characterization.l200_site.recommended_voltage_in_V - if isnothing(recV) - error("Metadata does not provide recommended voltage (null). Please provide operating voltage in settings") - end - recV - else - env.operating_voltage - end - - @info "Simulating at $(operation_voltage)V" - - # ----------------------------------------------------------- - # parameters common for all geometries - - siggen_geometry = Dict( - # -- crystal - "xtal_length" => meta.geometry.height_in_mm, - "xtal_radius" => meta.geometry.radius_in_mm, - # -- p+ contact - "pc_length" => meta.geometry.pp_contact.depth_in_mm, - "pc_radius" => meta.geometry.pp_contact.radius_in_mm, - # -- tapering - - # comment from David Radford: - - # The bottom taper is always at 45 degrees, but the top outer taper is not. - # It's width/angle is defined using either the outer_taper_width parameter or the taper_angle parameter. - # Likewise, the inner taper width/angle is defined using either the inner_taper_width parameter or the taper_angle parameter. - # The inner_taper_length can be anything from zero to the hole_length - - # size of 45-degree taper at bottom of ORTEC-type crystal (for r=z) - # -> for now, ignore bottom taper if not 45 - # -> works for 1) the ones with actual 45 taper, 2) the ones bulletized (saved as 45 deg taper with bullet radius) - # few that have proper non-45 taper - not possible. ToDo: put some warning - "bottom_taper_length" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm : 0, - # current metadata format: always angle - "bottom_taper_width" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm * tan(deg2rad(meta.geometry.taper.bottom.angle_in_deg)) : 0, - # z-length of outside taper for inverted-coax style -> Mariia: you mean top taper here? - "outer_taper_length" => meta.geometry.taper.top.height_in_mm, - "outer_taper_width" => meta.geometry.taper.top.height_in_mm * tan(deg2rad(meta.geometry.taper.top.angle_in_deg)), - # taper angle in degrees, for inner or outer taper -> how? why are they the same? how about borehole? not using this parameter... - # "taper_angle" => meta.geometry.taper.top.outer.angle_in_deg, - # "taper_angle" => meta.geometry.taper.top.inner.angle_in_deg, - # with this setting, it somehow takes taper angle: 45.000038 anyway -> force to zero? - # "taper_angle" => 0, - # depth of full-charge-collection boundary for Li contact - "Li_thickness" => dl_thickness_in_mm, - - # configuration for mjd_fieldgen (calculates electric fields & weighing potentials) - # detector bias for fieldgen, in Volts - "xtal_HV" => operation_voltage, - - # configuration for signal calculation - # crystal temperature in Kelvin - "xtal_temp" => env.crystal_temperature - ) - - # ----------------------------------------------------------- - # non common parameters - - # -- groove (for non PPCs) - if haskey(meta.geometry, :groove) - siggen_geometry["wrap_around_radius"] = meta.geometry.groove.radius_in_mm.outer - siggen_geometry["ditch_depth"] = meta.geometry.groove.depth_in_mm - siggen_geometry["ditch_thickness"] = meta.geometry.groove.radius_in_mm.outer - meta.geometry.groove.radius_in_mm.inner - end - - # -- borehole (for ICPC and Coax) - if haskey(meta.geometry, :borehole) - # length of hole, for inverted-coax style - siggen_geometry["hole_length"] = meta.geometry.borehole.depth_in_mm - # radius of hole, for inverted-coax style - siggen_geometry["hole_radius"] = meta.geometry.borehole.radius_in_mm - # z-length of inside (hole) taper for inverted-coax style - siggen_geometry["inner_taper_length"] = meta.geometry.taper.borehole.height_in_mm - siggen_geometry["inner_taper_width"] = meta.geometry.taper.borehole.height_in_mm * tan(deg2rad(meta.geometry.taper.borehole.angle_in_deg)) - end - - siggen_geometry -end diff --git a/ext/mjdsiggen_utils.jl b/ext/mjdsiggen_utils.jl deleted file mode 100644 index 9c1f07a..0000000 --- a/ext/mjdsiggen_utils.jl +++ /dev/null @@ -1,30 +0,0 @@ -using MJDSigGen: SigGenSetup -function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) - E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); - T = eltype(E_pot) - r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid - z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid - ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) - ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) - ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) - axs = (ax1, ax2, ax3) - grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) - data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); - data[:] = E_pot'[:]; - ElectricPotential(data, grid) -end - -function SolidStateDetectors.WeightingPotential(sim::SigGenSetup) - E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); - T = eltype(W_pot) - r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid - z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid - ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) - ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) - ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) - axs = (ax1, ax2, ax3) - grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) - data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); - data[:] = W_pot'[:]; - WeightingPotential(data, grid) -end \ No newline at end of file diff --git a/src/LegendGeSim.jl b/src/LegendGeSim.jl index e6e6933..3efb502 100644 --- a/src/LegendGeSim.jl +++ b/src/LegendGeSim.jl @@ -46,8 +46,6 @@ using LegendDataManagement using LegendDataTypes using LegendHDF5IO using LegendTextIO -using MJDSigGen - import LsqFit @@ -65,7 +63,6 @@ include("environment.jl") include("pss.jl") include("impurity.jl") include("legend_detector_to_ssd.jl") -include("legend_detector_to_siggen.jl") include("detector.jl") include("preamp.jl") @@ -86,7 +83,6 @@ include("pet_to_raw.jl") include("waveform_utils.jl") include("temp_utils.jl") -include("mjdsiggen_utils.jl") diff --git a/src/detector.jl b/src/detector.jl index 285561e..3e22c1a 100644 --- a/src/detector.jl +++ b/src/detector.jl @@ -28,60 +28,6 @@ end -#################################### -### FieldGen -#################################### - -""" - simulate_detector(det_meta, env, config_name, ps_simulator) - -PropDict, Environment, AbstractString, Siggen -> SigGenSetup - -Construct a fieldgen/siggen configuration file - according to geometry as given in LEGEND metadata - and environment settings given in . -Look up fieldgen generated electric field and weighting potential files - based on given , and if not found, call fieldgen. -""" -function simulate_detector(det_meta::PropDict, env::Environment, simulator::SiggenSimulator; - overwrite::Bool = false) - - # if no cached name given, force simulation from scratch (or else might read past "tmp" file) - if simulator.cached_name == "" - overwrite = true - end - - # returns the name of the resulting siggen config file - # and the name of (already or to be) generated weighting potential file - # ToDo: don't create if already present already at this stage -> check in siggen_config() if name exists and just return name - siggen_config_name, fieldgen_wp_name = siggen_config(det_meta, env, simulator; overwrite) - fieldgen_wp_name = joinpath("cache", fieldgen_wp_name) - - # if the WP file with such a name does not exist yet or want to overwrite it... - if !isfile(fieldgen_wp_name) || overwrite - imp_filename, offset = - if simulator.crystal_metadata_path != "" - # create impurity input on the fly - impurity_density_model(det_meta, simulator.crystal_metadata_path, simulator) - else - # user provided .dat file and corresponding offset - # (or if nothing is provided this will be "" and -1) - simulator.impurity_profile, simulator.offset_in_mm - end - #...call fieldgen -> will write the wp file - @info "_|~|_|~|_|~|_ Fieldgen simulation" - fieldgen(siggen_config_name; impurity_profile=imp_filename, offset_mm=offset) - @info "_|~|_|~|_|~|_ Fieldgen simulation complete" - else - #...do nothing, siggen will later read the files based on the generated conifg - @info "Cached simulation found. Reading cached fields from $fieldgen_wp_name" - end - - # a SigGenSetup object - SigGenSetup(siggen_config_name) -end - - #################################### ### SSD #################################### diff --git a/src/impurity.jl b/src/impurity.jl index 6fd577b..a4d67a3 100644 --- a/src/impurity.jl +++ b/src/impurity.jl @@ -142,52 +142,6 @@ function impurity_density_model(meta::PropDict, crystal_metadata::PropDict, fit_ end -## ---------- fieldgen - - -""" -Convert fit parameters from crystal metadata units (1e9 e/cm^3 VS mm) into Siggen units (1e10 e/cm^3 VS mm) -and create a .dat file (unformatted stream of Float32) to be later used when calling fieldgen(). -""" -function impurity_density_model(meta::PropDict, crystal_metadata::PropDict, fit_param::Vector{Float64}, ::SiggenSimulator) - - a,b,c,tau = fit_param - - T = Float32 - dist = T.(crystal_metadata.impurity_measurements.distance_from_seed_end_mm) - # sometimes might not have measurements at seed end? but last measurement can be taken as crystal length - L = dist[end] #mm - - # points in crystal axis from 0 to L with step of 1 mm - # always put 200 points to make all .dat files same length - imp = zeros(Float32, 200) # impurities in 10^10 cm^-3 - - # length is L+1 including both ends - Npoints = Int(ceil(L))+1 # step of 1 mm - zpoints = LinRange(0, L, Npoints) - for i in 1:Npoints - # divide by 10 to convert from 1e9 e/cm^3 to 1e10 e/cm^3 - # invert the order of the values because with siggen z=0 is tail (dirtiest part) - # i.e. z = L in our system (Mirion measurements) - imp[Npoints-i+1] = -(a + b * zpoints[i] + c * exp((zpoints[i] - L)/tau)) / 10. - end - - crystal_name = first(meta.name, length(meta.name)-1) - imp_filename = joinpath("cache", crystal_name*".dat") - open(imp_filename, "w") do io - write(io, imp) - end - @info "Impurity profile for siggen saved in $(imp_filename)" - - # return name of impurity profile file and offset - # offset from seed end - det_z0 = crystal_metadata.slices[Symbol(meta.production.slice)].detector_offset_in_mm - - # siggen offset is from tail i.e. at L offset = 0 - imp_filename, L - det_z0 -end - - #### OUTDATED, providing linear impurity profile with quadratic correction inside of siggen config # """ # fieldgen_impurity(meta) diff --git a/src/legend_detector_to_siggen.jl b/src/legend_detector_to_siggen.jl deleted file mode 100644 index 467f144..0000000 --- a/src/legend_detector_to_siggen.jl +++ /dev/null @@ -1,259 +0,0 @@ -""" - siggen_config(meta, env, siggen_sim, config_name) - -PropDict, Environment, SiggenSimulator, String -> String, String - -Construct Siggen/Fieldgen config based on - geometry as given in LEGEND metadata , - environment settings given in , - and fieldgen settings contained in . - -The resulting siggen config will be cached based on given -The function returns the path to the cached siggen config file, and - the relative path to the fieldgen WP file (to check if already exists later) -""" -function siggen_config(meta::PropDict, env::Environment, simulator::SiggenSimulator; overwrite::Bool=false) - # ---------------------------------------------------------------------------- - # set up & check if already present in cache - # ---------------------------------------------------------------------------- - - # if no cached name was given by user, make a temporaty name - needed for fieldgen to read from file - # also if no cached name was given means definitely doing from scratch - cached_name = simulator.cached_name - if simulator.cached_name == "" - cached_name = "tmp" - overwrite = true - end - - # append detector name - cached_name = meta.name * "_" * cached_name - - # filenames for fieldgen input/output - # ! no need to add cache/ folder because siggen does it by itself... - fieldgen_names = Dict( - "field_name" => cached_name * "_fieldgen_EV.dat", - "wp_name" => cached_name * "_fieldgen_WP.dat" - ) - - # siggen config filename - sig_conf_name = joinpath("cache", cached_name * "_siggen_config.txt") - - # if field exists, no need to construct config, just read cached simulation -> return - # if config exists, no need to construct again -> return -> maybe just the cached field? not sure - # of course unless asked to overwrite, then construct new config - if ( isfile(joinpath("cache", fieldgen_names["wp_name"])) || isfile(sig_conf_name) ) && !overwrite - return sig_conf_name, fieldgen_names["wp_name"] - end - - @info "Constructing Fieldgen/Siggen configuration file" - - if !ispath("cache") mkpath("cache") end - - - # ---------------------------------------------------------------------------- - # construct config - # ---------------------------------------------------------------------------- - - println("...detector geometry") - - # construct siggen geometry based on LEGEND metadata JSON + environment settings - siggen_geom = meta2siggen(meta, env) - - # collect geometry lines for siggen config file - detector_lines = Vector{String}() - push!(detector_lines, "# detector related parameters") - for (param, value) in siggen_geom - push!(detector_lines, "$param $value") - end - - println("...fieldgen configuration") - - # ---------------------------------------------------------------------------- - # drift velocity, field, wp - # ---------------------------------------------------------------------------- - - # drift velocity input (path to file given by user) - drift_file_path = simulator.drift_vel - # use default if not given by user - if simulator.drift_vel == "" - @warn "No drift velocity input given; using default" - drift_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "drift_vel_tcorr.tab") - end - - # copy to cache - cache_drift_file = cached_name * "_" * basename(drift_file_path) - cp(drift_file_path, joinpath("cache", cache_drift_file), force=true) # later think what to do here - - # add to field fieldgen_names - fieldgen_names["drift_name"] = cache_drift_file - - ## commenting this out: produced a bug - ## (empty file created, next time running empty file read and no field at point error) - # for name in ["field_name", "wp_name"] - # path = joinpath("cache", fieldgen_names[name]) - # if !isfile(path) - # touch(path) - # end - # end - - # add corresponding lines to existing detector geometry lines - push!(detector_lines, "") - push!(detector_lines, "# fieldgen filenames") - for (param, value) in fieldgen_names - push!(detector_lines, "$param $value") - end - - # ---------------------------------------------------------------------------- - # extra settings - # ---------------------------------------------------------------------------- - - # Note: this is similar but not identical to the part above because of the cache folder difference - # -> loop? copy-paste is not nice... - - # read lines from user input for fieldgen - config_file_path = simulator.fieldgen_config - if simulator.fieldgen_config == "" - @warn "No extra fieldgen settings given; using default" - config_file_path = joinpath(dirname(@__DIR__), "examples", "configs", "fieldgen_settings.txt") - end - # copy to cache -> maybe not needed, are appended in main settings? - cache_config_file = joinpath("cache", cached_name * "_" * basename(config_file_path)) - cp(config_file_path, cache_config_file, force=true) # later think what to do here - - fieldgen_lines = readlines(open(cache_config_file)) - fieldgen_lines = replace.(fieldgen_lines, "\t" => " ") - - # unite constructed detector geometry and fieldgen settings - total_lines = vcat(detector_lines, [""], fieldgen_lines) - - # ---------------------------------------------------------------------------- - # final siggen config file - # ---------------------------------------------------------------------------- - - # write a siggen/fieldgen config file - - writedlm(sig_conf_name, total_lines) - println("...fieldgen/siggen config written to $sig_conf_name") - - # return resulting fieldgen/siggen config name - # (currently kinda redundant but that's because we have no idea what's gonna happen later) - sig_conf_name, fieldgen_names["wp_name"] -end - - -""" - meta2siggen(meta, env) - -PropDict, Environment -> Dict - -Construct siggen type config in a form of Dict based on - geometry as given in LEGEND metadata - and environment settings given in . -""" -function meta2siggen(meta::PropDict, env::Environment) - # !!!! this is copy paste from legend_detector_to_SSD !! - # ToDo: rewrite to avoid copy-paste - - # ----------------------------------------------------------- - # DL thickness (quickfix) - - dl_thickness_in_mm = - if env.dl == "vendor" - dl_vendor = meta.characterization.manufacturer.dl_thickness_in_mm - if dl_vendor isa PropDicts.MissingProperty - throw(ArgumentError("No dead layer measurement provided by vendor! Please provide value or skip (default 0).")) - end - dl_vendor - else - env.dl - end - - @info "DL = $(dl_thickness_in_mm)mm" - - # ----------------------------------------------------------- - # operating voltage - - # sometimes recV is null -> in this case, complain and ask to provide (default env.opV is 0 if none given) - operation_voltage = - if env.operating_voltage == 0 - @warn "You did not provide operating voltage in environment settings -> taking recommended voltage from metadata" - recV = meta.characterization.l200_site.recommended_voltage_in_V - if isnothing(recV) - error("Metadata does not provide recommended voltage (null). Please provide operating voltage in settings") - end - recV - else - env.operating_voltage - end - - @info "Simulating at $(operation_voltage)V" - - # ----------------------------------------------------------- - # parameters common for all geometries - - siggen_geometry = Dict( - # -- crystal - "xtal_length" => meta.geometry.height_in_mm, - "xtal_radius" => meta.geometry.radius_in_mm, - # -- p+ contact - "pc_length" => meta.geometry.pp_contact.depth_in_mm, - "pc_radius" => meta.geometry.pp_contact.radius_in_mm, - # -- tapering - - # comment from David Radford: - - # The bottom taper is always at 45 degrees, but the top outer taper is not. - # It's width/angle is defined using either the outer_taper_width parameter or the taper_angle parameter. - # Likewise, the inner taper width/angle is defined using either the inner_taper_width parameter or the taper_angle parameter. - # The inner_taper_length can be anything from zero to the hole_length - - # size of 45-degree taper at bottom of ORTEC-type crystal (for r=z) - # -> for now, ignore bottom taper if not 45 - # -> works for 1) the ones with actual 45 taper, 2) the ones bulletized (saved as 45 deg taper with bullet radius) - # few that have proper non-45 taper - not possible. ToDo: put some warning - "bottom_taper_length" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm : 0, - # current metadata format: always angle - "bottom_taper_width" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm * tan(deg2rad(meta.geometry.taper.bottom.angle_in_deg)) : 0, - # z-length of outside taper for inverted-coax style -> Mariia: you mean top taper here? - "outer_taper_length" => meta.geometry.taper.top.height_in_mm, - "outer_taper_width" => meta.geometry.taper.top.height_in_mm * tan(deg2rad(meta.geometry.taper.top.angle_in_deg)), - # taper angle in degrees, for inner or outer taper -> how? why are they the same? how about borehole? not using this parameter... - # "taper_angle" => meta.geometry.taper.top.outer.angle_in_deg, - # "taper_angle" => meta.geometry.taper.top.inner.angle_in_deg, - # with this setting, it somehow takes taper angle: 45.000038 anyway -> force to zero? - # "taper_angle" => 0, - # depth of full-charge-collection boundary for Li contact - "Li_thickness" => dl_thickness_in_mm, - - # configuration for mjd_fieldgen (calculates electric fields & weighing potentials) - # detector bias for fieldgen, in Volts - "xtal_HV" => operation_voltage, - - # configuration for signal calculation - # crystal temperature in Kelvin - "xtal_temp" => env.crystal_temperature - ) - - # ----------------------------------------------------------- - # non common parameters - - # -- groove (for non PPCs) - if haskey(meta.geometry, :groove) - siggen_geometry["wrap_around_radius"] = meta.geometry.groove.radius_in_mm.outer - siggen_geometry["ditch_depth"] = meta.geometry.groove.depth_in_mm - siggen_geometry["ditch_thickness"] = meta.geometry.groove.radius_in_mm.outer - meta.geometry.groove.radius_in_mm.inner - end - - # -- borehole (for ICPC and Coax) - if haskey(meta.geometry, :borehole) - # length of hole, for inverted-coax style - siggen_geometry["hole_length"] = meta.geometry.borehole.depth_in_mm - # radius of hole, for inverted-coax style - siggen_geometry["hole_radius"] = meta.geometry.borehole.radius_in_mm - # z-length of inside (hole) taper for inverted-coax style - siggen_geometry["inner_taper_length"] = meta.geometry.taper.borehole.height_in_mm - siggen_geometry["inner_taper_width"] = meta.geometry.taper.borehole.height_in_mm * tan(deg2rad(meta.geometry.taper.borehole.angle_in_deg)) - end - - siggen_geometry -end diff --git a/src/mjdsiggen_utils.jl b/src/mjdsiggen_utils.jl deleted file mode 100644 index 2b0baf1..0000000 --- a/src/mjdsiggen_utils.jl +++ /dev/null @@ -1,29 +0,0 @@ -function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) - E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); - T = eltype(E_pot) - r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid - z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid - ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) - ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) - ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) - axs = (ax1, ax2, ax3) - grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) - data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); - data[:] = E_pot'[:]; - ElectricPotential(data, grid) -end - -function SolidStateDetectors.WeightingPotential(sim::SigGenSetup) - E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); - T = eltype(W_pot) - r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid - z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid - ax1 = SolidStateDetectors.DiscreteAxis{T, :r0, :fixed, ClosedInterval{T}}(r_axis[1]..r_axis[end], r_axis) - ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) - ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) - axs = (ax1, ax2, ax3) - grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) - data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); - data[:] = W_pot'[:]; - WeightingPotential(data, grid) -end \ No newline at end of file diff --git a/src/pss.jl b/src/pss.jl index 52e0e2c..fb4fd0d 100644 --- a/src/pss.jl +++ b/src/pss.jl @@ -42,6 +42,29 @@ Construct SSDSimulator instance based on simulation Currently SSDSimulator does not have any parameters """ +@with_kw struct SiggenSimulator <: PSSimulator + "Path to fieldgen settings" + fieldgen_config::AbstractString + + "Drift velocity correction (?)" + drift_vel::AbstractString + + ".dat/spe file with impurity profile" + impurity_profile::AbstractString="" + + "offset of detector Z=0 in crystal from seed start" + offset_in_mm::Real=-1 + + "Path to crystal jsons: alternative to .dat file" + crystal_metadata_path::AbstractString = "" + + "Name for cached simulation. Not caching if empty string" + cached_name::AbstractString = "" +end + + +SiggenSimulator(::Any) = throw(ErrorException("SiggenSimulator requires MJDSigGen to be loaded, e.g. via `import MJDSigGen`")) + # function SSDSimulator(sim_conf::LegendGeSimConfig) function SSDSimulator(simulation_settings::PropDict) if(!haskey(simulation_settings, :crystal_metadata_path)) @@ -75,81 +98,6 @@ function SSDSimulator(simulation_settings::PropDict) end -""" -Simulation method: siggen -""" -@with_kw struct SiggenSimulator <: PSSimulator - "Path to fieldgen settings" - fieldgen_config::AbstractString - - "Drift velocity correction (?)" - drift_vel::AbstractString - - ".dat/spe file with impurity profile" - impurity_profile::AbstractString="" - - "offset of detector Z=0 in crystal from seed start" - offset_in_mm::Real=-1 - - "Path to crystal jsons: alternative to .dat file" - crystal_metadata_path::AbstractString = "" - - "Name for cached simulation. Not caching if empty string" - cached_name::AbstractString = "" -end - - -""" - SiggeSimulator(sim_conf) - - LegendGeSimConfig -> SiggenSimulator - -Construct SiggeSimulator instance based on simulation - configuration given in . -""" -function SiggenSimulator(simulation_settings::PropDict) - # set defaults - inputs = Dict{String, Any}() - for k in (:fieldgen_config, :drift_vel, :impurity_profile, :crystal_metadata_path) - inputs[string(k)] = haskey(simulation_settings, k) ? simulation_settings[k] : "" - end - - # check if provided paths exist - for (param, path) in inputs - if (path != "") && !(isfile(path) || isdir(path)) - throw(ArgumentError("The input for $(param) that you provided does not exist: $(path)")) - end - end - - inputs["offset_in_mm"] = haskey(simulation_settings, :offset_in_mm) ? simulation_settings.offset_in_mm : -1 - - # throw error if both files are provided - if inputs["crystal_metadata_path"] != "" && inputs["impurity_profile"] != "" - throw(ArgumentError("You provided both an .spe/.dat file and path to crystal metadata! Provide one or the other as impurity input for siggen.")) - # warn if none are provided - will use constant impurity density - elseif inputs["crystal_metadata_path"] == "" && inputs["impurity_profile"] == "" - @warn "No impurity inputs given. Simulation with dummy constant impurity density." - # if we get here, one is provided - else - impinput = inputs["impurity_profile"] * inputs["crystal_metadata_path"] - @info "Impurity profile information based on $(impinput)" - end - - # check that offset is provided if .spe/.dat file is - if inputs["impurity_profile"] != "" && inputs["offset_in_mm"] == -1 - throw(ArgumentError("Provide offset_in_mm of this detector corresponding to impurity file $(inputs["impurity_profile"])!")) - end - - - SiggenSimulator( - inputs["fieldgen_config"], - inputs["drift_vel"], - inputs["impurity_profile"], - inputs["offset_in_mm"], - inputs["crystal_metadata_path"], - simulation_settings.cached_name - ) -end """ @@ -162,6 +110,7 @@ Construct a PSSSimulator supertype instance based on given simulation Returned type depends on the simulation method given in the config. """ + function PSSimulator(simulation_settings::PropDict) @info "Simulation method: $(simulation_settings.method)" @@ -179,7 +128,6 @@ function PSSimulator(simulation_settings::PropDict) error("This simulation method is not implemented!") end end - # ------------------------------------------------------------------- @@ -195,6 +143,7 @@ Constructs and returns a table with resulting pulses and a pss truth table (may be abolished in the future as unnecessary) """ function simulate_waveforms(stp_events::Table, detector::SolidStateDetectors.Simulation, simulator::SSDSimulator) + @info "here" @info("~.~.~ SolidStateDetectors") contact_charge_signals = SolidStateDetectors.simulate_waveforms( stp_events, @@ -236,96 +185,4 @@ function simulate_waveforms(stp_events::Table, detector::SolidStateDetectors.Sim ) pss_table, pss_truth -end - - -""" - simulate_waveforms(stp_events, detector) - -Table, AbstractString -> Table, Table - -Simulate pulses based on events given in - using Siggen with settings given in - -Constructs and returns a table with resulting pulses and a pss truth table - (may be abolished in the future as unnecessary) -""" -function simulate_waveforms(stp_events::Table, detector::SigGenSetup, ::SiggenSimulator) - T = Float32 # This should be somehow defined and be passed properly - @info("~.~.~ Siggen") - - nevents = length(stp_events) - wf_array = Array{RDWaveform}(undef, nevents) - for i in 1:nevents - # println("$i/$nevents") - if(i % 100 == 0) println("$i/$nevents") end - signal::Vector{Float32} = simulate_signal(stp_events[i].pos, stp_events[i].edep, detector) - time = range(T(0)u"ns", step = detector.step_time_out*u"ns", length = length(signal)) - - # push!(waveforms, signal) - wf_array[i] = RDWaveform(time, signal) - end - - ## Oliver said this should be the more efficient way of doing it - ## as opposed to an array of separate RDWaveform instances - # Δt = setup.step_time_out*u"ns" - # t_end = length(waveforms[1]) * Δt - # times = fill(0u"ns" : Δt : t_end, length(waveforms)) - # ArrayOfRDWaveforms((times, VectorOfSimilarVectors(waveforms))) - - waveforms = ArrayOfRDWaveforms(wf_array) - # why am I doing this? - waveforms = ArrayOfRDWaveforms((waveforms.time, VectorOfSimilarVectors(waveforms.signal))) - - pss_table = Table( - channel = [1 for idx in 1:length(waveforms)], # lists of ADCs that triggered, 1 for HADES all the time - ievt = stp_events.evtno, - waveform = waveforms - ) - - # using SSD we get mc truth from SSD output - # with siggen, let's just return the input - # maybe mc truth will not be propagated to (sim)raw anyway - pss_truth = Table( - detno = stp_events.detno, - edep = stp_events.edep, - ievt = stp_events.evtno, - pos = stp_events.pos, - thit = stp_events.thit - ) - - pss_table, pss_truth -end - - -""" - simulate_wf(pos, edep, siggen_setup) - -AbstractVector, AbstractVector, SigGenSetup -> Vector{Float32} - -Simulate a signal from energy depositions with corresponding positions - based on a given -""" -function simulate_signal(pos::AbstractVector, edep::AbstractVector, siggen_setup::SigGenSetup) - # this is not so nice, think about it - signal = zeros(Float32, siggen_setup.ntsteps_out) - - # round to siggen crystal grid precision - # a = string(siggen_setup.xtal_grid) # crystal grid precision e.g. 0.1 (mm) - # ndig = length(a[findfirst('.', a)+1:end]) # number of digits after . - - for i in 1:length(pos) - # pos_rounded = round.(ustrip.(pos[i]), digits=ndig) # strip of units and round to siggen precision - # in SSD the output is always in eV -> convert to eV - signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(uconvert(u"eV", edep[i])) - # signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(pos_rounded)) * ustrip(uconvert(u"eV", edep[i])) - - # somehow this doesn't work - # MJDSigGen.get_signal!(signal, siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(uconvert(u"eV", edep[i])) - end - - signal -end - - - +end \ No newline at end of file diff --git a/src/temp_utils.jl b/src/temp_utils.jl index 37d8be5..ac00cd3 100644 --- a/src/temp_utils.jl +++ b/src/temp_utils.jl @@ -123,9 +123,4 @@ function trigger_threshold(trigger::Trigger, preamp::GenericPreAmp, noise_model: trigger.threshold_keV == 0u"keV" ? std(noise_model.baseline_catalog.waveform[1].signal) * 3 : uconvert(u"eV", trigger.threshold_keV) / germanium_ionization_energy * preamp.gain end - -capacitance_matrix(sim::Simulation) = calculate_capacitance_matrix(sim) -function capacitance_matrix(sim::SigGenSetup) - c = sim.capacitance * u"pF" - [c -c; -c c] -end \ No newline at end of file +capacitance_matrix(sim::Simulation) = calculate_capacitance_matrix(sim) \ No newline at end of file From e1a8e946d511701f19aa5273f2316203cc4d282d Mon Sep 17 00:00:00 2001 From: Shailesh Giri Date: Mon, 18 Nov 2024 12:45:05 -0500 Subject: [PATCH 5/9] Update Project.toml --- Project.toml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 8850a82..a842556 100644 --- a/Project.toml +++ b/Project.toml @@ -50,7 +50,7 @@ ArgCheck = "1, 2" ArraysOfArrays = "0.5, 0.6" CurveFit = "0.3, 0.5, 0.6" DSP = "0.7" -DelimitedFiles = " 1" +DelimitedFiles = "1" Distributions = "0.23, 0.24, 0.25" ElasticArrays = "1" IntervalSets = "0.5, 0.6, 0.7" @@ -67,16 +67,16 @@ PropDicts = "0.2" RadiationDetectorDSP = "0.2.1" RadiationDetectorSignals = "0.3" RadiationSpectra = "0.5" -Random123 = "1.2, 1.7" # Compatible with 1.7? +Random123 = "1.2" RecipesBase = "1" Requires = "0.5, 1" SignalAnalysis = "0.4, 0.5, 0.6, 0.7, 0.8, 0.9" -SolidStateDetectors = "0.10.5" +SolidStateDetectors = "0.10.3" StaticArrays = "0.12, 1" Statistics = "1" StatsBase = "0.32, 0.33, 0.34" StructArrays = "0.5, 0.6" Tables = "1" -TypedTables = "1.2, 1.4" #Compatible with 1.4.6? +TypedTables = "1.2" Unitful = "1" -julia = "1.11" +julia = "1.10" From be0015ee7045b6263dd6672dc73d9a93a972ee0b Mon Sep 17 00:00:00 2001 From: Shailesh Giri Date: Tue, 19 Nov 2024 11:40:16 -0500 Subject: [PATCH 6/9] Resolving most of the comments --- ext/LegendGeSimExt.jl | 72 +++++++++++++++++++++---------------------- src/pss.jl | 1 - 2 files changed, 36 insertions(+), 37 deletions(-) diff --git a/ext/LegendGeSimExt.jl b/ext/LegendGeSimExt.jl index f21d147..734ec44 100644 --- a/ext/LegendGeSimExt.jl +++ b/ext/LegendGeSimExt.jl @@ -13,7 +13,7 @@ using DelimitedFiles using Distributions using DSP using ElasticArrays -using HDF5 +#using HDF5 using IntervalSets using JSON using LegendDataManagement @@ -50,7 +50,7 @@ using SolidStateDetectors: ConstructiveSolidGeometry as CSG ### FieldGen #################################### """ - SiggeSimulator(sim_conf::PropDict) + SiggenSimulator(sim_conf::PropDict) LegendGeSimConfig -> SiggenSimulator @@ -79,7 +79,7 @@ function SiggenSimulator(simulation_settings::PropDict) throw(ArgumentError("You provided both an .spe/.dat file and path to crystal metadata! Provide one or the other as impurity input for siggen.")) # warn if none are provided - will use constant impurity density elseif inputs["crystal_metadata_path"] == "" && inputs["impurity_profile"] == "" - @warn "No impurity inputs given. Simulation with dummy constant impurity density." + @warn "No impurity inputs given. Simulation with dummy constant impurity density of -0.9*10e10 e/cm3. See examples/fieldgen_settings.txt for more details." # if we get here, one is provided else impinput = inputs["impurity_profile"] * inputs["crystal_metadata_path"] @@ -117,38 +117,38 @@ Look up fieldgen generated electric field and weighting potential files function LegendGeSim.simulate_detector(det_meta::PropDict, env::Environment, simulator::SiggenSimulator; overwrite::Bool = false) -# if no cached name given, force simulation from scratch (or else might read past "tmp" file) -if simulator.cached_name == "" - overwrite = true -end + # if no cached name given, force simulation from scratch (or else might read past "tmp" file) + if simulator.cached_name == "" + overwrite = true + end -# returns the name of the resulting siggen config file -# and the name of (already or to be) generated weighting potential file -# ToDo: don't create if already present already at this stage -> check in siggen_config() if name exists and just return name -siggen_config_name, fieldgen_wp_name = siggen_config(det_meta, env, simulator; overwrite) -fieldgen_wp_name = joinpath("cache", fieldgen_wp_name) - -# if the WP file with such a name does not exist yet or want to overwrite it... -if !isfile(fieldgen_wp_name) || overwrite - imp_filename, offset = - if simulator.crystal_metadata_path != "" - # create impurity input on the fly - LegendGeSim.impurity_density_model(det_meta, simulator.crystal_metadata_path, simulator) - else - # user provided .dat file and corresponding offset - # (or if nothing is provided this will be "" and -1) - simulator.impurity_profile, simulator.offset_in_mm - end - #...call fieldgen -> will write the wp file - @info "_|~|_|~|_|~|_ Fieldgen simulation" - fieldgen(siggen_config_name; impurity_profile=imp_filename, offset_mm=offset) - @info "_|~|_|~|_|~|_ Fieldgen simulation complete" -else - #...do nothing, siggen will later read the files based on the generated conifg - @info "Cached simulation found. Reading cached fields from $fieldgen_wp_name" -end -# a SigGenSetup object -SigGenSetup(siggen_config_name) + # returns the name of the resulting siggen config file + # and the name of (already or to be) generated weighting potential file + # ToDo: don't create if already present already at this stage -> check in siggen_config() if name exists and just return name + siggen_config_name, fieldgen_wp_name = siggen_config(det_meta, env, simulator; overwrite) + fieldgen_wp_name = joinpath("cache", fieldgen_wp_name) + + # if the WP file with such a name does not exist yet or want to overwrite it... + if !isfile(fieldgen_wp_name) || overwrite + imp_filename, offset = + if simulator.crystal_metadata_path != "" + # create impurity input on the fly + LegendGeSim.impurity_density_model(det_meta, simulator.crystal_metadata_path, simulator) + else + # user provided .dat file and corresponding offset + # (or if nothing is provided this will be "" and -1) + simulator.impurity_profile, simulator.offset_in_mm + end + #...call fieldgen -> will write the wp file + @info "_|~|_|~|_|~|_ Fieldgen simulation" + fieldgen(siggen_config_name; impurity_profile=imp_filename, offset_mm=offset) + @info "_|~|_|~|_|~|_ Fieldgen simulation complete" + else + #...do nothing, siggen will later read the files based on the generated conifg + @info "Cached simulation found. Reading cached fields from $fieldgen_wp_name" + end + # a SigGenSetup object + SigGenSetup(siggen_config_name) end ## ---------- fieldgen @@ -535,7 +535,7 @@ end #mjdsiggen_utils.jl function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) - E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); + E_pot, W_pot, E_abs, E_r, E_z = MJDSigGen.read_fields(sim); T = eltype(E_pot) r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid @@ -550,7 +550,7 @@ function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) end function SolidStateDetectors.WeightingPotential(sim::SigGenSetup) - E_pot, W_pot, E_abs, E_r, E_z = LegendGeSim.MJDSigGen.read_fields(sim); + E_pot, W_pot, E_abs, E_r, E_z = MJDSigGen.read_fields(sim); T = eltype(W_pot) r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid z_axis = (0:(sim.zlen - 1)) * sim.xtal_grid diff --git a/src/pss.jl b/src/pss.jl index fb4fd0d..b7787dc 100644 --- a/src/pss.jl +++ b/src/pss.jl @@ -143,7 +143,6 @@ Constructs and returns a table with resulting pulses and a pss truth table (may be abolished in the future as unnecessary) """ function simulate_waveforms(stp_events::Table, detector::SolidStateDetectors.Simulation, simulator::SSDSimulator) - @info "here" @info("~.~.~ SolidStateDetectors") contact_charge_signals = SolidStateDetectors.simulate_waveforms( stp_events, From 616f7977d7932731e62f8fe76fc3fdb0f52d909a Mon Sep 17 00:00:00 2001 From: Shailesh Giri Date: Tue, 19 Nov 2024 12:13:18 -0500 Subject: [PATCH 7/9] Resolving more comments --- ext/LegendGeSimExt.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ext/LegendGeSimExt.jl b/ext/LegendGeSimExt.jl index 734ec44..783f1f7 100644 --- a/ext/LegendGeSimExt.jl +++ b/ext/LegendGeSimExt.jl @@ -58,7 +58,7 @@ Construct SiggeSimulator instance based on simulation configuration given in . """ -function SiggenSimulator(simulation_settings::PropDict) +function LegendGeSim.SiggenSimulator(simulation_settings::PropDict) # set defaults inputs = Dict{String, Any}() for k in (:fieldgen_config, :drift_vel, :impurity_profile, :crystal_metadata_path) @@ -79,7 +79,7 @@ function SiggenSimulator(simulation_settings::PropDict) throw(ArgumentError("You provided both an .spe/.dat file and path to crystal metadata! Provide one or the other as impurity input for siggen.")) # warn if none are provided - will use constant impurity density elseif inputs["crystal_metadata_path"] == "" && inputs["impurity_profile"] == "" - @warn "No impurity inputs given. Simulation with dummy constant impurity density of -0.9*10e10 e/cm3. See examples/fieldgen_settings.txt for more details." + @warn "No impurity inputs given. Simulation with dummy constant impurity density of -0.9e10 e/cm3. See examples/fieldgen_settings.txt for more details." # if we get here, one is provided else impinput = inputs["impurity_profile"] * inputs["crystal_metadata_path"] @@ -251,7 +251,7 @@ end """ - simulate_wf(pos, edep, siggen_setup) + simulate_signal(pos, edep, siggen_setup) AbstractVector, AbstractVector, SigGenSetup -> Vector{Float32} From b4fe12efa9fe71f8b4e76bf306ec610b6f6041d1 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 19 Nov 2024 18:16:10 +0100 Subject: [PATCH 8/9] Remove non-required packages --- ext/LegendGeSimExt.jl | 99 +++++++++++++++---------------------------- 1 file changed, 35 insertions(+), 64 deletions(-) diff --git a/ext/LegendGeSimExt.jl b/ext/LegendGeSimExt.jl index 783f1f7..90573c7 100644 --- a/ext/LegendGeSimExt.jl +++ b/ext/LegendGeSimExt.jl @@ -1,49 +1,19 @@ module LegendGeSimExt + using LegendGeSim + +using ArraysOfArrays: VectorOfSimilarVectors +using DelimitedFiles: writedlm +using IntervalSets: ClosedInterval, (..) using PropDicts -using LegendGeSim -using LegendGeSim: Environment -using MJDSigGen -using MJDSigGen: SigGenSetup -using LegendGeSim: SiggenSimulator -using ArgCheck -using ArraysOfArrays -using CurveFit -using DelimitedFiles -using Distributions -using DSP -using ElasticArrays -#using HDF5 -using IntervalSets -using JSON -using LegendDataManagement -using LegendDataTypes -using LegendHDF5IO -using LegendTextIO -using LinearAlgebra -using Parameters -using Polynomials -using PropDicts -using RadiationDetectorDSP -using RadiationDetectorSignals -using RadiationSpectra -using Random -using Random123 -using RecipesBase -using Requires -using SignalAnalysis -using SolidStateDetectors -using StaticArrays -using Statistics -using StatsBase -using StructArrays -using Tables -using TypedTables +using RadiationDetectorSignals: RDWaveform, ArrayOfRDWaveforms +using TypedTables: Table using Unitful -import LsqFit +import MJDSigGen +import SolidStateDetectors + -using SolidStateDetectors: ConstructiveSolidGeometry as CSG #################################### @@ -92,7 +62,7 @@ function LegendGeSim.SiggenSimulator(simulation_settings::PropDict) end - SiggenSimulator( + LegendGeSim.SiggenSimulator( inputs["fieldgen_config"], inputs["drift_vel"], inputs["impurity_profile"], @@ -114,7 +84,7 @@ Look up fieldgen generated electric field and weighting potential files based on given , and if not found, call fieldgen. """ -function LegendGeSim.simulate_detector(det_meta::PropDict, env::Environment, simulator::SiggenSimulator; +function LegendGeSim.simulate_detector(det_meta::PropDict, env::LegendGeSim.Environment, simulator::LegendGeSim.SiggenSimulator; overwrite::Bool = false) # if no cached name given, force simulation from scratch (or else might read past "tmp" file) @@ -141,23 +111,23 @@ function LegendGeSim.simulate_detector(det_meta::PropDict, env::Environment, sim end #...call fieldgen -> will write the wp file @info "_|~|_|~|_|~|_ Fieldgen simulation" - fieldgen(siggen_config_name; impurity_profile=imp_filename, offset_mm=offset) + MJDSigGen.fieldgen(siggen_config_name; impurity_profile=imp_filename, offset_mm=offset) @info "_|~|_|~|_|~|_ Fieldgen simulation complete" else #...do nothing, siggen will later read the files based on the generated conifg @info "Cached simulation found. Reading cached fields from $fieldgen_wp_name" end # a SigGenSetup object - SigGenSetup(siggen_config_name) + MJDSigGen.SigGenSetup(siggen_config_name) end ## ---------- fieldgen """ Convert fit parameters from crystal metadata units (1e9 e/cm^3 VS mm) into Siggen units (1e10 e/cm^3 VS mm) -and create a .dat file (unformatted stream of Float32) to be later used when calling fieldgen(). +and create a .dat file (unformatted stream of Float32) to be later used when calling `MJDSigGen.fieldgen()`. """ -function LegendGeSim.impurity_density_model(meta::PropDict, crystal_metadata::PropDict, fit_param::Vector{Float64}, ::SiggenSimulator) +function LegendGeSim.impurity_density_model(meta::PropDict, crystal_metadata::PropDict, fit_param::Vector{Float64}, ::LegendGeSim.SiggenSimulator) a,b,c,tau = fit_param @@ -201,7 +171,7 @@ Simulation method: siggen -function LegendGeSim.simulate_waveforms(stp_events::Table, detector::SigGenSetup, ::SiggenSimulator) +function LegendGeSim.simulate_waveforms(stp_events::Table, detector::MJDSigGen.SigGenSetup, ::LegendGeSim.SiggenSimulator) T = Float32 # This should be somehow defined and be passed properly @info("~.~.~ Siggen") @@ -259,7 +229,7 @@ Simulate a signal from energy depositions with corresponding positions

""" -function simulate_signal(pos::AbstractVector, edep::AbstractVector, siggen_setup::SigGenSetup) +function simulate_signal(pos::AbstractVector, edep::AbstractVector, siggen_setup::MJDSigGen.SigGenSetup) # this is not so nice, think about it signal = zeros(Float32, siggen_setup.ntsteps_out) @@ -270,24 +240,25 @@ function simulate_signal(pos::AbstractVector, edep::AbstractVector, siggen_setup for i in 1:length(pos) # pos_rounded = round.(ustrip.(pos[i]), digits=ndig) # strip of units and round to siggen precision # in SSD the output is always in eV -> convert to eV - signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(uconvert(u"eV", edep[i])) - # signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(pos_rounded)) * ustrip(uconvert(u"eV", edep[i])) + signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(u"eV", edep[i]) + # signal = signal .+ MJDSigGen.get_signal!(siggen_setup, Tuple(pos_rounded)) * ustrip(u"eV", edep[i]) # somehow this doesn't work - # MJDSigGen.get_signal!(signal, siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(uconvert(u"eV", edep[i])) + # MJDSigGen.get_signal!(signal, siggen_setup, Tuple(ustrip.(pos[i]))) * ustrip(u"eV", edep[i]) end signal end #Calculating capacitance matrix using Siggen -function LegendGeSim.capacitance_matrix(sim::SigGenSetup) +function LegendGeSim.capacitance_matrix(sim::MJDSigGen.SigGenSetup) c = sim.capacitance * u"pF" [c -c; -c c] end -#export capacitance_matrix + + #legend_detector_to_siggen -function siggen_config(meta::PropDict, env::Environment, simulator::SiggenSimulator; overwrite::Bool=false) +function siggen_config(meta::PropDict, env::LegendGeSim.Environment, simulator::LegendGeSim.SiggenSimulator; overwrite::Bool=false) # ---------------------------------------------------------------------------- # set up & check if already present in cache # ---------------------------------------------------------------------------- @@ -424,7 +395,7 @@ Construct siggen type config in a form of Dict based on geometry as given in LEGEND metadata and environment settings given in . """ -function meta2siggen(meta::PropDict, env::Environment) +function meta2siggen(meta::PropDict, env::LegendGeSim.Environment) # !!!! this is copy paste from legend_detector_to_SSD !! # ToDo: rewrite to avoid copy-paste @@ -487,10 +458,10 @@ function meta2siggen(meta::PropDict, env::Environment) # few that have proper non-45 taper - not possible. ToDo: put some warning "bottom_taper_length" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm : 0, # current metadata format: always angle - "bottom_taper_width" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm * tan(deg2rad(meta.geometry.taper.bottom.angle_in_deg)) : 0, + "bottom_taper_width" => meta.geometry.taper.bottom.angle_in_deg == 45 ? meta.geometry.taper.bottom.height_in_mm * tand(meta.geometry.taper.bottom.angle_in_deg) : 0, # z-length of outside taper for inverted-coax style -> Mariia: you mean top taper here? "outer_taper_length" => meta.geometry.taper.top.height_in_mm, - "outer_taper_width" => meta.geometry.taper.top.height_in_mm * tan(deg2rad(meta.geometry.taper.top.angle_in_deg)), + "outer_taper_width" => meta.geometry.taper.top.height_in_mm * tand(meta.geometry.taper.top.angle_in_deg), # taper angle in degrees, for inner or outer taper -> how? why are they the same? how about borehole? not using this parameter... # "taper_angle" => meta.geometry.taper.top.outer.angle_in_deg, # "taper_angle" => meta.geometry.taper.top.inner.angle_in_deg, @@ -526,7 +497,7 @@ function meta2siggen(meta::PropDict, env::Environment) siggen_geometry["hole_radius"] = meta.geometry.borehole.radius_in_mm # z-length of inside (hole) taper for inverted-coax style siggen_geometry["inner_taper_length"] = meta.geometry.taper.borehole.height_in_mm - siggen_geometry["inner_taper_width"] = meta.geometry.taper.borehole.height_in_mm * tan(deg2rad(meta.geometry.taper.borehole.angle_in_deg)) + siggen_geometry["inner_taper_width"] = meta.geometry.taper.borehole.height_in_mm * tand(meta.geometry.taper.borehole.angle_in_deg) end siggen_geometry @@ -534,7 +505,7 @@ end #mjdsiggen_utils.jl -function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) +function SolidStateDetectors.ElectricPotential(sim::MJDSigGen.SigGenSetup) E_pot, W_pot, E_abs, E_r, E_z = MJDSigGen.read_fields(sim); T = eltype(E_pot) r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid @@ -543,13 +514,13 @@ function SolidStateDetectors.ElectricPotential(sim::SigGenSetup) ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) axs = (ax1, ax2, ax3) - grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) + grid = SolidStateDetectors.Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); data[:] = E_pot'[:]; - ElectricPotential(data, grid) + SolidStateDetectors.ElectricPotential(data, grid) end -function SolidStateDetectors.WeightingPotential(sim::SigGenSetup) +function SolidStateDetectors.WeightingPotential(sim::MJDSigGen.SigGenSetup) E_pot, W_pot, E_abs, E_r, E_z = MJDSigGen.read_fields(sim); T = eltype(W_pot) r_axis = (0:(sim.rlen - 1)) * sim.xtal_grid @@ -558,10 +529,10 @@ function SolidStateDetectors.WeightingPotential(sim::SigGenSetup) ax2 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(zero(T)..zero(T), zeros(T, 1)) ax3 = SolidStateDetectors.DiscreteAxis{T, :reflecting, :reflecting, ClosedInterval{T}}(z_axis[1]..z_axis[end], z_axis) axs = (ax1, ax2, ax3) - grid = Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) + grid = SolidStateDetectors.Grid{T, 3, SolidStateDetectors.Cylindrical, typeof(axs)}( axs ) data = Array{T, 3}(undef, length(ax1), length(ax2), length(ax3)); data[:] = W_pot'[:]; - WeightingPotential(data, grid) + SolidStateDetectors.WeightingPotential(data, grid) end end #module LegendGeSimExt From c05ed170056a674279a720e075021190389d5e06 Mon Sep 17 00:00:00 2001 From: Felix Hagemann Date: Tue, 19 Nov 2024 18:24:57 +0100 Subject: [PATCH 9/9] Adapt documentation to MJDSigGen extension --- docs/Project.toml | 2 ++ docs/src/tutorials/simulate_fields_lit.jl | 3 ++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 5e8104d..b8d9bbb 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,8 +3,10 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" LegendHDF5IO = "c9265ca6-b027-5446-b1a4-febfa8dd10b0" LegendTestData = "33d6da08-6349-5f7c-b5a4-6ff4d8795aaf" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +MJDSigGen = "0b925020-e659-485c-98dd-189633cb4be4" Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" [compat] Documenter = "1" +MJDSigGen = "0.2" diff --git a/docs/src/tutorials/simulate_fields_lit.jl b/docs/src/tutorials/simulate_fields_lit.jl index 4166e13..98b2c03 100644 --- a/docs/src/tutorials/simulate_fields_lit.jl +++ b/docs/src/tutorials/simulate_fields_lit.jl @@ -89,8 +89,9 @@ plot( # ## Simulation with `Fieldgen` -# The same simulation settings can be used for siggen, changing only the `"method"` settings to `"siggen"` (or `"fieldgen"`) +# The same simulation settings can be used for siggen, changing only the `"method"` settings to `"siggen"` (or `"fieldgen"`) and loading the `MJDSigGen` package. +using MJDSigGen simulation_settings_siggen = deepcopy(simulation_settings) simulation_settings_siggen["method"] = "fieldgen"