Belle II Software light-2406-ragdoll
BFieldComponentKlm1.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/BFieldComponentKlm1.h>
10
11#include <framework/utilities/FileSystem.h>
12#include <framework/logging/Logger.h>
13
14#include <boost/iostreams/filtering_stream.hpp>
15#include <boost/iostreams/device/file.hpp>
16
17#include <cmath>
18
19using namespace std;
20using namespace Belle2;
21namespace io = boost::iostreams;
22
24{
25 if (m_mapFilename.empty()) {
26 B2ERROR("The filename for the magnetic field Component is empty !");
27 return;
28 }
29
30 string fullPath = FileSystem::findFile("/data/" + m_mapFilename);
31
32 if (!FileSystem::fileExists(fullPath)) {
33 B2ERROR("The magnetic field map file '" << m_mapFilename << "' could not be found !");
34 return;
35 }
36
37 //Load the map file
38 io::filtering_istream fieldMapFile;
39 fieldMapFile.push(io::file_source(fullPath));
40
41 //read the data from the file
42 B2DEBUG(10, "Loading the magnetic field from file '" << m_mapFilename << "' in to the memory...");
43
44 char dummy[100];
45
46 //store parameters for Barrel
47 fieldMapFile.getline(dummy, sizeof(dummy));
48 for (int layer = 0; layer < m_nBarrelLayers; ++layer) {
49 int l;
50 fieldMapFile >> l >> m_barrelZBreakpoint[layer] >> m_barrelRBreakpoint[layer]
51 >> m_barrelFieldZSlope1[layer]
53 >> m_barrelFieldZSlope2[layer]
55 >> m_barrelFieldRSlope1[layer]
57 >> m_barrelFieldRSlope2[layer]
59
60 if (layer != l)
61 B2ERROR("Barrel Parameter file is something wrong ! " << layer << " " << l);
62 }
63
64 //store parameters for Endcap
65 fieldMapFile.getline(dummy, sizeof(dummy));
66 fieldMapFile.getline(dummy, sizeof(dummy));
67 for (int GapIron = 0; GapIron < 2; GapIron++) {
68 for (int layer = 0; layer < m_nEndcapLayers + 1; layer++) {
69 for (int j = 0; j < 5; j++) {
70 int g, l, jj;
71 fieldMapFile >> g >> l >> jj
72 >> m_endcapZBreakpoint[GapIron][layer][j] >> m_endcapRBreakpoint[GapIron][layer][j]
73 >> m_endcapFieldZSlope[GapIron][layer][j]
74 >> m_endcapFieldZIntercept[GapIron][layer][j]
75 >> m_endcapFieldRSlope[GapIron][layer][j]
76 >> m_endcapFieldRIntercept[GapIron][layer][j];
77 if (GapIron != g || layer != l || jj != j)
78 B2ERROR("Endcap Parameter file is something wrong ! " << g << " " << l << " " << jj);
79 }
80 }
81 }
82
83 B2DEBUG(10, "... loaded ");
84}
85
86
87ROOT::Math::XYZVector BFieldComponentKlm1::calculate(const ROOT::Math::XYZVector& point) const
88{
89 //Get the r and z Component
90 double r = point.Rho();
91 double x = point.X();
92 double y = point.Y();
93 double z = point.Z() - m_mapOffset;
94 double absZ = std::abs(z);
95
96 ROOT::Math::XYZVector bField(0., 0., 0.);
97
98 // Barrel
99 if (absZ < m_barrelZMax) {
100
101 // This code assumes that we are in the barrel KLM. Use the field map
102 // inside an iron plate; use zero in the gap between the iron plates.
103
104 double cosphi = std::abs(x) / r;
105 double d = ((cosphi < m_cospi8) ? ((cosphi < m_cos3pi8) ? std::abs(y)
106 : (std::abs(x) + std::abs(y)) * m_cospi4)
107 : std::abs(x)) - m_barrelRMin - m_barrelGapHeightLayer0;
108
109 if (d >= 0.0) {
110 int layer = static_cast<int>(floor(d / m_dLayer));
111
112 // make sure layer is in a valid range
113 if (layer < 0 || layer > m_nBarrelLayers) return bField;
114 if (layer == m_nBarrelLayers) {
115 layer = m_nBarrelLayers - 1;
116 d = layer * m_dLayer; // for extra-thick outermost iron plate
117 }
118 if ((d - layer * m_dLayer) < m_barrelIronThickness) {
119 double br = ((absZ > m_barrelRBreakpoint[layer]) ?
120 m_barrelFieldRIntercept2[layer] + m_barrelFieldRSlope2[layer] * absZ :
121 m_barrelFieldRIntercept1[layer] + m_barrelFieldRSlope1[layer] * absZ);
122 if (z < 0.0) { br = -br; }
123 bField.SetXYZ(br * x / r, br * y / r,
124 ((absZ > m_barrelZBreakpoint[layer]) ?
125 m_barrelFieldZIntercept2[layer] + m_barrelFieldZSlope2[layer] * absZ :
126 m_barrelFieldZIntercept1[layer] + m_barrelFieldZSlope1[layer] * absZ));
127 }
128 }
129 } else if (absZ > m_endcapZMin) {
130 //Endcap
131 if (r <= m_endcapRMin) return bField;
132
133 // This code assumes that we are in the endcap KLM.
134 // The field in an inter-plate gap in the endcap is not zero.
135
136 double dz = absZ - (m_endcapZMin - m_endcapGapHeight);
137 int layer = static_cast<int>(floor(dz / m_dLayer));
138 // make sure layer is in a valid range
139 if (layer < 0 || layer >= m_nEndcapLayers) return bField;
140
141 int GapIron = 0; // in gap
142 if ((dz - m_dLayer * layer) > m_endcapGapHeight) GapIron = 1; // in Iron plate?
143 for (int j = 0; j < 5; j++) {
144 if (r < m_endcapRBreakpoint[GapIron][layer][j]) {
145 double br = m_endcapFieldRIntercept[GapIron][layer][j] + m_endcapFieldRSlope[GapIron][layer][j] * r;
146 if (z < 0.0) { br = -br; }
147 bField.SetX(br * x / r);
148 bField.SetY(br * y / r);
149 break;
150 }
151 }
152 for (int j = 0; j < 5; j++) {
153 if (r < m_endcapZBreakpoint[GapIron][layer][j]) {
154 bField.SetZ(m_endcapFieldZIntercept[GapIron][layer][j] + m_endcapFieldZSlope[GapIron][layer][j] * r);
155 break;
156 }
157 }
158 }
159
160 return bField;
161}
162
163
165{
166}
double m_endcapFieldZSlope[2][15][5]
Slopes of the linear approximation of Bz in the endcap.
double m_endcapGapHeight
Gap height of BKLM layer1-14.
double m_barrelFieldRIntercept1[15]
Intercept of Br before the beackpoint in the barrel.
int m_nBarrelLayers
The number of layers per 1 sector for BKLM.
virtual void initialize() override
Initializes the magnetic field Component.
double m_barrelFieldZSlope1[15]
Slope of Bz before the breakpoint in the barrel.
double m_endcapFieldRIntercept[2][15][5]
Intercepts of the linear approximation of Br in the endcap.
double m_barrelZMax
The maximum boundaries of BKLM region in r.
double m_barrelIronThickness
Thickness of Barrel iron plate.
double m_barrelFieldZIntercept1[15]
Intercept of Bz before the beackpoint in the barrel.
std::string m_mapFilename
The filename of the magnetic field map.
virtual void terminate() override
Terminates the magnetic field Component.
double m_mapOffset
Offset required because the accelerator group defines the Belle center as zero.
double m_endcapRMin
The minimum boundaries of EKLM region in r.
virtual ROOT::Math::XYZVector calculate(const ROOT::Math::XYZVector &point) const override
Calculates the magnetic field vector at the specified space point.
double m_barrelFieldRSlope2[15]
Slope of Br after the breakpoint in the barrel.
double m_endcapRBreakpoint[2][15][5]
r position of breakpoints between linear functions in the endcap First index indicates whether or not...
double m_barrelRBreakpoint[15]
z position of breakpoints between the two linear approximations of Br in the barrel.
double m_barrelGapHeightLayer0
Gap height of BKLM layer0.
double m_barrelRMin
The minimum boundaries of BKLM region in r.
double m_endcapFieldRSlope[2][15][5]
Slopes of the linear approximation of Br in the endcap.
int m_nEndcapLayers
The number of layers per 1 sector for EKLM.
double m_barrelFieldZIntercept2[15]
Intercept of Bz after the breakpoint in the barrel.
double m_barrelZBreakpoint[15]
z position of breakpoints between the two linear approximations of Bz in the barrel.
double m_barrelFieldZSlope2[15]
Slope of Bz after the breakpoint in the barrel.
double m_barrelFieldRSlope1[15]
Slope of Br before the breakpoint in the barrel.
double m_endcapZBreakpoint[2][15][5]
z position of breakpoints between linear functions in the endcap.
double m_endcapFieldZIntercept[2][15][5]
Intercepts of the linear approximation of Bz in the endcap.
double m_endcapZMin
The minimum boundaries of EKLM region in z.
double m_barrelFieldRIntercept2[15]
Intercept of Br after the beackpoint in the barrel.
double m_dLayer
depth of BKLM module?
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
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.