From bc4855d4026fca4bdf8bc9b4fb46cd60753705c0 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Wed, 3 May 2023 15:00:58 -0500 Subject: [PATCH 01/11] feat: pull CAD generated service meshes into simulation --- compact/tracking/cable_services.xml | 17 +++++++++++++++++ compact/tracking/get_meshes.xml | 19 +++++++++++++++++++ configurations/services_only.yml | 4 ++++ src/FileLoader.cpp | 2 +- 4 files changed, 41 insertions(+), 1 deletion(-) create mode 100644 compact/tracking/cable_services.xml create mode 100644 compact/tracking/get_meshes.xml create mode 100644 configurations/services_only.yml diff --git a/compact/tracking/cable_services.xml b/compact/tracking/cable_services.xml new file mode 100644 index 000000000..9cd66f32f --- /dev/null +++ b/compact/tracking/cable_services.xml @@ -0,0 +1,17 @@ + + + + + + + + + + + + + + + + + diff --git a/compact/tracking/get_meshes.xml b/compact/tracking/get_meshes.xml new file mode 100644 index 000000000..53aa27121 --- /dev/null +++ b/compact/tracking/get_meshes.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/configurations/services_only.yml b/configurations/services_only.yml new file mode 100644 index 000000000..6f6c5f3df --- /dev/null +++ b/configurations/services_only.yml @@ -0,0 +1,4 @@ +features: + tracking: + get_meshes: + cable_services: diff --git a/src/FileLoader.cpp b/src/FileLoader.cpp index 081b4edea..275196ef0 100644 --- a/src/FileLoader.cpp +++ b/src/FileLoader.cpp @@ -36,7 +36,7 @@ long load_file(Detector& /* desc */, int argc, char** argv) { // argument parsing std::string cache, file, url; - std::string cmd("curl --retry 5 -f {0} -o {1}"); + std::string cmd("curl --location --retry 5 -f {0} -o {1}"); for (int i = 0; i < argc && argv[i]; ++i) { if (0 == std::strncmp("cache:", argv[i], 6)) cache = (argv[i] + 6); From 5929df36060f07a04330c7b7e1c947d535d8fd83 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Wed, 3 May 2023 15:13:52 -0500 Subject: [PATCH 02/11] fix: keep the materials in epic --- compact/tracking/cable_services.xml | 21 ++++++++++++++++++--- compact/tracking/get_meshes.xml | 7 ------- 2 files changed, 18 insertions(+), 10 deletions(-) diff --git a/compact/tracking/cable_services.xml b/compact/tracking/cable_services.xml index 9cd66f32f..401dd4782 100644 --- a/compact/tracking/cable_services.xml +++ b/compact/tracking/cable_services.xml @@ -2,9 +2,24 @@ - - - + + + + + + + + + + + + + + + + + + diff --git a/compact/tracking/get_meshes.xml b/compact/tracking/get_meshes.xml index 53aa27121..010b7d875 100644 --- a/compact/tracking/get_meshes.xml +++ b/compact/tracking/get_meshes.xml @@ -9,11 +9,4 @@ - - - - - - - From 5d0161a5f3cfa933f8e4f50d4f4855cef12120ac Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Wed, 17 May 2023 13:20:17 -0500 Subject: [PATCH 03/11] update material scan script for testing the imported geometry --- scripts/subdetector_tests/material_scan.py | 97 ++++++++++++++++++++-- 1 file changed, 91 insertions(+), 6 deletions(-) diff --git a/scripts/subdetector_tests/material_scan.py b/scripts/subdetector_tests/material_scan.py index dd1afe3ee..148b35ed1 100644 --- a/scripts/subdetector_tests/material_scan.py +++ b/scripts/subdetector_tests/material_scan.py @@ -1,7 +1,7 @@ # SPDX-License-Identifier: LGPL-3.0-or-later # Copyright (C) 2023 Chao Peng ''' - A script to scan the materials thickness (X0 or Lambda) over eta + A script to scan the materials thickness (X0 or Lambda) by detectors over eta It uses dd4hep::rec::MaterialManager::placementsBetween to check the material layers, and then uses the detector alignment to transform world coordinates to local coordiantes, and then assigns the materials to a detector based @@ -34,6 +34,7 @@ COLORS = ['royalblue', 'indianred', 'forestgreen', 'gold', 'darkviolet', 'orange', 'darkturquoise'] # a specified color for Others, this should not be included in COLORS OTHERS_COLOR = 'silver' +MAX_N_MATS = 50 # FIXME: this is a work-around to an issue of dd4hep::rec::MaterialManager::placementsBetween @@ -156,6 +157,10 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= dest='compact', help='Top-level xml file of the detector description.' ) + parser.add_argument( + '-o', default='.', dest='outdir', + help='Output directory.' + ) parser.add_argument( '--path-r', type=float, default=400., # 120., help='R_xy (cm) where the scan stops.' @@ -188,6 +193,10 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= '--value-type', default='X0', help='Choose one in {}.'.format(list(ALLOWED_VALUE_TYPES.keys())) ) + parser.add_argument( + '--groupby-materials', action='store_true', + help='Enable it to group the results by materials instead of by detectors.' + ) parser.add_argument( '--detectors', # default='all', @@ -205,6 +214,8 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= print('Cannot find value {}, choose one in {}'.format(args.value_type, list(ALLOWED_VALUE_TYPES.keys()))) exit(-1) + os.makedirs(args.outdir, exist_ok=True) + # scan parameters path_r = args.path_r phi = args.phi @@ -246,12 +257,86 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= th_corr.scan_missing_thickness(desc, start_point, np.array(start_point) + np.array([100., 50., 0.]), dz=0.0001, epsilon=args.epsilon) + # value-type + vt = ALLOWED_VALUE_TYPES[args.value_type] + + # material-based scan + if args.groupby_materials: + vals = np.zeros(shape=(len(etas), MAX_N_MATS)) + mats = [] + for i, eta in enumerate(etas): + if i % PROGRESS_STEP == 0: + print('Scanned {:d}/{:d} for {:.2f} <= eta <= {:.2f}'.format(i, len(etas), etas[0], etas[-1]), + end='\r', flush=True) + # scan material layers + end_x = path_r*np.cos(phi/180.*np.pi) + end_y = path_r*np.sin(phi/180.*np.pi) + end_z = path_r*np.sinh(eta) + # print('({:.2f}, {:.2f}, {:.2f})'.format(end_x, end_y, end_z)) + dfr = material_scan(desc, start_point, (end_x, end_y, end_z), epsilon=args.epsilon, int_dets=dets, thickness_corrector=th_corr) + + # aggregated values for detectors + x0_vals = dfr.groupby('material')[vt[0]].sum().to_dict() + # update material list + for mat in x0_vals.keys(): + if mat not in mats: + mats.append(mat) + # save material scan results + for j, mat in enumerate(mats): + vals[i, j] = x0_vals.get(mat, 0.) + print('Scanned {:d}/{:d} lines for {:.2f} < eta < {:.2f}'.format(len(etas), len(etas), etas[0], etas[-1])) + + # aggregated data for plots + dfa = pd.DataFrame(data=vals[:, :len(mats)], columns=mats, index=etas) + # sort by total thickness + smats = [d for d in dfa.sum().sort_values(ascending=False).index] + dfa.to_csv(os.path.join(args.outdir, 'material_scan_agg.csv')) + + # plot material scan + fig, ax = plt.subplots(figsize=(16, 8), dpi=160, + gridspec_kw={'top': 0.9, 'bottom': 0.2, 'left': 0.08, 'right': 0.98}) + + # group some materials into "Others" if there're no enough colors + plot_mats = smats + if len(smats) > len(COLORS): + group_mats = smats[len(COLORS):] + dfa.loc[:, OTHERS_NAME] += dfa.loc[:, group_mats].sum(axis=1) + print('Materials {} are grouped into {} due to limited number of colors.'.format(group_mats, OTHERS_NAME)) + plot_mats = sdets[:len(COLORS)] + + bottom = np.zeros(len(dfa)) + width = np.mean(np.diff(dfa.index)) + for col, c in zip(plot_mats, COLORS): + ax.fill_between(dfa.index, bottom, dfa[col].values + bottom, label=col, step='mid', color=c) + bottom += dfa[col].values + # only plot Others if there are any + if OTHERS_NAME in dfa.columns and dfa[OTHERS_NAME].sum() > 0.: + ax.fill_between(dfa.index, bottom, dfa[OTHERS_NAME].values + bottom, label=OTHERS_NAME, step='mid', color=OTHERS_COLOR) + + # formatting + ax.tick_params(which='both', direction='in', labelsize=22) + ax.set_xlabel('$\eta$', fontsize=22) + ax.set_ylabel('{} ($R_{{xy}} \leq {:.1f}$ cm)'.format(vt[1], args.path_r), fontsize=22) + ax.xaxis.set_major_locator(MultipleLocator(0.5)) + ax.xaxis.set_minor_locator(MultipleLocator(0.1)) + ax.yaxis.set_major_locator(MultipleLocator(vt[2])) + ax.yaxis.set_minor_locator(MultipleLocator(vt[3])) + ax.grid(ls=':', which='both') + ax.set_axisbelow(False) + ax.set_xlim(args.eta_min, args.eta_max) + ax.set_ylim(0., ax.get_ylim()[1]*1.1) + ax.legend(bbox_to_anchor=(0.0, 0.9, 1.0, 0.1), ncol=2, loc="upper center", fontsize=22, + borderpad=0.2, labelspacing=0.2, columnspacing=0.6, borderaxespad=0.05, handletextpad=0.4) + # ax.set_yscale('log') + fig.savefig(os.path.join(args.outdir, 'material_scan.png')) + exit(0) + + # detector-based scan # array for detailed data # number of materials cannot be pre-determined, so just assign a large number to be safe - vals = np.zeros(shape=(len(etas), len(dets) + 1, 50)) + vals = np.zeros(shape=(len(etas), len(dets) + 1, MAX_N_MATS)) dets2 = dets + [OTHERS_NAME] det_mats = {d: [] for d in dets2} - vt = ALLOWED_VALUE_TYPES[args.value_type] for i, eta in enumerate(etas): if i % PROGRESS_STEP == 0: @@ -286,10 +371,10 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= if sort_dets: sdets = [d for d in dfa.sum().sort_values(ascending=False).index if d != OTHERS_NAME] dfa = dfa[sdets + [OTHERS_NAME]] - dfa.to_csv('material_scan_agg.csv') + dfa.to_csv(os.path.join(args.outdir, 'material_scan_agg.csv')) # plots - pdf = matplotlib.backends.backend_pdf.PdfPages('material_scan_details.pdf') + pdf = matplotlib.backends.backend_pdf.PdfPages(os.path.join(args.outdir, 'material_scan_details.pdf')) fig, ax = plt.subplots(figsize=(16, 8), dpi=160, gridspec_kw={'top': 0.9, 'bottom': 0.2, 'left': 0.08, 'right': 0.98}) @@ -326,7 +411,7 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= ax.legend(bbox_to_anchor=(0.0, 0.9, 1.0, 0.1), ncol=4, loc="upper center", fontsize=22, borderpad=0.2, labelspacing=0.2, columnspacing=0.6, borderaxespad=0.05, handletextpad=0.4) # ax.set_yscale('log') - fig.savefig('material_scan.png') + fig.savefig(os.path.join(args.outdir, 'material_scan.png')) pdf.savefig(fig) plt.close(fig) From 6b770e23f75dfdc0679e3629dab95a6a6e6dbf2d Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Sun, 23 Jul 2023 10:11:00 -0500 Subject: [PATCH 04/11] revert changes to material scan script --- scripts/subdetector_tests/material_scan.py | 97 ++-------------------- 1 file changed, 6 insertions(+), 91 deletions(-) diff --git a/scripts/subdetector_tests/material_scan.py b/scripts/subdetector_tests/material_scan.py index 148b35ed1..dd1afe3ee 100644 --- a/scripts/subdetector_tests/material_scan.py +++ b/scripts/subdetector_tests/material_scan.py @@ -1,7 +1,7 @@ # SPDX-License-Identifier: LGPL-3.0-or-later # Copyright (C) 2023 Chao Peng ''' - A script to scan the materials thickness (X0 or Lambda) by detectors over eta + A script to scan the materials thickness (X0 or Lambda) over eta It uses dd4hep::rec::MaterialManager::placementsBetween to check the material layers, and then uses the detector alignment to transform world coordinates to local coordiantes, and then assigns the materials to a detector based @@ -34,7 +34,6 @@ COLORS = ['royalblue', 'indianred', 'forestgreen', 'gold', 'darkviolet', 'orange', 'darkturquoise'] # a specified color for Others, this should not be included in COLORS OTHERS_COLOR = 'silver' -MAX_N_MATS = 50 # FIXME: this is a work-around to an issue of dd4hep::rec::MaterialManager::placementsBetween @@ -157,10 +156,6 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= dest='compact', help='Top-level xml file of the detector description.' ) - parser.add_argument( - '-o', default='.', dest='outdir', - help='Output directory.' - ) parser.add_argument( '--path-r', type=float, default=400., # 120., help='R_xy (cm) where the scan stops.' @@ -193,10 +188,6 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= '--value-type', default='X0', help='Choose one in {}.'.format(list(ALLOWED_VALUE_TYPES.keys())) ) - parser.add_argument( - '--groupby-materials', action='store_true', - help='Enable it to group the results by materials instead of by detectors.' - ) parser.add_argument( '--detectors', # default='all', @@ -214,8 +205,6 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= print('Cannot find value {}, choose one in {}'.format(args.value_type, list(ALLOWED_VALUE_TYPES.keys()))) exit(-1) - os.makedirs(args.outdir, exist_ok=True) - # scan parameters path_r = args.path_r phi = args.phi @@ -257,86 +246,12 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= th_corr.scan_missing_thickness(desc, start_point, np.array(start_point) + np.array([100., 50., 0.]), dz=0.0001, epsilon=args.epsilon) - # value-type - vt = ALLOWED_VALUE_TYPES[args.value_type] - - # material-based scan - if args.groupby_materials: - vals = np.zeros(shape=(len(etas), MAX_N_MATS)) - mats = [] - for i, eta in enumerate(etas): - if i % PROGRESS_STEP == 0: - print('Scanned {:d}/{:d} for {:.2f} <= eta <= {:.2f}'.format(i, len(etas), etas[0], etas[-1]), - end='\r', flush=True) - # scan material layers - end_x = path_r*np.cos(phi/180.*np.pi) - end_y = path_r*np.sin(phi/180.*np.pi) - end_z = path_r*np.sinh(eta) - # print('({:.2f}, {:.2f}, {:.2f})'.format(end_x, end_y, end_z)) - dfr = material_scan(desc, start_point, (end_x, end_y, end_z), epsilon=args.epsilon, int_dets=dets, thickness_corrector=th_corr) - - # aggregated values for detectors - x0_vals = dfr.groupby('material')[vt[0]].sum().to_dict() - # update material list - for mat in x0_vals.keys(): - if mat not in mats: - mats.append(mat) - # save material scan results - for j, mat in enumerate(mats): - vals[i, j] = x0_vals.get(mat, 0.) - print('Scanned {:d}/{:d} lines for {:.2f} < eta < {:.2f}'.format(len(etas), len(etas), etas[0], etas[-1])) - - # aggregated data for plots - dfa = pd.DataFrame(data=vals[:, :len(mats)], columns=mats, index=etas) - # sort by total thickness - smats = [d for d in dfa.sum().sort_values(ascending=False).index] - dfa.to_csv(os.path.join(args.outdir, 'material_scan_agg.csv')) - - # plot material scan - fig, ax = plt.subplots(figsize=(16, 8), dpi=160, - gridspec_kw={'top': 0.9, 'bottom': 0.2, 'left': 0.08, 'right': 0.98}) - - # group some materials into "Others" if there're no enough colors - plot_mats = smats - if len(smats) > len(COLORS): - group_mats = smats[len(COLORS):] - dfa.loc[:, OTHERS_NAME] += dfa.loc[:, group_mats].sum(axis=1) - print('Materials {} are grouped into {} due to limited number of colors.'.format(group_mats, OTHERS_NAME)) - plot_mats = sdets[:len(COLORS)] - - bottom = np.zeros(len(dfa)) - width = np.mean(np.diff(dfa.index)) - for col, c in zip(plot_mats, COLORS): - ax.fill_between(dfa.index, bottom, dfa[col].values + bottom, label=col, step='mid', color=c) - bottom += dfa[col].values - # only plot Others if there are any - if OTHERS_NAME in dfa.columns and dfa[OTHERS_NAME].sum() > 0.: - ax.fill_between(dfa.index, bottom, dfa[OTHERS_NAME].values + bottom, label=OTHERS_NAME, step='mid', color=OTHERS_COLOR) - - # formatting - ax.tick_params(which='both', direction='in', labelsize=22) - ax.set_xlabel('$\eta$', fontsize=22) - ax.set_ylabel('{} ($R_{{xy}} \leq {:.1f}$ cm)'.format(vt[1], args.path_r), fontsize=22) - ax.xaxis.set_major_locator(MultipleLocator(0.5)) - ax.xaxis.set_minor_locator(MultipleLocator(0.1)) - ax.yaxis.set_major_locator(MultipleLocator(vt[2])) - ax.yaxis.set_minor_locator(MultipleLocator(vt[3])) - ax.grid(ls=':', which='both') - ax.set_axisbelow(False) - ax.set_xlim(args.eta_min, args.eta_max) - ax.set_ylim(0., ax.get_ylim()[1]*1.1) - ax.legend(bbox_to_anchor=(0.0, 0.9, 1.0, 0.1), ncol=2, loc="upper center", fontsize=22, - borderpad=0.2, labelspacing=0.2, columnspacing=0.6, borderaxespad=0.05, handletextpad=0.4) - # ax.set_yscale('log') - fig.savefig(os.path.join(args.outdir, 'material_scan.png')) - exit(0) - - # detector-based scan # array for detailed data # number of materials cannot be pre-determined, so just assign a large number to be safe - vals = np.zeros(shape=(len(etas), len(dets) + 1, MAX_N_MATS)) + vals = np.zeros(shape=(len(etas), len(dets) + 1, 50)) dets2 = dets + [OTHERS_NAME] det_mats = {d: [] for d in dets2} + vt = ALLOWED_VALUE_TYPES[args.value_type] for i, eta in enumerate(etas): if i % PROGRESS_STEP == 0: @@ -371,10 +286,10 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= if sort_dets: sdets = [d for d in dfa.sum().sort_values(ascending=False).index if d != OTHERS_NAME] dfa = dfa[sdets + [OTHERS_NAME]] - dfa.to_csv(os.path.join(args.outdir, 'material_scan_agg.csv')) + dfa.to_csv('material_scan_agg.csv') # plots - pdf = matplotlib.backends.backend_pdf.PdfPages(os.path.join(args.outdir, 'material_scan_details.pdf')) + pdf = matplotlib.backends.backend_pdf.PdfPages('material_scan_details.pdf') fig, ax = plt.subplots(figsize=(16, 8), dpi=160, gridspec_kw={'top': 0.9, 'bottom': 0.2, 'left': 0.08, 'right': 0.98}) @@ -411,7 +326,7 @@ def material_scan(desc, start, end, epsilon, int_dets=None, thickness_corrector= ax.legend(bbox_to_anchor=(0.0, 0.9, 1.0, 0.1), ncol=4, loc="upper center", fontsize=22, borderpad=0.2, labelspacing=0.2, columnspacing=0.6, borderaxespad=0.05, handletextpad=0.4) # ax.set_yscale('log') - fig.savefig(os.path.join(args.outdir, 'material_scan.png')) + fig.savefig('material_scan.png') pdf.savefig(fig) plt.close(fig) From f214831967c41c351affc53403e6e5248196bdc0 Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Sun, 23 Jul 2023 10:57:52 -0500 Subject: [PATCH 05/11] add a script for collecting g4MaterialScan results --- scripts/subdetector_tests/g4_material_scan.py | 103 ++++++++++++++++++ 1 file changed, 103 insertions(+) create mode 100644 scripts/subdetector_tests/g4_material_scan.py diff --git a/scripts/subdetector_tests/g4_material_scan.py b/scripts/subdetector_tests/g4_material_scan.py new file mode 100644 index 000000000..352e2a6b5 --- /dev/null +++ b/scripts/subdetector_tests/g4_material_scan.py @@ -0,0 +1,103 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) 2023 Chao Peng +''' + A script to use g4MaterialScan and collect the results +''' + +import os +import math +import json +import fnmatch +import argparse +import subprocess +from io import StringIO +import pandas as pd +import numpy as np + +pd.set_option('display.max_rows', 1000) +PROGRESS_STEP = 10 + + +''' + A parser function to convert output from g4MaterialScan to a pandas dataframe + compact: is path to compact file + start_point: a list or tuple for 3D coordinates (e.g., 0,0,0) + direction: a lit or tuple for direction (e.g., 0.1,0.2,1.0) + return: a dataframe for materialScan results +''' +def g4_material_scan(compact, start_point, direction, timeout=200): + EXEC_NAME = 'g4MaterialScan' + cmd = '{} --compact={} --position=\"{}\" --direction=\"{}\"'\ + .format(EXEC_NAME, compact, start_point, direction) + output = os.popen(cmd).read() + # output = subprocess.check_output(cmd.split(' '), timeout=timeout).decode('UTF-8') + + # find material scan lines + lines = [] + add_line = False + + for l in output.split('\n'): + if add_line: + lines.append(l.strip()) + if l.strip().startswith('MaterialScan'): + add_line = True + + + # NOTE: the code below is for materialScan output as of 03/05/2023 + # it may need change if the output format is changed + scans = [] + first_digit = False + for i, l in enumerate(lines): + line = l.strip('| ') + if not first_digit and not line[:1].isdigit(): + continue + first_digit = True + if not line[:1].isdigit(): + break + # break coordinates for endpoint, which has a format of (x, y, z) + scans.append(line.strip('| ').translate({ord(i): None for i in '()'}).replace(',', ' ')) + + cols = [ + 'material', 'Z', 'A', 'density', + 'rad_length', 'int_length', 'thickness', 'path_length', + 'int_X0', 'int_lamda', 'end_x', 'end_y', 'end_z' + ] + + dft = pd.read_csv(StringIO('\n'.join(scans)), sep='\s+', header=None, index_col=0, names=cols) + return dft.astype({key: np.float64 for key in cols[1:]}) + + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument( + dest='compact', + help='Top-level xml file of the detector description.' + ) + parser.add_argument( + '--start-point', default='0,0,0', + help='Start point of the scan, use the format \"x,y,z\", unit is cm.' + ) + parser.add_argument( + '--eta-min', type=float, default=-2.0, + help='Minimum eta for the scan.' + ) + parser.add_argument( + '--eta-max', type=float, default=2.0, + help='Minimum eta for the scan.' + ) + parser.add_argument( + '--eta-nbins', type=int, default=401, + help='Number of bins for the eta scan.' + ) + parser.add_argument( + '--phi', type=float, default=20., + help='Phi angle of the scan, unit is degree.' + ) + args = parser.parse_args() + + if not os.path.exists(args.compact): + print('Cannot find {}'.format(args.compact)) + exit(-1) + + df = g4_material_scan(args.compact, args.start_point, '0.2,0.01,1.0') + print(df) From a31e124ce6e81bc3e164c6ffa8955e636d692ff2 Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Fri, 25 Aug 2023 13:52:24 -0500 Subject: [PATCH 06/11] update the scan script --- scripts/subdetector_tests/g4_material_scan.py | 41 ++++++++++++++----- 1 file changed, 30 insertions(+), 11 deletions(-) diff --git a/scripts/subdetector_tests/g4_material_scan.py b/scripts/subdetector_tests/g4_material_scan.py index 352e2a6b5..9a09f2999 100644 --- a/scripts/subdetector_tests/g4_material_scan.py +++ b/scripts/subdetector_tests/g4_material_scan.py @@ -5,14 +5,12 @@ ''' import os -import math -import json -import fnmatch import argparse import subprocess -from io import StringIO import pandas as pd import numpy as np +from io import StringIO +from collections import OrderedDict as odict pd.set_option('display.max_rows', 1000) PROGRESS_STEP = 10 @@ -27,8 +25,8 @@ ''' def g4_material_scan(compact, start_point, direction, timeout=200): EXEC_NAME = 'g4MaterialScan' - cmd = '{} --compact={} --position=\"{}\" --direction=\"{}\"'\ - .format(EXEC_NAME, compact, start_point, direction) + cmd = '{} --compact={} --position=\"{},{},{}\" --direction=\"{},{},{}\"'\ + .format(EXEC_NAME, compact, *np.hstack([start_point, direction])) output = os.popen(cmd).read() # output = subprocess.check_output(cmd.split(' '), timeout=timeout).decode('UTF-8') @@ -73,20 +71,24 @@ def g4_material_scan(compact, start_point, direction, timeout=200): dest='compact', help='Top-level xml file of the detector description.' ) + parser.add_argument( + '-o', '--output', default='g4_materials.csv', + help='A csv file contains the scan info.' + ) parser.add_argument( '--start-point', default='0,0,0', help='Start point of the scan, use the format \"x,y,z\", unit is cm.' ) parser.add_argument( - '--eta-min', type=float, default=-2.0, + '--eta-min', type=float, default=-4.0, help='Minimum eta for the scan.' ) parser.add_argument( - '--eta-max', type=float, default=2.0, + '--eta-max', type=float, default=4.0, help='Minimum eta for the scan.' ) parser.add_argument( - '--eta-nbins', type=int, default=401, + '--eta-nbins', type=int, default=801, help='Number of bins for the eta scan.' ) parser.add_argument( @@ -99,5 +101,22 @@ def g4_material_scan(compact, start_point, direction, timeout=200): print('Cannot find {}'.format(args.compact)) exit(-1) - df = g4_material_scan(args.compact, args.start_point, '0.2,0.01,1.0') - print(df) + start_point = np.array([float(v.strip()) for v in args.start_point.split(',')]) + etas = np.linspace(args.eta_min, args.eta_max, args.eta_nbins) + mats_indices = odict() + # a data buffer for the X0 values of (eta, material), 50 should be large enough + data = np.zeros(shape=(len(etas), 50)) + for i, eta in enumerate(etas): + if i % PROGRESS_STEP == 0: + print('Scanned {:d}/{:d} lines for {:.2f} < eta < {:.2f}'.format(i, len(etas), etas[0], etas[-1])) + direction = (np.cos(args.phi/180.*np.pi), np.sin(args.phi/180.*np.pi), np.sinh(eta)) + dfa = g4_material_scan(args.compact, start_point, direction) + dfa.loc[:, 'X0'] = dfa['int_X0'].diff(1).fillna(dfa['int_X0']) + for mat, xval in dfa.groupby('material')['X0'].sum().items(): + if mat not in mats_indices: + mats_indices[mat] = len(mats_indices) + j = mats_indices.get(mat) + data[i, j] = xval + + result = pd.DataFrame(columns=mats_indices.keys(), index=etas, data=data[:, :len(mats_indices)]) + result.to_csv(args.output) From f9466b17016d5b816f4098050773148a7ba0517e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 25 Aug 2023 18:52:45 +0000 Subject: [PATCH 07/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/subdetector_tests/g4_material_scan.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/subdetector_tests/g4_material_scan.py b/scripts/subdetector_tests/g4_material_scan.py index 9a09f2999..5df65cbbc 100644 --- a/scripts/subdetector_tests/g4_material_scan.py +++ b/scripts/subdetector_tests/g4_material_scan.py @@ -33,7 +33,7 @@ def g4_material_scan(compact, start_point, direction, timeout=200): # find material scan lines lines = [] add_line = False - + for l in output.split('\n'): if add_line: lines.append(l.strip()) From 86edaa4958324e8358b7ac7d9c3212cb581a2b38 Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Sat, 26 Aug 2023 14:04:55 -0500 Subject: [PATCH 08/11] add a feature to scan over a phi range (take the average value) --- scripts/subdetector_tests/g4_material_scan.py | 61 +++++++++++++------ 1 file changed, 43 insertions(+), 18 deletions(-) diff --git a/scripts/subdetector_tests/g4_material_scan.py b/scripts/subdetector_tests/g4_material_scan.py index 5df65cbbc..321fef634 100644 --- a/scripts/subdetector_tests/g4_material_scan.py +++ b/scripts/subdetector_tests/g4_material_scan.py @@ -65,6 +65,23 @@ def g4_material_scan(compact, start_point, direction, timeout=200): return dft.astype({key: np.float64 for key in cols[1:]}) +''' + A helper function to convert a string ([:][:]) to an array +''' +def args_array(arg, step=1, include_end=True): + vals = [float(x.strip()) for x in arg.split(':')] + # empty or only one value + if len(vals) < 2: + return np.array(vals) + # has step input + if len(vals) > 2: + step = vals[2] + # inclusion of the endpoint (max) + if include_end: + vals[1] += step + return np.arange(vals[0], vals[1], step) + + if __name__ == '__main__': parser = argparse.ArgumentParser() parser.add_argument( @@ -80,20 +97,12 @@ def g4_material_scan(compact, start_point, direction, timeout=200): help='Start point of the scan, use the format \"x,y,z\", unit is cm.' ) parser.add_argument( - '--eta-min', type=float, default=-4.0, - help='Minimum eta for the scan.' + '--eta', default='-4.0:4.0:0.1', + help='Eta range, in the format of \"[:][:]\".' ) parser.add_argument( - '--eta-max', type=float, default=4.0, - help='Minimum eta for the scan.' - ) - parser.add_argument( - '--eta-nbins', type=int, default=801, - help='Number of bins for the eta scan.' - ) - parser.add_argument( - '--phi', type=float, default=20., - help='Phi angle of the scan, unit is degree.' + '--phi', default='0:30:1', + help='Phi angle range, in the format of \"[:][:]\" (degree).' ) args = parser.parse_args() @@ -102,17 +111,33 @@ def g4_material_scan(compact, start_point, direction, timeout=200): exit(-1) start_point = np.array([float(v.strip()) for v in args.start_point.split(',')]) - etas = np.linspace(args.eta_min, args.eta_max, args.eta_nbins) + etas = args_array(args.eta) + phis = args_array(args.phi) + # sanity check + if not len(phis): + print('No phi values from the input {}, aborted!'.format(args.phi)) + exit(-1) mats_indices = odict() # a data buffer for the X0 values of (eta, material), 50 should be large enough data = np.zeros(shape=(len(etas), 50)) + # scan eta for i, eta in enumerate(etas): if i % PROGRESS_STEP == 0: - print('Scanned {:d}/{:d} lines for {:.2f} < eta < {:.2f}'.format(i, len(etas), etas[0], etas[-1])) - direction = (np.cos(args.phi/180.*np.pi), np.sin(args.phi/180.*np.pi), np.sinh(eta)) - dfa = g4_material_scan(args.compact, start_point, direction) - dfa.loc[:, 'X0'] = dfa['int_X0'].diff(1).fillna(dfa['int_X0']) - for mat, xval in dfa.groupby('material')['X0'].sum().items(): + print('Scanned {:d}/{:d} lines for {:.2f} < eta < {:.2f}, each line contains {:d} phi angles from {:.2f} to {:.2f}'\ + .format(i, len(etas), etas[0], etas[-1], len(phis), phis[0], phis[-1]) + ) + # average over phi + eta_scan = pd.DataFrame() + for phi in phis: + direction = (np.cos(phi/180.*np.pi), np.sin(phi/180.*np.pi), np.sinh(eta)) + dfa = g4_material_scan(args.compact, start_point, direction) + dfa.loc[:, 'X0'] = dfa['int_X0'].diff(1).fillna(dfa['int_X0']) + phi_scan = dfa.groupby('material')['X0'].sum().to_frame(name=phi) + # using pd.DataFrame.merge to combine results, using outer so new materials at specific phi angle can be added in + eta_scan = eta_scan.merge(phi_scan, how='outer', left_index=True, right_index=True) + # print(eta_scan) + # print(eta_scan.sum(axis=1)) + for mat, xval in eta_scan.sum(axis=1).items(): if mat not in mats_indices: mats_indices[mat] = len(mats_indices) j = mats_indices.get(mat) From b38b79b91837916cac0bf7531568dad486166ab2 Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Sat, 26 Aug 2023 14:31:25 -0500 Subject: [PATCH 09/11] fix an issue with calculating average material thickness --- scripts/subdetector_tests/g4_material_scan.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/subdetector_tests/g4_material_scan.py b/scripts/subdetector_tests/g4_material_scan.py index 321fef634..a48534a17 100644 --- a/scripts/subdetector_tests/g4_material_scan.py +++ b/scripts/subdetector_tests/g4_material_scan.py @@ -137,7 +137,8 @@ def args_array(arg, step=1, include_end=True): eta_scan = eta_scan.merge(phi_scan, how='outer', left_index=True, right_index=True) # print(eta_scan) # print(eta_scan.sum(axis=1)) - for mat, xval in eta_scan.sum(axis=1).items(): + # replace nan value with 0 and then calculate average values + for mat, xval in eta_scan.fillna(0.).mean(axis=1).items(): if mat not in mats_indices: mats_indices[mat] = len(mats_indices) j = mats_indices.get(mat) From e390f99b38ea378230b7a4c81fa7650cfcdabc4d Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Sat, 26 Aug 2023 14:33:42 -0500 Subject: [PATCH 10/11] fix a typo --- scripts/subdetector_tests/g4_material_scan.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/subdetector_tests/g4_material_scan.py b/scripts/subdetector_tests/g4_material_scan.py index a48534a17..949e879ea 100644 --- a/scripts/subdetector_tests/g4_material_scan.py +++ b/scripts/subdetector_tests/g4_material_scan.py @@ -98,11 +98,11 @@ def args_array(arg, step=1, include_end=True): ) parser.add_argument( '--eta', default='-4.0:4.0:0.1', - help='Eta range, in the format of \"[:][:]\".' + help='Eta range, in the format of \"[:[:]]\".' ) parser.add_argument( '--phi', default='0:30:1', - help='Phi angle range, in the format of \"[:][:]\" (degree).' + help='Phi angle range, in the format of \"[:[:]]\" (degree).' ) args = parser.parse_args() From ad4d2dac3144ad538b69b763392d85ed37d36e37 Mon Sep 17 00:00:00 2001 From: Chao Peng Date: Sat, 26 Aug 2023 14:35:06 -0500 Subject: [PATCH 11/11] remove g4_material_scan.py, it will be added in a separate PR --- scripts/subdetector_tests/g4_material_scan.py | 148 ------------------ 1 file changed, 148 deletions(-) delete mode 100644 scripts/subdetector_tests/g4_material_scan.py diff --git a/scripts/subdetector_tests/g4_material_scan.py b/scripts/subdetector_tests/g4_material_scan.py deleted file mode 100644 index 949e879ea..000000000 --- a/scripts/subdetector_tests/g4_material_scan.py +++ /dev/null @@ -1,148 +0,0 @@ -# SPDX-License-Identifier: LGPL-3.0-or-later -# Copyright (C) 2023 Chao Peng -''' - A script to use g4MaterialScan and collect the results -''' - -import os -import argparse -import subprocess -import pandas as pd -import numpy as np -from io import StringIO -from collections import OrderedDict as odict - -pd.set_option('display.max_rows', 1000) -PROGRESS_STEP = 10 - - -''' - A parser function to convert output from g4MaterialScan to a pandas dataframe - compact: is path to compact file - start_point: a list or tuple for 3D coordinates (e.g., 0,0,0) - direction: a lit or tuple for direction (e.g., 0.1,0.2,1.0) - return: a dataframe for materialScan results -''' -def g4_material_scan(compact, start_point, direction, timeout=200): - EXEC_NAME = 'g4MaterialScan' - cmd = '{} --compact={} --position=\"{},{},{}\" --direction=\"{},{},{}\"'\ - .format(EXEC_NAME, compact, *np.hstack([start_point, direction])) - output = os.popen(cmd).read() - # output = subprocess.check_output(cmd.split(' '), timeout=timeout).decode('UTF-8') - - # find material scan lines - lines = [] - add_line = False - - for l in output.split('\n'): - if add_line: - lines.append(l.strip()) - if l.strip().startswith('MaterialScan'): - add_line = True - - - # NOTE: the code below is for materialScan output as of 03/05/2023 - # it may need change if the output format is changed - scans = [] - first_digit = False - for i, l in enumerate(lines): - line = l.strip('| ') - if not first_digit and not line[:1].isdigit(): - continue - first_digit = True - if not line[:1].isdigit(): - break - # break coordinates for endpoint, which has a format of (x, y, z) - scans.append(line.strip('| ').translate({ord(i): None for i in '()'}).replace(',', ' ')) - - cols = [ - 'material', 'Z', 'A', 'density', - 'rad_length', 'int_length', 'thickness', 'path_length', - 'int_X0', 'int_lamda', 'end_x', 'end_y', 'end_z' - ] - - dft = pd.read_csv(StringIO('\n'.join(scans)), sep='\s+', header=None, index_col=0, names=cols) - return dft.astype({key: np.float64 for key in cols[1:]}) - - -''' - A helper function to convert a string ([:][:]) to an array -''' -def args_array(arg, step=1, include_end=True): - vals = [float(x.strip()) for x in arg.split(':')] - # empty or only one value - if len(vals) < 2: - return np.array(vals) - # has step input - if len(vals) > 2: - step = vals[2] - # inclusion of the endpoint (max) - if include_end: - vals[1] += step - return np.arange(vals[0], vals[1], step) - - -if __name__ == '__main__': - parser = argparse.ArgumentParser() - parser.add_argument( - dest='compact', - help='Top-level xml file of the detector description.' - ) - parser.add_argument( - '-o', '--output', default='g4_materials.csv', - help='A csv file contains the scan info.' - ) - parser.add_argument( - '--start-point', default='0,0,0', - help='Start point of the scan, use the format \"x,y,z\", unit is cm.' - ) - parser.add_argument( - '--eta', default='-4.0:4.0:0.1', - help='Eta range, in the format of \"[:[:]]\".' - ) - parser.add_argument( - '--phi', default='0:30:1', - help='Phi angle range, in the format of \"[:[:]]\" (degree).' - ) - args = parser.parse_args() - - if not os.path.exists(args.compact): - print('Cannot find {}'.format(args.compact)) - exit(-1) - - start_point = np.array([float(v.strip()) for v in args.start_point.split(',')]) - etas = args_array(args.eta) - phis = args_array(args.phi) - # sanity check - if not len(phis): - print('No phi values from the input {}, aborted!'.format(args.phi)) - exit(-1) - mats_indices = odict() - # a data buffer for the X0 values of (eta, material), 50 should be large enough - data = np.zeros(shape=(len(etas), 50)) - # scan eta - for i, eta in enumerate(etas): - if i % PROGRESS_STEP == 0: - print('Scanned {:d}/{:d} lines for {:.2f} < eta < {:.2f}, each line contains {:d} phi angles from {:.2f} to {:.2f}'\ - .format(i, len(etas), etas[0], etas[-1], len(phis), phis[0], phis[-1]) - ) - # average over phi - eta_scan = pd.DataFrame() - for phi in phis: - direction = (np.cos(phi/180.*np.pi), np.sin(phi/180.*np.pi), np.sinh(eta)) - dfa = g4_material_scan(args.compact, start_point, direction) - dfa.loc[:, 'X0'] = dfa['int_X0'].diff(1).fillna(dfa['int_X0']) - phi_scan = dfa.groupby('material')['X0'].sum().to_frame(name=phi) - # using pd.DataFrame.merge to combine results, using outer so new materials at specific phi angle can be added in - eta_scan = eta_scan.merge(phi_scan, how='outer', left_index=True, right_index=True) - # print(eta_scan) - # print(eta_scan.sum(axis=1)) - # replace nan value with 0 and then calculate average values - for mat, xval in eta_scan.fillna(0.).mean(axis=1).items(): - if mat not in mats_indices: - mats_indices[mat] = len(mats_indices) - j = mats_indices.get(mat) - data[i, j] = xval - - result = pd.DataFrame(columns=mats_indices.keys(), index=etas, data=data[:, :len(mats_indices)]) - result.to_csv(args.output)