Belle II Software  release-05-01-25
BFieldComponentRadial.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll, Hiroyuki Nakayama *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <geometry/bfieldmap/BFieldComponentRadial.h>
12 #include <geometry/bfieldmap/BFieldComponentBeamline.h>
13 
14 #include <framework/utilities/FileSystem.h>
15 #include <framework/logging/Logger.h>
16 
17 #include <boost/iostreams/filtering_stream.hpp>
18 #include <boost/iostreams/device/file.hpp>
19 #include <boost/iostreams/filter/gzip.hpp>
20 
21 #include <cmath>
22 
23 using namespace std;
24 using namespace Belle2;
25 namespace io = boost::iostreams;
26 
27 void BFieldComponentRadial::initialize()
28 {
29  if (m_mapFilename.empty()) {
30  B2ERROR("The filename for the radial magnetic field component is empty !");
31  return;
32  }
33 
34  string fullPath = FileSystem::findFile("/data/" + m_mapFilename);
35 
36  if (!FileSystem::fileExists(fullPath)) {
37  B2ERROR("The radial magnetic field map file '" << m_mapFilename << "' could not be found !");
38  return;
39  }
40 
41  //Load the map file
42  io::filtering_istream fieldMapFile;
43  fieldMapFile.push(io::gzip_decompressor());
44  fieldMapFile.push(io::file_source(fullPath));
45 
46  //Create the magnetic field map [r,z] and read the data from the file
47  B2DEBUG(10, "Loading the radial magnetic field from file '" << m_mapFilename << "' in to the memory...");
48  m_mapBuffer.clear();
49  m_mapBuffer.reserve(m_mapSize[0]*m_mapSize[1]);
50 
51  double r, z, Br, Bz;
52  for (int i = 0; i < m_mapSize[0]; ++i) {
53  for (int j = 0; j < m_mapSize[1]; j++) {
54  fieldMapFile >> r >> z >> Br >> Bz;
55  //Store the values
56  m_mapBuffer.push_back({Br, Bz});
57  }
58  }
59 
60  m_igridPitchR = 1 / m_gridPitchR;
61  m_igridPitchZ = 1 / m_gridPitchZ;
62  m_Layer = m_gapHeight + m_ironPlateThickness;
63  m_iLayer = 1 / m_Layer;
64 
65  B2DEBUG(10, "... loaded " << m_mapSize[0] << "x" << m_mapSize[1] << " (r,z) elements.");
66 }
67 
68 namespace Belle2 {
75  {
76  return {v.r * a, v.z * a};
77  }
78 
82  {
83  return {u.r + v.r, u.z + v.z};
84  }
86 }
87 
88 B2Vector3D BFieldComponentRadial::calculate(const B2Vector3D& point) const
89 {
90  B2Vector3D B;
91  // If both 'Radial' and 'Beamline' components are defined in xml file,
92  // '3d' component returns zero field where 'Beamline' component is defined.
93  // If no 'Beamline' component is defined in xml file, the following function will never be called.
94  if (BFieldComponentBeamline::Instance().isInRange(point)) return B;
95 
96  //Get the z component
97  double z = point.Z();
98  //Check if the point lies inside the magnetic field boundaries
99  if ((z < m_mapRegionZ[0]) || (z > m_mapRegionZ[1])) return B;
100 
101  //Get the r component
102  double r2 = point.Perp2();
103  //Check if the point lies inside the magnetic field boundaries
104  if ((r2 < m_mapRegionR[0]*m_mapRegionR[0]) || (r2 > m_mapRegionR[1]*m_mapRegionR[1])) return B;
105 
106  double r = sqrt(r2);
107  //Calculate the distance to the lower grid point
108  double wr = (r - m_mapRegionR[0]) * m_igridPitchR;
109  double wz = (z - m_mapRegionZ[0]) * m_igridPitchZ;
110  //Calculate the lower index of the point in the grid
111  int ir = static_cast<int>(wr);
112  int iz = static_cast<int>(wz);
113 
114  double za = z - m_mapOffset;
115  double dz_eklm = abs(za) - (m_endyokeZMin - m_gapHeight);
116  if (dz_eklm > 0 && r > m_slotRMin) {
117  double wl = dz_eklm * m_iLayer;
118  int il = static_cast<int>(wl);
119  if (il <= 16) {
120  double L = m_Layer * il;
121  double l = dz_eklm - L;
122  int ingap = l < m_gapHeight;
123  ir += ingap & (r < m_slotRMin + m_gridPitchR);
124 
125  double zlg = L + m_endyokeZMin, zl0 = zlg - m_gapHeight, zl1 = zlg + m_ironPlateThickness;
126  if (za < 0) {
127  zl0 = -zl0;
128  zlg = -zlg;
129  zl1 = -zl1;
130  }
131  double zoff = m_mapOffset - m_mapRegionZ[0];
132  int izl0 = static_cast<int>((zl0 + zoff) * m_igridPitchZ);
133  int izlg = static_cast<int>((zlg + zoff) * m_igridPitchZ);
134  int izl1 = static_cast<int>((zl1 + zoff) * m_igridPitchZ);
135  int ig = izlg == iz;
136  int idz = (izl0 == iz) - (izl1 == iz) + (ingap ? -ig : ig);
137  iz += (za > 0) ? idz : -idz;
138  }
139  }
140 
141  //Bring the index values within the range
142  ir = min(m_mapSize[0] - 2, ir);
143  iz = min(m_mapSize[1] - 2, iz);
144 
145  wr -= ir;
146  wz -= iz;
147 
148  double wz0 = 1 - wz, wr0 = 1 - wr;
149  //Calculate the linear approx. of the magnetic field vector
150  const unsigned int strideZ = m_mapSize[1];
151  const unsigned int i00 = ir * strideZ + iz, i01 = i00 + 1, i10 = i00 + strideZ, i11 = i01 + strideZ;
152  const BFieldPoint* H = m_mapBuffer.data();
153  double w00 = wz0 * wr0, w01 = wz * wr0, w10 = wz0 * wr, w11 = wz * wr;
154  BFieldPoint Brz = H[i00] * w00 + H[i01] * w01 + H[i10] * w10 + H[i11] * w11;
155 
156  if (r > 0.0) {
157  double norm = Brz.r / r;
158  B.SetXYZ(point.X()*norm, point.Y()*norm, Brz.z);
159  } else {
160  B.SetZ(Brz.z);
161  }
162 
163  return B;
164 }
165 
166 void BFieldComponentRadial::terminate()
167 {
168 }
Belle2::BFieldComponentRadial::BFieldPoint::r
double r
Magnetic field in r direction.
Definition: BFieldComponentRadial.h:50
Belle2::B2Vector3::Z
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:434
Belle2::B2Vector3< double >
Belle2::BFieldComponentRadial::BFieldPoint
Magnetic field data structure.
Definition: BFieldComponentRadial.h:49
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::operator*
BFieldComponentRadial::BFieldPoint operator*(const BFieldComponentRadial::BFieldPoint &v, double a)
multiply a radial bfield point by a real number
Definition: BFieldComponentRadial.cc:74
Belle2::B2Vector3::Perp2
DataType Perp2() const
The transverse component squared (R^2 in cylindrical coordinate system).
Definition: B2Vector3.h:195
Belle2::operator+
BFieldComponentRadial::BFieldPoint operator+(const BFieldComponentRadial::BFieldPoint &u, const BFieldComponentRadial::BFieldPoint &v)
Add two radial bfield points together.
Definition: BFieldComponentRadial.cc:80
Belle2::B2Vector3::X
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:430
Belle2::BFieldComponentRadial::BFieldPoint::z
double z
Magnetic field in z direction.
Definition: BFieldComponentRadial.h:51
Belle2::B2Vector3::Y
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:432