Commit c088bc12 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

Added ACTS Bfield wrapper for dd4hep fields.

- Use the field map and beamline magnetic fields as defined via compact
detector description.
- There is probably a better implementation but this is the most direct.
parent 76a682b5
......@@ -17,7 +17,7 @@ gaudi_add_library(JugBase
src/Utilities/Paths.cpp
src/Utilities/Options.cpp
src/Plugins/BFieldOptions.cpp
#src/Plugins/BFieldScalor.cpp
src/Plugins/DD4hepBField.cpp
src/Plugins/BFieldUtils.cpp
LINK
Gaudi::GaudiKernel Gaudi::GaudiAlgLib
......
// This file is part of the Acts project.
//
// Copyright (C) 2017 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
#pragma once
#include "Acts/Definitions/Algebra.hpp"
......@@ -35,33 +27,28 @@ class ScalableBField;
} // namespace Jug
using InterpolatedMapper2D = Acts::InterpolatedBFieldMapper<
Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
Acts::detail::EquidistantAxis>>;
Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
using InterpolatedMapper3D = Acts::InterpolatedBFieldMapper<Acts::detail::Grid<
Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
Acts::detail::EquidistantAxis>>;
using InterpolatedMapper3D =
Acts::InterpolatedBFieldMapper<Acts::detail::Grid<Acts::Vector3, Acts::detail::EquidistantAxis,
Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
using InterpolatedBFieldMap2D =
Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
using InterpolatedBFieldMap3D =
Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
using InterpolatedBFieldMap2D = Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
using InterpolatedBFieldMap3D = Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
namespace Jug {
namespace Options {
namespace Options {
using BFieldVariant =
std::variant<std::shared_ptr<InterpolatedBFieldMap2D>,
std::shared_ptr<InterpolatedBFieldMap3D>,
std::shared_ptr<Acts::ConstantBField>,
std::shared_ptr<Jug::BField::ScalableBField>>;
using BFieldVariant =
std::variant<std::shared_ptr<InterpolatedBFieldMap2D>, std::shared_ptr<InterpolatedBFieldMap3D>,
std::shared_ptr<Acts::ConstantBField>, std::shared_ptr<Jug::BField::ScalableBField>>;
// common bfield options, with a bf prefix
void addBFieldOptions(boost::program_options::options_description& opt);
// common bfield options, with a bf prefix
void addBFieldOptions(boost::program_options::options_description& opt);
// create the bfield maps
BFieldVariant readBField(const boost::program_options::variables_map& vm);
// create the bfield maps
BFieldVariant readBField(const boost::program_options::variables_map& vm);
} // namespace Options
} // namespace Jug
} // namespace Options
} // namespace Jug
// This file is part of the Acts project.
//
// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.
#pragma once
#include <variant>
#include "Acts/Definitions/Algebra.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/MagneticField/MagneticFieldProvider.hpp"
#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
namespace Jug::BField {
///// The Context to be handed around
//struct ScalableBFieldContext {
// double scalor = 1.;
//};
/** Use the dd4hep magnetic field in acts.
*
* \ingroup magnets
* \ingroup magsvc
*/
class DD4hepBField final : public Acts::MagneticFieldProvider {
public:
std::shared_ptr<dd4hep::Detector> m_det;
public:
struct Cache {
Cache(const Acts::MagneticFieldContext& /*mcfg*/) { }
};
Acts::MagneticFieldProvider::Cache makeCache(const Acts::MagneticFieldContext& mctx) const override
{
return Acts::MagneticFieldProvider::Cache::make<Cache>(mctx);
}
/** construct constant magnetic field from field vector.
*
* @param [in] DD4hep detector instance
*/
explicit DD4hepBField(dd4hep::Detector* det) : m_det(std::shared_ptr<dd4hep::Detector>(det)) {}
/** Retrieve magnetic field value.
*
* @param [in] position global position
* @return magnetic field vector
*/
Acts::Vector3 getField(const Acts::Vector3& position) const override;
/** retrieve magnetic field value.
*
* @param [in] position global position
* @param [in] cache Cache object (is ignored)
* @return magnetic field vector
*
* @note The @p position is ignored and only kept as argument to provide
* a consistent interface with other magnetic field services.
*/
Acts::Vector3 getField(const Acts::Vector3& position, Acts::MagneticFieldProvider::Cache& cache) const override;
/** @brief retrieve magnetic field value & its gradient
*
* @param [in] position global position
* @param [out] derivative gradient of magnetic field vector as (3x3)
* matrix
* @return magnetic field vector
*
* @note The @p position is ignored and only kept as argument to provide
* a consistent interface with other magnetic field services.
* @note currently the derivative is not calculated
* @todo return derivative
*/
Acts::Vector3 getFieldGradient(const Acts::Vector3& position,
Acts::ActsMatrix<3, 3>& /*derivative*/) const override;
/** @brief retrieve magnetic field value & its gradient
*
* @param [in] position global position
* @param [out] derivative gradient of magnetic field vector as (3x3)
* matrix
* @param [in] cache Cache object (is ignored)
* @return magnetic field vector
*
* @note The @p position is ignored and only kept as argument to provide
* a consistent interface with other magnetic field services.
* @note currently the derivative is not calculated
* @todo return derivative
*/
Acts::Vector3 getFieldGradient(const Acts::Vector3& position, Acts::ActsMatrix<3, 3>& /*derivative*/,
Acts::MagneticFieldProvider::Cache& cache) const override;
};
using BFieldVariant = std::variant<std::shared_ptr<const DD4hepBField>>;
} // namespace Jug::BField
......@@ -15,6 +15,7 @@ namespace dd4hep {
namespace Acts {
class TrackingGeometry;
class Surface;
class MagneticFieldProvider;
}
namespace genfit {
......@@ -41,7 +42,9 @@ public:
virtual std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> cellIDPositionConverter() const = 0;
virtual std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry() const = 0;
virtual std::shared_ptr<const Acts::MagneticFieldProvider> getFieldProvider() const = 0;
virtual double centralMagneticField() const = 0;
// receive Geant4 Geometry
//virtual G4VUserDetectorConstruction* getGeant4Geo() = 0;
virtual const VolumeSurfaceMap& surfaceMap() const = 0;
......
#ifndef JUGBASE_VectorHelpers_HH
#define JUGBASE_VectorHelpers_HH
#include "Acts/Definitions/Units.hpp"
#include "Acts/Definitions/Algebra.hpp"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/Objects.h"
namespace Jug::Helpers {
// Need some time to think about these...
//template<typename T>
//struct VectorToActs : public T {
// VectorToActs(const T& v) : T(v) {}
// operator Acts::Vector3() const {return {this->x(), this->y(), this->z()};}
//};
//template<typename T>
//struct ArrayToRoot : public T {
// ArrayToRoot(const T& v) : T(v) {}
// operator ROOT::Math::XYZVector() const {return {(*this)[0], (*this)[1], (*this)[1]}; }
//};
}
#endif
#include "JugBase/BField/DD4hepBField.h"
#include <cmath>
#include "Acts/Definitions/Units.hpp"
#include "Acts/Definitions/Algebra.hpp"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/Objects.h"
//#include "JugBase/VectorHelpers.hpp"
//using Vec = Jug::Helpers::VectorToActs<ROOT::Math::XYZVector>;
//using Vec2DD4hep = Jug::Helpers::ArrayToRoot<Acts::Vector3>;
namespace Jug::BField {
Acts::Vector3 DD4hepBField::getField(const Acts::Vector3& position) const
{
dd4hep::Position pos(position[0]/10.0,position[1]/10.0,position[2]/10.0);
auto field = m_det->field().magneticField(pos) * (Acts::UnitConstants::T / dd4hep::tesla);
return {field.x(), field.y(),field.z()};
}
Acts::Vector3 DD4hepBField::getField(const Acts::Vector3& position, Acts::MagneticFieldProvider::Cache& cache) const
{
dd4hep::Position pos(position[0]/10.0,position[1]/10.0,position[2]/10.0);
auto field = m_det->field().magneticField(pos) * (Acts::UnitConstants::T / dd4hep::tesla);
return {field.x(), field.y(),field.z()};
}
Acts::Vector3 DD4hepBField::getFieldGradient(const Acts::Vector3& position,
Acts::ActsMatrix<3, 3>& /*derivative*/) const
{
return this->getField(position);
}
Acts::Vector3 DD4hepBField::getFieldGradient(const Acts::Vector3& position, Acts::ActsMatrix<3, 3>& /*derivative*/,
Acts::MagneticFieldProvider::Cache& cache) const
{
return this->getField(position, cache);
}
} // namespace Jug::BField
......@@ -193,6 +193,12 @@ StatusCode GeoSvc::initialize() {
});
}
m_magneticField = std::make_shared<const Jug::BField::DD4hepBField>(m_dd4hepGeo);
for (int z : {0, 1000, 2000, 4000}) {
auto b = m_magneticField->getField({0.0, 0.0, double(z)});
debug() << "B(z=" << z << " mm) = " << b.transpose() << " T" << endmsg;
}
return StatusCode::SUCCESS;
}
......
......@@ -31,6 +31,7 @@
#include "DDRec/Surface.h"
#include "DD4hep/DD4hepUnits.h"
#include "JugBase/BField/DD4hepBField.h"
/** Draw the surfaces and save to obj file.
......@@ -69,6 +70,9 @@ private:
*/
std::shared_ptr<const dd4hep::rec::CellIDPositionConverter> m_cellid_converter = nullptr;
/// Acts magnetic field
std::shared_ptr<const Jug::BField::DD4hepBField> m_magneticField = nullptr;
/// XML-files with the detector description
Gaudi::Property<std::vector<std::string>> m_xmlFileNames{
this, "detectors", {}, "Detector descriptions XML-files"};
......@@ -114,6 +118,8 @@ public:
*/
virtual std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry() const;
virtual std::shared_ptr<const Acts::MagneticFieldProvider> getFieldProvider() const override { return m_magneticField; }
virtual double centralMagneticField() const
{
return m_dd4hepGeo->field().magneticField({0, 0, 0}).z() * (Acts::UnitConstants::T / dd4hep::tesla);
......
......@@ -38,8 +38,10 @@ gaudi_add_module(JugTrackPlugins
)
target_include_directories(JugTrackPlugins PUBLIC
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/JugBase>
$<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
)
target_compile_options(JugTrackPlugins PRIVATE -Wno-suggest-override)
......
#ifndef Jug_BFIELD_HH
#define Jug_BFIELD_HH 1
#include "Acts//Definitions/Units.hpp"
#include "Acts/Utilities/detail/AxisFwd.hpp"
#include "Acts/Utilities/detail/GridFwd.hpp"
#include <memory>
#include <tuple>
#include <variant>
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Definitions/Common.hpp"
// Forward declarations
namespace Acts {
template <typename G>
struct InterpolatedBFieldMapper;
template <typename M>
class InterpolatedBFieldMap;
class ConstantBField;
} // namespace Acts
namespace Jug {
// namespace BField {
// class ScalableBField;
//}
namespace BField {
/// The Context to be handed around
struct ScalableBFieldContext {
double scalor = 1.;
};
/// @ingroup MagneticField
///
/// @brief returns a given constant field value at every point
///
/// This class is based on the constant magnetic field class
/// but allows a event based context
class ScalableBField final {
public:
struct Cache {
double scalor = 1.;
/// @brief constructor with context
Cache(const Acts::MagneticFieldContext& mcfg)
{
scalor = std::any_cast<const ScalableBFieldContext>(mcfg).scalor;
}
};
/// @brief construct constant magnetic field from field vector
///
/// @param [in] B magnetic field vector in global coordinate system
explicit ScalableBField(Acts::Vector3 B) : m_BField(std::move(B)) {}
/// @brief construct constant magnetic field from components
///
/// @param [in] Bx magnetic field component in global x-direction
/// @param [in] By magnetic field component in global y-direction
/// @param [in] Bz magnetic field component in global z-direction
ScalableBField(double Bx = 0., double By = 0., double Bz = 0.) : m_BField(Bx, By, Bz) {}
/// @brief retrieve magnetic field value
///
/// @param [in] position global position
/// @return magnetic field vector
///
/// @note The @p position is ignored and only kept as argument to provide
/// a consistent interface with other magnetic field services.
Acts::Vector3 getField(const Acts::Vector3& /*position*/) const { return m_BField; }
/// @brief retrieve magnetic field value
///
/// @param [in] position global position
/// @param [in] cache Cache object (is ignored)
/// @return magnetic field vector
///
/// @note The @p position is ignored and only kept as argument to provide
/// a consistent interface with other magnetic field services.
Acts::Vector3 getField(const Acts::Vector3& /*position*/, Cache& cache) const
{
return m_BField * cache.scalor;
}
/// @brief retrieve magnetic field value & its gradient
///
/// @param [in] position global position
/// @param [out] derivative gradient of magnetic field vector as (3x3)
/// matrix
/// @return magnetic field vector
///
/// @note The @p position is ignored and only kept as argument to provide
/// a consistent interface with other magnetic field services.
/// @note currently the derivative is not calculated
/// @todo return derivative
Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/, Acts::ActsMatrix<3, 3>& /*derivative*/) const
{
return m_BField;
}
/// @brief retrieve magnetic field value & its gradient
///
/// @param [in] position global position
/// @param [out] derivative gradient of magnetic field vector as (3x3)
/// matrix
/// @param [in] cache Cache object (is ignored)
/// @return magnetic field vector
///
/// @note The @p position is ignored and only kept as argument to provide
/// a consistent interface with other magnetic field services.
/// @note currently the derivative is not calculated
/// @todo return derivative
Acts::Vector3 getFieldGradient(const Acts::Vector3& /*position*/, Acts::ActsMatrix<3, 3>& /*derivative*/,
Cache& cache) const
{
return m_BField * cache.scalor;
}
/// @brief check whether given 3D position is inside look-up domain
///
/// @param [in] position global 3D position
/// @return @c true if position is inside the defined look-up grid,
/// otherwise @c false
/// @note The method will always return true for the constant B-Field
bool isInside(const Acts::Vector3& /*position*/) const { return true; }
/// @brief update magnetic field vector from components
///
/// @param [in] Bx magnetic field component in global x-direction
/// @param [in] By magnetic field component in global y-direction
/// @param [in] Bz magnetic field component in global z-direction
void setField(double Bx, double By, double Bz) { m_BField << Bx, By, Bz; }
/// @brief update magnetic field vector
///
/// @param [in] B magnetic field vector in global coordinate system
void setField(const Acts::Vector3& B) { m_BField = B; }
private:
/// magnetic field vector
Acts::Vector3 m_BField;
};
} // namespace BField
} // namespace Jug
using InterpolatedMapper2D = Acts::InterpolatedBFieldMapper<
Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
using InterpolatedMapper3D =
Acts::InterpolatedBFieldMapper<Acts::detail::Grid<Acts::Vector3, Acts::detail::EquidistantAxis,
Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis>>;
using InterpolatedBFieldMap2D = Acts::InterpolatedBFieldMap<InterpolatedMapper2D>;
using InterpolatedBFieldMap3D = Acts::InterpolatedBFieldMap<InterpolatedMapper3D>;
using BFieldVariant = std::variant<std::shared_ptr<InterpolatedBFieldMap2D>, std::shared_ptr<InterpolatedBFieldMap3D>,
std::shared_ptr<Acts::ConstantBField>, std::shared_ptr<Jug::BField::ScalableBField>>;
#endif
......@@ -13,7 +13,7 @@
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugTrack/BField.h"
#include "JugBase/BField/DD4hepBField.h"
#include "JugTrack/GeometryContainers.hpp"
#include "JugTrack/IndexSourceLink.hpp"
#include "JugTrack/Measurement.hpp"
......
......@@ -15,9 +15,10 @@
#include <stdexcept>
#include <vector>
#include "JugBase/BField/DD4hepBField.h"
#include "JugTrack/IndexSourceLink.hpp"
#include "JugTrack/Track.hpp"
#include "JugTrack/BField.h"
#include "JugTrack/Measurement.hpp"
#include "JugTrack/Trajectories.hpp"
#include "JugTrack/ProtoTrack.hpp"
......
......@@ -32,12 +32,14 @@
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugBase/BField/DD4hepBField.h"
#include "JugTrack/GeometryContainers.hpp"
#include "JugTrack/Measurement.hpp"
#include "JugTrack/Index.hpp"
#include "JugTrack/IndexSourceLink.hpp"
#include "JugTrack/Track.hpp"
#include "JugTrack/BField.h"
#include "eicd/TrackerHitCollection.h"
......@@ -80,8 +82,18 @@ namespace Jug::Reco {
<< "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
return StatusCode::FAILURE;
}
m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3{0.0, 0.0, m_geoSvc->centralMagneticField()});
m_fieldctx = BFieldVariant(m_BField);
//m_BField = m_geoSvc->getFieldProvider();//std::make_shared<Acts::ConstantBField>(Acts::Vector3{0.0, 0.0, m_geoSvc->centralMagneticField()});
m_BField = std::dynamic_pointer_cast<const Jug::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
m_fieldctx = Jug::BField::BFieldVariant(m_BField);
for(int z : {0,1000,2000,4000,5899}){
auto b = m_BField->getField({0.0,0.0,double(z)})/(Acts::UnitConstants::T);
debug() << "B(z=" << z << " mm) = " << b.transpose() << " T" << endmsg;
}
// chi2 and #sourclinks per surface cutoffs
m_sourcelinkSelectorCfg = {
{Acts::GeometryIdentifier(), {15, 10}},
......
......@@ -16,7 +16,8 @@
#include <stdexcept>
#include <vector>
#include "JugTrack/BField.h"
#include "JugBase/BField/DD4hepBField.h"
#include "JugTrack/Index.hpp"
#include "JugTrack/IndexSourceLink.hpp"
#include "JugTrack/Measurement.hpp"
......@@ -28,7 +29,6 @@
//#include "Acts/Surfaces/PerigeeSurface.hpp"
#include "Acts/Definitions/Common.hpp"
#include "Acts/Geometry/TrackingGeometry.hpp"
#include "Acts/MagneticField/MagneticFieldProvider.hpp"
#include "Acts/TrackFinding/CombinatorialKalmanFilter.hpp"
#include "Acts/TrackFinding/MeasurementSelector.hpp"
......@@ -69,7 +69,7 @@ namespace Jug::Reco {
TrackFinderFunction m_trackFinderFunc;
SmartIF<IGeoSvc> m_geoSvc;
std::shared_ptr<Acts::ConstantBField> m_BField = nullptr;
std::shared_ptr<const Jug::BField::DD4hepBField> m_BField = nullptr;
Acts::GeometryContext m_geoctx;
Acts::CalibrationContext m_calibctx;
Acts::MagneticFieldContext m_fieldctx;
......
......@@ -18,7 +18,8 @@
#include "TrackFindingAlgorithm.h"
#include "JugTrack/BField.h"
#include "JugBase/BField/DD4hepBField.h"
#include <random>
#include <stdexcept>
......
......@@ -29,10 +29,11 @@
#include "JugBase/DataHandle.h"
#include "JugBase/IGeoSvc.h"
#include "JugBase/BField/DD4hepBField.h"
#include "JugTrack/GeometryContainers.hpp"
#include "JugTrack/IndexSourceLink.hpp"
#include "JugTrack/Track.hpp"
#include "JugTrack/BField.h"
#include "JugTrack/Measurement.hpp"
#include "eicd/TrackerHitCollection.h"
......@@ -68,8 +69,10 @@ namespace Jug::Reco {
<< "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
return StatusCode::FAILURE;
}
m_BField = std::make_shared<Acts::ConstantBField>(Acts::Vector3{0.0, 0.0, m_geoSvc->centralMagneticField()});
m_fieldctx = BFieldVariant(m_BField);
//m_BField = std::shared_ptr<const Jug::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
m_BField = std::dynamic_pointer_cast<const Jug::BField::DD4hepBField>(m_geoSvc->getFieldProvider());
//std::make_shared<Acts::ConstantBField>(Acts::Vector3{0.0, 0.0, m_geoSvc->centralMagneticField()});
m_fieldctx = Jug::BField::BFieldVariant(m_BField);