python/impl/plot_material_scan.py

0001 # This file is part of the ACTS project.
0002 #
0003 # Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 #
0005 # This Source Code Form is subject to the terms of the Mozilla Public
0006 # License, v. 2.0. If a copy of the MPL was not distributed with this
0007 # file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 # detray includes
0010 import plotting
0011
0012 # python includes
0013 import math
0014 import numpy as np
0015 import os
0016 import pandas as pd
0017
0018 # Common plot labels
0019 label_eta = r"$\eta$"
0020 label_phi = r"$\phi\,\mathrm{[rad]}$"
0021 label_thickness_x0 = r"thickness / $X_0$"
0022 label_path_x0 = r"path length / $X_0$"
0023 label_thickness_l0 = r"thickness / $\Lambda_0$"
0024 label_path_l0 = r"path length / $\Lambda_0$"
0025
0026 # Common options
0027 ldg_loc = "upper center"
0028
0029 """ Read the material scan data from file and prepare data frame """
0030
0031
0032 def read_material_data(inputdir, logging, det_name, read_cuda):
0033
0034     # Input data directory
0035     data_dir = os.fsencode(inputdir)
0036
0037     material_scan_file = cpu_material_trace_file = cuda_material_trace_file = ""
0038
0039     # Find the data files by naming convention
0040     for file in os.listdir(data_dir):
0041         filename = os.fsdecode(file)
0042
0043         if filename.find(det_name + "_material_scan") != -1:
0044             material_scan_file = inputdir + "/" + filename
0045         elif filename.find(det_name + "_navigation_material_trace_cpu") != -1:
0046             cpu_material_trace_file = inputdir + "/" + filename
0047         elif (
0048             read_cuda
0049             and filename.find(det_name + "_navigation_material_trace_cuda") != -1
0050         ):
0051             cuda_material_trace_file = inputdir + "/" + filename
0052
0053     df_scan = pd.read_csv(material_scan_file, float_precision="round_trip")
0054     df_cpu_trace = pd.read_csv(cpu_material_trace_file, float_precision="round_trip")
0055     df_cuda_trace = pd.DataFrame({})
0056     if read_cuda:
0057         df_cuda_trace = pd.read_csv(
0058             cuda_material_trace_file, float_precision="round_trip"
0059         )
0060
0061     return df_scan, df_cpu_trace, df_cuda_trace
0062
0063
0064 """ Calculate edges of bins to plot the material data """
0065
0066
0067 def get_n_bins(df):
0068     # Find the number of ray directions
0069     row_count = df.groupby(df["eta"]).count()
0070     y_bins = row_count["phi"].max()
0071     x_bins = int(len(df["eta"]) / y_bins)
0072     assert (
0073         len(df["eta"]) == x_bins * y_bins
0074     ), "Could not infer the number of rays correctly"
0075
0076     # Get the axis spacing
0077     x_range = np.max(df["eta"]) - np.min(df["eta"])
0078     x_binning = np.linspace(
0079         np.min(df["eta"]) - 0.5 * x_range / x_bins,
0080         np.max(df["eta"]) + 0.5 * x_range / x_bins,
0081         x_bins + 1,
0082     )
0083
0084     y_range = np.max(df["phi"]) - np.min(df["phi"])
0085     y_binning = np.linspace(
0086         np.min(df["phi"]) - 0.5 * y_range / y_bins,
0087         np.max(df["phi"]) + 0.5 * y_range / y_bins,
0088         y_bins + 1,
0089     )
0090
0091     return x_binning, y_binning
0092
0093
0094 """ Calculate the binwise errors: Standard Error on the Mean """
0095
0096
0097 def get_errors(df, n, name):
0098     # Number of entries per bin
0099     errors = []
0100     # Iterate over data per bin
0101     for i in range(0, n):
0102         # Project the next n rows of the data frame (the bin content is the
0103         # mean of the material along phi)
0104         bin_data = df.iloc[i * n : (i + 1) * n]
0105         # Calculate the error on the mean: stddev/sqrt(n)
0106         errors.append(np.std(bin_data[name], axis=0) / math.sqrt(n))
0107
0108     return errors
0109
0110
0111 """ Plot the material thickenss vs phi and eta in units of X_0 """
0112
0113
0114 def X0_vs_eta_phi(df, label, detector, plot_factory, out_format="pdf"):
0115
0116     # Histogram bin edges
0117     x_binning, y_binning = get_n_bins(df)
0118
0119     # Plot the thickness of every material slab in units of X_0
0120     hist_data = plot_factory.hist2D(
0121         x=df["eta"],
0122         y=df["phi"],
0123         z=df["mat_tX0"],
0124         label=detector,
0125         x_axis=plotting.axis_options(label=label_eta),
0126         y_axis=plotting.axis_options(label=label_phi),
0127         z_axis=plotting.axis_options(label=label_thickness_x0),
0128         x_bins=x_binning,
0129         y_bins=y_binning,
0130         figsize=(9, 7),
0131         show_stats=False,
0132     )
0133
0134     plot_factory.write_plot(hist_data, detector + "_" + label + "_t_X0_map", out_format)
0135
0136     # Plot path length through material of the respective ray in units of X_0
0137     hist_data = plot_factory.hist2D(
0138         x=df["eta"],
0139         y=df["phi"],
0140         z=df["mat_sX0"],
0141         label=detector,
0142         x_axis=plotting.axis_options(label=label_eta),
0143         y_axis=plotting.axis_options(label=label_phi),
0144         z_axis=plotting.axis_options(label=label_path_x0),
0145         x_bins=x_binning,
0146         y_bins=y_binning,
0147         figsize=(9, 7),
0148         show_stats=False,
0149     )
0150
0151     plot_factory.write_plot(hist_data, detector + "_" + label + "_s_X0_map", out_format)
0152
0153
0154 """ Plot the material thickenss vs phi and eta in units of L_0 """
0155
0156
0157 def L0_vs_eta_phi(df, label, detector, plot_factory, out_format="pdf"):
0158
0159     # Histogram bin edges
0160     x_binning, y_binning = get_n_bins(df)
0161
0162     # Plot the thickness of every material slab in units of L_0
0163     hist_data = plot_factory.hist2D(
0164         x=df["eta"],
0165         y=df["phi"],
0166         z=df["mat_tL0"],
0167         label=detector,
0168         x_axis=plotting.axis_options(label=label_eta),
0169         y_axis=plotting.axis_options(label=label_phi),
0170         z_axis=plotting.axis_options(label=label_thickness_l0),
0171         x_bins=x_binning,
0172         y_bins=y_binning,
0173         figsize=(9, 7),
0174         show_stats=False,
0175     )
0176
0177     plot_factory.write_plot(hist_data, detector + "_" + label + "_t_L0_map", out_format)
0178
0179     # Plot path length through material of the respective ray in units of L_0
0180     hist_data = plot_factory.hist2D(
0181         x=df["eta"],
0182         y=df["phi"],
0183         z=df["mat_sL0"],
0184         label=detector,
0185         x_axis=plotting.axis_options(label=label_eta),
0186         y_axis=plotting.axis_options(label=label_phi),
0187         z_axis=plotting.axis_options(label=label_path_l0),
0188         x_bins=x_binning,
0189         y_bins=y_binning,
0190         figsize=(9, 7),
0191         show_stats=False,
0192     )
0193
0194     plot_factory.write_plot(hist_data, detector + "_" + label + "_s_L0_map", out_format)
0195
0196
0197 """ Plot the material thickness in units of X_0 vs eta """
0198
0199
0200 def X0_vs_eta(df, label, detector, plot_factory, out_format="pdf"):
0201     # Where to place the legend box
0202     lgd_ops = plotting.legend_options(
0203         loc=ldg_loc, horiz_anchor=0.725, vert_anchor=1.145
0204     )
0205
0206     # Histogram bin edges
0207     x_binning, y_binning = get_n_bins(df)
0208
0209     # Same number of entries in every bin as per uniform ray scan
0210     n_phi = len(y_binning) - 1
0211
0212     hist_data = plot_factory.hist1D(
0213         x=df["eta"],
0214         w=df["mat_tX0"] / n_phi,
0215         errors=get_errors(df, n_phi, "mat_tX0"),
0216         normalize=False,
0217         label=rf"{detector}",
0218         x_axis=plotting.axis_options(label=label_eta),
0219         y_axis=plotting.axis_options(label=label_thickness_x0),
0220         bins=x_binning,
0221         show_stats=False,
0222         figsize=(9, 7),
0223         lgd_ops=lgd_ops,
0224     )
0225
0226     plot_factory.write_plot(hist_data, detector + "_" + label + "_t_X0", out_format)
0227
0228     hist_data = plot_factory.hist1D(
0229         x=df["eta"],
0230         w=df["mat_sX0"] / n_phi,
0231         errors=get_errors(df, n_phi, "mat_sX0"),
0232         normalize=False,
0233         label=rf"{detector}",
0234         x_axis=plotting.axis_options(label=label_eta),
0235         y_axis=plotting.axis_options(label=label_path_x0),
0236         bins=x_binning,
0237         show_stats=False,
0238         figsize=(9, 7),
0239         lgd_ops=lgd_ops,
0240     )
0241
0242     plot_factory.write_plot(hist_data, detector + "_" + label + "_s_X0", out_format)
0243
0244
0245 """ Plot the material thickness in units of L_0 vs eta """
0246
0247
0248 def L0_vs_eta(df, label, detector, plot_factory, out_format="pdf"):
0249     # Where to place the legend box
0250     lgd_ops = plotting.legend_options(
0251         loc=ldg_loc, horiz_anchor=0.725, vert_anchor=1.145
0252     )
0253
0254     # Histogram bin edges
0255     x_binning, y_binning = get_n_bins(df)
0256
0257     # Same number of entries in every bin as per uniform ray scan
0258     n_phi = len(y_binning) - 1
0259
0260     hist_data = plot_factory.hist1D(
0261         x=df["eta"],
0262         w=df["mat_tL0"] / n_phi,
0263         errors=get_errors(df, n_phi, "mat_tL0"),
0264         normalize=False,
0265         label=rf"{detector}",
0266         x_axis=plotting.axis_options(label=label_eta),
0267         y_axis=plotting.axis_options(label=label_thickness_l0),
0268         bins=x_binning,
0269         show_stats=False,
0270         figsize=(9, 7),
0271         lgd_ops=lgd_ops,
0272     )
0273
0274     plot_factory.write_plot(hist_data, detector + "_" + label + "_t_L0", out_format)
0275
0276     hist_data = plot_factory.hist1D(
0277         x=df["eta"],
0278         w=df["mat_sL0"] / n_phi,
0279         errors=get_errors(df, n_phi, "mat_sL0"),
0280         normalize=False,
0281         label=rf"{detector}",
0282         x_axis=plotting.axis_options(label=label_eta),
0283         y_axis=plotting.axis_options(label=label_path_l0),
0284         bins=x_binning,
0285         show_stats=False,
0286         figsize=(9, 7),
0287         lgd_ops=lgd_ops,
0288     )
0289
0290     plot_factory.write_plot(hist_data, detector + "_" + label + "_s_L0", out_format)
0291
0292
0293 """ Compare two material distributions """
0294
0295
0296 def compare_mat(df_truth, df_rec, label, detector, plot_factory, out_format="pdf"):
0297     # Where to place the legend box
0298     lgd_ops = plotting.legend_options(loc=ldg_loc, horiz_anchor=0.5, vert_anchor=1.29)
0299
0300     # Histogram bin edges
0301     x_binning, y_binning = get_n_bins(df_truth)
0302
0303     # Same number of entries in every bin as per uniform ray scan
0304     n_phi = len(y_binning) - 1
0305
0306     truth_data = plot_factory.hist1D(
0307         x=df_truth["eta"],
0308         w=df_truth["mat_sX0"] / n_phi,
0309         errors=get_errors(df_truth, n_phi, "mat_sX0"),
0310         normalize=False,
0311         label=rf"{detector}: scan",
0312         x_axis=plotting.axis_options(label=label_eta),
0313         y_axis=plotting.axis_options(label=label_path_x0),
0314         bins=x_binning,
0315         show_stats=False,
0316         figsize=(10, 10),
0317         layout="tight",
0318         lgd_ops=lgd_ops,
0319     )
0320
0321     # Add recorded data for comparison
0322     rec_data = plot_factory.add_hist(
0323         old_hist=truth_data,
0324         x=df_rec["eta"].to_numpy(dtype=np.double),
0325         w=df_rec["mat_sX0"] / n_phi,
0326         errors=get_errors(df_rec, n_phi, "mat_sX0"),
0327         normalize=False,
0328         label=rf"{detector}: navigator",
0329     )
0330
0331     # Add a ratio plot to hist_data
0332     ratio_data = plot_factory.add_ratio(
0333         nom=truth_data,
0334         denom=rec_data,
0335         label="scan/navigation",
0336         set_log=False,
0337         show_error=True,
0338     )
0339
0340     plot_factory.write_plot(
0341         ratio_data, detector + "_" + label + "_comparison_s_X0", out_format
0342     )