Scripts/MaterialMapping/MaterialComposition.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 "ActsExamples/Utilities/Options.hpp"
0010
0011 #include <algorithm>
0012 #include <cstddef>
0013 #include <exception>
0014 #include <filesystem>
0015 #include <fstream>
0016 #include <iostream>
0017 #include <sstream>
0018 #include <string>
0019 #include <vector>
0020
0021 #include <TApplication.h>
0022 #include <boost/program_options.hpp>
0023 #include <boost/version.hpp>
0024 #include <nlohmann/json.hpp>
0025
0026 #define BOOST_AVAILABLE 1
0027 #if BOOST_VERSION < 107200
0028 // Boost <=1.71 and lower do not have progress_display.hpp as a replacement yet
0029 #include <boost/progress.hpp>
0030
0031 using progress_display = boost::progress_display;
0032 #else
0033 // Boost >=1.72 can use this as a replacement
0034 #include <boost/timer/progress_display.hpp>
0035
0036 using progress_display = boost::timer::progress_display;
0037 #endif
0038
0039 #include "materialComposition.C"
0040
0041 using namespace boost::program_options;
0042 using VariableReals = ActsExamples::Options::VariableReals;
0043
0044 int main(int argc, char** argv) {
0045   std::cout << "*** Material Composition plotting " << std::endl;
0046
0047   try {
0048     options_description description("*** Usage:");
0049
0050     // Add the program options
0051     auto ao = description.add_options();
0052     ao("help,h", "Display this help message");
0053     ao("silent,s", bool_switch(), "Silent mode (without X-window/display).");
0054     ao("input,i", value<std::string>()->default_value(""),
0055        "Input ROOT file containing the input TTree.");
0056     ao("tree,t", value<std::string>()->default_value("material-tracks"),
0057        "Input TTree name.");
0058     ao("output,o", value<std::string>()->default_value(""),
0059        "Output ROOT file with histograms");
0060     ao("bins,b", value<unsigned int>()->default_value(60),
0061        "Number of bins in eta/phi");
0062     ao("eta,e", value<float>()->default_value(4.), "Eta range.");
0063     ao("sub-names", value<std::vector<std::string>>()->multitoken(),
0064        "Subdetector names.");
0065     ao("sub-rmin", value<VariableReals>(), "Minimal radial restrictions.");
0066     ao("sub-rmax", value<VariableReals>(), "Maximal radial restrictions.");
0067     ao("sub-zmin", value<VariableReals>(), "Minimal z radial restrictions");
0068     ao("sub-zmax", value<VariableReals>(), "Maximal z radial restrictions.");
0069     ao("config,c", value<std::string>(), "Configuration file (json).");
0070
0071     // Set up the variables map
0072     variables_map vm;
0073     store(command_line_parser(argc, argv).options(description).run(), vm);
0074     notify(vm);
0075
0076     if (vm.contains("help")) {
0077       std::cout << description;
0078     }
0079
0080     // Parse the parameters
0081     auto iFile = vm["input"].as<std::string>();
0082     auto iTree = vm["tree"].as<std::string>();
0083     auto oFile = vm["output"].as<std::string>();
0084
0085     // Bins & eta range
0086     unsigned int bins = vm["bins"].as<unsigned int>();
0087     float eta = vm["eta"].as<float>();
0088
0089     // Subdetector configurations
0090     std::vector<Region> dRegion = {};
0091
0092     if (vm.contains("config")) {
0093       std::filesystem::path config = vm["config"].as<std::string>();
0094       std::cout << "Reading region configuration from JSON: " << config
0095                 << std::endl;
0096
0097       if (!std::filesystem::exists(config)) {
0098         std::cerr << "Configuration file does not exist." << std::endl;
0099         return 1;
0100       }
0101
0102       std::ifstream ifs(config.string().c_str());
0103       nlohmann::ordered_json j = nlohmann::ordered_json::parse(ifs);
0104
0105       for (const auto& [key, regions] : j.items()) {
0106         dRegion.push_back(Region{key, {}});
0107         auto& reg = dRegion.back();
0108         std::cout << "Region(" << key << ")" << std::endl;
0109         for (const auto& region : regions) {
0110           float rmin = region["rmin"].template get<float>();
0111           float rmax = region["rmax"].template get<float>();
0112           float zmin = region["zmin"].template get<float>();
0113           float zmax = region["zmax"].template get<float>();
0114
0115           reg.boxes.push_back({rmin, rmax, zmin, zmax});
0116           std::cout << "* " << key << " r/z: " << rmin << "/" << rmax << " "
0117                     << zmin << "/" << zmax << std::endl;
0118         }
0119       }
0120     } else {
0121       auto snames = vm["sub-names"].as<std::vector<std::string>>();
0122       auto rmins = vm["sub-rmin"].as<VariableReals>().values;
0123       auto rmaxs = vm["sub-rmax"].as<VariableReals>().values;
0124       auto zmins = vm["sub-zmin"].as<VariableReals>().values;
0125       auto zmaxs = vm["sub-zmax"].as<VariableReals>().values;
0126
0127       std::size_t subs = snames.size();
0128
0129       if (subs != rmins.size() || subs != rmaxs.size() ||
0130           subs != zmins.size() || subs != zmaxs.size()) {
0131         std::cerr << "Configuration problem." << std::endl;
0132         return 1;
0133       }
0134
0135       // Create the regions
0136       for (unsigned int is = 0; is < subs; ++is) {
0137         dRegion.push_back(Region{
0138             snames[is],
0139             {{static_cast<float>(rmins[is]), static_cast<float>(rmaxs[is]),
0140               static_cast<float>(zmins[is]), static_cast<float>(zmaxs[is])}}});
0141       }
0142     }
0143
0144     TApplication* tApp =
0145         vm["silent"].as<bool>()
0146             ? nullptr
0147             : new TApplication("ResidualAndPulls", nullptr, nullptr);
0148
0149     materialComposition(iFile, iTree, oFile, bins, eta, dRegion);
0150
0151     if (tApp != nullptr) {
0152       tApp->Run();
0153     }
0154
0155   } catch (std::exception& e) {
0156     std::cerr << e.what() << "\n";
0157   }
0158
0159   std::cout << "*** Done." << std::endl;
0160   return 0;
0161 }