material_validation_utils.hpp

// This file is part of the ACTS project.
//
// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
 
#pragma once
 
// Project include(s)
#include "detray/core/concepts.hpp"
#include "detray/definitions/algebra.hpp"
#include "detray/material/concepts.hpp"
#include "detray/material/detail/material_accessor.hpp"
#include "detray/navigation/caching_navigator.hpp"
#include "detray/navigation/detail/print_state.hpp"
#include "detray/propagator/actors.hpp"
#include "detray/propagator/line_stepper.hpp"
#include "detray/propagator/propagator.hpp"
#include "detray/utils/logging.hpp"
#include "detray/utils/type_list.hpp"
 
// Detray IO include(s)
#include "detray/io/utils/create_path.hpp"
#include "detray/io/utils/file_handle.hpp"
 
// Test include(s)
#include "detray/test/utils/data_record.hpp"
 
// Vecmem include(s)
#include <vecmem/memory/memory_resource.hpp>
 
// System include(s)
#include <filesystem>
 
namespace detray::material_validator {

struct track_material {
  double phi{detail::invalid_value<double>()};
  double eta{detail::invalid_value<double>()};
  double sX0{0.f};
  double tX0{0.f};
  double sL0{0.f};
  double tL0{0.f};
};

struct get_material_record {
  template <typename mat_group_t, typename index_t, concepts::point2D point2_t,
            concepts::scalar scalar_t>
  DETRAY_HOST_DEVICE auto operator()(
      [[maybe_unused]] const mat_group_t &mat_group,
      [[maybe_unused]] const index_t &index,
      [[maybe_unused]] const point2_t &loc,
      [[maybe_unused]] const scalar_t cos_inc_angle) const {
    using material_t = typename mat_group_t::value_type;
 
    constexpr auto inv{detail::invalid_value<double>()};
    constexpr geometry::identifier inv_sf{};
 
    // Access homogeneous surface material or material maps
    if constexpr (concepts::surface_material<material_t>) {
      // Slab or rod
      const auto mat = detail::material_accessor::get(mat_group, index, loc);
 
      // Empty material can occur in material maps, skip it
      if (!mat) {
        // Set the pathlength and thickness to zero so that they are not counted
        return material_record{inv_sf, 0.f, 0.f, inv, inv};
      }
 
      const double seg{mat.path_segment(cos_inc_angle, loc[0])};
      const double t{mat.thickness()};
      const double mat_X0{mat.get_material().X0()};
      const double mat_L0{mat.get_material().L0()};
 
      return material_record{inv_sf, seg, t, mat_X0, mat_L0};
    } else {
      return material_record{inv_sf, 0., 0., inv, inv};
    }
  }
};

template <template <typename...> class vector_t>
struct material_tracer : public detray::base_actor {
  using track_material_type = track_material;
  using material_record_type = material_record;

  struct state {
    friend struct material_tracer;

    DETRAY_HOST
    explicit state(vecmem::memory_resource &resource)
        : m_mat_steps(&resource) {}

    DETRAY_HOST_DEVICE
    explicit state(vector_t<material_record_type> &&steps)
        : m_mat_steps(std::move(steps)) {}

    DETRAY_HOST_DEVICE
    const auto &get_track_material() const { return m_track_material; }

    DETRAY_HOST
    auto &&release_track_material() && { return std::move(m_track_material); }

    DETRAY_HOST_DEVICE
    const auto &get_material_steps() const { return m_mat_steps; }

    DETRAY_HOST_DEVICE
    auto &&release_material_steps() && { return std::move(m_mat_steps); }
 
   private:
    track_material_type m_track_material{};

    vector_t<material_record_type> m_mat_steps{};
  };

  template <typename propagator_state_t, concepts::algebra algebra_t>
  DETRAY_HOST_DEVICE void operator()(
      state &tracer, const propagator_state_t &prop_state,
      const detray::actor::parameter_transporter_result<algebra_t> &res) const {
    record_track_dir(tracer, prop_state);
 
    // Only count material if navigator encountered it
    const auto &navigation = prop_state.navigation();
    if (!navigation.encountered_sf_material()) {
      return;
    }
 
    // For now use default context
    typename propagator_state_t::detector_type::geometry_context gctx{};
 
    const auto &track_param = res.destination_params();
    dpoint2D<algebra_t> loc_pos = track_param.bound_local();
 
    record_mat_step(tracer, gctx, navigation.current_surface(), loc_pos,
                    track_param.dir());
  }

  template <typename propagator_state_t>
  DETRAY_HOST_DEVICE void operator()(
      state &tracer, const propagator_state_t &prop_state) const {
    using algebra_t = typename propagator_state_t::detector_type::algebra_type;
 
    record_track_dir(tracer, prop_state);
 
    // Only count material if navigator encountered it
    const auto &navigation = prop_state.navigation();
    if (!navigation.encountered_sf_material()) {
      return;
    }
 
    // For now use default context
    typename propagator_state_t::detector_type::geometry_context gctx{};
 
    // Current surface
    const auto sf = navigation.current_surface();
 
    const auto &track_param = prop_state.stepping()();
    dvector3D<algebra_t> glob_dir = track_param.dir();
    dpoint2D<algebra_t> loc_pos =
        sf.global_to_bound(gctx, track_param.pos(), glob_dir);
 
    record_mat_step(tracer, gctx, sf, loc_pos, glob_dir);
  }
 
 private:
  template <typename propagator_state_t>
  DETRAY_HOST_DEVICE inline auto record_track_dir(
      state &tracer, const propagator_state_t &prop_state) const {
    using algebra_t = typename propagator_state_t::detector_type::algebra_type;
    using vector3_t = dvector3D<algebra_t>;
 
    // Record the initial track direction
    vector3_t glob_dir = prop_state.stepping()().dir();
    if (detray::detail::is_invalid_value(tracer.m_track_material.eta) &&
        detray::detail::is_invalid_value(tracer.m_track_material.phi)) {
      tracer.m_track_material.eta = vector::eta(glob_dir);
      tracer.m_track_material.phi = vector::phi(glob_dir);
    }
  }

  template <concepts::detector detector_t>
  DETRAY_HOST_DEVICE inline auto record_mat_step(
      state &tracer, const typename detector_t::geometry_context &gctx,
      const tracking_surface<detector_t> sf,
      const dpoint2D<typename detector_t::algebra_type> &loc_pos,
      const dvector3D<typename detector_t::algebra_type> &glob_dir) const {
    // Fetch the material parameters and pathlength through the material
    const auto mat_params = sf.template visit_material<get_material_record>(
        loc_pos, cos_angle(gctx, sf, glob_dir, loc_pos));
 
    const double seg{mat_params.path};
    const double t{mat_params.thickness};
    const double mx0{mat_params.mat_X0};
    const double ml0{mat_params.mat_L0};
 
    // Fill the material record
    if (mx0 > 0.) {
      tracer.m_track_material.sX0 += seg / mx0;
      tracer.m_track_material.tX0 += t / mx0;
    }
    if (ml0 > 0.) {
      tracer.m_track_material.sL0 += seg / ml0;
      tracer.m_track_material.tL0 += t / ml0;
    }
    if (t > 0.) {
      tracer.m_mat_steps.push_back({sf.identifier(), seg, t, mx0, ml0});
    }
  }
};

template <concepts::detector detector_t>
inline auto record_material(
    const typename detector_t::geometry_context /*gctx*/,
    vecmem::memory_resource *host_mr, const detector_t &det,
    const propagation::config &cfg,
    const free_track_parameters<typename detector_t::algebra_type> &track) {
  using algebra_t = typename detector_t::algebra_type;
  using scalar_t = dscalar<algebra_t>;
 
  using stepper_t = line_stepper<algebra_t>;
  using navigator_t = caching_navigator<detector_t>;
 
  // Propagator with pathlimit aborter
  using material_tracer_t = material_validator::material_tracer<vecmem::vector>;
  using pathlimit_aborter_t = actor::pathlimit_aborter<scalar_t>;
  using parameter_updater_t =
      actor::parameter_updater<algebra_t,
                               actor::pointwise_material_interactor<algebra_t>,
                               material_tracer_t>;
  using actor_chain_t = actor_chain<pathlimit_aborter_t, parameter_updater_t>;
  using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
 
  // Propagator
  propagator_t prop{cfg};
 
  // Build actor and propagator states
  typename pathlimit_aborter_t::state pathlimit_aborter_state{
      cfg.stepping.path_limit};
  actor::parameter_updater_state<algebra_t> updater_state{cfg, track};
  typename actor::pointwise_material_interactor<algebra_t>::state
      interactor_state{};
  typename material_tracer_t::state mat_tracer_state{*host_mr};
 
  auto actor_states = detray::tie(pathlimit_aborter_state, updater_state,
                                  interactor_state, mat_tracer_state);
 
  typename propagator_t::state propagation{track, det, cfg.context};
 
  // Run the propagation
  bool success = prop.propagate(propagation, actor_states);
 
  return std::make_tuple(success,
                         std::move(mat_tracer_state).release_track_material(),
                         std::move(mat_tracer_state).release_material_steps());
}

inline auto compare_traces(
    const dvector<material_record> &reference, const track_material &ref_record,
    const dvector<material_record> &mat_trace, const track_material &mat_record,
    std::size_t trk_i = detail::invalid_value<std::size_t>(),
    const double rel_tol = 0.01, const bool verbose = true) {
  // Material traces contain records of different surfaces/material
  bool is_bad_comp{reference.size() != mat_trace.size()};
  // Overall recorded material is different (in case small errors accumulate)
  bool is_diff_mat{false};
 
  std::stringstream debug_msg{};
  debug_msg << "Track No. " << trk_i << ":\n----------------" << std::endl;
 
  if (is_bad_comp) {
    debug_msg << "-> Different no. of surfaces: " << mat_trace.size()
              << " (ref.: " << reference.size() << ")\n"
              << std::endl;
  } else {
    for (std::size_t j = 0u; j < math::min(reference.size(), mat_trace.size());
         ++j) {
      if (reference[j].geo_id != mat_trace[j].geo_id) {
        is_bad_comp = true;
        debug_msg << "-> Surfaces don't match: " << mat_trace[j].geo_id
                  << " (ref.: " << reference[j].geo_id << ")" << std::endl;
        continue;
      }
 
      // TODO: Use approx_equal from algebra-plugins
      // Compare thickness of the surface material
      if ((reference[j].thickness - mat_trace[j].thickness) /
              reference[j].thickness >
          rel_tol) {
        is_bad_comp = true;
        debug_msg << "-> On surface: " << reference[j].geo_id << ":"
                  << std::endl;
        debug_msg << "-> thickness: " << mat_trace[j].thickness
                  << ", thickness ref.: " << reference[j].thickness
                  << std::endl;
      }
 
      // Compare radiation length of the surface material
      if ((reference[j].mat_X0 - mat_trace[j].mat_X0) / reference[j].mat_X0 >
          rel_tol) {
        is_bad_comp = true;
        debug_msg << "-> On surface: " << reference[j].geo_id << ":"
                  << std::endl;
        debug_msg << "-> X0: " << mat_trace[j].mat_X0
                  << ", X0 ref.: " << reference[j].mat_X0 << std::endl;
      }
 
      // Compare interaction length of the surface material
      if ((reference[j].mat_L0 - mat_trace[j].mat_L0) / reference[j].mat_L0 >
          rel_tol) {
        is_bad_comp = true;
        debug_msg << "-> On surface: " << reference[j].geo_id << ":"
                  << std::endl;
        debug_msg << "-> L0: " << mat_trace[j].mat_L0
                  << ", L0 ref.: " << reference[j].mat_L0 << std::endl;
      }
 
      // Compare path of the track through the surface material
      if ((reference[j].path - mat_trace[j].path) / reference[j].path >
          rel_tol) {
        is_bad_comp = true;
        debug_msg << "-> On surface: " << reference[j].geo_id << ":"
                  << std::endl;
        debug_msg << "-> Mat. path: " << mat_trace[j].path / unit<double>::mm
                  << " mm, mat. path ref.: "
                  << reference[j].path / unit<double>::mm << " mm" << std::endl;
      }
    }
  }
 
  // Compare the total accumulated material along the track
  const double rel_error{(ref_record.sX0 - mat_record.sX0) / ref_record.sX0};
  const bool small_mat{ref_record.sX0 < rel_tol && mat_record.sX0 < rel_tol};
 
  // If almost no material was collected, the relative error can be
  // large, but also has little consequence for tracking
  if (!(small_mat || (rel_error <= rel_tol))) {
    // Already know that material cannot add up: Only flag the additional
    // cases here
    if (!is_bad_comp) {
      is_diff_mat = true;
    }
    debug_msg << "\nTotal material discrepancy of " << 100. * rel_error << "%\n"
              << std::endl;
  }
 
  if (verbose && (is_bad_comp || is_diff_mat)) {
    DETRAY_INFO_HOST(debug_msg.str());
  }
 
  return std::make_tuple(is_bad_comp, is_diff_mat);
}

inline auto write_material(const std::string &mat_file_name,
                           const dvector<track_material> &mat_records) {
  const auto file_path = std::filesystem::path{mat_file_name};
  assert(file_path.extension() == ".csv");
 
  // Make sure path to file exists
  io::create_path(file_path.parent_path());
 
  detray::io::file_handle outfile{
      mat_file_name, std::ios::out | std::ios::binary | std::ios::trunc};
  *outfile << "eta,phi,mat_sX0,mat_sL0,mat_tX0,mat_tL0" << std::endl;
 
  for (const auto &rec : mat_records) {
    *outfile << rec.eta << "," << rec.phi << "," << rec.sX0 << "," << rec.sL0
             << "," << rec.tX0 << "," << rec.tL0 << std::endl;
  }
}
 
}  // namespace detray::material_validator