packages/PHField/PHFieldInterpolated.h

0001 #pragma once
0002
0003 #include "PHField.h"
0004
0005 #if defined(__GNUC__) && !defined(__clang__)
0006     #pragma GCC diagnostic push
0007     #pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
0008     #include <Eigen/Dense>
0009     #pragma GCC diagnostic pop
0010 #else
0011     #include <Eigen/Dense>
0012 #endif
0013
0014 #include <array>
0015 #include <deque>
0016 #include <map>
0017 #include <optional>
0018 #include <string>
0019
0020 //! Provides a best-fit cubic approximation of the Field
0021 //! In 3D cartesian coordinates
0022 class PHFieldInterpolated : public PHField
0023 {
0024 public:
0025     typedef Eigen::Vector3i Indices_t;
0026     typedef Eigen::Vector3f Point_t;
0027     typedef Eigen::Vector3f Field_t;
0028
0029     //! There is only one fieldmap, and all points lie on a perfect grid
0030     static int const GRID_COUNT = 111;
0031     static float constexpr GRID_STEP = 20;
0032     static float constexpr GRID_MIN = - GRID_STEP * (GRID_COUNT / 2);
0033     static float constexpr GRID_MAX = + GRID_STEP * (GRID_COUNT / 2);
0034
0035     //! constructor
0036     //! @param[in] path to calibration file
0037     //! @param[in] magnetic field scaling factor (default 1.0)
0038     explicit PHFieldInterpolated () = default;
0039
0040     //! destructor
0041     ~PHFieldInterpolated () override = default;
0042
0043     //! access field value
0044     //! Follow the convention of G4ElectroMagneticField
0045     //! @param[in]  pointer to double[4] representing (x, y, z, t), where field is to be evaluated (in Geant4/CLHEP units) (cm)
0046     //! @param[out] pointer to double[3] representing (b_x, b_y, b_z), mutated to contain the value of the field (in Geant4/CLHEP units) (T)
0047     void GetFieldValue (double const*, double*) const override;
0048
0049     //! Loads the fieldmap from a ROOT file, expects specific contents
0050     void load_fieldmap (
0051         std::string const& = "/cvmfs/sphenix.sdcc.bnl.gov/alma9.2-gcc-14.2.0/release/release_ana/ana.499/share/calibrations/Field/Map/sphenix3dtrackingmapxyz.root",
0052         const float& = 1.0
0053     );
0054
0055     //! Returns the 20 element row of a design matrix for position point about the stored value m_center
0056     Eigen::VectorXf get_design_vector (Point_t const&) const;
0057     //! Caches the interpolation for the cell containing given point
0058     void cache_interpolation (Point_t const&) const;
0059     //! Returns an O(3) best fit interpolation of the field at the point, updating the cached interpolated parameters as needed
0060     Field_t get_interpolated (Point_t const&) const;
0061
0062     //! Gets the 3D grid indices from a 1D deque index
0063     static Indices_t get_indices (std::size_t const&);
0064     //! Gets the 1D deque index from 3D grid indices
0065     static std::size_t get_index (Indices_t const&);
0066
0067     //! return the 3D grid indices of the left-down-back corner of the grid cell containing point
0068     static Indices_t get_indices (Point_t const&);
0069     //! return point at the left-down-back corner of the grid cell at the given 3D grid indices
0070     static Point_t get_point (Indices_t const&);
0071
0072     //! return field at given 1D deque index
0073     Field_t get_field (std::size_t const& index) const { return m_field.at(index); }
0074     //! return field at the left-down-back corner of the grid cell at given 3D grid indices
0075     Field_t get_field (Indices_t const& indices) const { return get_field(get_index(indices)); }
0076     //! return the field at the left-down-back corner of the grid cell containing point
0077     Field_t get_field (Point_t const& point) const { return get_field(get_indices(point)); }
0078
0079     //! throw if the 3D grid indices are not in the grid and would give an invalid index
0080     static void validate_indices (Indices_t const&);
0081     //! throw if the point lies outside the grid and would give an invalid index
0082     static void validate_point (Point_t const&);
0083
0084     //! prints the contents of the loaded map to std::cout
0085     void print_map() const;
0086
0087     //! prints the cached Taylor coefficients used for interpolation to std::cout
0088     void print_coefficients() const;
0089
0090 private:
0091     //! The entire fieldmap, loaded into memory (yikes)
0092     //! deque for O(1) access, without requiring contiguous allocation
0093     std::deque<Field_t> m_field;
0094
0095     // the following should be wrapped in a class
0096     // a set of these should be maintained to buffer more cells
0097
0098     //! The indices of the left-down-back corner of the cell that the interpolation is cached for
0099     mutable Indices_t m_buffered_indices = {-1, -1, -1};
0100     //! The point about which the field is currently interpolated
0101     mutable Point_t m_center;
0102     //! The 20 coefficients for an O(3) Taylor expansion of a function in 3D
0103     //! (one per component of the magnetic field)
0104     mutable std::array<Eigen::VectorXf, 3>  m_coefficients;
0105 };