feat: material composition analysis (#1080)

* added material composition analysis

* adding projection vs phi

* clang format fix

* fixing typo in material mapping

Examples/Scripts/CMakeLists.txt

@@ -1 +1,2 @@
+add_subdirectory(MaterialMapping)
 add_subdirectory(TrackingPerformance)

Examples/Scripts/MaterialMapping/CMakeLists.txt

@@ -0,0 +1,8 @@
+add_executable(ActsAnalysisMaterialComposition MaterialComposition.cpp)
+target_link_libraries(ActsAnalysisMaterialComposition ActsExamplesFramework ROOT::Core ROOT::Hist ROOT::Tree ROOT::TreePlayer Boost::program_options)
+
+install(
+  TARGETS
+  ActsAnalysisMaterialComposition
+  RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})
+

Examples/Scripts/MaterialMapping/MaterialComposition.cpp

@@ -0,0 +1,116 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2021 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "ActsExamples/Utilities/Options.hpp"
+
+#include <exception>
+#include <iostream>
+#include <string>
+
+#include <TApplication.h>
+#include <boost/program_options.hpp>
+
+#define BOOST_AVAILABLE 1
+#if ((BOOST_VERSION / 100) % 1000) <= 71
+// Boost <=1.71 and lower do not have progress_display.hpp as a replacement yet
+#include <boost/progress.hpp>
+using progress_display = boost::progress_display;
+#else
+// Boost >=1.72 can use this as a replacement
+#include <boost/timer/progress_display.hpp>
+using progress_display = boost::timer::progress_display;
+#endif
+
+#include "materialComposition.C"
+
+using namespace boost::program_options;
+using VariableReals = ActsExamples::Options::VariableReals;
+
+int main(int argc, char** argv) {
+  std::cout << "*** Material Composition plotting " << std::endl;
+
+  try {
+    options_description description("*** Usage:");
+
+    // Add the program options
+    auto ao = description.add_options();
+    ao("help,h", "Display this help message");
+    ao("silent,s", bool_switch(), "Silent mode (without X-window/display).");
+    ao("input,i", value<std::string>()->default_value(""),
+       "Input ROOT file containing the input TTree.");
+    ao("tree,t", value<std::string>()->default_value("trackstates"),
+       "Input TTree name.");
+    ao("output,o", value<std::string>()->default_value(""),
+       "Output ROOT file with histograms");
+    ao("bins,b", value<unsigned int>()->default_value(60),
+       "Number of bins in eta/phi");
+    ao("eta,e", value<float>()->default_value(4.), "Eta range.");
+    ao("sub-names", value<std::vector<std::string>>()->multitoken(),
+       "Subdetector names.");
+    ao("sub-rmin", value<VariableReals>(), "Minimal radial restrictions.");
+    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.");
+
+    // Set up the variables map
+    variables_map vm;
+    store(command_line_parser(argc, argv).options(description).run(), vm);
+    notify(vm);
+
+    if (vm.count("help")) {
+      std::cout << description;
+    }
+
+    // Parse the parameters
+    auto iFile = vm["input"].as<std::string>();
+    auto iTree = vm["tree"].as<std::string>();
+    auto oFile = vm["output"].as<std::string>();
+
+    // Bins & eta range
+    unsigned int bins = vm["bins"].as<unsigned int>();
+    float eta = vm["eta"].as<float>();
+
+    // 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]});
+    }
+
+    TApplication* tApp = vm["silent"].as<bool>()
+                             ? nullptr
+                             : new TApplication("ResidualAndPulls", 0, 0);
+
+    materialComposition(iFile, iTree, oFile, bins, eta, dRegion);
+
+    if (tApp != nullptr) {
+      tApp->Run();
+    }
+
+  } catch (std::exception& e) {
+    std::cerr << e.what() << "\n";
+  }
+
+  std::cout << "*** Done." << std::endl;
+  return 1;
+}

Examples/Scripts/MaterialMapping/materialComposition.C

@@ -0,0 +1,231 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2021 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include <iostream>
+#include <map>
+#include <string>
+#include <tuple>
+
+#include <TCanvas.h>
+#include <TDirectory.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TProfile.h>
+#include <TTree.h>
+
+struct MaterialHistograms {
+
+  TProfile* x0_vs_eta = nullptr;
+  TProfile* l0_vs_eta = nullptr;
+
+  TProfile* x0_vs_phi= nullptr;
+  TProfile* l0_vs_phi = nullptr;
+
+  float s_x0 = 0.;
+  float s_l0 = 0.;
+
+  /// Material histogram constructor
+  ///
+  /// @param iA the atomic number
+  MaterialHistograms(const std::string& name, unsigned int iA,
+                     unsigned int bins, float eta) {
+    std::string x0NameEta =
+        (iA == 0) ? name + std::string("_x0_vs_eta_all")
+                  : name + std::string("_x0_vs_eta_A") + std::to_string(iA);
+    std::string l0NameEta =
+        (iA == 0) ? name + std::string("_l0_vs_eta_all")
+                  : name + std::string("_l0_vs_eta_A") + std::to_string(iA);
+
+    x0_vs_eta =
+        new TProfile(x0NameEta.c_str(), "X_{0} vs. #eta", bins, -eta, eta, 0., 5.);
+    l0_vs_eta =
+        new TProfile(l0NameEta.c_str(), "L_{0} vs. #eta", bins, -eta, eta, 0., 5.);
+
+    std::string x0NamePhi =
+        (iA == 0) ? name + std::string("_x0_vs_phi_all")
+                  : name + std::string("_x0_vs_phi_A") + std::to_string(iA);
+    std::string l0NamePhi =
+        (iA == 0) ? name + std::string("_l0_vs_phi_all")
+                  : name + std::string("_l0_vs_phi_A") + std::to_string(iA);
+
+    x0_vs_phi =
+        new TProfile(x0NamePhi.c_str(), "X_{0} vs. #phi", bins, -M_PI, M_PI, 0., 5.);
+    l0_vs_phi =
+        new TProfile(l0NamePhi.c_str(), "L_{0} vs. #phi", bins, -M_PI, M_PI, 0., 5.);
+
+  }
+
+  MaterialHistograms() = default;
+
+  /// This fills the event into the histograms
+  /// and clears the cache accordingly
+  ///
+  /// @param eta the pseudorapidity value
+  /// @param phi the phi value
+  ///
+  void fillAndClear(float eta, float phi) {
+    x0_vs_eta->Fill(eta, s_x0);
+    l0_vs_eta->Fill(eta, s_l0);
+
+    x0_vs_phi->Fill(phi, s_x0);
+    l0_vs_phi->Fill(phi, s_l0);
+
+    s_x0 = 0.;
+    s_l0 = 0.;
+  }
+
+  /// Write out the histograms, the TDirectory needs
+  /// to be se before
+  void write() {
+    x0_vs_eta->Write();
+    l0_vs_eta->Write();
+
+    x0_vs_phi->Write();
+    l0_vs_phi->Write();
+  }
+};
+
+using Region = std::tuple<std::string, float, float, float, float>;
+
+/// Plot the material composition
+///
+/// @param inFile the input root file
+/// @param treeNAme the input tree name (default: 'trackstates)
+/// @param outFile the output root file
+/// @param bins the number of bins
+/// @param eta the eta range
+void materialComposition(const std::string& inFile, const std::string& treeName,
+                         const std::string& outFile, unsigned int bins,
+                         float eta, const std::vector<Region>& regions) {
+  // Open the input file & get the tree
+  auto inputFile = TFile::Open(inFile.c_str());
+  auto inputTree = dynamic_cast<TTree*>(inputFile->Get(treeName.c_str()));
+  if (inputTree != nullptr) {
+    // Get the different atomic numbers
+    TCanvas* materialCanvas =
+        new TCanvas("materialCanvas", "Materials", 100, 100, 620, 400);
+    // Draw all the atomic elements & get the histogram
+    inputTree->Draw("mat_A>>hA(100,0.5,100.5)");
+    TH1F* histA = dynamic_cast<TH1F*>(gDirectory->Get("hA"));
+    histA->Draw();
+
+    auto outputFile = TFile::Open(outFile.c_str(), "recreate");
+
+    float v_eta;
+    float v_phi;
+    std::vector<float>* stepLength = new std::vector<float>;
+    std::vector<float>* stepX0 = new std::vector<float>;
+    std::vector<float>* stepL0 = new std::vector<float>;
+    std::vector<float>* stepA = new std::vector<float>;
+
+    std::vector<float>* startX = new std::vector<float>;
+    std::vector<float>* startY = new std::vector<float>;
+    std::vector<float>* startZ = new std::vector<float>;
+
+    std::vector<float>* endX = new std::vector<float>;
+    std::vector<float>* endY = new std::vector<float>;
+    std::vector<float>* endZ = new std::vector<float>;
+
+    inputTree->SetBranchAddress("v_eta", &v_eta);
+    inputTree->SetBranchAddress("v_phi", &v_phi);
+    inputTree->SetBranchAddress("mat_step_length", &stepLength);
+
+    inputTree->SetBranchAddress("mat_X0", &stepX0);
+    inputTree->SetBranchAddress("mat_L0", &stepL0);
+    inputTree->SetBranchAddress("mat_A", &stepA);
+
+    inputTree->SetBranchAddress("mat_sx", &startX);
+    inputTree->SetBranchAddress("mat_sy", &startY);
+    inputTree->SetBranchAddress("mat_sz", &startZ);
+
+    inputTree->SetBranchAddress("mat_ex", &endX);
+    inputTree->SetBranchAddress("mat_ey", &endY);
+    inputTree->SetBranchAddress("mat_ez", &endZ);
+
+    // Loop over all entries ---------------
+    unsigned int entries = inputTree->GetEntries();
+
+#ifdef BOOST_AVAILABLE
+    std::cout << "*** Event Loop: " << std::endl;
+    progress_display event_loop_progress(entries * regions.size());
+#endif
+
+    // Loop of the regions
+    for (auto& r : regions) {
+      std::string rName = std::get<0>(r);
+
+      // The material histograms ordered by atomic mass
+      std::map<unsigned int, MaterialHistograms> mCache;
+
+      // The material histograms for all
+      mCache[0] = MaterialHistograms(rName, 0, bins, eta);
+      for (unsigned int ib = 1; ib <= 100; ++ib) {
+        if (histA->GetBinContent(ib) > 0.) {
+          mCache[ib] = MaterialHistograms(rName, ib, bins, eta);
+        }
+      }
+
+      for (unsigned int ie = 0; ie < entries; ++ie) {
+        // Get the entry
+        inputTree->GetEntry(ie);
+
+#ifdef BOOST_AVAILABLE
+        ++event_loop_progress;
+#endif
+
+        // Accumulate the material per track
+        size_t steps = stepLength->size();
+        for (unsigned int is = 0; is < steps; ++is) {
+          float sX = startX->at(is);
+          float sY = startY->at(is);
+          float sZ = startZ->at(is);
+          float sR = sqrt(sX * sX + sY * sY);
+
+          float eX = endX->at(is);
+          float eY = endY->at(is);
+          float eZ = endZ->at(is);
+          float eR = sqrt(eX * eX + eY * eY);
+
+          float minR = std::get<1>(r);
+          float maxR = std::get<2>(r);
+          float minZ = std::get<3>(r);
+          float maxZ = std::get<4>(r);
+
+          if (minR > sR or minZ > sZ or maxR < eR or maxZ < eZ) {
+            continue;
+          }
+
+          float step = stepLength->at(is);
+          float X0 = stepX0->at(is);
+          float L0 = stepL0->at(is);
+
+          // The integral one
+          auto& all = mCache[0];
+          all.s_x0 += step / X0;
+          all.s_l0 += step / L0;
+
+          unsigned int sA = histA->FindBin(stepA->at(is));
+          // The current one
+          auto& current = mCache[sA];
+          current.s_x0 += step / X0;
+          current.s_l0 += step / L0;
+        }
+        // Fill the profiles and clear the cache
+        for (auto& [key, cache] : mCache) {
+          cache.fillAndClear(v_eta, v_phi);
+        }
+      }
+      // -----------------------------------
+
+      for (auto [key, cache] : mCache) {
+        cache.write();
+      }
+    }
+    outputFile->Close();
+  }
+}