Add logic to load B-field from mesh file (9ca63005) · Commits · XRFLab109 / Allpix Squared

src/modules/MagneticFieldReader/MagneticFieldReaderModule.cpp

+78 −0

Original line number	Diff line number	Diff line
		@@ -62,5 +62,83 @@ void MagneticFieldReaderModule::initialize() {
		<< Units::display(detector->getMagneticField(detector->getLocalPosition(position)), {"T", "mT"});
		}
		LOG(INFO) << "Set constant magnetic field: " << Units::display(b_field, {"T", "mT"});
		} else if(field_model == MagneticField::MESH) {
		type = MagneticFieldType::CUSTOM;
		auto field_data = read_field();

		auto dims = field_data.getDimensionality(); // Vector dimensionality of the field (should be 3 for B-field)
		if(dims != 3) {
		throw std::invalid_argument("B-field vectors were not three dimensional.");
		}

		auto cellsize = field_data.getSize(); // Physical dimensions for each cell in the field mesh
		auto ncells = field_data.getDimensions(); // Number of cells in each direction of the field mesh
		auto field_mesh = field_data.getData(); // The field mesh B-field data as a flattened array

		if(ncells[0] * ncells[1] * ncells[2] * dims != field_mesh->size()) {
		throw std::invalid_argument("B-field does not span the given dimensions");
		}

		MagneticFieldFunction function = [field_mesh, cellsize, ncells, dims](const ROOT::Math::XYZPoint& coord) {
		// Find the nearest field mesh cell to the given input coordinate (assuming the mesh is centred about the origin)
		auto xind = std::round((coord.X() / cellsize[0]) + (static_cast<double>(ncells[0]) / 2.));
		auto yind = std::round((coord.Y() / cellsize[1]) + (static_cast<double>(ncells[1]) / 2.));
		auto zind = std::round((coord.Z() / cellsize[2]) + (static_cast<double>(ncells[2]) / 2.));

		if(xind < 0 \|\| xind >= static_cast<double>(ncells[0]) \|\| yind < 0 \|\| yind >= static_cast<double>(ncells[1]) \|\|
		zind < 0 \|\| zind >= static_cast<double>(ncells[2])) {
		return ROOT::Math::XYZVector(0, 0, 0); // Default to no field if mesh doesn't cover input coordinate
		} else {
		auto field = field_mesh.get();

		long unsigned int ix = static_cast<long unsigned int>(xind);
		long unsigned int iy = static_cast<long unsigned int>(yind);
		long unsigned int iz = static_cast<long unsigned int>(zind);

		auto bfieldx = field->at(ix * ncells[1] * ncells[2] * dims + iy * ncells[2] * dims + iz * dims);
		auto bfieldy = field->at(ix * ncells[1] * ncells[2] * dims + iy * ncells[2] * dims + iz * dims + 1);
		auto bfieldz = field->at(ix * ncells[1] * ncells[2] * dims + iy * ncells[2] * dims + iz * dims + 2);
		return ROOT::Math::XYZVector(bfieldx, bfieldy, bfieldz);
		}
		};

		geometryManager_->setMagneticFieldFunction(function, type);
		auto detectors = geometryManager_->getDetectors();
		for(auto& detector : detectors) {
		// TODO the magnetic field is calculated once for the center position of the detector. This could be extended to
		// a function enabling a gradient in the magnetic field inside the sensor
		auto position = detector->getPosition();
		detector->setMagneticField(detector->getOrientation().Inverse() * geometryManager_->getMagneticField(position));
		LOG(DEBUG) << "Magnetic field in detector " << detector->getName() << ": "
		<< Units::display(detector->getMagneticField(detector->getLocalPosition(position)), {"T", "mT"});
		}
		LOG(INFO) << "Set meshed magnetic field from file";
		}
		}

		/**
		* The field data read from files are shared between module instantiations using the static
		* FieldParser's getByFileName method.
		*/
		FieldParser<double> MagneticFieldReaderModule::field_parser_(FieldQuantity::VECTOR);
		FieldData<double> MagneticFieldReaderModule::read_field() {

		try {
		LOG(TRACE) << "Fetching magnetic field from mesh file";

		// Get field from file
		auto field_data = field_parser_.getByFileName(config_.getPath("file_name", true), "T");

		LOG(INFO) << "Set magnetic field with " << field_data.getDimensions().at(0) << "x"
		<< field_data.getDimensions().at(1) << "x" << field_data.getDimensions().at(2) << " cells";

		// Return the field data
		return field_data;
		} catch(std::invalid_argument& e) {
		throw InvalidValueError(config_, "file_name", e.what());
		} catch(std::runtime_error& e) {
		throw InvalidValueError(config_, "file_name", e.what());
		} catch(std::bad_alloc& e) {
		throw InvalidValueError(config_, "file_name", "file too large");
		}
		}

src/modules/MagneticFieldReader/MagneticFieldReaderModule.hpp

+9 −0

Original line number	Diff line number	Diff line
		@@ -17,6 +17,7 @@
		#include "core/config/Configuration.hpp"
		#include "core/geometry/GeometryManager.hpp"
		#include "core/messenger/Messenger.hpp"
		#include "tools/field_parser.h"

		#include "core/module/Module.hpp"

		@@ -34,6 +35,7 @@ namespace allpix {
		*/
		enum class MagneticField {
		CONSTANT, ///< Constant magnetic field
		MESH, ///< Magnetic field defined by a mesh
		};

		public:
		@@ -52,5 +54,12 @@ namespace allpix {

		private:
		GeometryManager* geometryManager_;

		/**
		* @brief Read field from a file in init or apf format
		* @return Data of the field read from file
		*/
		FieldData<double> read_field();
		static FieldParser<double> field_parser_;
		};
		} // namespace allpix