| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | #pragma once | ||
| 2 | |||
| 3 | #include <cmath> | ||
| 4 | |||
| 5 | #include "na64sw-config.h" | ||
| 6 | #include "str-fmt.hh" | ||
| 7 | #include "na64util/numerical/linalg.hh" | ||
| 8 | |||
| 9 | //#if defined(ROOT_FOUND) && ROOT_FOUND | ||
| 10 | //#include <TVector3.h> | ||
| 11 | //#endif | ||
| 12 | |||
| 13 | namespace na64dp { | ||
| 14 | namespace util { | ||
| 15 | |||
| 16 | /// Template adapter for `na64sw_rotation_matrix_*()` | ||
| 17 | template<typename FloatT> Matrix3T<FloatT> | ||
| 18 | rotation_matrix( na64sw_RotationOrder, const FloatT * angles ); | ||
| 19 | |||
| 20 | /// Forwards to `na64sw_rotation_matrix_f()` | ||
| 21 | template<> Matrix3T<float> | ||
| 22 | rotation_matrix<float>( na64sw_RotationOrder, const float * angles ); | ||
| 23 | |||
| 24 | /// Forwards to `na64sw_rotation_matrix_d()` | ||
| 25 | template<> Matrix3T<double> | ||
| 26 | rotation_matrix<double>( na64sw_RotationOrder, const double * angles ); | ||
| 27 | |||
| 28 | /**\brief Framework-wide standard rotation order | ||
| 29 | * | ||
| 30 | * This constant defines default axis order applied for the angular rotation | ||
| 31 | * triplets wherever appliable. | ||
| 32 | * */ | ||
| 33 | extern const na64sw_RotationOrder gStdRotationOrder; | ||
| 34 | |||
| 35 | /**\brief Converts wire units to DRS-normalized coordinate | ||
| 36 | * | ||
| 37 | * This is widely used formula that maps wire number \f$n\f$ of \f$N\f$N | ||
| 38 | * equidistant wires to range \f$(0, 1)\f$. Depending on wire coordinate type | ||
| 39 | * half-stride shift for position can be applied (applied for integer types). | ||
| 40 | * | ||
| 41 | * \warning This function assumes that conversion from discrete elements | ||
| 42 | * (slabs, tubes, etc) is done from integer types while cluster | ||
| 43 | * centers and other floating point position values already account | ||
| 44 | * for half-stride shift. | ||
| 45 | * */ | ||
| 46 | template< typename PositionT> typename std::enable_if< std::is_floating_point<PositionT>::value | ||
| 47 | , Float_t | ||
| 48 | >::type | ||
| 49 | wire_to_normalized_DRS(PositionT nWire, size_t nWires) { | ||
| 50 | return nWire/nWires - .5; | ||
| 51 | } | ||
| 52 | |||
| 53 | template< typename PositionT> typename std::enable_if< std::is_integral<PositionT>::value | ||
| 54 | , Float_t | ||
| 55 | >::type | ||
| 56 | wire_to_normalized_DRS(PositionT nWire, size_t nWires) { | ||
| 57 | return (nWire + .5)/nWires - .5; | ||
| 58 | } | ||
| 59 | |||
| 60 | /**\brief Combined transformation (rotation, translation, scaling and shear) | ||
| 61 | * | ||
| 62 | * This class defines complete set of affine transformation between two \f$\mathbb{R}_3\f$ | ||
| 63 | * spaces in terms of: | ||
| 64 | * | ||
| 65 | * * a translation vector \f$\vec{o}\f$ | ||
| 66 | * * a rotation matrix \f$R\f$ | ||
| 67 | * * a scale/shear matrix \f$W\f$ composed from local bases vectors | ||
| 68 | * \f$\vec{u}\f$, \f$\vec{v}\f$, \f$\vec{w}\f$. | ||
| 69 | * | ||
| 70 | * As a result, having a measurement vector \f$\vec{m}\f$ expressed in local | ||
| 71 | * (detector's) frame one can obtain global coordinates as: | ||
| 72 | * | ||
| 73 | * \f[ | ||
| 74 | * \vec{r} = T(\vec{m}) = \vec{o} + R (W \vec{m}) | ||
| 75 | * \f] | ||
| 76 | * | ||
| 77 | * This class can be constructed from (and can provide back) quantities | ||
| 78 | * that are more convenient to deal with: | ||
| 79 | * | ||
| 80 | * * rotation angles triplet in certain (user-defined) order | ||
| 81 | * * offset vector | ||
| 82 | * * local bases vectors of dimensions length (usually corresponding to | ||
| 83 | * detector's sensitive volume). | ||
| 84 | * | ||
| 85 | * \note Matrix of affine transform (rotation and scaling in local basis) can | ||
| 86 | * be used to "rotate" the covariance. | ||
| 87 | * | ||
| 88 | * \warning Dealing with angles, please note that due to mathematical | ||
| 89 | * ambiguity returned angular values of \f$> \pi/2\f$ can not always | ||
| 90 | * be spatially resolved matching the exact input. | ||
| 91 | * | ||
| 92 | * \todo Support or combination (merging) operations and transformations chaining | ||
| 93 | * \todo Full support for inversion | ||
| 94 | * \todo Method for covariance transformation | ||
| 95 | * | ||
| 96 | * \ingroup numerical-utils | ||
| 97 | * */ | ||
| 98 | struct Transformation { | ||
| 99 | private: | ||
| 100 | /// Set if affine and rotation matrix caches are valid | ||
| 101 | mutable bool _cachesValid; | ||
| 102 | /// (cache) spatial rotation matrix cache (order-invariant) | ||
| 103 | mutable Matrix3 _rotation; | ||
| 104 | /// (cache) shear and scale matrix (local basis) | ||
| 105 | mutable Matrix3 _shearAndScale; | ||
| 106 | /// (cache) full affine transformation matrix cache | ||
| 107 | mutable Matrix3 _affineMatrix; | ||
| 108 | protected: | ||
| 109 | /// Rotation order in use | ||
| 110 | na64sw_RotationOrder _rotationOrder; | ||
| 111 | /// True for inverted transformation | ||
| 112 | /// | ||
| 113 | /// Translation (offset) is not a linear transform, so to avoid homegenous | ||
| 114 | /// coordinates we simply use a special flag for inverted transformations | ||
| 115 | bool _invertOffset; | ||
| 116 | /// Transformation offset | ||
| 117 | Vec3 _offset; | ||
| 118 | /// Local basis (\f$\vec{u}, \vec{v}, \vec{w}\f$) as shear+scale product. | ||
| 119 | Vec3 _u, _v, _w; | ||
| 120 | /// Rotation angles | ||
| 121 | Float_t _angles[3]; | ||
| 122 | |||
| 123 | /// Recalculates affine transformation with respect to values set | ||
| 124 | void _recache_affine_transforms() const; | ||
| 125 | /// Resets "cache valid" flag (for descendant classes) | ||
| 126 | 13 | void _invalidate_caches() { _cachesValid = false; } | |
| 127 | public: | ||
| 128 | /// Creates identity transformation | ||
| 129 | Transformation(); | ||
| 130 | /// Creates transformation defined by values | ||
| 131 | Transformation( const Vec3 & offset_ | ||
| 132 | , const Vec3 & u_, const Vec3 & v_, const Vec3 & w_ | ||
| 133 | , Float_t alpha, Float_t beta, Float_t gamma | ||
| 134 | , na64sw_RotationOrder rotationOrder=gStdRotationOrder | ||
| 135 | , bool invertOffset=false | ||
| 136 | ); | ||
| 137 | /// Copy ctr | ||
| 138 | Transformation(const Transformation &); | ||
| 139 | |||
| 140 | |||
| 141 | |||
| 142 | // simple getters/setters | ||
| 143 | |||
| 144 | /// Returns `true` for reverse transformation | ||
| 145 | bool is_inverse() const { return _invertOffset; } | ||
| 146 | ///\brief Returns offset in global (lab) frame | ||
| 147 | /// | ||
| 148 | /// Offset vector set for the transformation. | ||
| 149 | /// | ||
| 150 | /// \note Returns same value of same meaning for direct and inverted | ||
| 151 | /// transform. | ||
| 152 | 18 | const Vec3 & o() const { return _offset; } | |
| 153 | /// Sets the offset vector (in lab frame) | ||
| 154 | /// | ||
| 155 | /// Does not invalidate the cache with affine transformations | ||
| 156 | 1 | void o(const Vec3 & o_) { _offset = o_; } | |
| 157 | |||
| 158 | /// Returns 1st basis vector | ||
| 159 | 3 | const Vec3 & u() const { return _u; } | |
| 160 | ///\brief Returns 1st basis vector rotated to the global frame | ||
| 161 | Vec3 gU() const { return rotation_matrix()*_u; } | ||
| 162 | ///\brief Sets 1st basis vector | ||
| 163 | /// | ||
| 164 | /// Invalidates the cache of affine transformations | ||
| 165 | 1 | void u(const Vec3 & u_) { _invalidate_caches(); _u = u_; } | |
| 166 | |||
| 167 | /// Returns 2nd basis vector | ||
| 168 | 3 | const Vec3 & v() const { return _v; } | |
| 169 | ///\brief Returns 2nd basis vector rotated to the global frame | ||
| 170 | Vec3 gV() const { return rotation_matrix()*_v; } | ||
| 171 | ///\brief Sets 2nd basis vector | ||
| 172 | /// | ||
| 173 | /// Invalidates the cache of affine transformations | ||
| 174 | 1 | void v(const Vec3 & v_) { _invalidate_caches(); _v = v_; } | |
| 175 | |||
| 176 | /// Returns 3rd basis vector | ||
| 177 | 3 | const Vec3 & w() const { return _w; } | |
| 178 | ///\brief Returns 3rd basis vector rotated to the global frame | ||
| 179 | Vec3 gW() const { return rotation_matrix()*_w; } | ||
| 180 | ///\brief Sets 3rd basis vector | ||
| 181 | /// | ||
| 182 | /// Invalidates the cache of affine transformations | ||
| 183 | 1 | void w(const Vec3 & w_) { _invalidate_caches(); _w = w_; } | |
| 184 | |||
| 185 | /// Returns rotation order | ||
| 186 | na64sw_RotationOrder rotation_order() const { return _rotationOrder; } | ||
| 187 | ///\brief Sets rotation order | ||
| 188 | /// | ||
| 189 | /// Invalidates cache of affine transformations | ||
| 190 | void rotation_order(na64sw_RotationOrder order_) | ||
| 191 | { _invalidate_caches(); _rotationOrder = order_; } | ||
| 192 | |||
| 193 | ///\brief retruns N-th rotation angle | ||
| 194 | /// | ||
| 195 | /// Returns one of three rotation angles (called in ctr as alpha, beta, | ||
| 196 | /// gamma) that can have different meaning depending on rotation order set. | ||
| 197 | /// | ||
| 198 | /// \throw Runtime error if index is out of range. | ||
| 199 | Float_t angle(int n) const; | ||
| 200 | ///\brief Sets N-th rotation angle | ||
| 201 | /// | ||
| 202 | /// Sets one of three rotation angles (called in ctr as alpha, beta, | ||
| 203 | /// gamma) that can have different meaning depending on rotation order set. | ||
| 204 | /// Invalidates cache of affine transformations. | ||
| 205 | /// | ||
| 206 | /// \throw Runtime error if index is out of range. | ||
| 207 | void angle(int n, Float_t); | ||
| 208 | |||
| 209 | |||
| 210 | // Complex getters and operations | ||
| 211 | |||
| 212 | /// Returns rotation matrix | ||
| 213 | const Matrix3 & rotation_matrix() const; | ||
| 214 | /// Returns shear-and-scale transformation matrix built on local basis | ||
| 215 | const Matrix3 & basis() const; | ||
| 216 | /// Returns full affine transform matrix (without global offset) | ||
| 217 | const Matrix3 & affine_matrix() const; | ||
| 218 | |||
| 219 | /// Applies spatial transformation to the vector, modifying it | ||
| 220 | 5 | void apply(Vec3 & m) const { m = this->operator()(m); } | |
| 221 | /// Shortcut for copying form of `apply()` | ||
| 222 |
2/2✓ Branch 1 taken 15 times.
✓ Branch 4 taken 15 times.
|
15 | Vec3 operator()(const Vec3 & m) const { return o() + affine_matrix()*m; } |
| 223 | |||
| 224 | /// Returns inverse affine transformation | ||
| 225 | Transformation inverted() const; | ||
| 226 | }; | ||
| 227 | |||
| 228 | #if 0 // TODO | ||
| 229 | /**\brief Multiple ordered chained transformation | ||
| 230 | * | ||
| 231 | * Having chained representation of affine trnaforms is useful for ceertain | ||
| 232 | * cases. | ||
| 233 | * */ | ||
| 234 | class TransformationChain : public std::vector<Transformation> { | ||
| 235 | public: | ||
| 236 | /// Applies spatial transformation to the vector, modifying it | ||
| 237 | void apply(Vec3 & m) const { m = this->operator()(m); } | ||
| 238 | /// Shortcut for copying form of `apply()` | ||
| 239 | Vec3 operator()(const Vec3 & m) const { return o() + affine_matrix()*m; } | ||
| 240 | |||
| 241 | /// Returns inverse spatial transformation | ||
| 242 | Transformation inverted() const; | ||
| 243 | }; | ||
| 244 | #endif | ||
| 245 | |||
| 246 | #if 0 | ||
| 247 | /**\brief Performs "standard" NA64sw rotation | ||
| 248 | * | ||
| 249 | * This function rotates point at given `r` by three angles provided in `rot` | ||
| 250 | * array. The `rot` array should contain axes-sorted angles (X, Y, Z), while | ||
| 251 | * the particular order of the rotations effectively applied is different | ||
| 252 | * and is the subject of convention. */ | ||
| 253 | void apply_std_rotation( Vec3 & r, const float rot[3]); | ||
| 254 | |||
| 255 | /**\brief Sets plane's cardinal vectors wrt given placements info | ||
| 256 | * | ||
| 257 | * Sets the cardinal vectors \f$o, u, v\f$ of the plane. \f$o\f$ is the | ||
| 258 | * center in a global (lab) reference system, \f$u\f$ is ort vector for | ||
| 259 | * measured coordinate and \f$v\f$ is the wire vector. | ||
| 260 | * | ||
| 261 | * Given length units will be preserved, angular measure is expected to be | ||
| 262 | * given in radians. | ||
| 263 | * | ||
| 264 | * Order of spatial vector components (placement center and size) is | ||
| 265 | * standard: x,y,z. | ||
| 266 | * | ||
| 267 | * Order of rotation angles is standard: x,y,z. Yet, note, that in NA64sw | ||
| 268 | * rotation order is defined as sequence of intrinsic rotations (Tait-Btyan | ||
| 269 | * angles) in order Z,Y,X, so Z rotation is applied first, then Y applied | ||
| 270 | * within the modified frame and, finally, X rotation is applied in modified | ||
| 271 | * frame. | ||
| 272 | * */ | ||
| 273 | void cardinal_vectors( const float center[3] | ||
| 274 | , const float size[3] | ||
| 275 | , const float rotation[3] | ||
| 276 | , Vec3 & o | ||
| 277 | , Vec3 & u, Vec3 & v, Vec3 & w | ||
| 278 | ); | ||
| 279 | #endif | ||
| 280 | |||
| 281 | #define NA64SW_FOR_EVERY_UNITS( MACRO, ... ) \ | ||
| 282 | /* length */ \ | ||
| 283 | MACRO( cm, 1, "length" __VA_ARGS__ ) \ | ||
| 284 | MACRO( m, 100, "length" __VA_ARGS__ ) \ | ||
| 285 | MACRO( mm, 0.1, "length" __VA_ARGS__ ) \ | ||
| 286 | /* magnetic field induction */ \ | ||
| 287 | MACRO( kGauss, 1, "mag.induction" __VA_ARGS__ ) \ | ||
| 288 | MACRO( Gauss, 0.001, "mag.induction" __VA_ARGS__ ) \ | ||
| 289 | MACRO( T, 10, "mag.induction" __VA_ARGS__ ) \ | ||
| 290 | /* energy */ \ | ||
| 291 | MACRO( GeV, 1, "energy" __VA_ARGS__ ) \ | ||
| 292 | MACRO( MeV, 0.001, "energy" __VA_ARGS__ ) \ | ||
| 293 | MACRO( keV, 1.0e-6, "energy" __VA_ARGS__ ) \ | ||
| 294 | /* angle */ \ | ||
| 295 | MACRO( rad, 1, "angle" __VA_ARGS__ ) \ | ||
| 296 | MACRO( deg, M_PI/180., "angle" __VA_ARGS__ ) \ | ||
| 297 | MACRO( mrad, 1e-3, "angle" __VA_ARGS__ ) \ | ||
| 298 | /*...*/ | ||
| 299 | |||
| 300 | struct Units { | ||
| 301 | #define DECLARE_CONSTEXPR(unitName, unitNum, unitType) static constexpr float unitName = unitNum; | ||
| 302 | NA64SW_FOR_EVERY_UNITS(DECLARE_CONSTEXPR) | ||
| 303 | #undef DECLARE_CONSTEXPR | ||
| 304 | |||
| 305 | /**\brief Returns converion factor for units in string form | ||
| 306 | * | ||
| 307 | * First argument must be one of the units to convert from. | ||
| 308 | * Second argument may be either units to convert to, or expected name of | ||
| 309 | * quantity to convert ("length", "energy", "mag.induction", etc), or empty | ||
| 310 | * string. Then returned value depends on: | ||
| 311 | * | ||
| 312 | * - if second argument is name of the units, a multiplication factor to | ||
| 313 | * convert from `strUnitsFrom` to `strUnitsToOrQuantity` is returned. If | ||
| 314 | * conversion is requested for units of different quantity, | ||
| 315 | * an `errors::QuantityMismatch` exception is thrown; | ||
| 316 | * - if second argument is name of the quantity, a multiplication factor | ||
| 317 | * to convert from `strUnitsFrom` to standard NA64sw units is returned, | ||
| 318 | * checking that quantity of `strUnitsFrom` is of what | ||
| 319 | * is `strUnitsToOrQuantity`. | ||
| 320 | * - if second argument is empty string, a multiplication factor | ||
| 321 | * to convert from `strUnitsFrom` to standard NA64sw units is returned. | ||
| 322 | * | ||
| 323 | * If `strUnitsFrom` does not match to any of units (of certain quantity or | ||
| 324 | * at all, an `errrors::UnknownUnit` or `errors::UnknownUnitOrQuantity` is | ||
| 325 | * thrown. | ||
| 326 | * */ | ||
| 327 | static float get_units_conversion_factor( const std::string & strUnitsFrom | ||
| 328 | , const std::string & strUnitsToOrQuantity="" ); | ||
| 329 | }; | ||
| 330 | |||
| 331 | } // namespace ::na64dp::util | ||
| 332 | |||
| 333 | |||
| 334 | // | ||
| 335 | // Errors | ||
| 336 | |||
| 337 | namespace errors { | ||
| 338 | |||
| 339 | /// Quantity unit name is unknown for the framework | ||
| 340 | class UnknownUnit : public errors::GenericRuntimeError { | ||
| 341 | public: | ||
| 342 | const std::string unitName; | ||
| 343 | const std::string unitOrQuantity; | ||
| 344 | ✗ | UnknownUnit( const std::string & unitName_ ) | |
| 345 | ✗ | : errors::GenericRuntimeError( util::format("Unknown unit name: \"%s\"", unitName_.c_str()).c_str()) | |
| 346 | ✗ | , unitName(unitName_) | |
| 347 | ✗ | {} | |
| 348 | ✗ | UnknownUnit( const std::string & unitName_ | |
| 349 | , const std::string & unitOrQuantity_ ) | ||
| 350 | ✗ | : errors::GenericRuntimeError( util::format("Unknown unit name: \"%s\"" | |
| 351 | " for conversion to/quantity \"%s\"" | ||
| 352 | , unitName_.c_str() | ||
| 353 | , unitOrQuantity_.c_str()).c_str() | ||
| 354 | ) | ||
| 355 | ✗ | , unitName(unitName_) | |
| 356 | ✗ | , unitOrQuantity(unitOrQuantity_) | |
| 357 | ✗ | {} | |
| 358 | }; | ||
| 359 | |||
| 360 | /// Quantity mismatch during unit conversion | ||
| 361 | class QuantityMismatch : public errors::GenericRuntimeError { | ||
| 362 | public: | ||
| 363 | const std::string unitName; | ||
| 364 | const std::string expectedQuantity; | ||
| 365 | ✗ | QuantityMismatch( const std::string & unitName_ | |
| 366 | , const std::string & expectedQuantity_ ) | ||
| 367 | ✗ | : errors::GenericRuntimeError( util::format("Quantity mismatch;" | |
| 368 | " got \"%s\" units, expected %s quantity." | ||
| 369 | , unitName_.c_str() | ||
| 370 | , expectedQuantity_.c_str() | ||
| 371 | ).c_str()) | ||
| 372 | ✗ | , expectedQuantity(expectedQuantity_) | |
| 373 | ✗ | {} | |
| 374 | }; | ||
| 375 | |||
| 376 | } // namespace ::na64dp::errors | ||
| 377 | |||
| 378 | } // namespace na64dp | ||
| 379 | |||
| 380 |