src/MagneticField/SolenoidBField.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/MagneticField/SolenoidBField.hpp"
0010
0011 #include "Acts/Utilities/VectorHelpers.hpp"
0012
0013 #include <algorithm>
0014 #include <cmath>
0015
0016 #define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
0017
0018 #include <boost/exception/exception.hpp>
0019 #include <boost/math/special_functions/ellint_1.hpp>
0020 #include <boost/math/special_functions/ellint_2.hpp>
0021
0022 Acts::SolenoidBField::SolenoidBField(Config config) : m_cfg(config) {
0023   m_dz = m_cfg.length / m_cfg.nCoils;
0024   m_R2 = m_cfg.radius * m_cfg.radius;
0025   // we need to scale so we reproduce the expected B field strength
0026   // at the center of the solenoid
0027   Vector2 field = multiCoilField({0, 0}, 1.);  // scale = 1
0028   m_scale = m_cfg.bMagCenter / field.norm();
0029 }
0030
0031 Acts::MagneticFieldProvider::Cache Acts::SolenoidBField::makeCache(
0032     const MagneticFieldContext& mctx) const {
0033   return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
0034 }
0035
0036 Acts::Vector3 Acts::SolenoidBField::getField(const Vector3& position) const {
0037   using VectorHelpers::perp;
0038   Vector2 rzPos(perp(position), position.z());
0039   Vector2 rzField = multiCoilField(rzPos, m_scale);
0040   Vector3 xyzField(0, 0, rzField[1]);
0041
0042   if (rzPos[0] != 0.) {
0043     // add xy field component, radially symmetric
0044     Vector3 rDir = Vector3(position.x(), position.y(), 0).normalized();
0045     xyzField += rDir * rzField[0];
0046   }
0047
0048   return xyzField;
0049 }
0050
0051 Acts::Result<Acts::Vector3> Acts::SolenoidBField::getField(
0052     const Vector3& position, MagneticFieldProvider::Cache& /*cache*/) const {
0053   return Result<Vector3>::success(getField(position));
0054 }
0055
0056 Acts::Vector2 Acts::SolenoidBField::getField(const Vector2& position) const {
0057   return multiCoilField(position, m_scale);
0058 }
0059
0060 Acts::Result<Acts::Vector3> Acts::SolenoidBField::getFieldGradient(
0061     const Vector3& position, ActsMatrix<3, 3>& /*derivative*/,
0062     MagneticFieldProvider::Cache& /*cache*/) const {
0063   return Result<Vector3>::success(getField(position));
0064 }
0065
0066 Acts::Vector2 Acts::SolenoidBField::multiCoilField(const Vector2& pos,
0067                                                    double scale) const {
0068   // iterate over all coils
0069   Vector2 resultField(0, 0);
0070   for (std::size_t coil = 0; coil < m_cfg.nCoils; coil++) {
0071     Vector2 shiftedPos =
0072         Vector2(pos[0], pos[1] + m_cfg.length * 0.5 - m_dz * (coil + 0.5));
0073     resultField += singleCoilField(shiftedPos, scale);
0074   }
0075
0076   return resultField;
0077 }
0078
0079 Acts::Vector2 Acts::SolenoidBField::singleCoilField(const Vector2& pos,
0080                                                     double scale) const {
0081   return {B_r(pos, scale), B_z(pos, scale)};
0082 }
0083
0084 double Acts::SolenoidBField::B_r(const Vector2& pos, double scale) const {
0085   //              _
0086   //     2       /  pi / 2          2    2          - 1 / 2
0087   // E (k )  =   |         ( 1  -  k  sin {theta} )         dtheta
0088   //  1         _/  0
0089   using boost::math::ellint_1;
0090   //              _          ____________________
0091   //     2       /  pi / 2| /       2    2
0092   // E (k )  =   |        |/ 1  -  k  sin {theta} dtheta
0093   //  2         _/  0
0094   using boost::math::ellint_2;
0095
0096   double r = std::abs(pos[0]);
0097   double z = pos[1];
0098
0099   if (r == 0) {
0100     return 0.;
0101   }
0102
0103   //                            _                             _
0104   //              mu  I        |  /     2 \                    |
0105   //                0     kz   |  |2 - k  |    2          2    |
0106   // B (r, z)  =  ----- ------ |  |-------|E (k )  -  E (k )   |
0107   //  r            4pi     ___ |  |      2| 2          1       |
0108   //                    | /  3 |_ \2 - 2k /                   _|
0109   //                    |/ Rr
0110   double k_2 = k2(r, z);
0111   double k = std::sqrt(k_2);
0112   double constant =
0113       scale * k * z / (4 * M_PI * std::sqrt(m_cfg.radius * r * r * r));
0114
0115   double B = (2. - k_2) / (2. - 2. * k_2) * ellint_2(k_2) - ellint_1(k_2);
0116
0117   // pos[0] is still signed!
0118   return r / pos[0] * constant * B;
0119 }
0120
0121 double Acts::SolenoidBField::B_z(const Vector2& pos, double scale) const {
0122   //              _
0123   //     2       /  pi / 2          2    2          - 1 / 2
0124   // E (k )  =   |         ( 1  -  k  sin {theta} )         dtheta
0125   //  1         _/  0
0126   using boost::math::ellint_1;
0127   //              _          ____________________
0128   //     2       /  pi / 2| /       2    2
0129   // E (k )  =   |        |/ 1  -  k  sin {theta} dtheta
0130   //  2         _/  0
0131   using boost::math::ellint_2;
0132
0133   double r = std::abs(pos[0]);
0134   double z = pos[1];
0135
0136   //                         _                                       _
0137   //             mu  I      |  /         2      \                     |
0138   //               0     k  |  | (R + r)k  - 2r |     2          2    |
0139   // B (r,z)  =  ----- ---- |  | -------------- | E (k )  +  E (k )   |
0140   //  z           4pi    __ |  |           2    |  2          1       |
0141   //                   |/Rr |_ \   2r(1 - k )   /                    _|
0142
0143   if (r == 0) {
0144     double res = scale / 2. * m_R2 / (std::sqrt(m_R2 + z * z) * (m_R2 + z * z));
0145     return res;
0146   }
0147
0148   double k_2 = k2(r, z);
0149   double k = std::sqrt(k_2);
0150   double constant = scale * k / (4 * M_PI * std::sqrt(m_cfg.radius * r));
0151   double B = ((m_cfg.radius + r) * k_2 - 2. * r) / (2. * r * (1. - k_2)) *
0152                  ellint_2(k_2) +
0153              ellint_1(k_2);
0154
0155   return constant * B;
0156 }
0157
0158 double Acts::SolenoidBField::k2(double r, double z) const {
0159   //  2           4Rr
0160   // k   =  ---------------
0161   //               2      2
0162   //        (R + r)   +  z
0163   return 4 * m_cfg.radius * r /
0164          ((m_cfg.radius + r) * (m_cfg.radius + r) + z * z);
0165 }