From b49dc4a63bed2faa5591f144fc0625165949ddca Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Mon, 16 Oct 2023 13:35:28 +0200 Subject: [PATCH] Add functionality to output the positions of the ray surface intersections, 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. --- .../unit_tests/io/io_json_detector_reader.cpp | 2 +- .../detray/validation/detector_helix_scan.hpp | 27 ++++ .../detray/validation/detector_ray_scan.hpp | 30 +++- tests/validation/python/plot_helpers.py | 30 ++++ tests/validation/python/plot_material_scan.py | 74 ++++----- tests/validation/python/plot_ray_scan.py | 147 ++++++++++++++++++ tests/validation/python/pyplot_factory.py | 89 +++++++++-- .../python/run_material_validation.py | 13 +- .../python/run_ray_scan_validation.py | 137 ++++++++++++++++ tests/validation/src/detector_validation.cpp | 50 ++++-- tests/validation/src/material_validation.cpp | 33 +++- 11 files changed, 567 insertions(+), 65 deletions(-) create mode 100644 tests/validation/python/plot_helpers.py create mode 100644 tests/validation/python/plot_ray_scan.py create mode 100644 tests/validation/python/run_ray_scan_validation.py diff --git a/tests/unit_tests/io/io_json_detector_reader.cpp b/tests/unit_tests/io/io_json_detector_reader.cpp index eb953c7175..2e674432c8 100644 --- a/tests/unit_tests/io/io_json_detector_reader.cpp +++ b/tests/unit_tests/io/io_json_detector_reader.cpp @@ -87,7 +87,7 @@ TEST(io, json_telescope_detector_reader) { tel_det_config> 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); diff --git a/tests/validation/include/detray/validation/detector_helix_scan.hpp b/tests/validation/include/detray/validation/detector_helix_scan.hpp index 4988ae1d92..661f3f99a1 100644 --- a/tests/validation/include/detray/validation/detector_helix_scan.hpp +++ b/tests/validation/include/detray/validation/detector_helix_scan.hpp @@ -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" @@ -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()}; + typename detector_t::geometry_context gctx{}; // Index of the volume that the helix origin lies in dindex start_index{0u}; @@ -96,6 +99,14 @@ class helix_scan : public test::fixture_base<> { uniform_track_generator( 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; @@ -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::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(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] = diff --git a/tests/validation/include/detray/validation/detector_ray_scan.hpp b/tests/validation/include/detray/validation/detector_ray_scan.hpp index a64130c011..2cdccaace5 100644 --- a/tests/validation/include/detray/validation/detector_ray_scan.hpp +++ b/tests/validation/include/detray/validation/detector_ray_scan.hpp @@ -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" @@ -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()}; + typename detector_t::geometry_context gctx{}; // Get the volume adjaceny matrix from ray scan volume_graph graph(m_det); @@ -128,16 +131,41 @@ class ray_scan : public test::fixture_base<> { auto ray_generator = uniform_track_generator(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(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] = diff --git a/tests/validation/python/plot_helpers.py b/tests/validation/python/plot_helpers.py new file mode 100644 index 0000000000..960eccca7f --- /dev/null +++ b/tests/validation/python/plot_helpers.py @@ -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) + +#------------------------------------------------------------------------------- diff --git a/tests/validation/python/plot_material_scan.py b/tests/validation/python/plot_material_scan.py index 2007456f9f..f057375d9b 100644 --- a/tests/validation/python/plot_material_scan.py +++ b/tests/validation/python/plot_material_scan.py @@ -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 @@ -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) @@ -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) @@ -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) diff --git a/tests/validation/python/plot_ray_scan.py b/tests/validation/python/plot_ray_scan.py new file mode 100644 index 0000000000..7359591af8 --- /dev/null +++ b/tests/validation/python/plot_ray_scan.py @@ -0,0 +1,147 @@ +# 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 + +import plot_helpers +from pyplot_factory import legend_options + +# python includes +import numpy as np + +""" Plot the intersection points of the detector with the rays - xy view """ +def intersection_points_xy(opts, df, detector, scan_type, plotFactory, out_format = "png"): + + n_rays = np.max(df['index']) + 1 + tracks = "rays" if scan_type == "ray" else "helices" + + # Reduce data to the requested z-range (50mm tolerance) + min_z = opts.z_range[0] + max_z = opts.z_range[1] + assert min_z < max_z, "xy plotting range: min z must be smaller that max z" + sensitive_range = lambda data: ((data['z'] > min_z) & (data['z'] < max_z) & (data['type'] == 1)) + portal_range = lambda data: ((data['z'] > min_z) & (data['z'] < max_z) & (data['type'] == 0)) + passive_range = lambda data: ((data['z'] > min_z) & (data['z'] < max_z) & (data['type'] == 2)) + + senstive_x, senstive_y = plot_helpers.filter_data( + data = df, + filter = sensitive_range, + variables = ['x', 'y']) + + # Plot the xy coordinates of the filtered intersections points + lgd_ops = legend_options('upper center', 4, 0.4, 0.005) + hist_data = plotFactory.scatter( + figsize = (8, 8), + x = senstive_x, + y = senstive_y, + xLabel = r'$x\,\mathrm{[mm]}$', + yLabel = r'$y\,\mathrm{[mm]}$', + label = "sensitives", + color = 'C5', + showStats = lambda x, y: f"{n_rays} {tracks}", + lgd_ops = lgd_ops) + + # Portal surfaces + if not opts.hide_portals: + portal_x, portal_y = plot_helpers.filter_data( + data = df, + filter = portal_range, + variables = ['x', 'y']) + + plotFactory.highlight_region(hist_data, portal_x, portal_y, 'C0', \ + "portals") + + # Passive surfaces + if not opts.hide_passives: + passive_x, passive_y = plot_helpers.filter_data( + data = df, + filter = passive_range, + variables = ['x', 'y']) + + plotFactory.highlight_region(hist_data, passive_x, passive_y, 'C2', \ + "passives") + + # Refine legend + hist_data.lgd.legendHandles[0].set_visible(False) + for handle in hist_data.lgd.legendHandles[1:]: + handle.set_sizes([40]) + + # For this plot, move the legend ouside + hist_data.lgd.set_bbox_to_anchor((0.5, 1.11)) + + # Adjust spacing in box + for vpack in hist_data.lgd._legend_handle_box.get_children()[:1]: + for hpack in vpack.get_children(): + hpack.get_children()[0].set_width(0) + + detector_name = detector.replace(' ', '_') + plotFactory.write_plot(hist_data, f"{detector_name}_{scan_type}_scan_xy", out_format) + + +""" Plot the intersection points of the detector with the rays - rz view """ +def intersection_points_rz(opts, df, detector, scan_type, plotFactory, out_format = "png"): + + n_rays = np.max(df['index']) + 1 + tracks = "rays" if scan_type == "ray" else "helices" + + # Reduce data to the requested z-range + sensitive_range = lambda data: (data['type'] == 1) + portal_range = lambda data: (data['type'] == 0) + passive_range = lambda data: (data['type'] == 2) + + sensitive_x, sensitive_y, sensitive_z = plot_helpers.filter_data( + data = df, + filter = sensitive_range, + variables = ['x', 'y', 'z']) + + # Plot the xy coordinates of the filtered intersections points + lgd_ops = legend_options('upper center', 4, 0.8, 0.005) + hist_data = plotFactory.scatter( + figsize = (12, 6), + x = sensitive_z, + y = np.hypot(sensitive_x, sensitive_y), + xLabel = r'$z\,\mathrm{[mm]}$', + yLabel = r'$r\,\mathrm{[mm]}$', + label = "sensitives", + color = 'C5', + showStats = lambda x, y: f"{n_rays} {tracks}", + lgd_ops = lgd_ops) + + # Portal surfaces + if not opts.hide_portals: + portal_x, portal_y, portal_z = plot_helpers.filter_data( + data = df, + filter = portal_range, + variables = ['x', 'y', 'z']) + + plotFactory.highlight_region(hist_data, portal_z, \ + np.hypot(portal_x, portal_y), \ + 'C0', "portals") + + # Passive surfaces + if not opts.hide_passives: + passive_x, passive_y, passive_z = plot_helpers.filter_data( + data = df, + filter = passive_range, + variables = ['x', 'y', 'z']) + + plotFactory.highlight_region(hist_data, passive_z, \ + np.hypot(passive_x, passive_y), \ + 'C2', "passives") + + # Refine legend + hist_data.lgd.legendHandles[0].set_visible(False) + for handle in hist_data.lgd.legendHandles[1:]: + handle.set_sizes([40]) + + # For this plot, move the legend ouside + hist_data.lgd.set_bbox_to_anchor((0.5, 1.15)) + + # Adjust spacing in box + for vpack in hist_data.lgd._legend_handle_box.get_children()[:1]: + for hpack in vpack.get_children(): + hpack.get_children()[0].set_width(0) + + detector_name = detector.replace(' ', '_') + plotFactory.write_plot(hist_data, f"{detector_name}_{scan_type}_scan_rz", out_format) diff --git a/tests/validation/python/pyplot_factory.py b/tests/validation/python/pyplot_factory.py index 1784bb7316..966a6d6fd7 100644 --- a/tests/validation/python/pyplot_factory.py +++ b/tests/validation/python/pyplot_factory.py @@ -4,8 +4,6 @@ # # Mozilla Public License Version 2.0 -# -*- coding: utf-8 -*- - # python includes from collections import namedtuple import math @@ -15,6 +13,16 @@ import matplotlib.pyplot as plt import matplotlib.colors as mcolors +import matplotlib.style as style +style.use('tableau-colorblind10') +#style.use('seaborn-colorblind') + +plt.rcParams.update({ + 'text.usetex': True, + 'font.size': 20, + 'font.family': 'serif', +}) + #------------------------------------------------------------------------------- # Global identifiers #------------------------------------------------------------------------------- @@ -22,6 +30,13 @@ """ Pass plotting data between functions """ plt_data = namedtuple('plt_data', 'fig ax lgd data bins mu rms errors') +""" Wrap the configuration for a legend """ +legend_options = namedtuple('legend_options', 'loc ncol colspacing handletextpad') + +""" Conveniently get the legend options """ +def get_legend_options(): + return legend_options('upper right', 1, 1, 1) + #------------------------------------------------------------------------------- # Data Plotting #------------------------------------------------------------------------------- @@ -40,8 +55,11 @@ def __init__(self, outDir, logger, atlas_badge = ""): """ Add legend to a plot. Labbels must be defined. """ - def add_legend(self, ax): - return ax.legend(loc="upper right") + def add_legend(self, ax, options = get_legend_options()): + return ax.legend(loc = options.loc, + ncol = options.ncol, + columnspacing = options.colspacing, + handletextpad = options.handletextpad) """ @@ -56,7 +74,8 @@ def hist1D(self, x, w = None, setLog = False, normalize = False, showError = False, - showStats = True): + showStats = True, + lgd_ops = get_legend_options()): # Create fresh plot fig = plt.figure(figsize = (8, 6)) @@ -108,7 +127,7 @@ def hist1D(self, x, w = None, ax.grid(True, alpha = 0.25) # Add legend - lgd = self.add_legend(ax) + lgd = self.add_legend(ax, lgd_ops) # Calculate the bin error binCenters = 0.5 * (bins[1:] + bins[:-1]) @@ -167,11 +186,12 @@ def hist2D(self, x, y, z = None, # Fill data data, xbins, ybins, hist = ax.hist2d( x, y, weights = z, - range = [(xMin, xMax), (yMin, yMax)], - bins = (xBins, yBins), + range = [(xMin, xMax), (yMin, yMax)], + bins = (xBins, yBins), label=f"{label} ({len(x)*len(y)} entries)", - facecolor = mcolors.to_rgba(color, alpha), - edgecolor = None) + facecolor = mcolors.to_rgba(color, alpha), + edgecolor = None, + rasterized = True) # Add some additional information if showStats: @@ -198,6 +218,55 @@ def hist2D(self, x, y, z = None, return plt_data(fig, ax, None, data, None, None, None, None) + """ Create a 2D scatter plot """ + def scatter(self, x, y, + xLabel = "", yLabel = "", title = "", label = "", + color = 'tab:blue', alpha = 1, + figsize = (8, 6), + showStats = lambda x, _: f"{len(x)} entries", + lgd_ops = get_legend_options()): + + fig = plt.figure(figsize = figsize) + ax = fig.add_subplot(1, 1, 1) + + # Create empty plot with blank marker containing the extra label + ax.plot([], [], ' ', label=showStats(x, y)) + scatter = ax.scatter(x, y, + label = label, + c = color, + s = 0.1, + alpha = alpha, + rasterized=True) + + # Refine plot + ax.set_title(title) + ax.set_xlabel(xLabel) + ax.set_ylabel(yLabel) + ax.grid(True, alpha = 0.25) + + # Add legend + lgd = self.add_legend(ax, lgd_ops) + + return plt_data(fig, ax, lgd, scatter, None, None, None, None) + + + """ Add new data in a different color to a scatter plot """ + def highlight_region(self, plotData, x, y, color, label = ""): + + if label == "": + plotData.ax.scatter(x, y, c = color, s = 0.1) + else: + plotData.ax.scatter(x, y, c = color, s = 0.1, label=label) + + # Update legend + lgd = plotData.lgd + handles, labels = lgd.axes.get_legend_handles_labels() + lgd._legend_box = None + lgd._init_legend_box(handles, labels) + lgd._set_loc(lgd._loc) + lgd.set_title(lgd.get_title().get_text()) + + """ Safe a plot to disk """ def write_plot(self, plot_data, name = "plot", file_format = "svg", outPrefix = ""): diff --git a/tests/validation/python/run_material_validation.py b/tests/validation/python/run_material_validation.py index e93a5b0f8e..069a881c7b 100644 --- a/tests/validation/python/run_material_validation.py +++ b/tests/validation/python/run_material_validation.py @@ -99,6 +99,11 @@ def __main__(): #----------------------------------------------------------------prepare data + # Get detector name + detector_name = mat_scan_file.removeprefix('material_scan_') + detector_name = detector_name.removesuffix('.csv') + detector_name = detector_name.replace('_', ' ') + df = pd.read_csv(mat_scan_file) plot_factory = pyplot_factory(outdir + "material_", logging) @@ -107,10 +112,10 @@ def __main__(): # The histograms are not re-weighted (if the rays are not evenly distributed # the material in some bins might be artificially high)! - plot_material_scan.X0_vs_eta_phi(df, "toy detector", plot_factory, out_format) - plot_material_scan.L0_vs_eta_phi(df, "toy detector", plot_factory, out_format) - plot_material_scan.X0_vs_eta(df, "toy detector", plot_factory, out_format) - plot_material_scan.L0_vs_eta(df, "toy detector", plot_factory, out_format) + plot_material_scan.X0_vs_eta_phi(df, detector_name, plot_factory, out_format) + plot_material_scan.L0_vs_eta_phi(df, detector_name, plot_factory, out_format) + plot_material_scan.X0_vs_eta(df, detector_name, plot_factory, out_format) + plot_material_scan.L0_vs_eta(df, detector_name, plot_factory, out_format) #------------------------------------------------------------------------------- diff --git a/tests/validation/python/run_ray_scan_validation.py b/tests/validation/python/run_ray_scan_validation.py new file mode 100644 index 0000000000..ca10cc894a --- /dev/null +++ b/tests/validation/python/run_ray_scan_validation.py @@ -0,0 +1,137 @@ +# 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 + +import plot_ray_scan +from pyplot_factory import pyplot_factory + +# python includes +import argparse +import logging +import numpy as np +import pandas as pd +import os +import sys +from datetime import datetime +import matplotlib.pyplot as plt + + +plt.rc('text', usetex=True) +plt.rc('text.latex', preamble=r'\usepackage{amsmath}') + + +def __main__(): + +#----------------------------------------------------------------arg parsing + + parser = argparse.ArgumentParser(description = "Detray Ray Scan") + parser.add_argument("--debug", "-d", + help=("Enables debug logging"), + action="store_true") + parser.add_argument("--logfile", + help=("Write log in file"), + default = "", type=str) + parser.add_argument("--input", "-i", + help=("Input ray scan data file."), + default = "", type=str, required=True) + parser.add_argument("--outdir", "-o", + help=("Output directory for plots."), + default = "./geometry_plots/", type=str) + parser.add_argument("--output_format", "-of", + help=("Format of the plot files (svg|png|pdf)."), + default = "png", type=str) + parser.add_argument("--z_range", "-zrng", nargs=2, + help=("z range for the xy-view."), + default = [-50, 50], type=float) + parser.add_argument("--hide_portals", + action="store_true", default=False) + parser.add_argument("--hide_passives", + action="store_true", default=False) + + args = parser.parse_args() + +#---------------------------------------------------------------------config + + # Check output path + if not os.path.isdir(args.outdir): + os.mkdir(args.outdir, 0o755) + outdir = args.outdir + + # Set log level + logLevel = logging.INFO + if args.debug: + logLevel = logging.DEBUG + + # Check logfile path + if args.logfile != "": + logDirName = os.path.dirname(args.logfile) + + if logDirName != "" and not os.path.isdir(logDirName): + os.mkdir(logDirName, 0o755) + + if not os.path.isfile(args.logfile): + with open(args.logfile, 'x'): pass + + # Write log in logfile + logging.basicConfig(filename=args.logfile, + format=("%(levelname)s (%(module)s):" + " %(message)s"), level=logLevel) + else: + # Write log to terminal + logging.basicConfig(format=("%(levelname)s (%(module)s):" + " %(message)s"), level=logLevel) + + logging.info("\n--------------------------------------------------------\n" + "Running ray scan validation "+\ + str(datetime.now().strftime("%d/%m/%Y %H:%M"))+\ + "\n--------------------------------------------------------\n") + + # Check input data files from material scan + if args.input == "": + logging.error(f"Please specify an input data file!") + sys.exit(1) + + if not os.path.isfile(args.input): + logging.error(f"Data file does not exist! ({args.input})") + sys.exit(1) + + if not args.output_format in ["svg", "png", "pdf"]: + logging.error(f"Unknown output file format: {out_format}") + sys.exit(1) + + ray_scan_file = args.input + out_format = args.output_format + +#----------------------------------------------------------------prepare data + + # Get detector name + if ray_scan_file.find('ray_scan_') != -1: + detector_name = ray_scan_file.removeprefix('ray_scan_') + scan_type = "ray" + elif ray_scan_file.find('helix_scan_') != -1: + detector_name = ray_scan_file.removeprefix('helix_scan_') + scan_type = "helix" + else: + logging.error('Input filename needs to contain \'ray_scan\' prefix') + detector_name = detector_name.removesuffix('.csv') + detector_name = detector_name.replace('_', ' ') + + df = pd.read_csv(ray_scan_file) + + plot_factory = pyplot_factory(outdir + "geometry_", logging) + +#------------------------------------------------------------------------run + + plot_ray_scan.intersection_points_xy(args, df, detector_name, + scan_type, plot_factory, out_format) + plot_ray_scan.intersection_points_rz(args, df, detector_name, scan_type, + plot_factory, out_format) + +#------------------------------------------------------------------------------- + +if __name__ == "__main__": + __main__() + +#------------------------------------------------------------------------------- diff --git a/tests/validation/src/detector_validation.cpp b/tests/validation/src/detector_validation.cpp index e1eff3caa6..ab1bc0ad17 100644 --- a/tests/validation/src/detector_validation.cpp +++ b/tests/validation/src/detector_validation.cpp @@ -43,24 +43,31 @@ int main(int argc, char **argv) { ::testing::InitGoogleTest(&argc, argv); // Options parsing - po::options_description desc("\ndetray validation options"); + po::options_description desc("\ndetray detector validation options"); desc.add_options()("help", "produce help message")( "write_volume_graph", "writes the volume graph to file")( - "write_ray_scan", "writes the ray scan intersections to file")( + "write_scan_data", "writes the ray/helix scan intersections to file")( "geometry_file", po::value(), "geometry input file")( "material_file", po::value(), "material input file")( "phi_steps", po::value()->default_value(50u), "# phi steps for particle gun")( "theta_steps", po::value()->default_value(50u), "# theta steps for particle gun")( + "theta_range", po::value>()->multitoken(), + "min, max range of theta values for particle gun")( + "origin", po::value>()->multitoken(), + "coordintates for particle gun origin position")( "p_mag", po::value()->default_value(10.f), "absolute momentum of the test particle [GeV]")( "overstep_tol", po::value()->default_value(-100.f), "overstepping tolerance [um] NOTE: Must be negative!"); po::variables_map vm; - po::store(po::parse_command_line(argc, argv, desc), vm); + po::store(parse_command_line(argc, argv, desc, + po::command_line_style::unix_style ^ + po::command_line_style::allow_short), + vm); po::notify(vm); // Help message @@ -82,10 +89,9 @@ int main(int argc, char **argv) { con_chk_cfg.write_graph(true); throw std::invalid_argument("Writing of voume graph not implemented"); } - if (vm.count("write_ray_scan")) { + if (vm.count("write_scan_data")) { ray_scan_cfg.write_intersections(true); hel_scan_cfg.write_intersections(true); - throw std::invalid_argument("Writing of ray scan not implemented"); } // Input files @@ -106,11 +112,36 @@ int main(int argc, char **argv) { const std::size_t phi_steps{vm["phi_steps"].as()}; ray_scan_cfg.track_generator().phi_steps(phi_steps); + hel_nav_cfg.track_generator().phi_steps(phi_steps); } if (vm.count("theta_steps")) { const std::size_t theta_steps{vm["theta_steps"].as()}; ray_scan_cfg.track_generator().theta_steps(theta_steps); + hel_nav_cfg.track_generator().theta_steps(theta_steps); + } + if (vm.count("theta_range")) { + const auto theta_range = vm["theta_range"].as>(); + if (theta_range.size() == 2u) { + ray_scan_cfg.track_generator().theta_range(theta_range[0], + theta_range[1]); + hel_nav_cfg.track_generator().theta_range(theta_range[0], + theta_range[1]); + } else { + throw std::invalid_argument("Theta range needs two arguments"); + } + } + if (vm.count("origin")) { + const auto origin = vm["origin"].as>(); + if (origin.size() == 3u) { + ray_scan_cfg.track_generator().origin( + {origin[0], origin[1], origin[2]}); + hel_nav_cfg.track_generator().origin( + {origin[0], origin[1], origin[2]}); + } else { + throw std::invalid_argument( + "Particle gun origin needs three arguments"); + } } if (vm.count("p_mag")) { const scalar_t p_mag{vm["p_mag"].as()}; @@ -140,22 +171,21 @@ int main(int argc, char **argv) { con_chk_cfg); // Navigation link consistency, discovered by ray intersection + ray_scan_cfg.name("ray_scan_" + names.at(0)); detray::detail::register_checks(det, names, ray_scan_cfg); // Navigation link consistency, discovered by helix intersection - // Full number of helices only works in double precision - // hel_nav_cfg.track_generator().phi_steps(100u).theta_steps(100u); - hel_scan_cfg.track_generator().phi_steps(30u).theta_steps(30u); + hel_scan_cfg.name("helix_scan_" + names.at(0)); detray::detail::register_checks(det, names, hel_scan_cfg); // Comparision of straight line navigation with ray scan + str_nav_cfg.name("straight_line_navigation_" + names.at(0)); detray::detail::register_checks( det, names, str_nav_cfg); // Comparision of navigation in a constant B-field with helix - // For now, run with reduced number of helices, until helix test is fixed - hel_nav_cfg.track_generator().phi_steps(30u).theta_steps(30u); + hel_nav_cfg.name("helix_navigation_" + names.at(0)); detray::detail::register_checks(det, names, hel_nav_cfg); diff --git a/tests/validation/src/material_validation.cpp b/tests/validation/src/material_validation.cpp index 9606204e0f..06f5a8faa9 100644 --- a/tests/validation/src/material_validation.cpp +++ b/tests/validation/src/material_validation.cpp @@ -33,12 +33,13 @@ int main(int argc, char **argv) { // Use the most general type to be able to read in all detector files using detector_t = detray::detector<>; + using scalar_t = detector_t::scalar_type; // Filter out the google test flags ::testing::InitGoogleTest(&argc, argv); // Options parsing - po::options_description desc("\ndetray validation options"); + po::options_description desc("\ndetray material validation options"); desc.add_options()("help", "produce help message")( "geometry_file", po::value(), "geometry input file")( @@ -46,10 +47,17 @@ int main(int argc, char **argv) { "phi_steps", po::value()->default_value(50u), "# phi steps for particle gun")( "eta_steps", po::value()->default_value(50u), - "# eta steps for particle gun"); + "# eta steps for particle gun")( + "eta_range", po::value>()->multitoken(), + "min, max range of eta values for particle gun")( + "origin", po::value>()->multitoken(), + "coordintates for particle gun origin position"); po::variables_map vm; - po::store(po::parse_command_line(argc, argv, desc), vm); + po::store(parse_command_line(argc, argv, desc, + po::command_line_style::unix_style ^ + po::command_line_style::allow_short), + vm); po::notify(vm); // Help message @@ -93,6 +101,25 @@ int main(int argc, char **argv) { mat_scan_cfg.track_generator().eta_steps(eta_steps); mat_scan_cfg.track_generator().eta_range(-4.f, 4.f); } + if (vm.count("eta_range")) { + const auto eta_range = vm["eta_range"].as>(); + if (eta_range.size() == 2u) { + mat_scan_cfg.track_generator().eta_range(eta_range[0], + eta_range[1]); + } else { + throw std::invalid_argument("Eta range needs two arguments"); + } + } + if (vm.count("origin")) { + const auto origin = vm["origin"].as>(); + if (origin.size() == 3u) { + mat_scan_cfg.track_generator().origin( + {origin[0], origin[1], -origin[2]}); + } else { + throw std::invalid_argument( + "Particle gun origin needs three arguments"); + } + } vecmem::host_memory_resource host_mr;