fix: Skipping un-physical steps in material mapping (#1882) (#1968)

This PR fixes #1882.

The reason why the outer boundary was not seen is due to the accumulation of invalid material on it. Skipping un-physical steps (step length = 0) restores a sane behaviour of the mapping procedure.

I have also added a check on existence of the `event_id` branch in the `TChain`. If the branch doesn't exist, the batch size is evaluated to be 1.

@Corentin-Allaire @goetzgaycken

Examples/Io/Root/src/RootMaterialTrackReader.cpp

@@ -15,16 +15,33 @@
 #include <TChain.h>
 #include <TFile.h>
 #include <TMath.h>
+#include <TTree.h>
 
 ActsExamples::RootMaterialTrackReader::RootMaterialTrackReader(
     const Config& config, Acts::Logging::Level level)
     : ActsExamples::IReader(),
       m_logger{Acts::getDefaultLogger(name(), level)},
       m_cfg(config) {
+  if (m_cfg.fileList.empty()) {
+    throw std::invalid_argument{"No input files given"};
+  }
+
   m_inputChain = new TChain(m_cfg.treeName.c_str());
 
-  // Set the branches
-  m_inputChain->SetBranchAddress("event_id", &m_eventId);
+  // loop over the input files
+  for (const auto& inputFile : m_cfg.fileList) {
+    // add file to the input chain
+    m_inputChain->Add(inputFile.c_str());
+    ACTS_DEBUG("Adding File " << inputFile << " to tree '" << m_cfg.treeName
+                              << "'.");
+  }
+
+  // get the number of entries, which also loads the tree
+  size_t nentries = m_inputChain->GetEntries();
+
+  bool eventIdPresent =
+      (TTree::kMatch == m_inputChain->SetBranchAddress("event_id", &m_eventId));
+
   m_inputChain->SetBranchAddress("v_x", &m_v_x);
   m_inputChain->SetBranchAddress("v_y", &m_v_y);
   m_inputChain->SetBranchAddress("v_z", &m_v_z);
@@ -54,20 +71,10 @@ ActsExamples::RootMaterialTrackReader::RootMaterialTrackReader(
     m_inputChain->SetBranchAddress("sur_z", &m_sur_z);
     m_inputChain->SetBranchAddress("sur_pathCorrection", &m_sur_pathCorrection);
   }
-  if (m_cfg.fileList.empty()) {
-    throw std::invalid_argument{"No input files given"};
-  }
 
-  // loop over the input files
-  for (const auto& inputFile : m_cfg.fileList) {
-    // add file to the input chain
-    m_inputChain->Add(inputFile.c_str());
-    ACTS_DEBUG("Adding File " << inputFile << " to tree '" << m_cfg.treeName
-                              << "'.");
-  }
-
-  m_events = static_cast<size_t>(m_inputChain->GetMaximum("event_id") + 1);
-  size_t nentries = m_inputChain->GetEntries();
+  m_events = eventIdPresent
+                 ? static_cast<size_t>(m_inputChain->GetMaximum("event_id") + 1)
+                 : nentries;
   m_batchSize = nentries / m_events;
   ACTS_DEBUG("The full chain has "
              << nentries << " entries for " << m_events
@@ -78,6 +85,11 @@ ActsExamples::RootMaterialTrackReader::RootMaterialTrackReader(
 
   // If the events are not in order, get the entry numbers for ordered events
   if (not m_cfg.orderedEvents) {
+    if (not eventIdPresent) {
+      throw std::invalid_argument{
+          "'event_id' branch is missing in your tree. This is not compatible "
+          "with unordered events."};
+    }
     m_entryNumbers.resize(nentries);
     m_inputChain->Draw("event_id", "", "goff");
     // Sort to get the entry numbers of the ordered events
@@ -147,6 +159,9 @@ ActsExamples::ProcessCode ActsExamples::RootMaterialTrackReader::read(
       rmTrack.first.first = Acts::Vector3(m_v_x, m_v_y, m_v_z);
       rmTrack.first.second = Acts::Vector3(m_v_px, m_v_py, m_v_pz);
 
+      ACTS_VERBOSE("Track vertex:  " << rmTrack.first.first);
+      ACTS_VERBOSE("Track momentum:" << rmTrack.first.second);
+
       // Fill the individual steps
       size_t msteps = m_step_length->size();
       ACTS_VERBOSE("Reading " << msteps << " material steps.");
@@ -155,21 +170,40 @@ ActsExamples::ProcessCode ActsExamples::RootMaterialTrackReader::read(
       rmTrack.second.materialInL0 = 0.;
 
       for (size_t is = 0; is < msteps; ++is) {
+        ACTS_VERBOSE("====================");
+        ACTS_VERBOSE("[" << is + 1 << "/" << msteps << "] STEP INFORMATION: ");
+
+        double s = (*m_step_length)[is];
+        if (s == 0) {
+          ACTS_VERBOSE("invalid step length... skipping!");
+          continue;
+        }
+
         double mX0 = (*m_step_X0)[is];
         double mL0 = (*m_step_L0)[is];
-        double s = (*m_step_length)[is];
+
         rmTrack.second.materialInX0 += s / mX0;
         rmTrack.second.materialInL0 += s / mL0;
         /// Fill the position & the material
         Acts::MaterialInteraction mInteraction;
         mInteraction.position =
             Acts::Vector3((*m_step_x)[is], (*m_step_y)[is], (*m_step_z)[is]);
+        ACTS_VERBOSE("POSITION : " << (*m_step_x)[is] << ", " << (*m_step_y)[is]
+                                   << ", " << (*m_step_z)[is]);
         mInteraction.direction =
             Acts::Vector3((*m_step_dx)[is], (*m_step_dy)[is], (*m_step_dz)[is]);
+        ACTS_VERBOSE("DIRECTION: " << (*m_step_dx)[is] << ", "
+                                   << (*m_step_dy)[is] << ", "
+                                   << (*m_step_dz)[is]);
         mInteraction.materialSlab = Acts::MaterialSlab(
             Acts::Material::fromMassDensity(mX0, mL0, (*m_step_A)[is],
                                             (*m_step_Z)[is], (*m_step_rho)[is]),
             s);
+        ACTS_VERBOSE("MATERIAL: " << mX0 << ", " << mL0 << ", "
+                                  << (*m_step_A)[is] << ", " << (*m_step_Z)[is]
+                                  << ", " << (*m_step_rho)[is]);
+        ACTS_VERBOSE("====================");
+
         if (m_cfg.readCachedSurfaceInformation) {
           // add the surface information to the interaction this allows the
           // mapping to be speed up