Merge branch 'meshed-magnetic-field' into 'master' (d381c5ec) · Commits · XRFLab109 / Allpix Squared

src/modules/DepositionGeant4/DepositionGeant4Module.cpp

+14 −12

Original line number	Diff line number	Diff line
		@@ -31,7 +31,6 @@

		#include "G4FieldManager.hh"
		#include "G4TransportationManager.hh"
		#include "G4UniformMagField.hh"

		#include "core/config/exceptions.h"
		#include "core/geometry/GeometryManager.hpp"
		@@ -469,17 +468,10 @@ void DepositionGeant4Module::finalizeThread() {

		void DepositionGeant4Module::construct_sensitive_detectors_and_fields() {
		if(geo_manager_->hasMagneticField()) {
		MagneticFieldType magnetic_field_type_ = geo_manager_->getMagneticFieldType();

		if(magnetic_field_type_ == MagneticFieldType::CONSTANT) {
		ROOT::Math::XYZVector b_field = geo_manager_->getMagneticField(ROOT::Math::XYZPoint(0., 0., 0.));
		G4MagneticField* magField = new G4UniformMagField(G4ThreeVector(b_field.x(), b_field.y(), b_field.z()));
		G4FieldManager* globalFieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
		globalFieldMgr->SetDetectorField(magField);
		globalFieldMgr->CreateChordFinder(magField);
		} else {
		throw ModuleError("Magnetic field enabled, but not constant. This can't be handled by this module yet.");
		}
		auto* magneticField = new MagneticField(geo_manager_);
		auto* globalFieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
		globalFieldMgr->SetDetectorField(magneticField);
		globalFieldMgr->CreateChordFinder(magneticField);
		}

		// Loop through all detectors and set the sensitive detector action that handles the particle passage
		@@ -612,3 +604,13 @@ void DepositionGeant4Module::record_module_statistics() {
		total_charges_ += sensor->getTotalDepositedCharge();
		}
		}

		MagneticField::MagneticField(GeometryManager* geometry_manager) : geometry_manager_(geometry_manager){};

		void MagneticField::GetFieldValue(const double Point[4], double* Bfield) const { // NOLINT
		const auto point = ROOT::Math::XYZPoint(Point[0], Point[1], Point[2]);
		ROOT::Math::XYZVector bfield_vector = geometry_manager_->getMagneticField(point);
		Bfield[0] = bfield_vector.x();
		Bfield[1] = bfield_vector.y();
		Bfield[2] = bfield_vector.z();
		}

src/modules/DepositionGeant4/DepositionGeant4Module.hpp

+12 −0

Original line number	Diff line number	Diff line
		@@ -16,6 +16,7 @@
		#include <memory>
		#include <string>

		#include <G4MagneticField.hh>
		#include <G4UserLimits.hh>

		#include "core/config/Configuration.hpp"
		@@ -133,6 +134,17 @@ namespace allpix {
		// Mutex used for the construction of histograms
		std::mutex histogram_mutex_;
		};

		class MagneticField : public G4MagneticField {
		protected:
		GeometryManager* geometry_manager_;

		public:
		explicit MagneticField(GeometryManager* geometry_manager);
		~MagneticField() override = default;
		// The Geant4 API expects a const double Point[4], not the std::array<> the linter suggests
		virtual void GetFieldValue(const double Point[4], double* Bfield) const override; // NOLINT
		};
		} // namespace allpix

		#endif /* ALLPIX_SIMPLE_DEPOSITION_MODULE_H */

src/modules/MagneticFieldReader/MagneticFieldReaderModule.cpp

+77 −9

Original line number	Diff line number	Diff line
		@@ -52,6 +52,49 @@ void MagneticFieldReaderModule::initialize() {
		MagneticFieldFunction function = [b_field](const ROOT::Math::XYZPoint&) { return b_field; };

		geometryManager_->setMagneticFieldFunction(function, type);
		LOG(INFO) << "Set constant magnetic field: " << Units::display(b_field, {"T", "mT"});
		} else if(field_model == MagneticField::MESH) {
		type = MagneticFieldType::CUSTOM;

		auto offset = config_.get<ROOT::Math::XYZVector>("offset", ROOT::Math::XYZVector());

		auto fallback_field = config_.get<ROOT::Math::XYZVector>("magnetic_field", ROOT::Math::XYZVector());
		auto field_data = read_field();

		auto fieldsize = field_data.getSize(); // Physical dimensions of the field
		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

		MagneticFieldFunction function =
		[fallback_field, field_mesh, fieldsize, ncells, offset](const ROOT::Math::XYZPoint& coord) {
		// Find the nearest field mesh cell to the given input coordinate (assuming the mesh is centered about the
		// origin)
		auto xind = std::round(((coord.X() + offset.X()) * static_cast<double>(ncells[0]) / fieldsize[0]) +
		(static_cast<double>(ncells[0]) / 2.));
		auto yind = std::round(((coord.Y() + offset.Y()) * static_cast<double>(ncells[1]) / fieldsize[1]) +
		(static_cast<double>(ncells[1]) / 2.));
		auto zind = std::round(((coord.Z() + offset.Z()) * static_cast<double>(ncells[2]) / fieldsize[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 fallback_field;
		}
		auto ix = static_cast<uint64_t>(xind);
		auto iy = static_cast<uint64_t>(yind);
		auto iz = static_cast<uint64_t>(zind);

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

		geometryManager_->setMagneticFieldFunction(function, type);
		LOG(INFO) << "Set meshed magnetic field from file";
		}

		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
		@@ -61,6 +104,31 @@ void MagneticFieldReaderModule::initialize() {
		LOG(DEBUG) << "Magnetic field in detector " << detector->getName() << ": "
		<< Units::display(detector->getMagneticField(detector->getLocalPosition(position)), {"T", "mT"});
		}
		LOG(INFO) << "Set constant magnetic field: " << Units::display(b_field, {"T", "mT"});
		}

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

src/modules/MagneticFieldReader/README.md

+28 −5

Original line number	Diff line number	Diff line
		@@ -9,21 +9,44 @@ module_maintainers: ["Paul Schuetze (<paul.schuetze@desy.de>)"]

		## Description

		Unique module, adds a magnetic field to the full volume, including the active sensors. By default, the magnetic field is turned off.
		Unique module, adds a magnetic field to the full volume of the simulation, including the active sensors.
		By default, the magnetic field is turned off.

		The magnetic field reader only provides constant magnetic fields, read in as a three-dimensional vector. The magnetic field is forwarded to the GeometryManager, enabling the magnetic field for the particle propagation via Geant4, as well as to all detectors for enabling a Lorentz drift during the charge propagation.
		The magnetic field reader provides the possibility of using a simple constant magnetic field permeating the entire simulated setup,
		or a meshed field which is centered around the origin of the global coordinate system. The magnetic field is forwarded to the
		GeometryManager, enabling the magnetic field for the particle propagation via Geant4, as well as to all detectors for
		enabling a Lorentz drift during the charge propagation.

		For the constant model, the field is set as a three-dimensional vector.

		For the mesh model, the field needs to be provided in form of an APF or INIT file which provides total size of the field,
		the bin size of the mesh as well as the actual field data. In addition, an offset of this field from the origin of the global
		coordinate system can be provided via the `offset` parameter. Here, the value set via `magnetic_field` is sued as a fallback
		field value used outside the volume provided by the field file.

		## Parameters

		* `model` : Type of the magnetic field model, currently only constant possible.
		* `magnetic_field` : Vector describing the magnetic field.
		* `model` : Type of the magnetic field model, possible values are constant and mesh.
		* `magnetic_field` : Vector describing the magnetic field. In the model mesh this is used as the fallback field outside the meshed region.
		* `file_name` : Path to the APF or INIT file containing the magnetic field to be used. Only used in the mesh model.
		* `offset` : Offset of the meshed magnetic field center from the global center of origin of the simulation. Defaults to `0, 0, 0` and is only used in the mesh model.

		## Usage

		An example is given below
		An example for a constant magnetic field is

		```ini
		[MagneticFieldReader]
		model = "constant"
		magnetic_field = 500mT 3.8T 0T
		```

		The configuration for a meshed field may look like the following:

		```ini
		[MagneticFieldReader]
		model = "mesh"
		magnetic_field = 500mT 0T 0T
		file_name = "path/to/magnetic_field.apf"
		offset = 5cm, 6cm, 4mm
		```