Merge branch 'p-detfield-impr-2' into 'master' (c9614be3) · Commits · XRFLab109 / Allpix Squared

src/core/geometry/DetectorField.hpp

+33 −3

Original line number	Diff line number	Diff line
		@@ -174,6 +174,23 @@ namespace allpix {
		*/
		void setModel(const std::shared_ptr<DetectorModel>& model) { model_ = model; }

		protected:
		/**
		* @brief Helper to calculate field size normalization factors and configure them
		* @param bins The bins of the flat field array
		* @param size Physical extent of the field
		* @param mapping Specification of the mapping of the field onto the pixel plane
		* @param scales Scaling factors for the field size, given in fractions of the field size in x and y
		* @param offset Offset of the field from the pixel center, given in fractions of the field size in x and y
		* @param thickness_domain Domain in local coordinates in the thickness direction where the field holds
		*/
		void set_grid_parameters(std::array<size_t, 3> bins,
		std::array<double, 3> size,
		FieldMapping mapping,
		std::array<double, 2> scales,
		std::array<double, 2> offset,
		std::pair<double, double> thickness_domain);

		private:
		/**
		* @brief Helper function to retrieve the return type from a calculated index of the field data vector
		@@ -183,14 +200,27 @@ namespace allpix {
		template <std::size_t... I> inline auto get_impl(size_t offset, std::index_sequence<I...>) const noexcept;

		/**
		* @brief Helper function to calculate the field index based on the distance from its center and to return the values
		* @brief Helper function to calculate the field index based on the distance from its center
		* @param index Absolute index in the field grid
		* @param x Distance in local-coordinate x from the center of the field to obtain the values for
		* @param y Distance in local-coordinate y from the center of the field to obtain the values for
		* @param z Distance in local-coordinate z from the center of the field to obtain the values for
		* @param extrapolate_z Flag whether we should extrapolate
		* @return Value(s) of the field at the queried point
		*/
		T get_field_from_grid(const double x, const double y, const double z, const bool extrapolate_z) const noexcept;
		inline bool get_grid_index(
		size_t& index, const double x, const double y, const double z, const bool extrapolate_z) const noexcept;

		/**
		* @brief Map x and y coordinates of a position and a reference point onto a pixel given the chosen mapping.
		*
		* @param pos Position to calculate coordinates for
		* @param ref Reference position to calculate relative position to
		*
		* @return Tuple with relative x and y coordinates, mapped into the chosen area, and booleans indicating whether
		* flipping of vector components is necessary
		*/
		inline std::tuple<double, double, bool, bool> map_coordinates(const ROOT::Math::XYZPoint& pos,
		const ROOT::Math::XYPoint& ref) const;

		/**
		* @brief Fast floor-to-int implementation without overflow protection as std::floor

src/core/geometry/DetectorField.tpp

+95 −60

Original line number	Diff line number	Diff line
		@@ -81,7 +81,14 @@ namespace allpix {
		T ret_val;
		// Compute using the grid or a function depending on the setting
		if(type_ == FieldType::GRID) {
		ret_val = get_field_from_grid(x * normalization_[0] + 0.5, y * normalization_[1] + 0.5, z, extrapolate_z);
		// Calculate the linearized index of the bin in the field vector
		size_t index;
		if(!get_grid_index(index, x * normalization_[0] + 0.5, y * normalization_[1] + 0.5, z, extrapolate_z)) {
		return {};
		}

		// Fetch the field value from the given index
		ret_val = get_impl(index, std::make_index_sequence<N>{});
		} else {
		// Calculate the field from the configured function:
		ret_val = function_(ROOT::Math::XYZPoint(x, y, z));
		@@ -112,12 +119,44 @@ namespace allpix {
		return {};
		}

		T ret_val;
		if(type_ == FieldType::GRID) {
		// Map the coordinates onto the chosen pixel fraction
		auto [px, py, flip_x, flip_y] = map_coordinates(pos, ref);

		// Intentionally do floating-point equality comparison to avoid us landing on the edge of the field
		px -= (px == 1.0 ? std::numeric_limits<double>::epsilon() : 0.);
		py -= (py == 1.0 ? std::numeric_limits<double>::epsilon() : 0.);

		// Calculate the linearized index of the bin in the field vector
		size_t index;
		if(!get_grid_index(index, px, py, z, extrapolate_z)) {
		return {};
		}

		// Fetch the field value from the given index
		ret_val = get_impl(index, std::make_index_sequence<N>{});

		// Flip vector if necessary
		flip_vector_components(ret_val, flip_x, flip_y);
		} else {
		// Calculate the coordinates relative to the reference point:
		auto x = pos.x() - ref.x() + offset_[0];
		auto y = pos.y() - ref.y() + offset_[1];

		T ret_val;
		if(type_ == FieldType::GRID) {
		// Calculate the field from the configured function:
		ret_val = function_(ROOT::Math::XYZPoint(x, y, z));
		}

		return ret_val;
		}

		template <typename T, size_t N>
		std::tuple<double, double, bool, bool> DetectorField<T, N>::map_coordinates(const ROOT::Math::XYZPoint& pos,
		const ROOT::Math::XYPoint& ref) const {
		// Calculate the coordinates relative to the reference point:
		auto x = pos.x() - ref.x() + offset_[0];
		auto y = pos.y() - ref.y() + offset_[1];

		// Do we need to flip the position vector components?
		auto flip_x =
		@@ -157,41 +196,26 @@ namespace allpix {
		py += (y >= 0 ? 0. : 1.0);
		}

		// Intentionally do floating-point equality comparison to avoid us landing on the edge of the field
		px -= (px == 1.0 ? std::numeric_limits<double>::epsilon() : 0.);
		py -= (py == 1.0 ? std::numeric_limits<double>::epsilon() : 0.);

		ret_val = get_field_from_grid(px, py, z, extrapolate_z);

		// Flip vector if necessary
		flip_vector_components(ret_val, flip_x, flip_y);
		} else {
		// Calculate the field from the configured function:
		ret_val = function_(ROOT::Math::XYZPoint(x, y, z));
		}

		return ret_val;
		return {px, py, flip_x, flip_y};
		}

		// Maps the field indices onto the range of -d/2 < x < d/2, where d is the scale of the field in coordinate x.
		// This means, {x,y,z} = (0,0,0) is in the center of the field.
		template <typename T, size_t N>
		T DetectorField<T, N>::get_field_from_grid(const double x,
		const double y,
		const double z,
		const bool extrapolate_z) const noexcept {
		bool DetectorField<T, N>::get_grid_index(
		size_t& index, const double x, const double y, const double z, const bool extrapolate_z) const noexcept {

		// Compute indices
		// If the number of bins in x or y is 1, the field is assumed to be 2-dimensional and the respective index
		// is forced to zero. This circumvents that the field size in the respective dimension would otherwise be zero
		auto x_ind = (bins_[0] == 1 ? 0 : int_floor(x * static_cast<double>(bins_[0])));
		if(x_ind < 0 \|\| x_ind >= static_cast<int>(bins_[0])) {
		return {};
		return false;
		}

		auto y_ind = (bins_[1] == 1 ? 0 : int_floor(y * static_cast<double>(bins_[1])));
		if(y_ind < 0 \|\| y_ind >= static_cast<int>(bins_[1])) {
		return {};
		return false;
		}

		auto z_ind = int_floor(static_cast<double>(bins_[2]) * (z - thickness_domain_.first) /
		@@ -199,15 +223,13 @@ namespace allpix {
		// Clamp to field indices if required - we do this here (again) to not be affected by floating-point rounding:
		z_ind = (extrapolate_z ? std::clamp(z_ind, 0, static_cast<int>(bins_[2]) - 1) : z_ind);
		if(z_ind < 0 \|\| z_ind >= static_cast<int>(bins_[2])) {
		return {};
		return false;
		}

		// Compute total index
		size_t tot_ind = static_cast<size_t>(x_ind) * bins_[1] * bins_[2] * N + static_cast<size_t>(y_ind) * bins_[2] * N +
		index = static_cast<size_t>(x_ind) * bins_[1] * bins_[2] * N + static_cast<size_t>(y_ind) * bins_[2] * N +
		static_cast<size_t>(z_ind) * N;

		// Retrieve field
		return get_impl(tot_ind, std::make_index_sequence<N>{});
		return true;
		}

		/**
		@@ -232,12 +254,26 @@ namespace allpix {
		std::array<double, 2> scales,
		std::array<double, 2> offset,
		std::pair<double, double> thickness_domain) {
		if(model_ == nullptr) {
		throw std::invalid_argument("field not initialized with detector model parameters");
		}
		set_grid_parameters(bins, size, mapping, scales, offset, std::move(thickness_domain));

		if(bins[0] * bins[1] * bins[2] * N != field->size()) {
		throw std::invalid_argument("field does not match the given dimensions");
		}

		// Store the field
		field_ = std::move(field);
		};

		template <typename T, size_t N>
		void DetectorField<T, N>::set_grid_parameters(std::array<size_t, 3> bins,
		std::array<double, 3> size,
		FieldMapping mapping,
		std::array<double, 2> scales,
		std::array<double, 2> offset,
		std::pair<double, double> thickness_domain) {
		if(model_ == nullptr) {
		throw std::invalid_argument("field not initialized with detector model parameters");
		}
		if(thickness_domain.first + 1e-9 < model_->getSensorCenter().z() - model_->getSensorSize().z() / 2.0 \|\|
		model_->getSensorCenter().z() + model_->getSensorSize().z() / 2.0 < thickness_domain.second - 1e-9) {
		throw std::invalid_argument("thickness domain is outside sensor dimensions");
		@@ -246,7 +282,6 @@ namespace allpix {
		throw std::invalid_argument("end of thickness domain is before begin");
		}

		field_ = std::move(field);
		bins_ = bins;
		mapping_ = mapping;