InteractionList.hpp

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

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

#pragma once

#include "Acts/Material/MaterialSlab.hpp"
#include "Acts/Utilities/TypeTraits.hpp"
#include "ActsFatras/EventData/Particle.hpp"

#include <bitset>
#include <tuple>
#include <type_traits>
#include <utility>

namespace ActsFatras {
namespace detail {

template <class T, class Tuple>
struct TupleIndexOf;
template <class T, class... Types>
struct TupleIndexOf<T, std::tuple<T, Types...>> {
  static constexpr std::size_t value = 0u;
};
template <class T, class U, class... Types>
struct TupleIndexOf<T, std::tuple<U, Types...>> {
  static constexpr std::size_t value =
      1u + TupleIndexOf<T, std::tuple<Types...>>::value;
};

// Construct an index sequence for a subset of the tuple elements.
//
// Whether an element is part of the subset is defined by the predicate
// template type. It must take the element type as its only template parameter
// and must provide a static `value` member value. If the value evaluates to
// `true`, then the corresponding index will be part of the index sequence.
//
// Example: The tuple contains four elements, where all but the third one (i=2)
// should be selected. This leads to the following recursive expansion
// where the index sequence of the subset is filled from the front.
//
//        TupleFilterImpl<..., kCounter=4>          // select index=3
//     -> TupleFilterImpl<..., kCounter=3, 3>       // skip   index=2
//     -> TupleFilterImpl<..., kCounter=2, 3>       // select index=1
//     -> TupleFilterImpl<..., kCounter=1, 1, 3>    // select index=0
//     -> TupleFilterImpl<..., kCounter=0, 0, 1, 3> // terminate
//
template <template <typename> typename predicate_t, typename tuple_t,
          size_t kCounter, size_t... kIndices>
struct TupleFilterImpl {
  static constexpr auto kIndex = kCounter - 1u;
  static constexpr bool kElementSelection =
      predicate_t<std::tuple_element_t<kIndex, tuple_t>>::value;
  // recursive type if the element would be selected
  using SelectElement = typename TupleFilterImpl<predicate_t, tuple_t, kIndex,
                                                 kIndex, kIndices...>::Type;
  // recursive type if the element would be skipped
  using SkipElement =
      typename TupleFilterImpl<predicate_t, tuple_t, kIndex, kIndices...>::Type;
  // select recursive type based on the selector decision
  using Type =
      std::conditional_t<kElementSelection, SelectElement, SkipElement>;
};
template <template <typename> typename predicate_t, typename tuple_t,
          size_t... kIndices>
struct TupleFilterImpl<predicate_t, tuple_t, 0u, kIndices...> {
  using Type = std::index_sequence<kIndices...>;
};
template <template <typename> typename predicate_t, typename tuple_t>
using TupleFilter = typename TupleFilterImpl<predicate_t, tuple_t,
                                             std::tuple_size_v<tuple_t>>::Type;

template <typename process_t>
class IsPointLikeProcess {
  struct MockUniformRandomBitGenerator {
    using result_type = unsigned int;

    static constexpr result_type min() { return 0u; }
    static constexpr result_type max() { return 1u << 15u; }
    constexpr result_type operator()() { return 0u; }
  };

  METHOD_TRAIT(generatePathLimits_method_t, generatePathLimits);

  using scalar_pair_t = std::pair<Particle::Scalar, Particle::Scalar>;

 public:
  static constexpr bool value = Acts::Concepts::has_method<
      const process_t, scalar_pair_t, generatePathLimits_method_t,
      MockUniformRandomBitGenerator&, const Particle&>;
};

template <typename process_t>
struct IsContinuousProcess {
  static constexpr bool value = not IsPointLikeProcess<process_t>::value;
};

template <typename processes_t>
using ContinuousIndices = TupleFilter<IsContinuousProcess, processes_t>;
template <typename processes_t>
using PointLikeIndices = TupleFilter<IsPointLikeProcess, processes_t>;

}  // namespace detail

template <typename... processes_t>
class InteractionList {
  using Mask = std::bitset<sizeof...(processes_t)>;
  using Processes = std::tuple<processes_t...>;
  using ContinuousIndices = detail::ContinuousIndices<Processes>;
  using PointLikeIndices = detail::PointLikeIndices<Processes>;

 public:
  struct Selection {
    Particle::Scalar x0Limit =
        std::numeric_limits<Particle::Scalar>::infinity();
    Particle::Scalar l0Limit =
        std::numeric_limits<Particle::Scalar>::infinity();
    size_t x0Process = SIZE_MAX;
    size_t l0Process = SIZE_MAX;
  };

  void disable(size_t process) { m_mask.set(process); }
  template <typename process_t>
  void disable() {
    m_mask.set(detail::TupleIndexOf<process_t, Processes>::value);
  }

  template <size_t kProcess>
  std::tuple_element_t<kProcess, Processes>& get() {
    return std::get<kProcess>(m_processes);
  }
  template <typename process_t>
  process_t& get() {
    return std::get<process_t>(m_processes);
  }

  template <typename generator_t>
  bool runContinuous(generator_t& rng, const Acts::MaterialSlab& slab,
                     Particle& particle,
                     std::vector<Particle>& generated) const {
    return runContinuousImpl(rng, slab, particle, generated,
                             ContinuousIndices());
  }

  template <typename generator_t>
  Selection armPointLike(generator_t& rng, const Particle& particle) const {
    Selection selection;
    armPointLikeImpl(rng, particle, selection, PointLikeIndices());
    return selection;
  }

  template <typename generator_t>
  bool runPointLike(generator_t& rng, size_t processIndex, Particle& particle,
                    std::vector<Particle>& generated) const {
    return runPointLikeImpl(rng, processIndex, particle, generated,
                            PointLikeIndices());
  }

 private:
  // allow processes to be masked. defaults to zeros -> no masked processes
  Mask m_mask;
  Processes m_processes;

  // for the `runContinuous` call, we need to iterate over all available
  // processes and apply the ones that implement the continuous process
  // interface. this is done using an index-based compile-time recursive call.
  template <typename generator_t, std::size_t kI0, std::size_t... kIs>
  bool runContinuousImpl(generator_t& rng, const Acts::MaterialSlab& slab,
                         Particle& particle, std::vector<Particle>& generated,
                         std::index_sequence<kI0, kIs...> /*indices*/) const {
    const auto& process = std::get<kI0>(m_processes);
    // only call process if it is not masked
    if (not m_mask[kI0] and process(rng, slab, particle, generated)) {
      // exit early in case the process signals an abort
      return true;
    }
    return runContinuousImpl(rng, slab, particle, generated,
                             std::index_sequence<kIs...>());
  }
  template <typename generator_t>
  bool runContinuousImpl(generator_t& /*rng*/,
                         const Acts::MaterialSlab& /*slab*/,
                         Particle& /*particle*/,
                         std::vector<Particle>& /*generated*/,
                         std::index_sequence<> /*indices*/) const {
    return false;
  }

  // for the `armPointLike` call, we need to iterate over all available
  // processes and select the ones that generate the smallest limits. this is
  // done using an index-based compile-time recursive call.
  template <typename generator_t, std::size_t kI0, std::size_t... kIs>
  void armPointLikeImpl(generator_t& rng, const Particle& particle,
                        Selection& selection,
                        std::index_sequence<kI0, kIs...> /*indices*/) const {
    // only arm the process if it is not masked
    if (not m_mask[kI0]) {
      auto [x0Limit, l0Limit] =
          std::get<kI0>(m_processes).generatePathLimits(rng, particle);
      if (x0Limit < selection.x0Limit) {
        selection.x0Limit = x0Limit;
        selection.x0Process = kI0;
      }
      if (l0Limit < selection.l0Limit) {
        selection.l0Limit = l0Limit;
        selection.l0Process = kI0;
      }
    }
    // continue with the remaining processes
    armPointLikeImpl(rng, particle, selection, std::index_sequence<kIs...>());
  }
  template <typename generator_t>
  void armPointLikeImpl(generator_t& /*rng*/, const Particle& /*particle*/,
                        Selection& /*selection*/,
                        std::index_sequence<> /*indices*/) const {}

  // for the `runPointLike` call we need to call just one process. since we can
  // not select a tuple element with a runtime index, we need to iterate over
  // all processes with a compile-time recursive function until we reach the
  // requested one.
  template <typename generator_t, size_t kI0, size_t... kIs>
  bool runPointLikeImpl(generator_t& rng, size_t processIndex,
                        Particle& particle, std::vector<Particle>& generated,
                        std::index_sequence<kI0, kIs...> /*indices*/) const {
    if (kI0 == processIndex) {
      if (m_mask[kI0]) {
        // the selected process is masked. since nothing is executed the
        // particle continues to be alive; not a break condition.
        return false;
      }
      return std::get<kI0>(m_processes).run(rng, particle, generated);
    }
    // continue the iteration with the remaining processes
    return runPointLikeImpl(rng, processIndex, particle, generated,
                            std::index_sequence<kIs...>());
  }
  template <typename generator_t>
  bool runPointLikeImpl(generator_t& /*rng*/, size_t /*processIndex*/,
                        Particle& /*particle*/,
                        std::vector<Particle>& /*generated*/,
                        std::index_sequence<> /*indices*/) const {
    // the requested process index is outside the possible range. **do not**
    // treat this as an error to simplify the case of an empty physics lists or
    // a default process index.
    return false;
  }
};

}  // namespace ActsFatras