feat: MaterialComposition plots with python config (#2531)

Changes the material composition script to optionally accept a region config in the form of a JSON config file like

```json
{
  "beampipe": [
    {
      "rmin": 0,
      "rmax": 24.4,
      "zmin": -4000,
      "zmax": 4000
    }
  ],
  "pixel": [
    {
      "rmin": 25,
      "rmax": 200,
      "zmin": -2400,
      "zmax": 2400
    }
  ],
  "sstrips": [
    {
      "rmin": 202,
      "rmax": 720,
      "zmin": -3150,
      "zmax": 3150
    }
  ],
  "lstrips": [
    {
      "rmin": 720,
      "rmax": 1140,
      "zmin": -3150,
      "zmax": 3150
    }
  ],
  "solenoid": [
    {
      "rmin": 1160,
      "rmax": 1200,
      "zmin": -3000,
      "zmax": 3000
    }
  ],
  "ecalbarrel": [
    {
      "rmin": 1250,
      "rmax": 1500,
      "zmin": -3050,
      "zmax": 3050
    }
  ],
  "ecalendcap": [
    {
      "rmin": 315,
      "rmax": 1500,
      "zmin": 3200,
      "zmax": 3450
    },
    {
      "rmin": 315,
      "rmax": 1500,
      "zmin": -3450,
      "zmax": -3200
    }
  ]
}
```

for the example of the ODD. This then also allows to have multiple $rz$ regions make up a combined grouping.

Examples/Scripts/MaterialMapping/CMakeLists.txt

@@ -1,5 +1,5 @@
 add_executable(ActsAnalysisMaterialComposition MaterialComposition.cpp)
-target_link_libraries(ActsAnalysisMaterialComposition ActsExamplesFramework ROOT::Core ROOT::Hist ROOT::Tree ROOT::TreePlayer Boost::program_options)
+target_link_libraries(ActsAnalysisMaterialComposition ActsExamplesFramework ROOT::Core ROOT::Hist ROOT::Tree ROOT::TreePlayer Boost::program_options nlohmann_json::nlohmann_json)
 
 install(
   TARGETS

Examples/Scripts/MaterialMapping/MaterialComposition.cpp

@@ -11,13 +11,16 @@
 #include <algorithm>
 #include <cstddef>
 #include <exception>
+#include <filesystem>
+#include <fstream>
 #include <iostream>
 #include <sstream>
 #include <string>
 #include <vector>
 
 #include <TApplication.h>
 #include <boost/program_options.hpp>
+#include <nlohmann/json.hpp>
 
 #define BOOST_AVAILABLE 1
 #if ((BOOST_VERSION / 100) % 1000) <= 71
@@ -62,6 +65,7 @@ int main(int argc, char** argv) {
     ao("sub-rmax", value<VariableReals>(), "Maximal radial restrictions.");
     ao("sub-zmin", value<VariableReals>(), "Minimal z radial restrictions");
     ao("sub-zmax", value<VariableReals>(), "Maximal z radial restrictions.");
+    ao("config,c", value<std::string>(), "Configuration file (json).");
 
     // Set up the variables map
     variables_map vm;
@@ -83,24 +87,57 @@ int main(int argc, char** argv) {
 
     // Subdetector configurations
     std::vector<Region> dRegion = {};
-    auto snames = vm["sub-names"].as<std::vector<std::string>>();
-    auto rmins = vm["sub-rmin"].as<VariableReals>().values;
-    auto rmaxs = vm["sub-rmax"].as<VariableReals>().values;
-    auto zmins = vm["sub-zmin"].as<VariableReals>().values;
-    auto zmaxs = vm["sub-zmax"].as<VariableReals>().values;
-
-    size_t subs = snames.size();
-
-    if (subs != rmins.size() or subs != rmaxs.size() or subs != zmins.size() or
-        subs != zmaxs.size()) {
-      std::cerr << "Configuration problem." << std::endl;
-      return 1;
-    }
 
-    // Create the regions
-    for (unsigned int is = 0; is < subs; ++is) {
-      dRegion.push_back(
-          {snames[is], rmins[is], rmaxs[is], zmins[is], zmaxs[is]});
+    if (vm.count("config") > 0) {
+      std::filesystem::path config = vm["config"].as<std::string>();
+      std::cout << "Reading region configuration from JSON: " << config
+                << std::endl;
+
+      if (!std::filesystem::exists(config)) {
+        std::cerr << "Configuration file does not exist." << std::endl;
+        return 1;
+      }
+
+      std::ifstream ifs(config.string().c_str());
+      nlohmann::ordered_json j = nlohmann::ordered_json::parse(ifs);
+
+      for (const auto& [key, regions] : j.items()) {
+        dRegion.push_back(Region{key, {}});
+        auto& reg = dRegion.back();
+        std::cout << "Region(" << key << ")" << std::endl;
+        for (const auto& region : regions) {
+          float rmin = region["rmin"].template get<float>();
+          float rmax = region["rmax"].template get<float>();
+          float zmin = region["zmin"].template get<float>();
+          float zmax = region["zmax"].template get<float>();
+
+          reg.boxes.push_back({rmin, rmax, zmin, zmax});
+          std::cout << "* " << key << " r/z: " << rmin << "/" << rmax << " "
+                    << zmin << "/" << zmax << std::endl;
+        }
+      }
+    } else {
+      auto snames = vm["sub-names"].as<std::vector<std::string>>();
+      auto rmins = vm["sub-rmin"].as<VariableReals>().values;
+      auto rmaxs = vm["sub-rmax"].as<VariableReals>().values;
+      auto zmins = vm["sub-zmin"].as<VariableReals>().values;
+      auto zmaxs = vm["sub-zmax"].as<VariableReals>().values;
+
+      size_t subs = snames.size();
+
+      if (subs != rmins.size() or subs != rmaxs.size() or
+          subs != zmins.size() or subs != zmaxs.size()) {
+        std::cerr << "Configuration problem." << std::endl;
+        return 1;
+      }
+
+      // Create the regions
+      for (unsigned int is = 0; is < subs; ++is) {
+        dRegion.push_back(Region{
+            snames[is],
+            {{static_cast<float>(rmins[is]), static_cast<float>(rmaxs[is]),
+              static_cast<float>(zmins[is]), static_cast<float>(zmaxs[is])}}});
+      }
     }
 
     TApplication* tApp =

Examples/Scripts/MaterialMapping/materialComposition.C

@@ -94,7 +94,19 @@ struct MaterialHistograms {
   }
 };
 
-using Region = std::tuple<std::string, float, float, float, float>;
+struct Region {
+  std::string name;
+  std::vector<std::tuple<float, float, float, float>> boxes;
+
+  bool inside(float r, float z) const {
+    for (const auto& [minR, maxR, minZ, maxZ] : boxes) {
+      if (minR <= r && r < maxR && minZ <= z && z < maxZ) {
+        return true;
+      }
+    }
+    return false;
+  }
+};
 
 /// Plot the material composition
 ///
@@ -157,7 +169,7 @@ void materialComposition(const std::string& inFile, const std::string& treeName,
 
   // Loop of the regions
   for (auto& region : regions) {
-    const auto [rName, minR, maxR, minZ, maxZ] = region;
+    const auto rName = region.name;
 
     // The material histograms ordered by atomic mass
     std::map<unsigned int, MaterialHistograms> mCache;
@@ -186,7 +198,7 @@ void materialComposition(const std::string& inFile, const std::string& treeName,
         float z = stepZ->at(is);
         float r = std::hypot(x, y);
 
-        if (minR > r or minZ > z or maxR < r or maxZ < z) {
+        if (!region.inside(r, z)) {
           continue;
         }
 
@@ -203,7 +215,8 @@ void materialComposition(const std::string& inFile, const std::string& treeName,
         // The current one
         auto currentIt = mCache.find(sA);
         if (currentIt == mCache.end()) {
-          throw std::runtime_error{"Unknown atomic number " +std::to_string(sA)};
+          throw std::runtime_error{"Unknown atomic number " +
+                                   std::to_string(sA)};
         }
         auto& current = currentIt->second;
         current.s_x0 += step / X0;