Merge branch 'perf_detfield' into 'master' (46b3ef24) · Commits · XRFLab109 / Allpix Squared

src/core/geometry/DetectorField.hpp

+17 −10

Original line number	Diff line number	Diff line
		@@ -35,6 +35,7 @@ namespace allpix {
		CONSTANT, ///< Constant field
		LINEAR, ///< Linear field (linearity determined by function)
		GRID, ///< Field supplied through a regularized grid
		CUSTOM1D, ///< Custom field function, dependent only on z
		CUSTOM, ///< Custom field function
		};

		@@ -90,7 +91,7 @@ namespace allpix {
		/**
		* @brief Field instance of a detector
		*
		* Contains the a pointer to the field dat along with the field sizes, binning and potential field distortions such as
		* Contains the a pointer to the field data along with the field sizes, binning and potential field distortions such as
		* scaling or offset parameters.
		*/
		template <typename T, size_t N = 3> class DetectorField {
		@@ -173,20 +174,26 @@ namespace allpix {
		* @param offset The calculated global index to start from
		* @note The index sequence is expanded to the number of elements requested, depending on the template instance
		*/
		template <std::size_t... I> auto get_impl(size_t offset, std::index_sequence<I...>) const;
		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
		* @param dist Distance from the center of the field to obtain the values for, given in local coordinates
		* @param extrapolate_z Switch to either extrapolate the field along z when outside the grid or return zero
		* @param flip_x Flip vector component x of resulting field vector
		* @param flip_y Flip vector component y of resulting field vector
		* @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
		* @return Value(s) of the field at the queried point
		*/
		T get_field_from_grid(const ROOT::Math::XYZPoint& dist,
		const bool extrapolate_z = false,
		const bool flip_x = false,
		const bool flip_y = false) const;
		T get_field_from_grid(const double x, const double y, const double z) const noexcept;

		/**
		* @brief Fast floor-to-int implementation without overflow protection as std::floor
		* @param x Double-precision floating point value
		* @return Integer floored towards negative infinity
		* */
		static inline int int_floor(double x) noexcept {
		auto i = static_cast<int>(x);
		return i - static_cast<int>(static_cast<double>(i) > x);
		};

		/**
		* Field properties

src/core/geometry/DetectorField.tpp

+64 −67

Original line number	Diff line number	Diff line
		@@ -28,6 +28,20 @@ namespace allpix {
		return {};
		}

		// 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 {

		// For per-pixel fields, resort to getRelativeTo with current pixel as reference:
		if(mapping_ != FieldMapping::SENSOR) {
		// Calculate center of current pixel from index as reference point:
		@@ -41,14 +55,13 @@ namespace allpix {
		// Shift the coordinates by the offset configured for the field:
		auto x = pos.x() + offset_[0];
		auto y = pos.y() + offset_[1];
		auto z = pos.z();

		auto pitch = model_->getPixelSize();

		// Compute corresponding field replica coordinates:
		// WARNING This relies on the origin of the local coordinate system
		auto replica_x = static_cast<int>(std::floor((x + 0.5 * pitch.x()) * normalization_[0]));
		auto replica_y = static_cast<int>(std::floor((y + 0.5 * pitch.y()) * normalization_[1]));
		auto replica_x = int_floor((x + 0.5 * pitch.x()) * normalization_[0]);
		auto replica_y = int_floor((y + 0.5 * pitch.y()) * normalization_[1]);

		// Convert to the replica frame:
		x -= (replica_x + 0.5) / normalization_[0] - 0.5 * pitch.x();
		@@ -58,19 +71,11 @@ namespace allpix {
		x *= ((replica_x % 2) == 1 ? -1 : 1);
		y *= ((replica_y % 2) == 1 ? -1 : 1);

		// Compute using the grid or a function depending on the setting
		T ret_val;
		// Compute using the grid or a function depending on the setting
		if(type_ == FieldType::GRID) {
		ret_val = get_field_from_grid(
		ROOT::Math::XYZPoint(x * normalization_[0] + 0.5, y * normalization_[1] + 0.5, pos.z()), extrapolate_z);
		ret_val = get_field_from_grid(x * normalization_[0] + 0.5, y * normalization_[1] + 0.5, pos.z());
		} else {
		// Check if we need to extrapolate along the z axis or if is inside thickness domain:
		if(extrapolate_z) {
		z = std::clamp(z, thickness_domain_.first, thickness_domain_.second);
		} else if(z < thickness_domain_.first \|\| thickness_domain_.second < z) {
		return {};
		}

		// Calculate the field from the configured function:
		ret_val = function_(ROOT::Math::XYZPoint(x, y, z));
		}
		@@ -79,6 +84,7 @@ namespace allpix {
		flip_vector_components(ret_val, replica_x % 2, replica_y % 2);
		return ret_val;
		}
		}

		/**
		* Get a value from the field assigned to a specific pixel. This means, we cannot wrap around at the pixel edges and
		@@ -93,10 +99,15 @@ namespace allpix {
		return {};
		}

		// 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 {};
		}

		// Calculate the coordinates relative to the reference point:
		auto x = pos.x() - ref.x() + offset_[0];
		auto y = pos.y() - ref.y() + offset_[1];
		auto z = pos.z();

		T ret_val;
		if(type_ == FieldType::GRID) {
		@@ -139,15 +150,11 @@ namespace allpix {
		py += (y >= 0 ? 0. : 1.0);
		}

		ret_val = get_field_from_grid(ROOT::Math::XYZPoint(px, py, z), extrapolate_z, flip_x, flip_y);
		} else {
		// Check if we need to extrapolate along the z axis or if is inside thickness domain:
		if(extrapolate_z) {
		z = std::clamp(z, thickness_domain_.first, thickness_domain_.second);
		} else if(z < thickness_domain_.first \|\| thickness_domain_.second < z) {
		return {};
		}
		ret_val = get_field_from_grid(px, py, 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));
		}
		@@ -158,31 +165,24 @@ namespace allpix {
		// 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 ROOT::Math::XYZPoint& dist,
		const bool extrapolate_z,
		const bool flip_x,
		const bool flip_y) const {
		T DetectorField<T, N>::get_field_from_grid(const double x, const double y, const double 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 : static_cast<int>(std::floor(dist.x() * static_cast<double>(bins_[0]))));
		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 {};
		}

		auto y_ind = (bins_[1] == 1 ? 0 : static_cast<int>(std::floor(dist.y() * static_cast<double>(bins_[1]))));
		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 {};
		}

		auto z_ind = static_cast<int>(std::floor(static_cast<double>(bins_[2]) * (dist.z() - thickness_domain_.first) /
		(thickness_domain_.second - thickness_domain_.first)));

		// Check if we need to extrapolate along the z axis:
		if(extrapolate_z) {
		z_ind = std::clamp(z_ind, 0, static_cast<int>(bins_[2]) - 1);
		} else if(z_ind < 0 \|\| z_ind >= static_cast<int>(bins_[2])) {
		auto z_ind = int_floor(static_cast<double>(bins_[2]) * (z - thickness_domain_.first) /
		(thickness_domain_.second - thickness_domain_.first));
		if(z_ind < 0 \|\| z_ind >= static_cast<int>(bins_[2])) {
		return {};
		}

		@@ -191,10 +191,7 @@ namespace allpix {
		static_cast<size_t>(z_ind) * N;

		// Retrieve field
		auto field_vector = get_impl(tot_ind, std::make_index_sequence<N>{});
		// Flip sign of vector components if necessary
		flip_vector_components(field_vector, flip_x, flip_y);
		return field_vector;
		return get_impl(tot_ind, std::make_index_sequence<N>{});
		}

		/**
		@@ -204,7 +201,7 @@ namespace allpix {
		*/
		template <typename T, size_t N>
		template <std::size_t... I>
		auto DetectorField<T, N>::get_impl(size_t offset, std::index_sequence<I...>) const {
		auto DetectorField<T, N>::get_impl(size_t offset, std::index_sequence<I...>) const noexcept {
		return T{(*field_)[offset + I]...};
		}

src/core/geometry/DetectorModel.cpp

+2 −3

Original line number	Diff line number	Diff line
		@@ -288,15 +288,14 @@ std::vector<SupportLayer> DetectorModel::getSupportLayers() const {
		std::optional<DetectorModel::Implant> DetectorModel::isWithinImplant(const ROOT::Math::XYZPoint& local_pos) const {

		// Bail out if we have no implants - no need to transform coordinates:
		auto implants = getImplants();
		if(implants.empty()) {
		if(implants_.empty()) {
		return std::nullopt;
		}

		auto [xpixel, ypixel] = getPixelIndex(local_pos);
		auto inPixelPos = local_pos - getPixelCenter(xpixel, ypixel);

		for(const auto& implant : implants) {
		for(const auto& implant : implants_) {
		if(implant.contains(inPixelPos)) {
		return implant;
		}

src/modules/DopingProfileReader/DopingProfileReaderModule.cpp

+2 −5

Original line number	Diff line number	Diff line
		@@ -27,7 +27,6 @@ DopingProfileReaderModule::DopingProfileReaderModule(Configuration& config, Mess
		}

		void DopingProfileReaderModule::initialize() {
		FieldType type = FieldType::GRID;

		// Check field strength
		auto field_model = config_.get<DopingProfile>("model");
		@@ -80,16 +79,14 @@ void DopingProfileReaderModule::initialize() {

		} else if(field_model == DopingProfile::CONSTANT) {
		LOG(TRACE) << "Adding constant doping concentration";
		type = FieldType::CONSTANT;

		auto concentration = config_.get<double>("doping_concentration");
		LOG(INFO) << "Set constant doping concentration of " << Units::display(concentration, {"/cm/cm/cm"});
		FieldFunction<double> function = [concentration](const ROOT::Math::XYZPoint&) noexcept { return concentration; };

		detector_->setDopingProfileFunction(function, type);
		detector_->setDopingProfileFunction(function, FieldType::CONSTANT);
		} else if(field_model == DopingProfile::REGIONS) {
		LOG(TRACE) << "Adding doping concentration depending on sensor region";
		type = FieldType::CUSTOM;

		auto concentration = config_.getMatrix<double>("doping_concentration");
		std::map<double, double> concentration_map;
		@@ -115,7 +112,7 @@ void DopingProfileReaderModule::initialize() {
		}
		};

		detector_->setDopingProfileFunction(function, type);
		detector_->setDopingProfileFunction(function, FieldType::CUSTOM1D);
		}

		// Produce doping_concentration_histograms if needed

src/modules/ElectricFieldReader/ElectricFieldReaderModule.cpp

+14 −10

Original line number	Diff line number	Diff line
		@@ -128,10 +128,11 @@ void ElectricFieldReaderModule::initialize() {
		<< Units::display(config_.get<double>("minimum_position"), {"um", "mm"}) << " and maximum field "
		<< Units::display(config_.get<double>("maximum_field"), "V/cm") << " at electrode";
		detector_->setElectricFieldFunction(
		get_parabolic_field_function(thickness_domain), thickness_domain, FieldType::CUSTOM);
		get_parabolic_field_function(thickness_domain), thickness_domain, FieldType::CUSTOM1D);
		} else if(field_model == ElectricField::CUSTOM) {
		LOG(TRACE) << "Adding custom electric field";
		detector_->setElectricFieldFunction(get_custom_field_function(), thickness_domain, FieldType::CUSTOM);
		auto [field_function, field_type] = get_custom_field_function();
		detector_->setElectricFieldFunction(field_function, thickness_domain, field_type);
		}

		// Produce histograms if needed
		@@ -195,7 +196,7 @@ ElectricFieldReaderModule::get_parabolic_field_function(std::pair<double, double
		};
		}

		FieldFunction<ROOT::Math::XYZVector> ElectricFieldReaderModule::get_custom_field_function() {
		std::pair<FieldFunction<ROOT::Math::XYZVector>, FieldType> ElectricFieldReaderModule::get_custom_field_function() {

		auto field_functions = config_.getArray<std::string>("field_function");
		auto field_parameters = config_.getArray<double>("field_parameters");
		@@ -219,9 +220,10 @@ FieldFunction<ROOT::Math::XYZVector> ElectricFieldReaderModule::get_custom_field
		}

		LOG(DEBUG) << "Value of custom field at pixel center: " << Units::display(z->Eval(0., 0., 0.), "V/cm");
		return [z = std::move(z)](const ROOT::Math::XYZPoint& pos) {
		return {[z = std::move(z)](const ROOT::Math::XYZPoint& pos) {
		return ROOT::Math::XYZVector(0, 0, z->Eval(pos.x(), pos.y(), pos.z()));
		};
		},
		FieldType::CUSTOM1D};
		} else if(field_functions.size() == 3) {
		LOG(DEBUG) << "Found definition of 3D custom field, applying to three Cartesian axes";
		auto x = std::make_shared<TFormula>("ex", field_functions.at(0).c_str(), false);
		@@ -250,10 +252,12 @@ FieldFunction<ROOT::Math::XYZVector> ElectricFieldReaderModule::get_custom_field
		LOG(DEBUG) << "Value of custom field at pixel center: "
		<< Units::display(ROOT::Math::XYZVector(x->Eval(0., 0., 0.), y->Eval(0., 0., 0.), z->Eval(0., 0., 0.)),
		{"V/cm"});
		return [x = std::move(x), y = std::move(y), z = std::move(z)](const ROOT::Math::XYZPoint& pos) {
		return ROOT::Math::XYZVector(
		x->Eval(pos.x(), pos.y(), pos.z()), y->Eval(pos.x(), pos.y(), pos.z()), z->Eval(pos.x(), pos.y(), pos.z()));
		};
		return {[x = std::move(x), y = std::move(y), z = std::move(z)](const ROOT::Math::XYZPoint& pos) {
		return ROOT::Math::XYZVector(x->Eval(pos.x(), pos.y(), pos.z()),
		y->Eval(pos.x(), pos.y(), pos.z()),
		z->Eval(pos.x(), pos.y(), pos.z()));
		},
		FieldType::CUSTOM};
		} else {
		throw InvalidValueError(config_,
		"field_function",