Commit 72790006 authored by Sam Wood's avatar Sam Wood Committed by Simon Spannagel
Browse files

Add constant fallback B-field config option for non-mesh regions

parent 1af19424

src/modules/MagneticFieldReader/MagneticFieldReaderModule.cpp

+26 −22
Original line number Diff line number Diff line
@@ -64,6 +64,8 @@ void MagneticFieldReaderModule::initialize() { 
        LOG(INFO) << "Set constant magnetic field: " << Units::display(b_field, {"T", "mT"});

    } else if(field_model == MagneticField::MESH) {

        type = MagneticFieldType::CUSTOM;



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

        auto field_data = read_field();



        auto dims = field_data.getDimensionality(); // Vector dimensionality of the field (should be 3 for B-field)
@@ -79,15 +81,17 @@ void MagneticFieldReaderModule::initialize() { 
            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)

        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

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

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