EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SolenoidBField.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SolenoidBField.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
12 
13 #include <boost/math/special_functions/ellint_1.hpp>
14 #include <boost/math/special_functions/ellint_2.hpp>
15 
16 Acts::SolenoidBField::SolenoidBField(Config config) : m_cfg(std::move(config)) {
19  // we need to scale so we reproduce the expected B field strength
20  // at the center of the solenoid
21  Vector2D field = multiCoilField({0, 0}, 1.); // scale = 1
22  m_scale = m_cfg.bMagCenter / field.norm();
23 }
24 
26  using VectorHelpers::perp;
27  Vector2D rzPos(perp(position), position.z());
28  Vector2D rzField = multiCoilField(rzPos, m_scale);
29  Vector3D xyzField(0, 0, rzField[1]);
30 
31  if (rzPos[0] != 0.) {
32  // add xy field component, radially symmetric
33  Vector3D rDir = Vector3D(position.x(), position.y(), 0).normalized();
34  xyzField += rDir * rzField[0];
35  }
36 
37  return xyzField;
38 }
39 
41  Cache& /*cache*/) const {
42  return getField(position);
43 }
44 
46  return multiCoilField(position, m_scale);
47 }
48 
50  const Vector3D& position, ActsMatrixD<3, 3>& /*derivative*/) const {
51  return getField(position);
52 }
53 
55  const Vector3D& position, ActsMatrixD<3, 3>& /*derivative*/,
56  Cache& /*cache*/) const {
57  return getField(position);
58 }
59 
61  double scale) const {
62  // iterate over all coils
63  Vector2D resultField(0, 0);
64  for (size_t coil = 0; coil < m_cfg.nCoils; coil++) {
65  Vector2D shiftedPos =
66  Vector2D(pos[0], pos[1] + m_cfg.length * 0.5 - m_dz * (coil + 0.5));
67  resultField += singleCoilField(shiftedPos, scale);
68  }
69 
70  return resultField;
71 }
72 
74  double scale) const {
75  return {B_r(pos, scale), B_z(pos, scale)};
76 }
77 
78 double Acts::SolenoidBField::B_r(const Vector2D& pos, double scale) const {
79  // _
80  // 2 / pi / 2 2 2 - 1 / 2
81  // E (k ) = | ( 1 - k sin {theta} ) dtheta
82  // 1 _/ 0
83  using boost::math::ellint_1;
84  // _ ____________________
85  // 2 / pi / 2| / 2 2
86  // E (k ) = | |/ 1 - k sin {theta} dtheta
87  // 2 _/ 0
88  using boost::math::ellint_2;
89 
90  double r = std::abs(pos[0]);
91  double z = pos[1];
92 
93  if (r == 0) {
94  return 0.;
95  }
96 
97  // _ _
98  // mu I | / 2 \ |
99  // 0 kz | |2 - k | 2 2 |
100  // B (r, z) = ----- ------ | |-------|E (k ) - E (k ) |
101  // r 4pi ___ | | 2| 2 1 |
102  // | / 3 |_ \2 - 2k / _|
103  // |/ Rr
104  double k_2 = k2(r, z);
105  double k = std::sqrt(k_2);
106  double constant =
107  scale * k * z / (4 * M_PI * std::sqrt(m_cfg.radius * r * r * r));
108 
109  double B = (2. - k_2) / (2. - 2. * k_2) * ellint_2(k_2) - ellint_1(k_2);
110 
111  // pos[0] is still signed!
112  return r / pos[0] * constant * B;
113 }
114 
115 double Acts::SolenoidBField::B_z(const Vector2D& pos, double scale) const {
116  // _
117  // 2 / pi / 2 2 2 - 1 / 2
118  // E (k ) = | ( 1 - k sin {theta} ) dtheta
119  // 1 _/ 0
120  using boost::math::ellint_1;
121  // _ ____________________
122  // 2 / pi / 2| / 2 2
123  // E (k ) = | |/ 1 - k sin {theta} dtheta
124  // 2 _/ 0
125  using boost::math::ellint_2;
126 
127  double r = std::abs(pos[0]);
128  double z = pos[1];
129 
130  // _ _
131  // mu I | / 2 \ |
132  // 0 k | | (R + r)k - 2r | 2 2 |
133  // B (r,z) = ----- ---- | | -------------- | E (k ) + E (k ) |
134  // z 4pi __ | | 2 | 2 1 |
135  // |/Rr |_ \ 2r(1 - k ) / _|
136 
137  if (r == 0) {
138  double res = scale / 2. * m_R2 / (std::sqrt(m_R2 + z * z) * (m_R2 + z * z));
139  return res;
140  }
141 
142  double k_2 = k2(r, z);
143  double k = std::sqrt(k_2);
144  double constant = scale * k / (4 * M_PI * std::sqrt(m_cfg.radius * r));
145  double B = ((m_cfg.radius + r) * k_2 - 2. * r) / (2. * r * (1. - k_2)) *
146  ellint_2(k_2) +
147  ellint_1(k_2);
148 
149  return constant * B;
150 }
151 
152 double Acts::SolenoidBField::k2(double r, double z) const {
153  // 2 4Rr
154  // k = ---------------
155  // 2 2
156  // (R + r) + z
157  return 4 * m_cfg.radius * r /
158  ((m_cfg.radius + r) * (m_cfg.radius + r) + z * z);
159 }