Commit c8254314 authored by Simon Spannagel's avatar Simon Spannagel
Browse files

MagneticFieldReader: small cleanup, addition of offset parameter, documentation

parent 1ba50be8

src/modules/DepositionGeant4/DepositionGeant4Module.cpp

+5 −8
Original line number Diff line number Diff line
@@ -468,8 +468,8 @@ void DepositionGeant4Module::finalizeThread() { 


void DepositionGeant4Module::construct_sensitive_detectors_and_fields() {

    if(geo_manager_->hasMagneticField()) {

        G4MagneticField* magneticField = new MagneticField(geo_manager_);

        G4FieldManager* globalFieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();

        auto* magneticField = new MagneticField(geo_manager_);

        auto* globalFieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();

        globalFieldMgr->SetDetectorField(magneticField);

        globalFieldMgr->CreateChordFinder(magneticField);

    }
@@ -607,12 +607,9 @@ void DepositionGeant4Module::record_module_statistics() { 


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



// The Geant4 API expects a const double Point[4], not the std::array<> the linter suggests

void MagneticField::GetFieldValue(const double Point[4], double* Bfield) const { // NOLINT

    G4cout << "Getting magnetic field from geometry manager" << G4endl;

    ROOT::Math::XYZVector bfield_vector =

        geometry_manager_->getMagneticField(ROOT::Math::XYZPoint(Point[0], Point[1], Point[2]));

    G4cout << bfield_vector << G4endl;

    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/MagneticFieldReader/MagneticFieldReaderModule.cpp

+22 −20
Original line number Diff line number Diff line
@@ -56,37 +56,39 @@ void MagneticFieldReaderModule::initialize() { 
    } 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 dims = field_data.getDimensionality(); // Vector dimensionality of the field (will be 3 for B-field)

        auto cellsize = field_data.getSize();       // Physical dimensions for each cell in the field mesh

        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, 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

            [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() / 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.));

                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;

                } else {

                    auto* field = field_mesh.get();



                }

                auto ix = static_cast<uint64_t>(xind);

                auto iy = static_cast<uint64_t>(yind);

                auto iz = static_cast<uint64_t>(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);

                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);

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

```