Skip to content

Commit

Permalink
Add functionality to output the positions of the ray surface intersec…
Browse files Browse the repository at this point in the history 
…tions,

as well as python tools to plot them in an xy - and an rz view.
The data can both be plotted for the ray and the helix scan,
which needs to be turned on by the --write_scan_data option.
  • Loading branch information
@niermann999
niermann999 committed Oct 18, 2023
1 parent 8b3e0d5 commit b49dc4a
Show file tree
Hide file tree
Showing 11 changed files with 567 additions and 65 deletions.
2 changes: 1 addition & 1 deletion tests/unit_tests/io/io_json_detector_reader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,7 @@ TEST(io, json_telescope_detector_reader) {
tel_det_config<rectangle2D<>> tel_cfg{rec2};
tel_cfg.positions(positions);

// Wire chamber
// Telescope detector
vecmem::host_memory_resource host_mr;
auto [telescope_det, telescope_names] =
create_telescope_detector(host_mr, tel_cfg);
Expand Down
27 changes: 27 additions & 0 deletions tests/validation/include/detray/validation/detector_helix_scan.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,9 @@
#pragma once

// Project include(s)
#include "detray/geometry/surface.hpp"
#include "detray/intersection/detail/trajectories.hpp"
#include "detray/io/common/detail/file_handle.hpp"
#include "detray/simulation/event_generator/track_generators.hpp"
#include "detray/test/types.hpp"
#include "tests/common/test_base/fixture_base.hpp"
Expand Down Expand Up @@ -81,6 +83,7 @@ class helix_scan : public test::fixture_base<> {
using nav_link_t = typename detector_t::surface_type::navigation_link;

constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
typename detector_t::geometry_context gctx{};

// Index of the volume that the helix origin lies in
dindex start_index{0u};
Expand All @@ -96,6 +99,14 @@ class helix_scan : public test::fixture_base<> {
uniform_track_generator<free_track_parameters_t>(
m_cfg.track_generator());

detray::io::detail::file_handle outfile{
m_cfg.name(), ".csv",
std::ios::out | std::ios::binary | std::ios::trunc};

if (m_cfg.write_intersections()) {
*outfile << "index,type,x,y,z," << std::endl;
}

std::cout << "INFO: Running helix scan on: " << m_names.at(0) << "\n("
<< trk_state_generator.size() << " helices) ...\n"
<< std::endl;
Expand All @@ -113,6 +124,22 @@ class helix_scan : public test::fixture_base<> {
const auto intersection_record = particle_gun::shoot_particle(
m_det, helix, 14.9f * unit<scalar_t>::um);

// Csv output
if (m_cfg.write_intersections()) {
for (const auto &single_ir : intersection_record) {
const auto &intersection = single_ir.second;
const auto sf =
detray::surface{m_det, intersection.sf_desc};
auto glob_pos = sf.local_to_global(gctx, intersection.local,
helix.dir());
*outfile
<< n_tracks << ","
<< static_cast<int>(intersection.sf_desc.barcode().id())
<< "," << glob_pos[0] << "," << glob_pos[1] << ","
<< glob_pos[2] << "," << std::endl;
}
}

// Create a trace of the volume indices that were encountered
// and check that portal intersections are connected
auto [portal_trace, surface_trace] =
Expand Down
30 changes: 29 additions & 1 deletion tests/validation/include/detray/validation/detector_ray_scan.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,10 @@
#pragma once

// Project include(s)
#include "detray/geometry/surface.hpp"
#include "detray/geometry/volume_graph.hpp"
#include "detray/intersection/detail/trajectories.hpp"
#include "detray/io/common/detail/file_handle.hpp"
#include "detray/simulation/event_generator/track_generators.hpp"
#include "detray/test/types.hpp"
#include "tests/common/test_base/fixture_base.hpp"
Expand Down Expand Up @@ -112,6 +114,7 @@ class ray_scan : public test::fixture_base<> {
using nav_link_t = typename detector_t::surface_type::navigation_link;

constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
typename detector_t::geometry_context gctx{};

// Get the volume adjaceny matrix from ray scan
volume_graph graph(m_det);
Expand All @@ -128,16 +131,41 @@ class ray_scan : public test::fixture_base<> {
auto ray_generator =
uniform_track_generator<ray_t>(m_cfg.track_generator());

// Csv output file
detray::io::detail::file_handle outfile{
m_cfg.name(), ".csv",
std::ios::out | std::ios::binary | std::ios::trunc};

if (m_cfg.write_intersections()) {
*outfile << "index,type,x,y,z," << std::endl;
}

std::cout << "INFO: Running ray scan on: " << m_names.at(0) << "\n("
<< ray_generator.size() << " rays) ...\n"
<< std::endl;

for (const auto ray : ray_generator) {
for (const auto &ray : ray_generator) {

// Record all intersections and surfaces along the ray
const auto intersection_record =
particle_gun::shoot_particle(m_det, ray);

// Csv output
if (m_cfg.write_intersections()) {
for (const auto &single_ir : intersection_record) {
const auto &intersection = single_ir.second;
const auto sf =
detray::surface{m_det, intersection.sf_desc};
auto glob_pos =
sf.local_to_global(gctx, intersection.local, ray.dir());
*outfile
<< n_tracks << ","
<< static_cast<int>(intersection.sf_desc.barcode().id())
<< "," << glob_pos[0] << "," << glob_pos[1] << ","
<< glob_pos[2] << "," << std::endl;
}
}

// Create a trace of the volume indices that were encountered
// and check that portal intersections are connected
auto [portal_trace, surface_trace] =
Expand Down
30 changes: 30 additions & 0 deletions tests/validation/python/plot_helpers.py
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
# Detray library, part of the ACTS project (R&D line)
#
# (c) 2023 CERN for the benefit of the ACTS project
#
# Mozilla Public License Version 2.0

from collections import namedtuple

#-------------------------------------------------------------------------------
# Common helpers for plotting measurement data
#-------------------------------------------------------------------------------

""" Filter the data in a data frame by a given prescription """
def filter_data(data, filter = lambda df: [], variables = []):
dataColl = []

# Get global data
if len(filter(data)) == 0:
for var in variables:
dataColl.append(data[var].to_numpy())

# Filtered data
else:
filtered = data.loc[filter]
for var in variables:
dataColl.append(filtered[var].to_numpy())

return tuple(dataColl)

#-------------------------------------------------------------------------------
74 changes: 38 additions & 36 deletions tests/validation/python/plot_material_scan.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@
#
# Mozilla Public License Version 2.0

from pyplot_factory import get_legend_options

# python includes
import numpy as np
import math
Expand All @@ -29,7 +31,7 @@ def get_n_bins(df):


""" Plot the material thickenss vs phi and eta in units of X_0 """
def X0_vs_eta_phi(df, detector, plotFactory, format = "svg"):
def X0_vs_eta_phi(df, detector, plotFactory, out_format = "pdf"):

# Histogram bin edges
xBinning, yBinning = get_n_bins(df)
Expand All @@ -39,34 +41,32 @@ def X0_vs_eta_phi(df, detector, plotFactory, format = "svg"):
x = df["eta"],
y = df["phi"],
z = df['mat_tX0'],
title = r'\textbf{Material Scan - thickness / }$\mathbf{X_0}$',
label = detector,
xLabel = r'$\mathbf{\eta}$',
yLabel = r'$\mathbf{\phi}\,\mathrm{[rad]}$',
xLabel = r'$\eta$',
yLabel = r'$\phi\,\mathrm{[rad]}$',
zLabel = r'thickness / $X_0$',
xBins = xBinning, yBins = yBinning,
showStats = False)

plotFactory.write_plot(hist_data, "t_X0_map", format)
plotFactory.write_plot(hist_data, "t_X0_map", out_format)

# Plot path length through material of the respective ray in units of X_0
hist_data = plotFactory.hist2D(
x = df["eta"],
y = df["phi"],
z = df['mat_sX0'],
title = r'\textbf{Material Scan - path length / }$\mathbf{X_0}$',
label = detector,
xLabel = r'$\mathbf{\eta}$',
yLabel = r'$\mathbf{\phi}\,\mathrm{[rad]}$',
xLabel = r'$\eta$',
yLabel = r'$\phi\,\mathrm{[rad]}$',
zLabel = r'path length / $X_0$',
xBins = xBinning, yBins = yBinning,
showStats = False)

plotFactory.write_plot(hist_data, "s_X0_map", format)
plotFactory.write_plot(hist_data, "s_X0_map", out_format)


""" Plot the material thickenss vs phi and eta in units of L_0 """
def L0_vs_eta_phi(df, detector, plotFactory, format = "svg"):
def L0_vs_eta_phi(df, detector, plotFactory, out_format = "pdf"):

# Histogram bin edges
xBinning, yBinning = get_n_bins(df)
Expand All @@ -76,93 +76,95 @@ def L0_vs_eta_phi(df, detector, plotFactory, format = "svg"):
x = df["eta"],
y = df["phi"],
z = df['mat_tL0'],
title = r'\textbf{Material Scan - thickness /} $\mathbf{\Lambda_0}$',
label = detector,
xLabel = r'$\mathbf{\eta}$',
yLabel = r'$\mathbf{\phi}\,\mathrm{[rad]}$',
xLabel = r'$\eta$',
yLabel = r'$\phi\,\mathrm{[rad]}$',
zLabel = r'thickness / $\Lambda_0$',
xBins = xBinning, yBins = yBinning,
showStats = False)

plotFactory.write_plot(hist_data, "t_L0_map", format)
plotFactory.write_plot(hist_data, "t_L0_map", out_format)

# Plot path length through material of the respective ray in units of L_0
hist_data = plotFactory.hist2D(
x = df["eta"],
y = df["phi"],
z = df['mat_sL0'],
title = r'\textbf{Material Scan - path length /} $\mathbf{\Lambda_0}$',
label = detector,
xLabel = r'$\mathbf{\eta}$',
yLabel = r'$\mathbf{\phi}\,\mathrm{[rad]}$',
xLabel = r'$\eta$',
yLabel = r'$\phi\,\mathrm{[rad]}$',
zLabel = r'path length / $\Lambda_0$',
xBins = xBinning, yBins = yBinning,
showStats = False)

plotFactory.write_plot(hist_data, "s_L0_map", format)
plotFactory.write_plot(hist_data, "s_L0_map", out_format)


""" Plot the material thickness in units of X_0 vs eta """
def X0_vs_eta(df, detector, plotFactory, format = "svg"):
def X0_vs_eta(df, detector, plotFactory, out_format = "pdf"):

# Histogram bin edges
xBinning, _ = get_n_bins(df)
lgd_ops = get_legend_options()
lgd_ops._replace(loc = 'upper center')

hist_data = plotFactory.hist1D(
x = df['eta'],
w = df['mat_tX0'],
normalize = False,
title = r'\textbf{Material Scan - thickness / }$\mathbf{X_0}$',
label = rf'{detector}',
xLabel = r'$\mathbf{\eta}$',
xLabel = r'$\eta$',
yLabel = r'thickness / $X_0$',
bins = xBinning,
showStats = False)
showStats = False,
lgd_ops = lgd_ops)

plotFactory.write_plot(hist_data, "t_X0", format)
plotFactory.write_plot(hist_data, "t_X0", out_format)

hist_data = plotFactory.hist1D(
x = df['eta'],
w = df['mat_sX0'],
normalize = False,
title = r'\textbf{Material Scan - path length / }$\mathbf{X_0}$',
label = rf'{detector}',
xLabel = r'$\mathbf{\eta}$',
xLabel = r'$\eta$',
yLabel = r'path length / $X_0$',
bins = xBinning,
showStats = False)
showStats = False,
lgd_ops = lgd_ops)

plotFactory.write_plot(hist_data, "s_X0", format)
plotFactory.write_plot(hist_data, "s_X0", out_format)


""" Plot the material thickness in units of L_0 vs eta """
def L0_vs_eta(df, detector, plotFactory, format = "svg"):
def L0_vs_eta(df, detector, plotFactory, out_format = "pdf"):

# Histogram bin edges
xBinning, _ = get_n_bins(df)
lgd_ops = get_legend_options()
lgd_ops._replace(loc = 'upper center')

hist_data = plotFactory.hist1D(
x = df['eta'],
w = df['mat_tL0'],
normalize = False,
title = r'\textbf{Material Scan - thickness / }$\mathbf{\Lambda_0}$',
label = rf'{detector}',
xLabel = r'$\mathbf{\eta}$',
xLabel = r'$\eta$',
yLabel = r'thickness / $\Lambda_0$',
bins = xBinning,
showStats = False)
showStats = False,
lgd_ops = lgd_ops)

plotFactory.write_plot(hist_data, "t_L0", format)
plotFactory.write_plot(hist_data, "t_L0", out_format)

hist_data = plotFactory.hist1D(
x = df['eta'],
w = df['mat_sL0'],
normalize = False,
title = r'\textbf{Material Scan - path length / }$\mathbf{\Lambda_0}$',
label = rf'{detector}',
xLabel = r'$\mathbf{\eta}$',
xLabel = r'$\eta$',
yLabel = r'path length / $\Lambda_0$',
bins = xBinning,
showStats = False)
showStats = False,
lgd_ops = lgd_ops)

plotFactory.write_plot(hist_data, "s_L0", format)
plotFactory.write_plot(hist_data, "s_L0", out_format)
Loading

0 comments on commit b49dc4a

Please sign in to comment.