HitsPrinter.cpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HitsPrinter.cpp

// This file is part of the Acts project.
//
// Copyright (C) 2019 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 "HitsPrinter.hpp"

#include "Acts/Digitization/PlanarModuleCluster.hpp"
#include "Acts/Geometry/GeometryIdentifier.hpp"
#include "Acts/Utilities/Logger.hpp"
#include "Acts/Utilities/MultiIndex.hpp"
#include "ActsExamples/EventData/GeometryContainers.hpp"
#include "ActsExamples/Framework/AlgorithmContext.hpp"
#include "ActsExamples/Utilities/Range.hpp"

#include <algorithm>
#include <ostream>
#include <stdexcept>
#include <utility>
#include <vector>

ActsExamples::HitsPrinter::HitsPrinter(
    const ActsExamples::HitsPrinter::Config& cfg, Acts::Logging::Level level)
    : IAlgorithm("HitsPrinter", level), m_cfg(cfg) {
  if (m_cfg.inputClusters.empty()) {
    throw std::invalid_argument("Input clusters collection is not configured");
  }
  if (m_cfg.inputMeasurementParticlesMap.empty()) {
    throw std::invalid_argument(
        "Input hit-particles map collection is not configured");
  }
  if (m_cfg.inputHitIds.empty()) {
    throw std::invalid_argument("Input hit ids collection is not configured");
  }

  m_inputClusters.initialize(m_cfg.inputClusters);
  m_inputMeasurementParticlesMap.initialize(m_cfg.inputMeasurementParticlesMap);
  m_inputHitIds.initialize(m_cfg.inputHitIds);
}

ActsExamples::ProcessCode ActsExamples::HitsPrinter::execute(
    const ActsExamples::AlgorithmContext& ctx) const {
  const auto& clusters = m_inputClusters(ctx);
  const auto& hitParticlesMap = m_inputMeasurementParticlesMap(ctx);
  const auto& hitIds = m_inputHitIds(ctx);

  if (clusters.size() != hitIds.size()) {
    ACTS_ERROR(
        "event "
        << ctx.eventNumber
        << " input clusters and hit ids collections have inconsistent size");
    return ProcessCode::ABORT;
  }
  ACTS_INFO("event " << ctx.eventNumber << " collection '"
                     << m_cfg.inputClusters << "' contains " << clusters.size()
                     << " hits");

  // print hits selected by index
  if (0 < m_cfg.selectIndexLength) {
    size_t ihit = m_cfg.selectIndexStart;
    // Saturated addition that does not overflow and exceed SIZE_MAX.
    // From http://locklessinc.com/articles/sat_arithmetic/
    size_t nend = ihit + m_cfg.selectIndexLength;
    nend |= -static_cast<int>(nend < ihit);
    // restrict to available hits
    nend = std::min(clusters.size(), nend);

    if (nend <= ihit) {
      ACTS_WARNING("event "
                   << ctx.eventNumber << " collection '" << m_cfg.inputClusters
                   << " hit index selection is outside the available range");
    } else {
      ACTS_INFO("event " << ctx.eventNumber << " collection '"
                         << m_cfg.inputClusters << "' contains "
                         << (nend - ihit) << " hits in index selection ["
                         << m_cfg.selectIndexStart << ", "
                         << m_cfg.selectIndexLength << ")");

      Acts::GeometryIdentifier prevGeoId;
      for (; ihit < nend; ++ihit) {
        // ihit is already known to be within the available cluster range.
        // the `ic` iterator does not need to be checked for validity.
        auto ic = clusters.nth(ihit);
        Acts::GeometryIdentifier geoId = ic->first;
        const Acts::PlanarModuleCluster& c = ic->second;
        auto hitId = hitIds.at(ihit);

        if (geoId != prevGeoId) {
          ACTS_INFO("on geometry id " << geoId);
          prevGeoId = geoId;
        }
        ACTS_INFO("  cluster " << ihit << " hit " << hitId << " size "
                               << c.digitizationCells().size());
        // get all contributing particles
        for (auto [hid, pid] : makeRange(hitParticlesMap.equal_range(ihit))) {
          ACTS_INFO("    generated by particle " << pid);
        }
      }
    }
  }

  // print hits within geometry selection
  auto geoSelection = makeRange(clusters.begin(), clusters.begin());
  if (m_cfg.selectModule != 0u) {
    geoSelection = selectModule(clusters, m_cfg.selectVolume, m_cfg.selectLayer,
                                m_cfg.selectModule);
  } else if (m_cfg.selectLayer != 0u) {
    geoSelection = selectLayer(clusters, m_cfg.selectVolume, m_cfg.selectLayer);
  } else if (m_cfg.selectVolume != 0u) {
    geoSelection = selectVolume(clusters, m_cfg.selectVolume);
  }
  if (not geoSelection.empty()) {
    ACTS_INFO("event " << ctx.eventNumber << " collection '"
                       << m_cfg.inputClusters << "' contains "
                       << geoSelection.size() << " hits in volume "
                       << m_cfg.selectVolume << " layer " << m_cfg.selectLayer
                       << " module " << m_cfg.selectModule);
    // we could also use for (const Acts::PlanarModuleCluster& c : rangeModule)
    // for simplicity, but then we could not get the hit index.
    Acts::GeometryIdentifier prevGeoId;
    for (auto ic = geoSelection.begin(); ic != geoSelection.end(); ++ic) {
      auto ihit = clusters.index_of(ic);
      auto hitId = hitIds[ihit];

      Acts::GeometryIdentifier geoId = ic->first;
      const Acts::PlanarModuleCluster& c = ic->second;
      if (geoId != prevGeoId) {
        ACTS_INFO("on geometry id " << geoId);
        prevGeoId = geoId;
      }
      ACTS_INFO("  cluster " << ihit << " hit " << hitId << " size "
                             << c.digitizationCells().size());
      // get all contributing particles
      for (auto [hid, pid] : makeRange(hitParticlesMap.equal_range(ihit))) {
        ACTS_INFO("    generated by particle " << pid);
      }
    }
  }

  return ProcessCode::SUCCESS;
}