Scripts/MaterialMapping/materialPlotHelper.cpp

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 #include "materialPlotHelper.hpp"
0010
0011 #include <iomanip>
0012 #include <ostream>
0013 #include <string>
0014 #include <unordered_map>
0015
0016 #include <nlohmann/json.hpp>
0017 using json = nlohmann::json;
0018
0019 /// Information on a given surface.
0020
0021 struct sinfo {
0022   std::string name;
0023   std::string idname;
0024   std::string id;
0025   int type;
0026   float pos;
0027   float range_min;
0028   float range_max;
0029 };
0030
0031 std::ostream& Acts::operator<<(std::ostream& os, Acts::GeometryIdentifier id) {
0032   os << "[ " << std::setw(3) << id.volume();
0033   os << " | " << std::setw(3) << id.boundary();
0034   os << " | " << std::setw(3) << id.layer();
0035   os << " | " << std::setw(3) << id.approach();
0036   os << " | " << std::setw(4) << id.sensitive() << " ]";
0037   return os;
0038 }
0039
0040 std::unordered_map<std::uint64_t, json> load_geometry_file(
0041     std::string geometry_file) {
0042   json geom;
0043   {
0044     std::ifstream gj(geometry_file);
0045     if (!gj.good()) {
0046       std::cerr << "WARNING: " << geometry_file << " not found." << std::endl;
0047     } else {
0048       try {
0049         gj >> geom;
0050       } catch (...) {
0051         std::cerr << "WARNING: Failed to parse " << geometry_file << "."
0052                   << std::endl;
0053       }
0054     }
0055   }
0056   std::unordered_map<std::uint64_t, json> surface_bounds;
0057   const auto& entries = geom["Surfaces"]["entries"];
0058   for (const auto& entry : entries) {
0059     // Handle both old (std::uint64_t) and new (object with component fields)
0060     // formats
0061     std::uint64_t gid;
0062     if (entry["value"]["geo_id"].is_number()) {
0063       // Specification: geo_id is a simple std::uint64_t
0064       gid = entry["value"]["geo_id"].get<std::uint64_t>();
0065     } else if (entry["value"]["geo_id"].is_object()) {
0066       // Specification: geo_id is an object with component fields
0067       // Reconstruct the std::uint64_t from individual components
0068       // Layout: volume (bits 56-63), boundary (bits 48-55), layer (bits 36-47),
0069       //         approach (bits 28-35), sensitive (bits 0-27)
0070       std::uint64_t gid_value = 0;
0071       const auto& geo_id_obj = entry["value"]["geo_id"];
0072
0073       if (geo_id_obj.contains("volume") && !geo_id_obj["volume"].is_null()) {
0074         std::uint32_t volume = static_cast<std::uint32_t>(
0075             geo_id_obj["volume"].get<std::uint32_t>());
0076         gid_value |= (static_cast<std::uint64_t>(volume) & 0xFF) << 56;
0077       }
0078       if (geo_id_obj.contains("boundary") &&
0079           !geo_id_obj["boundary"].is_null()) {
0080         std::uint32_t boundary = static_cast<std::uint32_t>(
0081             geo_id_obj["boundary"].get<std::uint32_t>());
0082         gid_value |= (static_cast<std::uint64_t>(boundary) & 0xFF) << 48;
0083       }
0084       if (geo_id_obj.contains("layer") && !geo_id_obj["layer"].is_null()) {
0085         std::uint32_t layer = static_cast<std::uint32_t>(
0086             geo_id_obj["layer"].get<std::uint32_t>());
0087         gid_value |= (static_cast<std::uint64_t>(layer) & 0xFFF) << 36;
0088       }
0089       if (geo_id_obj.contains("approach") &&
0090           !geo_id_obj["approach"].is_null()) {
0091         std::uint32_t approach = static_cast<std::uint32_t>(
0092             geo_id_obj["approach"].get<std::uint32_t>());
0093         gid_value |= (static_cast<std::uint64_t>(approach) & 0xFF) << 28;
0094       }
0095       if (geo_id_obj.contains("sensitive") &&
0096           !geo_id_obj["sensitive"].is_null()) {
0097         std::uint32_t sensitive = static_cast<std::uint32_t>(
0098             geo_id_obj["sensitive"].get<std::uint32_t>());
0099         gid_value |= (static_cast<std::uint64_t>(sensitive) & 0xFFFFF) << 8;
0100       }
0101       if (geo_id_obj.contains("extra") && !geo_id_obj["extra"].is_null()) {
0102         std::uint64_t extra = geo_id_obj["extra"].get<std::uint64_t>();
0103         gid_value |= (extra & 0xFF);
0104       }
0105       gid = gid_value;
0106     } else {
0107       std::cerr << "WARNING: Unknown geo_id format in geometry file."
0108                 << std::endl;
0109       continue;
0110     }
0111     surface_bounds[gid] = entry["value"]["bounds"];
0112   }
0113   return surface_bounds;
0114 }
0115
0116 /// Initialise the information on each surface.
0117
0118 void Initialise_info(sinfo& surface_info,
0119                      const std::map<std::string, std::string>& surfaceName,
0120                      const std::uint64_t& id, const int& type, const float& pos,
0121                      const float& range_min, const float& range_max) {
0122   Acts::GeometryIdentifier ID(id);
0123   std::ostringstream layerID;
0124   layerID << ID;
0125   std::string surface_id = layerID.str();
0126
0127   std::string Id_temp = surface_id;
0128   std::string delimiter = " | ";
0129   std::size_t del_pos = 0;
0130   std::vector<std::string> Ids;
0131   while ((del_pos = Id_temp.find(delimiter)) != std::string::npos) {
0132     Ids.push_back(Id_temp.substr(0, del_pos));
0133     Id_temp.erase(0, del_pos + delimiter.length());
0134   }
0135   Ids.push_back(Id_temp);
0136
0137   for (int tag = 0; tag < 5; tag++) {
0138     Ids[tag].erase(std::remove(Ids[tag].begin(), Ids[tag].end(), ' '),
0139                    Ids[tag].end());
0140     Ids[tag].erase(std::remove(Ids[tag].begin(), Ids[tag].end(), '['),
0141                    Ids[tag].end());
0142     Ids[tag].erase(std::remove(Ids[tag].begin(), Ids[tag].end(), ']'),
0143                    Ids[tag].end());
0144   }
0145
0146   surface_info.idname = "v" + Ids[0] + "_b" + Ids[1] + "_l" + Ids[2] + "_a" +
0147                         Ids[3] + "_s" + Ids[4];
0148   surface_info.type = type;
0149
0150   if (surfaceName.contains(surface_id)) {
0151     surface_info.name = surfaceName.at(surface_id);
0152   } else {
0153     surface_info.name = "";
0154   }
0155
0156   surface_info.id = surface_id;
0157   surface_info.pos = pos;
0158   surface_info.range_min = range_min;
0159   surface_info.range_max = range_max;
0160 }