Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-03 07:39:02

0001 import json
0002 
0003 
0004 def dumper(obj):
0005     try:
0006         return obj.toJSON()
0007     except:
0008         return obj.__dict__
0009 
0010 
0011 class index_info(object):
0012     def __init__(self, index):
0013         self.index = index
0014         self.boundaries = []
0015         self.layers_no_approach = []
0016         self.layers_with_approach = []
0017         self.approaches = []
0018         self.passiveDiscBinningR = 0
0019         self.passiveDiscBinningPhi = 0
0020         self.passiveCylinderBinningZ = 0
0021         self.passiveCylinderBinningPhi = 0
0022         self.activeBinningRorZ = 0
0023         self.activeBinningPhi = 0
0024 
0025     def __repr__(self):
0026         return repr(
0027             (
0028                 self.index,
0029                 self.boundaries,
0030                 self.layers_no_approach,
0031                 self.layers_with_approach,
0032                 self.approaches,
0033                 self.passiveDiscBinningR,
0034                 self.passiveDiscBinningPhi,
0035                 self.passiveCylinderBinningZ,
0036                 self.passiveCylinderBinningPhi,
0037                 self.activeBinningRorZ,
0038                 self.activeBinningPhi,
0039             )
0040         )
0041 
0042 
0043 def append_index_if_missing(dictionary, name, index):
0044     if name not in dictionary:
0045         dictionary[name] = index_info(index)
0046 
0047 
0048 def extract_coords(coords, is_disc):
0049     x_values = [coords[1], coords[2]]
0050     y_values = [coords[0], coords[0]]
0051     if is_disc:
0052         x_values = [coords[2], coords[2]]
0053         y_values = [coords[1], coords[0]]
0054     return x_values, y_values
0055 
0056 
0057 def dump_geo(filename, plot, output_folder, dump_steering, steering_file):
0058     f = open(filename)
0059     data = json.load(f)
0060 
0061     # First you want to read the Volume entries and save the relevant quantities
0062     # - names
0063     index_to_names = []
0064     for entry in data["Volumes"]["entries"]:
0065         index_to_names.append(entry["value"]["NAME"])
0066 
0067     # In case volume names are not found check in the volume names with dummy values
0068     if not index_to_names:
0069         for entry in data["Surfaces"]["entries"]:
0070             if "volume" in entry:
0071                 vol = entry["volume"]
0072                 if "volume" + str(vol) not in index_to_names:
0073                     index_to_names.append("volume" + str(vol))
0074 
0075     # Stores the information to easily decide on which layers you want the material to be mapped
0076     steering_map = {}
0077 
0078     # Once you have read the volumes extensions, you read the surfaces representing layers and boundaries.
0079     # All these surfaces belong to a volume, they have therefore a volume index.
0080     # You are interested only in:
0081     # - surfaces with layer index (NO approach index):
0082     #    -- even layer indices represent "active" layers, i.e. the ones corresponding to sensitive layers
0083     #    -- odd event indices represent navigation layers
0084     # - surfaces with layer index AND approach index:
0085     #    -- indicate approach layers to "active" layers:
0086     #       e.g. it can happen that for cylinders: 1 sits in front of the active layer,
0087     #                                              2 sits in the middle of the layer,
0088     #                                              3 sits behind the layer.
0089     #                                              ...
0090     # - surfaces with boundary index (no layer index in this case):
0091     #    -- they indicate boundaries between volumes. You are interested in boundaries between volumes
0092     #       containing at least an active layer.
0093 
0094     index_to_extends_layers_bounds_cylinders = [[] for _ in range(len(index_to_names))]
0095     index_to_extends_layers_bounds_discs = [[] for _ in range(len(index_to_names))]
0096 
0097     index_to_extends_layers_cylinders = [[] for _ in range(len(index_to_names))]
0098     index_to_extends_layers_discs = [[] for _ in range(len(index_to_names))]
0099 
0100     for entry in data["Surfaces"]["entries"]:
0101         if "layer" in entry:
0102             extends = []
0103             vol = entry["volume"]
0104 
0105             if "sensitive" in entry:
0106                 continue
0107 
0108             # Filling extends with the following quantities:
0109             # for cylinders:
0110             #     radius, z min, z max, approach index, layer index
0111             # for discs:
0112             #     inner radius, outer radius, z position, approach index, layer index
0113 
0114             z_shift = 0.0
0115             if entry["value"]["transform"]["translation"] != None:
0116                 z_shift = entry["value"]["transform"]["translation"][2]
0117 
0118             approach_index = -1
0119             if "approach" in entry:
0120                 approach_index = entry["approach"]
0121 
0122             if entry["value"]["type"] == "CylinderSurface":
0123                 # radius, z min, z max, approach index, layer index
0124                 extends = [
0125                     entry["value"]["bounds"]["values"][0],
0126                     z_shift - entry["value"]["bounds"]["values"][1],
0127                     z_shift + entry["value"]["bounds"]["values"][1],
0128                     approach_index,
0129                     entry["layer"],
0130                 ]
0131                 index_to_extends_layers_cylinders[vol - 1].append(extends)
0132             elif entry["value"]["type"] == "DiscSurface":
0133                 # inner radius, outer radius, z position, approach index, layer index
0134                 extends = [
0135                     entry["value"]["bounds"]["values"][0],
0136                     entry["value"]["bounds"]["values"][1],
0137                     z_shift,
0138                     approach_index,
0139                     entry["layer"],
0140                 ]
0141                 index_to_extends_layers_discs[vol - 1].append(extends)
0142             else:
0143                 surface_type = entry["value"]["type"]
0144                 print(
0145                     f"WARNING: Processing surface with unknown type '{surface_type}'. Only CylinderSurface and DiscSurface are considered."
0146                 )
0147 
0148         if "boundary" in entry:
0149             extends = []
0150             vol = entry["volume"]
0151 
0152             # Filling extends with the following quantities:
0153             # for cylinders:
0154             #     radius, z min, z max, boundary index
0155             # for discs:
0156             #     inner radius, outer radius, z position, boundary index
0157 
0158             z_shift = 0.0
0159             if entry["value"]["transform"]["translation"] != None:
0160                 z_shift = entry["value"]["transform"]["translation"][2]
0161 
0162             if entry["value"]["type"] == "CylinderSurface":
0163                 extends = [
0164                     entry["value"]["bounds"]["values"][0],
0165                     z_shift - entry["value"]["bounds"]["values"][1],
0166                     z_shift + entry["value"]["bounds"]["values"][1],
0167                     entry["boundary"],
0168                 ]
0169                 index_to_extends_layers_bounds_cylinders[vol - 1].append(extends)
0170             elif entry["value"]["type"] == "DiscSurface":
0171                 extends = [
0172                     entry["value"]["bounds"]["values"][0],
0173                     entry["value"]["bounds"]["values"][1],
0174                     z_shift,
0175                     entry["boundary"],
0176                 ]
0177                 index_to_extends_layers_bounds_discs[vol - 1].append(extends)
0178             else:
0179                 surface_type = entry["value"]["type"]
0180                 print(
0181                     f"WARNING: Processing surface with unknown type '{surface_type}'. Only CylinderSurface and DiscSurface are considered."
0182                 )
0183 
0184     # Steering the information and collect it into an output file if needed
0185     from itertools import chain
0186 
0187     interesting_volumes = []
0188     v_index = 0
0189     is_disc = False
0190     for elements in chain(
0191         index_to_extends_layers_cylinders, index_to_extends_layers_discs
0192     ):
0193         for coords in elements:
0194             if coords[3] > 0:
0195                 continue
0196             if v_index not in interesting_volumes:
0197                 interesting_volumes.append(v_index)
0198             append_index_if_missing(steering_map, index_to_names[v_index], v_index + 1)
0199             steering_map[index_to_names[v_index]].layers_no_approach.append(coords[4])
0200         v_index = v_index + 1
0201         if not is_disc and v_index == len(index_to_extends_layers_cylinders):
0202             v_index = 0
0203             is_disc = True
0204 
0205     v_index = 0
0206     is_disc = False
0207     for elements in chain(
0208         index_to_extends_layers_bounds_cylinders, index_to_extends_layers_bounds_discs
0209     ):
0210         for coords in elements:
0211             if v_index in interesting_volumes:
0212                 append_index_if_missing(
0213                     steering_map, index_to_names[v_index], v_index + 1
0214                 )
0215                 steering_map[index_to_names[v_index]].boundaries.append(coords[3])
0216         v_index = v_index + 1
0217         if not is_disc and v_index == len(index_to_extends_layers_bounds_cylinders):
0218             v_index = 0
0219             is_disc = True
0220 
0221     v_index = 0
0222     is_disc = False
0223     for elements in chain(
0224         index_to_extends_layers_cylinders, index_to_extends_layers_discs
0225     ):
0226         for coords in elements:
0227             if coords[3] == -1:
0228                 continue
0229             if (
0230                 coords[4]
0231                 not in steering_map[index_to_names[v_index]].layers_with_approach
0232             ):
0233                 steering_map[index_to_names[v_index]].layers_with_approach.append(
0234                     coords[4]
0235                 )
0236             if coords[3] not in steering_map[index_to_names[v_index]].approaches:
0237                 steering_map[index_to_names[v_index]].approaches.append(coords[3])
0238         v_index = v_index + 1
0239         if not is_disc and v_index == len(index_to_extends_layers_bounds_cylinders):
0240             v_index = 0
0241             is_disc = True
0242 
0243     if dump_steering:
0244         output_map = {"SteeringField": steering_map}
0245         with open(steering_file, "w", encoding="utf-8") as f:
0246             json.dump(output_map, f, default=dumper, ensure_ascii=False, indent=4)
0247 
0248     # Once you have all the data in hands, you make a few nice plots to show volumes and surfaces
0249     if plot:
0250         import matplotlib.pyplot as plt
0251 
0252         plt.rcParams.update({"figure.max_open_warning": 0})
0253         from matplotlib.pyplot import cm
0254         import numpy as np
0255 
0256         color = cm.rainbow(np.linspace(0, 1, len(index_to_extends_layers_cylinders)))
0257 
0258         is_in_legend = []
0259 
0260         plt.figure(figsize=(20, 10))
0261 
0262         # Plot all layers (w/o approach layers) + boundaries
0263         v_index = 0
0264         is_disc = False
0265         for elements in chain(
0266             index_to_extends_layers_cylinders, index_to_extends_layers_discs
0267         ):
0268             for coords in elements:
0269                 # skip approach layers
0270                 if coords[3] > 0:
0271                     continue
0272                 x_values, y_values = extract_coords(coords, is_disc)
0273                 if index_to_names[v_index] not in is_in_legend:
0274                     plt.plot(
0275                         x_values,
0276                         y_values,
0277                         c=color[v_index],
0278                         label="v: " + str(v_index + 1) + ", " + index_to_names[v_index],
0279                     )
0280                     is_in_legend.append(index_to_names[v_index])
0281                 else:
0282                     plt.plot(x_values, y_values, c=color[v_index])
0283             v_index = v_index + 1
0284             if not is_disc and v_index == len(index_to_extends_layers_cylinders):
0285                 v_index = 0
0286                 is_disc = True
0287 
0288         v_index = 0
0289         is_disc = False
0290         for elements in chain(
0291             index_to_extends_layers_bounds_cylinders,
0292             index_to_extends_layers_bounds_discs,
0293         ):
0294             for coords in elements:
0295                 if v_index in interesting_volumes:
0296                     x_values, y_values = extract_coords(coords, is_disc)
0297                     plt.plot(x_values, y_values, c=color[v_index])
0298             v_index = v_index + 1
0299             if not is_disc and v_index == len(index_to_extends_layers_bounds_cylinders):
0300                 v_index = 0
0301                 is_disc = True
0302 
0303         plt.xlabel("z [mm]")
0304         plt.ylabel("R [mm]")
0305         plt.title("Volumes and Layers (no approach layers)")
0306         plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
0307         plt.savefig(output_folder + "/volumes_and_layers.png", bbox_inches="tight")
0308 
0309         # Plot each volume: layers + approach layers
0310         v_index = 0
0311         is_disc = False
0312         approach_colors = ["black", "blue", "red", "green", "orange", "purple", "pink"]
0313         for elements in chain(
0314             index_to_extends_layers_cylinders, index_to_extends_layers_discs
0315         ):
0316             l_index = 0
0317             if not elements:
0318                 v_index = v_index + 1
0319                 continue
0320             plt.figure(figsize=(20, 10))
0321             color_layers = cm.rainbow(np.linspace(0, 1, len(elements)))
0322             has_elements = False
0323             for coords in elements:
0324                 if coords[3] > 0:
0325                     continue
0326                 has_elements = True
0327                 x_values, y_values = extract_coords(coords, is_disc)
0328                 plt.plot(
0329                     x_values,
0330                     y_values,
0331                     c=color_layers[l_index],
0332                     label="l: " + str(coords[4]),
0333                 )
0334                 l_index = l_index + 1
0335 
0336                 a_is_disc = False
0337                 count = 0
0338                 for a_coords in chain(
0339                     index_to_extends_layers_cylinders[v_index],
0340                     index_to_extends_layers_discs[v_index],
0341                 ):
0342                     if a_coords[4] == coords[4] and a_coords[3] > 0:
0343                         ax_values, ay_values = extract_coords(a_coords, a_is_disc)
0344                         plt.plot(
0345                             ax_values,
0346                             ay_values,
0347                             linestyle=(0, (5, 10)),
0348                             c=approach_colors[a_coords[3]],
0349                             label="l: " + str(coords[4]) + ", a: " + str(a_coords[3]),
0350                         )
0351                     count = count + 1
0352                     if count == len(index_to_extends_layers_cylinders[v_index]):
0353                         a_is_disc = True
0354 
0355             v_index = v_index + 1
0356             if not is_disc and v_index == len(index_to_extends_layers_cylinders):
0357                 v_index = 0
0358                 is_disc = True
0359 
0360             if not has_elements:
0361                 continue
0362 
0363             plt.xlabel("z [mm]")
0364             plt.ylabel("R [mm]")
0365             plt.title(index_to_names[v_index - 1])
0366             plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
0367             plt.savefig(output_folder + "/layers_for_volume_" + str(v_index) + ".png")
0368 
0369         plt.figure(figsize=(20, 10))
0370 
0371         # Plot boundaries
0372         v_index = 0
0373         is_disc = False
0374         for elements in chain(
0375             index_to_extends_layers_bounds_cylinders,
0376             index_to_extends_layers_bounds_discs,
0377         ):
0378             for coords in elements:
0379                 x_values, y_values = extract_coords(coords, is_disc)
0380                 if v_index in interesting_volumes:
0381                     plt.plot(
0382                         x_values,
0383                         y_values,
0384                         linestyle="--",
0385                         c=color[v_index],
0386                         label="v: " + str(v_index + 1) + ", b: " + str(coords[3]),
0387                     )
0388             v_index = v_index + 1
0389             if not is_disc and v_index == len(index_to_extends_layers_bounds_cylinders):
0390                 v_index = 0
0391                 is_disc = True
0392 
0393         plt.xlabel("z [mm]")
0394         plt.ylabel("R [mm]")
0395         plt.title("Boundary surfaces")
0396         plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
0397         plt.savefig(output_folder + "/boundaries.png")
0398 
0399         plt.figure(figsize=(20, 10))
0400 
0401         # Plot all approach layers
0402         v_index = 0
0403         is_disc = False
0404         add_to_legend = []
0405         for elements in chain(
0406             index_to_extends_layers_cylinders, index_to_extends_layers_discs
0407         ):
0408             if not elements:
0409                 v_index = v_index + 1
0410                 continue
0411             for coords in elements:
0412                 if coords[3] == -1:
0413                     continue
0414                 x_values, y_values = extract_coords(coords, is_disc)
0415                 if coords[3] not in add_to_legend:
0416                     plt.plot(
0417                         x_values,
0418                         y_values,
0419                         c=approach_colors[coords[3]],
0420                         linestyle="--",
0421                         label="approach index = " + str(coords[3]),
0422                     )
0423                     add_to_legend.append(coords[3])
0424                 else:
0425                     plt.plot(
0426                         x_values, y_values, c=approach_colors[coords[3]], linestyle="--"
0427                     )
0428             v_index = v_index + 1
0429             if not is_disc and v_index == len(index_to_extends_layers_bounds_cylinders):
0430                 v_index = 0
0431                 is_disc = True
0432 
0433         plt.xlabel("z [mm]")
0434         plt.ylabel("R [mm]")
0435         plt.title("Approach layers")
0436         plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
0437         plt.savefig(output_folder + "/approach_layers.png")
0438 
0439 
0440 def read_and_modify(filename, plot, output_folder, steering_file, output_file):
0441     with open(filename) as f:
0442         data = json.load(f)
0443 
0444     with open(steering_file) as f:
0445         full_data = json.load(f)
0446 
0447     layer_data = full_data["SteeringField"]
0448 
0449     index_to_names = []
0450 
0451     # Allows to dump the binning defined for the material layers
0452     dump_binning_for_material = False
0453 
0454     # Allows to check which layers are configured to carry material
0455     check_material_layers = True
0456 
0457     # First you want to read the Volume entries and save the relevant quantities
0458     # - names
0459     for entry in data["Volumes"]["entries"]:
0460         index_to_names.append(entry["value"]["NAME"])
0461 
0462     # In case volume names are not found check in the volume names with dummy values
0463     if not index_to_names:
0464         for entry in data["Surfaces"]["entries"]:
0465             if "volume" in entry:
0466                 vol = entry["volume"]
0467                 if "volume" + str(vol) not in index_to_names:
0468                     index_to_names.append("volume" + str(vol))
0469 
0470     for entry in data["Surfaces"]["entries"]:
0471         if "layer" in entry:
0472             vol = entry["volume"]
0473 
0474             if index_to_names[vol - 1] in layer_data:
0475                 if "approach" in entry:
0476                     if (
0477                         entry["layer"]
0478                         in layer_data[index_to_names[vol - 1]]["layers_with_approach"]
0479                         and entry["approach"]
0480                         in layer_data[index_to_names[vol - 1]]["approaches"]
0481                     ):
0482                         entry["value"]["material"]["mapMaterial"] = True
0483                 else:
0484                     if (
0485                         entry["layer"]
0486                         in layer_data[index_to_names[vol - 1]]["layers_no_approach"]
0487                     ):
0488                         entry["value"]["material"]["mapMaterial"] = True
0489 
0490                 if entry["value"]["material"]["mapMaterial"]:
0491                     for val in entry["value"]["material"]["binUtility"]["binningdata"]:
0492                         if val["value"] == "binZ" or val["value"] == "binR":
0493                             val["bins"] = layer_data[index_to_names[vol - 1]][
0494                                 "activeBinningRorZ"
0495                             ]
0496                         else:
0497                             val["bins"] = layer_data[index_to_names[vol - 1]][
0498                                 "activeBinningPhi"
0499                             ]
0500                         if val["bins"] == 0:
0501                             print(
0502                                 "ERROR!!! Using binning value == 0! Check you input for",
0503                                 index_to_names[vol - 1],
0504                             )
0505                             return
0506 
0507             approach_index = "None"
0508             if "approach" in entry:
0509                 approach_index = entry["approach"]
0510 
0511             if dump_binning_for_material and entry["value"]["material"]["mapMaterial"]:
0512                 print(
0513                     "Volume: ",
0514                     entry["volume"],
0515                     index_to_names[vol - 1],
0516                     " - Layer: ",
0517                     entry["layer"],
0518                     " - Approach:",
0519                     approach_index,
0520                 )
0521                 for val in entry["value"]["material"]["binUtility"]["binningdata"]:
0522                     print("-->", val["value"], ": ", val["bins"])
0523 
0524         if "boundary" in entry:
0525             extends = []
0526             vol = entry["volume"]
0527 
0528             if (
0529                 index_to_names[vol - 1] in layer_data
0530                 and entry["boundary"]
0531                 in layer_data[index_to_names[vol - 1]]["boundaries"]
0532             ):
0533                 entry["value"]["material"]["mapMaterial"] = True
0534                 for val in entry["value"]["material"]["binUtility"]["binningdata"]:
0535                     if entry["value"]["type"] == "CylinderSurface":
0536                         if val["value"] == "binZ":
0537                             val["bins"] = layer_data[index_to_names[vol - 1]][
0538                                 "passiveCylinderBinningZ"
0539                             ]
0540                         else:
0541                             val["bins"] = layer_data[index_to_names[vol - 1]][
0542                                 "passiveCylinderBinningPhi"
0543                             ]
0544                     elif entry["value"]["type"] == "DiscSurface":
0545                         if val["value"] == "binR":
0546                             val["bins"] = layer_data[index_to_names[vol - 1]][
0547                                 "passiveDiscBinningR"
0548                             ]
0549                         else:
0550                             val["bins"] = layer_data[index_to_names[vol - 1]][
0551                                 "passiveDiscBinningPhi"
0552                             ]
0553                     else:
0554                         surface_type = entry["value"]["type"]
0555                         print(
0556                             f"WARNING: Processing surface with unknown type '{surface_type}'. Only CylinderSurface and DiscSurface are considered."
0557                         )
0558                     if val["bins"] == 0:
0559                         print(
0560                             "ERROR!!! Using binning value == 0! Check you input for",
0561                             index_to_names[vol - 1],
0562                         )
0563                         return
0564 
0565             if dump_binning_for_material and entry["value"]["material"]["mapMaterial"]:
0566                 print(
0567                     "Volume: ",
0568                     entry["volume"],
0569                     index_to_names[vol - 1],
0570                     " - Boundary:",
0571                     entry["boundary"],
0572                 )
0573                 for val in entry["value"]["material"]["binUtility"]["binningdata"]:
0574                     print("-->", val["value"], ": ", val["bins"])
0575 
0576     # Once you have all the data in hands, you make a few nice plots to show volumes and surfaces
0577 
0578     with open(output_file, "w", encoding="utf-8") as f:
0579         json.dump(data, f, ensure_ascii=False, indent=4)
0580 
0581     if plot and check_material_layers:
0582         import matplotlib.pyplot as plt
0583         from matplotlib.pyplot import cm
0584         import numpy as np
0585 
0586         plt.figure(figsize=(20, 10))
0587 
0588         material_layer_cylinders = [[] for _ in range(len(index_to_names))]
0589         material_layer_discs = [[] for _ in range(len(index_to_names))]
0590 
0591         material_approach_cylinders = [[] for _ in range(len(index_to_names))]
0592         material_approach_discs = [[] for _ in range(len(index_to_names))]
0593 
0594         material_boundary_cylinders = [[] for _ in range(len(index_to_names))]
0595         material_boundary_discs = [[] for _ in range(len(index_to_names))]
0596 
0597         for entry in data["Surfaces"]["entries"]:
0598             if not entry["value"]["material"]["mapMaterial"]:
0599                 continue
0600 
0601             z_shift = 0.0
0602             if entry["value"]["transform"]["translation"] != None:
0603                 z_shift = entry["value"]["transform"]["translation"][2]
0604 
0605             if "layer" in entry:
0606                 extends = []
0607                 vol = entry["volume"]
0608 
0609                 if entry["value"]["type"] == "CylinderSurface":
0610                     extends = [
0611                         entry["value"]["bounds"]["values"][0],
0612                         z_shift - entry["value"]["bounds"]["values"][1],
0613                         z_shift + entry["value"]["bounds"]["values"][1],
0614                     ]
0615                     if "approach" in entry:
0616                         material_approach_cylinders[vol - 1].append(extends)
0617                     else:
0618                         material_layer_cylinders[vol - 1].append(extends)
0619 
0620                 elif entry["value"]["type"] == "DiscSurface":
0621                     extends = [
0622                         entry["value"]["bounds"]["values"][0],
0623                         entry["value"]["bounds"]["values"][1],
0624                         z_shift,
0625                     ]
0626                     if "approach" in entry:
0627                         material_approach_discs[vol - 1].append(extends)
0628                     else:
0629                         material_layer_discs[vol - 1].append(extends)
0630                 else:
0631                     surface_type = entry["value"]["type"]
0632                     print(
0633                         f"WARNING: Processing surface with unknown type '{surface_type}'. Only CylinderSurface and DiscSurface are considered."
0634                     )
0635 
0636             if "boundary" in entry:
0637                 extends = []
0638                 vol = entry["volume"]
0639 
0640                 if entry["value"]["type"] == "CylinderSurface":
0641                     extends = [
0642                         entry["value"]["bounds"]["values"][0],
0643                         z_shift - entry["value"]["bounds"]["values"][1],
0644                         z_shift + entry["value"]["bounds"]["values"][1],
0645                     ]
0646                     material_boundary_cylinders[vol - 1].append(extends)
0647 
0648                 elif entry["value"]["type"] == "DiscSurface":
0649                     extends = [
0650                         entry["value"]["bounds"]["values"][0],
0651                         entry["value"]["bounds"]["values"][1],
0652                         z_shift,
0653                     ]
0654                     material_boundary_discs[vol - 1].append(extends)
0655                 else:
0656                     print(
0657                         "WARNING: Processing surface with unknown type. Only CylinderSurface and DiscSurface are considered."
0658                     )
0659 
0660         from itertools import chain
0661 
0662         v_index = 0
0663         is_first = True
0664         is_disc = False
0665         for elements in chain(material_layer_cylinders, material_layer_discs):
0666             l_index = 0
0667             for coords in elements:
0668                 x_values, y_values = extract_coords(coords, is_disc)
0669                 if is_first:
0670                     plt.plot(x_values, y_values, c="black", label="layer")
0671                     is_first = False
0672                 else:
0673                     plt.plot(x_values, y_values, c="black")
0674                 l_index = l_index + 1
0675             v_index = v_index + 1
0676             if not is_disc and v_index == len(material_layer_cylinders):
0677                 is_disc = True
0678                 v_index = 0
0679 
0680         v_index = 0
0681         is_first = True
0682         is_disc = False
0683         for elements in chain(material_approach_cylinders, material_approach_discs):
0684             l_index = 0
0685             for coords in elements:
0686                 x_values, y_values = extract_coords(coords, is_disc)
0687                 if is_first:
0688                     plt.plot(x_values, y_values, c="red", label="approach")
0689                     is_first = False
0690                 else:
0691                     plt.plot(x_values, y_values, c="red")
0692                 l_index = l_index + 1
0693             v_index = v_index + 1
0694             if not is_disc and v_index == len(material_approach_cylinders):
0695                 is_disc = True
0696                 v_index = 0
0697 
0698         v_index = 0
0699         is_first = True
0700         is_disc = False
0701         for elements in chain(material_boundary_cylinders, material_boundary_discs):
0702             l_index = 0
0703             for coords in elements:
0704                 x_values, y_values = extract_coords(coords, is_disc)
0705                 if is_first:
0706                     plt.plot(x_values, y_values, c="blue", label="boundary")
0707                     is_first = False
0708                 else:
0709                     plt.plot(x_values, y_values, c="blue")
0710                 l_index = l_index + 1
0711             v_index = v_index + 1
0712             if not is_disc and v_index == len(material_boundary_cylinders):
0713                 is_disc = True
0714                 v_index = 0
0715 
0716         plt.xlabel("z [mm]")
0717         plt.ylabel("R [mm]")
0718         plt.title("Layers with material")
0719         plt.legend(loc="center left", bbox_to_anchor=(1, 0.5))
0720         plt.savefig(output_folder + "/material_layers.png")
0721 
0722 
0723 import argparse
0724 import os
0725 
0726 # Initialize parser
0727 parser = argparse.ArgumentParser()
0728 parser.add_argument(
0729     "--geometry",
0730     default="",
0731     type=str,
0732     help="Specify the input file for geometry visualisation",
0733 )
0734 parser.add_argument(
0735     "--plot",
0736     default=True,
0737     action="store_true",
0738     help="Enable plot creation for geometry visualisation (Default : True)",
0739 )
0740 parser.add_argument(
0741     "--output_folder",
0742     default="plots",
0743     type=str,
0744     help="Specify the output folder for plots (Default : plots)",
0745 )
0746 parser.add_argument(
0747     "--dump_steering",
0748     default=False,
0749     action="store_true",
0750     help="Enable production of steering file for material mapping (Default : False)",
0751 )
0752 parser.add_argument(
0753     "--edit",
0754     default=False,
0755     action="store_true",
0756     help="Enable editing of input file for creation of json for material mapping (Default : False)",
0757 )
0758 parser.add_argument(
0759     "--steering_file",
0760     default="",
0761     type=str,
0762     help="Specify the steering file to guide editing of json for material mapping",
0763 )
0764 parser.add_argument(
0765     "--output_map",
0766     default="",
0767     type=str,
0768     help="Specify the output json for material mapping",
0769 )
0770 args = parser.parse_args()
0771 
0772 print(" --- Geometry visualisation and material map fie creation --- ")
0773 print()
0774 
0775 if not args.geometry:
0776     print("Error: Missing input geometry file. Please specify --geometry.")
0777     exit()
0778 
0779 if not os.path.exists(args.geometry):
0780     print("Error: Invalid file path/name in --geometry. Please check your input!")
0781     exit()
0782 
0783 if not args.geometry.endswith(".json"):
0784     print("Error: Invalid file format in --geometry. Please check your input!")
0785     exit()
0786 
0787 print("\t parsing file : ", args.geometry)
0788 if args.plot:
0789     print("\t job configured to produce plots in output folder: ", args.output_folder)
0790     if not os.path.exists(args.output_folder):
0791         os.mkdir(args.output_folder)
0792 
0793 if args.dump_steering and args.edit:
0794     print("Error: Wrong job configuration. --dump_steering and --edit can't be \
0795         both true at the same time.")
0796     print(
0797         "\t Decide if you want to dump the steering file OR to read an existing file for editing the geometry file."
0798     )
0799     exit()
0800 
0801 if args.dump_steering:
0802     if not args.steering_file:
0803         print("Error: Missing output steering file. Please specify --steering_file.")
0804         exit()
0805     if not args.steering_file.endswith(".json"):
0806         print("Error: Invalid file format in --steering_file. It must end with .json!")
0807         exit()
0808     print(
0809         "\t job configured to produce steering file for material mapping with name: ",
0810         args.steering_file,
0811     )
0812 
0813 if args.edit:
0814     print("\t job configured to edit the input geometry file following a steering file")
0815     if not args.steering_file:
0816         print("Error: Missing input steering file. Please specify --steering_file.")
0817         exit()
0818     if not os.path.exists(args.steering_file):
0819         print(
0820             "Error: Invalid file path/name in --steering_file. Please check your input!"
0821         )
0822         exit()
0823     if not args.steering_file.endswith(".json"):
0824         print("Error: Invalid file format in --steering_file. Please check your input!")
0825         exit()
0826     if not args.output_map:
0827         print("Error: Missing output map file. Please specify --output_map.")
0828         exit()
0829     if not args.output_map.endswith(".json"):
0830         print("Error: Invalid file format in --output_map. Please check your input!")
0831         exit()
0832 
0833 
0834 if args.plot or args.dump_steering:
0835     dump_geo(
0836         args.geometry,
0837         args.plot,
0838         args.output_folder,
0839         args.dump_steering,
0840         args.steering_file,
0841     )
0842 
0843 if args.edit:
0844     read_and_modify(
0845         args.geometry,
0846         args.plot,
0847         args.output_folder,
0848         args.steering_file,
0849         args.output_map,
0850     )