Belle II Software light-2406-ragdoll
BFieldComponentRadial.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <geometry/bfieldmap/BFieldComponentRadial.h>
10#include <geometry/bfieldmap/BFieldComponentBeamline.h>
11
12#include <framework/utilities/FileSystem.h>
13#include <framework/logging/Logger.h>
14
15#include <boost/iostreams/filtering_stream.hpp>
16#include <boost/iostreams/device/file.hpp>
17#include <boost/iostreams/filter/gzip.hpp>
18
19#include <cmath>
20
21using namespace std;
22using namespace Belle2;
23namespace io = boost::iostreams;
24
26{
27 if (m_mapFilename.empty()) {
28 B2ERROR("The filename for the radial magnetic field component is empty !");
29 return;
30 }
31
32 string fullPath = FileSystem::findFile("/data/" + m_mapFilename);
33
34 if (!FileSystem::fileExists(fullPath)) {
35 B2ERROR("The radial magnetic field map file '" << m_mapFilename << "' could not be found !");
36 return;
37 }
38
39 //Load the map file
40 io::filtering_istream fieldMapFile;
41 fieldMapFile.push(io::gzip_decompressor());
42 fieldMapFile.push(io::file_source(fullPath));
43
44 //Create the magnetic field map [r,z] and read the data from the file
45 B2DEBUG(10, "Loading the radial magnetic field from file '" << m_mapFilename << "' in to the memory...");
46 m_mapBuffer.clear();
47 m_mapBuffer.reserve(m_mapSize[0]*m_mapSize[1]);
48
49 double r, z, Br, Bz;
50 for (int i = 0; i < m_mapSize[0]; ++i) {
51 for (int j = 0; j < m_mapSize[1]; j++) {
52 fieldMapFile >> r >> z >> Br >> Bz;
53 //Store the values
54 m_mapBuffer.push_back({Br, Bz});
55 }
56 }
57
61 m_iLayer = 1 / m_Layer;
62
63 B2DEBUG(10, "... loaded " << m_mapSize[0] << "x" << m_mapSize[1] << " (r,z) elements.");
64}
65
66namespace Belle2 {
73 {
74 return {v.r * a, v.z * a};
75 }
76
80 {
81 return {u.r + v.r, u.z + v.z};
82 }
84}
85
86ROOT::Math::XYZVector BFieldComponentRadial::calculate(const ROOT::Math::XYZVector& point) const
87{
88 ROOT::Math::XYZVector B;
89 // If both 'Radial' and 'Beamline' components are defined in xml file,
90 // '3d' component returns zero field where 'Beamline' component is defined.
91 // If no 'Beamline' component is defined in xml file, the following function will never be called.
92 if (BFieldComponentBeamline::Instance().isInRange(point)) return B;
93
94 //Get the z component
95 double z = point.Z();
96 //Check if the point lies inside the magnetic field boundaries
97 if ((z < m_mapRegionZ[0]) || (z > m_mapRegionZ[1])) return B;
98
99 //Get the r component
100 double r2 = point.Perp2();
101 //Check if the point lies inside the magnetic field boundaries
102 if ((r2 < m_mapRegionR[0]*m_mapRegionR[0]) || (r2 > m_mapRegionR[1]*m_mapRegionR[1])) return B;
103
104 double r = sqrt(r2);
105 //Calculate the distance to the lower grid point
106 double wr = (r - m_mapRegionR[0]) * m_igridPitchR;
107 double wz = (z - m_mapRegionZ[0]) * m_igridPitchZ;
108 //Calculate the lower index of the point in the grid
109 int ir = static_cast<int>(wr);
110 int iz = static_cast<int>(wz);
111
112 double za = z - m_mapOffset;
113 double dz_eklm = abs(za) - (m_endyokeZMin - m_gapHeight);
114 if (dz_eklm > 0 && r > m_slotRMin) {
115 double wl = dz_eklm * m_iLayer;
116 int il = static_cast<int>(wl);
117 if (il <= 16) {
118 double L = m_Layer * il;
119 double l = dz_eklm - L;
120 int ingap = l < m_gapHeight;
121 ir += ingap & (r < m_slotRMin + m_gridPitchR);
122
123 double zlg = L + m_endyokeZMin, zl0 = zlg - m_gapHeight, zl1 = zlg + m_ironPlateThickness;
124 if (za < 0) {
125 zl0 = -zl0;
126 zlg = -zlg;
127 zl1 = -zl1;
128 }
129 double zoff = m_mapOffset - m_mapRegionZ[0];
130 int izl0 = static_cast<int>((zl0 + zoff) * m_igridPitchZ);
131 int izlg = static_cast<int>((zlg + zoff) * m_igridPitchZ);
132 int izl1 = static_cast<int>((zl1 + zoff) * m_igridPitchZ);
133 int ig = izlg == iz;
134 int idz = (izl0 == iz) - (izl1 == iz) + (ingap ? -ig : ig);
135 iz += (za > 0) ? idz : -idz;
136 }
137 }
138
139 //Bring the index values within the range
140 ir = min(m_mapSize[0] - 2, ir);
141 iz = min(m_mapSize[1] - 2, iz);
142
143 wr -= ir;
144 wz -= iz;
145
146 double wz0 = 1 - wz, wr0 = 1 - wr;
147 //Calculate the linear approx. of the magnetic field vector
148 const unsigned int strideZ = m_mapSize[1];
149 const unsigned int i00 = ir * strideZ + iz, i01 = i00 + 1, i10 = i00 + strideZ, i11 = i01 + strideZ;
150 const BFieldPoint* H = m_mapBuffer.data();
151 double w00 = wz0 * wr0, w01 = wz * wr0, w10 = wz0 * wr, w11 = wz * wr;
152 BFieldPoint Brz = H[i00] * w00 + H[i01] * w01 + H[i10] * w10 + H[i11] * w11;
153
154 if (r > 0.0) {
155 double norm = Brz.r / r;
156 B.SetXYZ(point.X()*norm, point.Y()*norm, Brz.z);
157 } else {
158 B.SetZ(Brz.z);
159 }
160
161 return B;
162}
163
165{
166}
double m_iLayer
reciprocal of height of the layer (gap + iron plate) in endyoke
double m_slotRMin
minimum radius for the gap in endyoke
virtual void initialize() override
Initializes the magnetic field component.
double m_gridPitchZ
The grid pitch in z.
std::string m_mapFilename
The filename of the magnetic field map.
double m_Layer
height of the layer (gap + iron plate) in endyoke
virtual void terminate() override
Terminates the magnetic field component.
double m_mapRegionZ[2]
The min and max boundaries of the map region in z.
double m_mapOffset
Offset required because the accelerator group defines the Belle center as zero.
virtual ROOT::Math::XYZVector calculate(const ROOT::Math::XYZVector &point) const override
Calculates the magnetic field vector at the specified space point.
double m_ironPlateThickness
thickness of iron plate in endyoke
double m_gridPitchR
The grid pitch in r.
double m_mapRegionR[2]
The min and max boundaries of the map region in r.
double m_endyokeZMin
minimum Z of endyoke
std::vector< BFieldPoint > m_mapBuffer
The memory buffer for the magnetic field map.
double m_gapHeight
height of the gap in endyoke
double m_igridPitchZ
The reciprocal of grid pitch in z.
int m_mapSize[2]
The size of the map in r and z.
double m_igridPitchR
The reciprocal of grid pitch in r.
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
Definition: FileSystem.cc:151
static bool fileExists(const std::string &filename)
Check if the file with given filename exists.
Definition: FileSystem.cc:32
B2Vector3< DataType > operator*(DataType a, const B2Vector3< DataType > &p)
non-memberfunction Scaling of 3-vectors with a real number
Definition: B2Vector3.h:537
B2Vector3< DataType > operator+(const TVector3 &a, const B2Vector3< DataType > &b)
non-memberfunction for adding a TVector3 to a B2Vector3
Definition: B2Vector3.h:544
static BFieldComponentBeamline & Instance()
BFieldComponentBeamline instance.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
double r
Magnetic field in r direction.
double z
Magnetic field in z direction.