Scripts/MaterialMapping/MaterialComposition.cpp

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