Simulate realistic waveforms
You can also download this tutorial as a Jupyter notebook and a plain Julia source file.
Here we will use an example pet file csv format from LegendTestData corresponding to the Public Inverted Coax.
using LegendGeSim
using PlotsGet inputs from legend-test-data
using LegendTestDataDetector metadata
ldsim_path = joinpath(legend_test_data_path(), "data", "ldsim")
detector_name = "invcoax-metadata"
detector_metadata_filename = joinpath(ldsim_path, detector_name*".json");Alternatively, enter your own path to a real LEGEND detector JSON
detector_metadata_filename = "path/to/V04545A.json"PET input file
path_to_pet_file = joinpath(ldsim_path, "single-invcoax-th228-geant4.csv");Settings
See manual on Field Simulation and Ideal Pulse Simulation for a detailed explanation of environment and simulation settings, as well as the noise model settings
environment_settings = Dict(
"crystal_temperature_in_K" => 77,
"medium" => "vacuum",
);simple settings for point charge simulation with dummy constant impurity
simulation_settings = Dict(
"method" => "SSD",
"cached_name" => "", # a non-empty string will cache the simulation results
);daq_settings = Dict(
"preamp" => Dict(
"type" => "generic",
"t_decay_in_us" => 43, # from V04545A HADES data
"t_rise_in_ns" => 100, # by eye
"gain_ADC_eV" => 0.0138, # by eye from V04545A HADES data FEP @ 36100 ADC
"offset_in_ADC" => 11900, # from V04545A HADES data mean() of baseline
"noise_sigma_in_keV" => 2 # by eye
),
"fadc" => Dict(
"type" => "generic",
"sampling_interval" => 16 # ns, from HADES data
),
"trigger" => Dict(
"type" => "trapezoidal",
"window_lengths" => [250,250,250],
"threshold" => 9 # keV
),
"daq" => Dict(
"type" => "generic",
"nsamples" => 3748, # from HADES data
"baseline_length" => 1770 # by eye from data
)
);noise_model = Dict(
"type" => "sim"
);Simulate from scratch (pet -> raw)
raw_table = LegendGeSim.simulate_raw(detector_metadata_filename, path_to_pet_file, environment_settings, simulation_settings, daq_settings, noise_model; n_waveforms=10)Table with 11 columns and 20 rows:
baseline channel energy ievt numtraces packet_id ⋯
┌──────────────────────────────────────────────────────────────
1 │ 11921.9 1 466.734 keV 2 1.0 0.0 ⋯
2 │ 11914.8 1 226.265 keV 7 1.0 0.0 ⋯
3 │ 11916.1 1 82.8155 keV 12 1.0 0.0 ⋯
4 │ 11911.6 1 484.73 keV 13 1.0 0.0 ⋯
5 │ 11918.3 1 553.003 keV 17 1.0 0.0 ⋯
6 │ 11915.7 1 226.476 keV 40 1.0 0.0 ⋯
7 │ 11917.8 1 2478.94 keV 44 1.0 0.0 ⋯
8 │ 11912.8 1 338.107 keV 55 1.0 0.0 ⋯
9 │ 11915.0 1 90.4273 keV 62 1.0 0.0 ⋯
10 │ 11919.1 1 226.851 keV 64 1.0 0.0 ⋯
11 │ 11917.3 2 467.897 keV 2 1.0 0.0 ⋯
12 │ 11916.0 2 224.927 keV 7 1.0 0.0 ⋯
13 │ 11916.1 2 83.7165 keV 12 1.0 0.0 ⋯
14 │ 11914.5 2 485.214 keV 13 1.0 0.0 ⋯
15 │ 11917.0 2 552.744 keV 17 1.0 0.0 ⋯
16 │ 11914.9 2 225.988 keV 40 1.0 0.0 ⋯
17 │ 11914.6 2 2478.6 keV 44 1.0 0.0 ⋯
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱plot(raw_table.waveform)
The file contains double the amount of waveforms (20) compared to what we asked to simulate (n_waveforms = 10). That's because at the pulse simulation level, n+ contact pulses are kept as well. The first 10 entries are p+ contact waveforms.
length(raw_table)20plot(raw_table.waveform[1:10], legend=false, title="only p+ contact waveforms")
You can save simulated waveforms in a file that can be used later
using LegendHDF5IOraw_name = "cache/test_100wfs_raw.hdf5"
lh5open(raw_name, "w") do f
LegendHDF5IO.writedata(f.data_store, "raw", raw_table[1:10])
endSimulate from pss table in code
Rather that simulating from scratch pet->raw, you may also input an already ready pss table with ideal pulses, and simulate only the pss->raw step, i.e. the DAQ chain simulation
pss_table, pss_truth = LegendGeSim.simulate_pulses(detector_metadata_filename, path_to_pet_file, environment_settings, simulation_settings, noise_model; n_waveforms=10);[ Info: ---------------------- pet -> stp (stepping info)
Processing file: /home/runner/.julia/artifacts/11b2a7ab47d6f84b449b8f16112d0f3489ec697d/legend-exp-legend-testdata-cd70a8d/data/ldsim/single-invcoax-th228-geant4.csv for detector Public Inverted Coax
[ Info: ...clustering
1411 hits before clustering
0.001162 seconds (16.68 k allocations: 1.497 MiB)
244 hits after clustering
[ Info: ...removing hits with no energy deposits
[ Info: ...removing events outside of the detector
[ Info: Simulation method: SSD
┌ Warning: No crystal metadata path given. Simulation with dummy constant impurity density.
└ @ LegendGeSim ~/work/LegendGeSim.jl/LegendGeSim.jl/src/pss.jl:74
[ Info: ---------------------- stp -> pss (ideal pulses)
[ Info: _||_||_||_ Simulate detector
[ Info: Simulating with SSD from scratch for given settings
[ Info: DL = 0mm
┌ Warning: You did not provide operating voltage in environment settings -> taking recommended voltage from metadata
└ @ LegendGeSim ~/work/LegendGeSim.jl/LegendGeSim.jl/src/legend_detector_to_ssd.jl:45
[ Info: Simulating at 3800.0V
...electric potential
...electric field
...weighting potential 1
...weighting potential 2
[ Info: //\//\//\ Adding fano noise
[ Info: DL = 0mm
┌ Warning: You did not provide operating voltage in environment settings -> taking recommended voltage from metadata
└ @ LegendGeSim ~/work/LegendGeSim.jl/LegendGeSim.jl/src/legend_detector_to_ssd.jl:45
[ Info: Simulating at 3800.0V
[ Info: ~.~.~.~.~ Simulate charge pulses
[ Info: ~.~.~ SolidStateDetectors
[ Info: Detector has 2 contacts
[ Info: Table has 10 physics events (23 single charge depositions).
[ Info: Generating waveforms...plot(pss_table.waveform[1:10])
raw_table1 = LegendGeSim.pss_to_raw(pss_table, pss_truth, daq_settings, noise_model)Table with 11 columns and 20 rows:
baseline channel energy ievt numtraces packet_id ⋯
┌──────────────────────────────────────────────────────────────
1 │ 11915.8 1 467.8 keV 2 1.0 0.0 ⋯
2 │ 11913.3 1 227.318 keV 7 1.0 0.0 ⋯
3 │ 11913.6 1 83.8068 keV 12 1.0 0.0 ⋯
4 │ 11915.7 1 484.61 keV 13 1.0 0.0 ⋯
5 │ 11916.6 1 554.221 keV 17 1.0 0.0 ⋯
6 │ 11918.0 1 227.34 keV 40 1.0 0.0 ⋯
7 │ 11912.3 1 2479.98 keV 44 1.0 0.0 ⋯
8 │ 11921.8 1 337.837 keV 55 1.0 0.0 ⋯
9 │ 11917.5 1 89.2925 keV 62 1.0 0.0 ⋯
10 │ 11914.7 1 226.288 keV 64 1.0 0.0 ⋯
11 │ 11916.7 2 467.631 keV 2 1.0 0.0 ⋯
12 │ 11914.6 2 226.594 keV 7 1.0 0.0 ⋯
13 │ 11916.6 2 82.8852 keV 12 1.0 0.0 ⋯
14 │ 11913.0 2 485.072 keV 13 1.0 0.0 ⋯
15 │ 11918.8 2 553.952 keV 17 1.0 0.0 ⋯
16 │ 11916.4 2 226.341 keV 40 1.0 0.0 ⋯
17 │ 11920.0 2 2479.71 keV 44 1.0 0.0 ⋯
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱plot(
plot(pss_table.waveform[1:10]),
plot(raw_table1.waveform[1:10]),
size=(800,400)
)
Save the pss file for the next section
using LegendHDF5IOpss_name = "cache/test_100wfs_pss.hdf5"
lh5open(pss_name, "w") do f
LegendHDF5IO.writedata(f.data_store, "pss/pss", pss_table[1:10])
LegendHDF5IO.writedata(f.data_store, "pss/truth", pss_truth[1:10])
endSimulate from pre-saved pss hdf5 file
Rather than simulating form scratch pet->raw, you may input the name of a pre-saved pss file containing pss and pss truth information
raw_table2 = LegendGeSim.pss_to_raw(pss_name, daq_settings, noise_model)Table with 11 columns and 10 rows:
baseline channel energy ievt numtraces packet_id ⋯
┌──────────────────────────────────────────────────────────────
1 │ 11912.6 1 468.364 keV 2 1.0 0.0 ⋯
2 │ 11916.8 1 225.726 keV 7 1.0 0.0 ⋯
3 │ 11918.8 1 82.3564 keV 12 1.0 0.0 ⋯
4 │ 11918.2 1 484.489 keV 13 1.0 0.0 ⋯
5 │ 11917.6 1 553.34 keV 17 1.0 0.0 ⋯
6 │ 11916.6 1 227.08 keV 40 1.0 0.0 ⋯
7 │ 11915.6 1 2480.69 keV 44 1.0 0.0 ⋯
8 │ 11919.7 1 338.041 keV 55 1.0 0.0 ⋯
9 │ 11913.0 1 91.4586 keV 62 1.0 0.0 ⋯
10 │ 11919.3 1 225.429 keV 64 1.0 0.0 ⋯This page was generated using Literate.jl.