Merge branch 'p-detfield-improvements' into 'master' (e32c55e7) · Commits · XRFLab109 / Allpix Squared

src/core/geometry/Detector.cpp

+3 −3

Original line number	Diff line number	Diff line
		@@ -55,9 +55,9 @@ void Detector::set_model(std::shared_ptr<DetectorModel> model) {
		model_ = std::move(model);

		// Initialize the detector fields with the model:
		electric_field_.set_model(model_);
		weighting_potential_.set_model(model_);
		doping_profile_.set_model(model_);
		electric_field_.setModel(model_);
		weighting_potential_.setModel(model_);
		doping_profile_.setModel(model_);

		build_transform();
		}

src/core/geometry/DetectorField.hpp

+22 −16

Original line number	Diff line number	Diff line
		@@ -67,6 +67,26 @@ namespace allpix {
		*/
		template <typename T = ROOT::Math::XYZVector> using FieldFunction = std::function<T(const ROOT::Math::XYZPoint& pos)>;

		/**
		* @brief FieldTable is a linearized 5x5 matrix
		*/
		class FieldTable : public std::array<double, 25> {
		public:
		/**
		* @brief Helper function to translate an iterator of the array into a coordinate of the 5x5 matrix.
		*
		* The central pixel has coordinates 0,0, the others around positive or negative values, respectively. This allows
		* to directly add these coordinates to any pixel index of the sensor.
		*
		* @param it Iterator to the array
		* @return Pair of x and y coordinates.
		*/
		std::pair<int, int> getCoordinates(const FieldTable::iterator& it) {
		const auto i = std::distance(this->begin(), it);
		return {static_cast<int>(i % 5) - 2, static_cast<int>(i / 5) - 2};
		}
		};

		/**
		* @brief Helper function to invert the field vector when flipping the field direction at pixel/field boundaries
		* @param field Field value, templated to support vector fields and scalar fields
		@@ -75,19 +95,6 @@ namespace allpix {
		*/
		template <typename T> void flip_vector_components(T& field, bool x, bool y);

		/*
		* Vector field template specialization of helper function for field flipping
		*/
		template <> inline void flip_vector_components<ROOT::Math::XYZVector>(ROOT::Math::XYZVector& vec, bool x, bool y) {
		vec.SetXYZ((x ? -vec.x() : vec.x()), (y ? -vec.y() : vec.y()), vec.z());
		}

		/*
		* Scalar field template specialization of helper function for field flipping
		* Here, no inversion of the field components is required
		*/
		template <> inline void flip_vector_components<double>(double&, bool, bool) {}

		/**
		* @brief Field instance of a detector
		*
		@@ -95,7 +102,6 @@ namespace allpix {
		* scaling or offset parameters.
		*/
		template <typename T, size_t N = 3> class DetectorField {
		friend class Detector;

		public:
		/**
		@@ -162,13 +168,13 @@ namespace allpix {
		std::pair<double, double> thickness_domain,
		FieldType type = FieldType::CUSTOM);

		private:
		/**
		* @brief Set the detector model this field is used for
		* @param model The detector model
		*/
		void set_model(const std::shared_ptr<DetectorModel>& model) { model_ = model; }
		void setModel(const std::shared_ptr<DetectorModel>& model) { model_ = model; }

		private:
		/**
		* @brief Helper function to retrieve the return type from a calculated index of the field data vector
		* @param offset The calculated global index to start from

src/core/geometry/DetectorField.tpp

+37 −6

Original line number	Diff line number	Diff line
		@@ -28,16 +28,17 @@ namespace allpix {
		return {};
		}

		if(type_ == FieldType::CONSTANT) {
		// Constant field - return value:
		return function_({});
		} else if(type_ == FieldType::LINEAR \|\| type_ == FieldType::CUSTOM1D) {

		// Check if we need to extrapolate along the z axis or if is inside thickness domain:
		auto z = (extrapolate_z ? std::clamp(pos.z(), thickness_domain_.first, thickness_domain_.second) : pos.z());
		if(z < thickness_domain_.first \|\| thickness_domain_.second < z) {
		return {};
		}

		if(type_ == FieldType::CONSTANT) {
		// Constant field - return value:
		return function_({});
		} else if(type_ == FieldType::LINEAR \|\| type_ == FieldType::CUSTOM1D) {
		// Linear field or custom field function with z dependency only - calculate value from configured function:
		return function_(ROOT::Math::XYZPoint(0, 0, z));
		} else {
		@@ -52,6 +53,12 @@ namespace allpix {
		return getRelativeTo(pos, ref, extrapolate_z);
		}

		// Check if we need to extrapolate along the z axis or if is inside thickness domain:
		auto z = (extrapolate_z ? std::clamp(pos.z(), thickness_domain_.first, thickness_domain_.second) : pos.z());
		if(z < thickness_domain_.first \|\| thickness_domain_.second < z) {
		return {};
		}

		// Shift the coordinates by the offset configured for the field:
		auto x = pos.x() + offset_[0];
		auto y = pos.y() + offset_[1];
		@@ -150,6 +157,10 @@ 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
		@@ -256,4 +267,24 @@ namespace allpix {
		function_ = std::move(function);
		type_ = type;
		}

		/*
		* Vector field template specialization of helper function for field flipping
		*/
		template <> inline void flip_vector_components<ROOT::Math::XYZVector>(ROOT::Math::XYZVector& vec, bool x, bool y) {
		vec.SetXYZ((x ? -vec.x() : vec.x()), (y ? -vec.y() : vec.y()), vec.z());
		}

		/*
		* Map field template specialization of helper function for field flipping.
		* Here, no inversion of the field components is required since the map does not rotate.
		*/
		template <> inline void flip_vector_components<FieldTable>(FieldTable&, bool, bool) {}

		/*
		* Scalar field template specialization of helper function for field flipping
		* Here, no inversion of the field components is required
		*/
		template <> inline void flip_vector_components<double>(double&, bool, bool) {}

		} // namespace allpix

src/tools/field_parser.h

+39 −13

Original line number	Diff line number	Diff line
		@@ -31,7 +31,7 @@
		#include <utility>

		// Mime type version for APF files
		#define APF_MIME_TYPE_VERSION 1
		#define APF_MIME_TYPE_VERSION 2

		namespace allpix {

		@@ -42,6 +42,7 @@ namespace allpix {
		UNKNOWN = 0, ///< Unknown field quantity
		SCALAR = 1, ///< Scalar field, i.e. one entry per field position
		VECTOR = 3, ///< Vector field, i.e. three entries per field position
		MAP = 25, ///< Field with a 5x5 map for each entry
		};

		/**
		@@ -72,15 +73,17 @@ namespace allpix {
		* @param dimensions Number of bins of the field in each coordinate
		* @param size Physical extent of the field in each dimension, given in internal units
		* @param data Shared pointer to the flat field data
		* @param norm Normalization factor of field values, defaults to 1.0.
		*/
		FieldData(std::string header,
		std::array<size_t, 3> dimensions,
		std::array<T, 3> size,
		std::shared_ptr<std::vector<T>> data)
		: header_(std::move(header)), dimensions_(dimensions), size_(size), data_(std::move(data)){};
		std::shared_ptr<std::vector<T>> data,
		double norm = 1.)
		: header_(std::move(header)), dimensions_(dimensions), size_(size), data_(std::move(data)), norm_(norm){};

		/**
		* @brief Function to obtain the header (human readbale content description) of the field data
		* @brief Function to obtain the header (human readable content description) of the field data
		* @return header string
		*/
		std::string getHeader() const { return header_; }
		@@ -114,18 +117,26 @@ namespace allpix {
		return dim;
		}

		/**
		* @brief Get the norm value read from the file header. In some fields this may represent a normalization factor for
		* field values. Returns 1.0 if not read present in file.
		* @return [description]
		*/
		double getNorm() const { return norm_; }

		private:
		std::string header_;
		std::array<size_t, 3> dimensions_{};
		std::array<T, 3> size_{};
		std::shared_ptr<std::vector<T>> data_;
		double norm_{1.};

		friend class cereal::access;

		// Versioned serialization function:
		template <class Archive> void serialize(Archive& archive, std::uint32_t const version) {
		// For now, we only know one version of this file type:
		if(version != 1) {
		if(version < 1 \|\| version > 2) {
		throw std::runtime_error("unknown format version " + std::to_string(version));
		}

		@@ -133,6 +144,10 @@ namespace allpix {
		archive(header_);
		archive(dimensions_);
		archive(size_);
		// We store the norm factor only in file version 2
		if(version == 2) {
		archive(norm_);
		}
		archive(data_);
		}
		};
		@@ -145,7 +160,7 @@ namespace cereal::detail {
		static std::uint32_t registerVersion() {
		::cereal::detail::StaticObject<Versions>::getInstance().mapping.emplace(
		std::type_index(typeid(allpix::FieldData<T>)).hash_code(), APF_MIME_TYPE_VERSION);
		return 3;
		return APF_MIME_TYPE_VERSION;
		}
		static void unused() { (void)version; } // NOLINT
		}; /* end Version */
		@@ -322,7 +337,14 @@ namespace allpix {
		// WARNING the usage of this field as storage for the field units differs from the original INIT format!
		file >> tmp;
		check_unit_match(allpix::trim(tmp), units);
		file >> tmp; // ignore cluster length
		// Attempt to read the next value as sum value
		file >> tmp;
		double norm = 1.;
		try {
		norm = std::stod(tmp);
		} catch(...) {
		LOG(DEBUG) << "No norm value present in file.";
		}
		file >> tmp >> tmp >> tmp; // ignore the incident pion direction
		file >> tmp >> tmp >> tmp; // ignore the magnetic field (specify separately)
		double thickness = NAN, xpixsz = NAN, ypixsz = NAN;
		@@ -374,8 +396,11 @@ namespace allpix {
		}
		LOG_PROGRESS(INFO, "read_init") << "Reading field data: finished.";

		return FieldData<T>(
		header, std::array<size_t, 3>{{xsize, ysize, zsize}}, std::array<T, 3>{{xpixsz, ypixsz, thickness}}, field);
		return FieldData<T>(header,
		std::array<size_t, 3>{{xsize, ysize, zsize}},
		std::array<T, 3>{{xpixsz, ypixsz, thickness}},
		field,
		norm);
		}

		size_t N_;
		@@ -474,7 +499,8 @@ namespace allpix {

		// Write INIT file header
		file << field_data.getHeader() << std::endl; // Header line
		file << (!units.empty() ? units : "internal") << " ##EVENTS##" << std::endl; // Use placeholder for units
		file << (!units.empty() ? units : "internal") << " " << field_data.getNorm()
		<< std::endl; // Use placeholder for units
		file << "##TURN## ##TILT## 1.0" << std::endl; // Unused
		file << "0.0 0.0 0.0" << std::endl; // Magnetic field (unused)