Belle II Software development
GeoMagneticField.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/GeoMagneticField.h>
10#include <geometry/bfieldmap/BFieldMap.h>
11
12#include <geometry/bfieldmap/BFieldComponentConstant.h>
13#include <geometry/bfieldmap/BFieldComponentRadial.h>
14#include <geometry/bfieldmap/BFieldComponentQuad.h>
15#include <geometry/bfieldmap/BFieldComponentBeamline.h>
16#include <geometry/bfieldmap/BFieldComponentKlm1.h>
17#include <geometry/bfieldmap/BFieldComponent3d.h>
18#include <geometry/dbobjects/MagneticFieldComponent3D.h>
19#include <geometry/CreatorFactory.h>
20#include <framework/database/DBImportObjPtr.h>
21#include <framework/dbobjects/MagneticFieldComponentConstant.h>
22
23#include <framework/logging/Logger.h>
24#include <framework/gearbox/GearDir.h>
25#include <framework/gearbox/Unit.h>
26#include <framework/utilities/FileSystem.h>
27
28#include <boost/iostreams/filtering_stream.hpp>
29#include <boost/iostreams/device/file.hpp>
30#include <boost/iostreams/filter/gzip.hpp>
31#include <boost/bind/bind.hpp>
32
33using namespace std;
34using namespace Belle2;
35
36//-----------------------------------------------------------------
37// Register the Creator
38//-----------------------------------------------------------------
39
40geometry::CreatorFactory<GeoMagneticField> GeoMagneticFieldFactory("GeoMagneticField");
41
42//-----------------------------------------------------------------
43// Implementation
44//-----------------------------------------------------------------
45
46GeoMagneticField::GeoMagneticField() : CreatorBase()
47{
48 //Add the function pointers called for reading the components to the map
49 m_componentTypeMap.insert(make_pair("Constant", boost::bind(&GeoMagneticField::readConstantBField, this, boost::placeholders::_1)));
50 m_componentTypeMap.insert(make_pair("Radial", boost::bind(&GeoMagneticField::readRadialBField, this, boost::placeholders::_1)));
51 m_componentTypeMap.insert(make_pair("Quad", boost::bind(&GeoMagneticField::readQuadBField, this, boost::placeholders::_1)));
52 m_componentTypeMap.insert(make_pair("Beamline", boost::bind(&GeoMagneticField::readBeamlineBField, this, boost::placeholders::_1)));
53 m_componentTypeMap.insert(make_pair("Klm1", boost::bind(&GeoMagneticField::readKlm1BField, this, boost::placeholders::_1)));
54 m_componentTypeMap.insert(make_pair("3d", boost::bind(&GeoMagneticField::read3dBField, this, boost::placeholders::_1)));
55}
56
57
58GeoMagneticField::~GeoMagneticField() = default;
59
60MagneticField GeoMagneticField::createConfiguration(const GearDir& content)
61{
62 MagneticField fieldmap;
63 //Read the magnetic field components
64 for (const GearDir& component : content.getNodes("Components/Component")) {
65 //Get the type of the magnetic field and call the appropriate function
66 string compType = component.getString("attribute::type");
67 B2DEBUG(10, "GeoMagneticField creator: Loading the parameters for the component type'" << compType << "'");
68 if (compType == "3d") add3dBField(component, fieldmap);
69 else if (compType == "Constant") addConstantBField(component, fieldmap);
70 else B2ERROR("The magnetic field type " << compType << " can not yet be stored in the database");
71 }
72 return fieldmap;
73}
74
75void GeoMagneticField::addConstantBField(const GearDir& component, MagneticField& fieldmap)
76{
77 double xValue = component.getWithUnit("X");
78 double yValue = component.getWithUnit("Y");
79 double zValue = component.getWithUnit("Z");
80 double RminValue = component.getLength("MinR", 0); // stored in cm
81 double RmaxValue = component.getLength("MaxR"); // stored in cm
82 double ZminValue = component.getLength("MinZ"); // stored in cm
83 double ZmaxValue = component.getLength("MaxZ"); // stored in cm
84 auto field = new MagneticFieldComponentConstant(ROOT::Math::XYZVector(xValue, yValue, zValue),
85 RminValue, RmaxValue, ZminValue, ZmaxValue);
86 fieldmap.addComponent(field);
87}
88
89void GeoMagneticField::add3dBField(const GearDir& component, MagneticField& fieldmap)
90{
91 string mapFilename = component.getString("MapFilename");
92 int mapSizeR = component.getInt("NumberGridPointsR");
93 int mapSizeZ = component.getInt("NumberGridPointsZ");
94 int mapSizePhi = component.getInt("NumberGridPointsPhi");
95
96 double minZ = component.getLength("ZMin");
97 double maxZ = component.getLength("ZMax");
98 double minR = component.getLength("RadiusMin");
99 double maxR = component.getLength("RadiusMax");
100
101 const string fullPath = FileSystem::findFile("/data/" + mapFilename);
102 if (!FileSystem::fileExists(fullPath)) {
103 B2ERROR("The 3d magnetic field map file '" << mapFilename << "' could not be found !");
104 return;
105 }
106
107 std::vector<ROOT::Math::XYZVector> bmap;
108 bmap.reserve(mapSizeR * mapSizeZ * mapSizePhi);
109 // Load B-field map file
110 boost::iostreams::filtering_istream fieldMapFile;
111 fieldMapFile.push(boost::iostreams::gzip_decompressor());
112 fieldMapFile.push(boost::iostreams::file_source(fullPath));
113
114 char tmp[256];
115 for (int k = 0; k < mapSizeZ; k++) { // z
116 for (int i = 0; i < mapSizeR; i++) { // r
117 for (int j = 0; j < mapSizePhi; j++) { // phi
118 double Br, Bz, Bphi;
119 //r[m] phi[deg] z[m] Br[T] Bphi[T] Bz[T]
120 fieldMapFile.getline(tmp, 256);
121 // sscanf(tmp+33,"%lf %lf %lf",&Br,&Bphi,&Bz);
122 char* next;
123 Br = strtod(tmp + 33, &next);
124 Bphi = strtod(next, &next);
125 Bz = strtod(next, nullptr);
126 bmap.emplace_back(-Br * Unit::T, -Bphi * Unit::T, -Bz * Unit::T);
127 }
128 }
129 }
130 auto* field = new MagneticFieldComponent3D(
131 minR, maxR, minZ, maxZ,
132 mapSizeR, mapSizePhi, mapSizeZ,
133 std::move(bmap)
134 );
135 fieldmap.addComponent(field);
136}
137
138void GeoMagneticField::createPayloads(const GearDir& content, const IntervalOfValidity& iov)
139{
140 DBImportObjPtr<MagneticField> importObj;
141 importObj.construct(createConfiguration(content));
142 importObj.import(iov);
143}
144
145void GeoMagneticField::create(const GearDir& content, G4LogicalVolume& /*topVolume*/, geometry::GeometryTypes)
146{
147 // clear any existing BField
148 BFieldMap::Instance().clear();
149
150 //Loop over all components of the magnetic field
151 CompTypeMap::iterator findIter;
152 for (const GearDir& component : content.getNodes("Components/Component")) {
153 //Get the type of the magnetic field and call the appropriate function
154 string compType = component.getString("attribute::type");
155 B2DEBUG(10, "GeoMagneticField creator: Loading the parameters for the component type'" << compType << "'");
156
157 findIter = m_componentTypeMap.find(compType);
158 if (findIter != m_componentTypeMap.end()) {
159 findIter->second(component);
160 } else {
161 B2ERROR("The magnetic field component type '" << compType << "' is unknown !");
162 }
163 }
164
165 BFieldMap::Instance().initialize();
166}
167
168
169//============================================================
170// Protected methods
171//============================================================
172
173void GeoMagneticField::readConstantBField(const GearDir& component)
174{
175 double xValue = component.getDouble("X");
176 double yValue = component.getDouble("Y");
177 double zValue = component.getDouble("Z");
178 double RmaxValue = component.getLength("MaxR"); // stored in cm
179 double ZminValue = component.getLength("MinZ"); // stored in cm
180 double ZmaxValue = component.getLength("MaxZ"); // stored in cm
181
182 auto& bComp = BFieldMap::Instance().addBFieldComponent<BFieldComponentConstant>();
183 bComp.setMagneticFieldValues(xValue, yValue, zValue, RmaxValue, ZminValue, ZmaxValue);
184}
185
186
187void GeoMagneticField::readRadialBField(const GearDir& component)
188{
189 string mapFilename = component.getString("MapFilename");
190 int mapSizeR = component.getInt("NumberGridPointsR");
191 int mapSizeZ = component.getInt("NumberGridPointsZ");
192
193 double mapRegionMinZ = component.getLength("ZMin");
194 double mapRegionMaxZ = component.getLength("ZMax");
195 double mapOffset = component.getLength("ZOffset");
196
197 double mapRegionMinR = component.getLength("RadiusMin");
198 double mapRegionMaxR = component.getLength("RadiusMax");
199
200 double gridPitchR = component.getLength("GridPitchR");
201 double gridPitchZ = component.getLength("GridPitchZ");
202
203 double slotRMin = component.getLength("SlotRMin");
204 double endyokeZMin = component.getLength("EndyokeZMin");
205 double gapHeight = component.getLength("GapHeight");
206 double ironThickness = component.getLength("IronPlateThickness");
207
208 auto& bComp = BFieldMap::Instance().addBFieldComponent<BFieldComponentRadial>();
209 bComp.setMapFilename(mapFilename);
210 bComp.setMapSize(mapSizeR, mapSizeZ);
211 bComp.setMapRegionZ(mapRegionMinZ, mapRegionMaxZ, mapOffset);
212 bComp.setMapRegionR(mapRegionMinR, mapRegionMaxR);
213 bComp.setGridPitch(gridPitchR, gridPitchZ);
214 bComp.setKlmParameters(slotRMin, endyokeZMin, gapHeight, ironThickness);
215}
216
217void GeoMagneticField::readQuadBField(const GearDir& component)
218{
219 string mapFilenameHER = component.getString("MapFilenameHER");
220 string mapFilenameLER = component.getString("MapFilenameLER");
221 string mapFilenameHERleak = component.getString("MapFilenameHERleak");
222 string apertFilenameHER = component.getString("ApertFilenameHER");
223 string apertFilenameLER = component.getString("ApertFilenameLER");
224
225 int mapSizeHER = component.getInt("MapSizeHER");
226 int mapSizeLER = component.getInt("MapSizeLER");
227 int mapSizeHERleak = component.getInt("MapSizeHERleak");
228 int apertSizeHER = component.getInt("ApertSizeHER");
229 int apertSizeLER = component.getInt("ApertSizeLER");
230
231 auto& bComp = BFieldMap::Instance().addBFieldComponent<BFieldComponentQuad>();
232 bComp.setMapFilename(mapFilenameHER, mapFilenameLER, mapFilenameHERleak);
233 bComp.setApertFilename(apertFilenameHER, apertFilenameLER);
234 bComp.setMapSize(mapSizeHER, mapSizeLER, mapSizeHERleak);
235 bComp.setApertSize(apertSizeHER, apertSizeLER);
236}
237
238void GeoMagneticField::readBeamlineBField(const GearDir& component)
239{
240 string mapFilenameHER = component.getString("MapFilenameHER");
241 string mapFilenameLER = component.getString("MapFilenameLER");
242 string interFilenameHER = component.getString("InterFilenameHER");
243 string interFilenameLER = component.getString("InterFilenameLER");
244
245 double mapRegionMinZ = component.getLength("ZMin");
246 double mapRegionMaxZ = component.getLength("ZMax");
247
248 double mapRegionMinR = component.getLength("RadiusMin");
249 double mapRegionMaxR = component.getLength("RadiusMax");
250
251 double beamAngle = component.getLength("BeamAngle");
252
253 auto& bComp = BFieldMap::Instance().addBFieldComponent<BFieldComponentBeamline>();
254 bComp.setMapFilename(mapFilenameHER, mapFilenameLER);
255 bComp.setInterpolateFilename(interFilenameHER, interFilenameLER);
256 bComp.setMapRegionZ(mapRegionMinZ, mapRegionMaxZ);
257 bComp.setMapRegionR(mapRegionMinR, mapRegionMaxR);
258 bComp.setBeamAngle(beamAngle);
259}
260
261void GeoMagneticField::readKlm1BField(const GearDir& component)
262{
263 string mapFilename = component.getString("MapFilename");
264
265 int nBarrelLayers = component.getInt("NumberBarrelLayers");
266 int nEndcapLayers = component.getInt("NumberEndcapLayers");
267
268 double barrelRMin = component.getLength("BarrelRadiusMin");
269 double barrelZMax = component.getLength("BarrelZMax");
270 double mapOffset = component.getLength("ZOffset");
271
272 double endcapRMin = component.getLength("EdncapRadiusMin");
273 double endcapZMin = component.getLength("EndcapZMin");
274
275 double barrelGapHeightLayer0 = component.getLength("BarrelGapHeightLayer0");
276 double barrelIronThickness = component.getLength("BarrelIronThickness");
277 double endcapGapHeight = component.getLength("EndcapGapHeight");
278 double dLayer = component.getLength("DLayer");
279
280
281 auto& bComp = BFieldMap::Instance().addBFieldComponent<BFieldComponentKlm1>();
282 bComp.setMapFilename(mapFilename);
283 bComp.setNLayers(nBarrelLayers, nEndcapLayers);
284 bComp.setBarrelRegion(barrelRMin, barrelZMax, mapOffset);
285 bComp.setEndcapRegion(endcapRMin, endcapZMin);
286 bComp.setLayerParam(barrelGapHeightLayer0, barrelIronThickness, endcapGapHeight, dLayer);
287}
288
289void GeoMagneticField::read3dBField(const GearDir& component)
290{
291 string mapFilename = component.getString("MapFilename");
292 int mapSizeR = component.getInt("NumberGridPointsR");
293 int mapSizeZ = component.getInt("NumberGridPointsZ");
294 int mapSizePhi = component.getInt("NumberGridPointsPhi");
295
296 double mapRegionMinZ = component.getLength("ZMin");
297 double mapRegionMaxZ = component.getLength("ZMax");
298 double mapOffset = component.getLength("ZOffset");
299
300 double mapRegionMinR = component.getLength("RadiusMin");
301 double mapRegionMaxR = component.getLength("RadiusMax");
302
303 double gridPitchR = component.getLength("GridPitchR");
304 double gridPitchZ = component.getLength("GridPitchZ");
305 double gridPitchPhi = component.getLength("GridPitchPhi");
306 bool mirrorPhi = bool(component.getInt("MirrorPhi"));
307
308 double excludeRegionMinZ = component.getLength("ExcludeZMin");
309 double excludeRegionMaxZ = component.getLength("ExcludeZMax");
310 double excludeRegionMinR = component.getLength("ExcludeRadiusMin");
311 double excludeRegionMaxR = component.getLength("ExcludeRadiusMax");
312
313 double errorRegionMinR = component.getLength("BiasRadiusMin");
314 double errorRegionMaxR = component.getLength("BiasRadiusMax");
315 double errorBr = component.getDouble("BiasBr");
316 double errorBphi = component.getDouble("BiasBphi");
317 double errorBz = component.getDouble("BiasBz");
318
319 bool doInterpolation = bool(component.getInt("enableInterpolation"));
320 string enableCoordinate = component.getString("EnableCoordinate");
321
322 auto& bComp = BFieldMap::Instance().addBFieldComponent<BFieldComponent3d>();
323 bComp.setMapFilename(mapFilename);
324 bComp.setMapSize(mapSizeR, mapSizePhi, mapSizeZ);
325 bComp.setMapRegionZ(mapRegionMinZ, mapRegionMaxZ, mapOffset);
326 bComp.setMapRegionR(mapRegionMinR, mapRegionMaxR);
327 bComp.setGridPitch(gridPitchR, gridPitchPhi, gridPitchZ);
328 bComp.setExcludeRegionZ(excludeRegionMinZ, excludeRegionMaxZ);
329 bComp.setExcludeRegionR(excludeRegionMinR, excludeRegionMaxR);
330 bComp.setErrorRegionR(errorRegionMinR, errorRegionMaxR, errorBr, errorBphi, errorBz);
331 bComp.mirrorPhi(mirrorPhi);
332 bComp.doInterpolation(doInterpolation);
333 bComp.enableCoordinate(enableCoordinate);
334
335}
Belle2::BFieldComponent3d
The BFieldComponent3d class.
Definition: BFieldComponent3d.h:36
Belle2::BFieldComponent3d::setMapFilename
void setMapFilename(const std::string &filename)
Sets the filename of the magnetic field map.
Definition: BFieldComponent3d.h:71
Belle2::BFieldComponentBeamline
The BFieldComponentBeamline class.
Definition: BFieldComponentBeamline.h:31
Belle2::BFieldComponentBeamline::setMapFilename
void setMapFilename(const std::string &filename_her, const std::string &filename_ler)
Sets the filename of the magnetic field map.
Definition: BFieldComponentBeamline.h:83
Belle2::BFieldComponentConstant
The BFieldComponentConstant class.
Definition: BFieldComponentConstant.h:24
Belle2::BFieldComponentConstant::setMagneticFieldValues
void setMagneticFieldValues(double x, double y, double z, double rmax, double zmin, double zmax)
Sets the values for the homogeneous magnetic field vector.
Definition: BFieldComponentConstant.cc:22
Belle2::BFieldComponentKlm1
The Bfieldcomponentklm1 class.
Definition: BFieldComponentKlm1.h:33
Belle2::BFieldComponentKlm1::setMapFilename
void setMapFilename(const std::string &filename)
Sets the filename of the magnetic field map.
Definition: BFieldComponentKlm1.h:74
Belle2::BFieldComponentQuad
The BFieldComponentQuad class.
Definition: BFieldComponentQuad.h:35
Belle2::BFieldComponentQuad::setMapFilename
void setMapFilename(const std::string &filenameHER, const std::string &filenameLER, const std::string &filenameHERleak)
Sets the filename of the magnetic field map.
Definition: BFieldComponentQuad.h:146
Belle2::BFieldComponentRadial
The BFieldComponentRadial class.
Definition: BFieldComponentRadial.h:34
Belle2::BFieldComponentRadial::setMapFilename
void setMapFilename(const std::string &filename)
Sets the filename of the magnetic field map.
Definition: BFieldComponentRadial.h:74
Belle2::BFieldMap::initialize
void initialize()
Initialize the magnetic field after adding all components.
Definition: BFieldMap.cc:21
Belle2::BFieldMap::clear
void clear()
Clear the existing components.
Definition: BFieldMap.cc:33
Belle2::BFieldMap::Instance
static BFieldMap & Instance()
Static method to get a reference to the BFieldMap instance.
Definition: BFieldMap.cc:15
Belle2::DBImportBase::import
bool import(const IntervalOfValidity &iov)
Import the object to database.
Definition: DBImportBase.cc:36
Belle2::DBImportObjPtr
Class for importing a single object to the database.
Definition: DBImportObjPtr.h:23
Belle2::DBImportObjPtr::construct
void construct(Args &&... params)
Construct an object of type T in this DBImportObjPtr using the provided constructor arguments.
Definition: DBImportObjPtr.h:47
Belle2::FileSystem::findFile
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
Belle2::FileSystem::fileExists
static bool fileExists(const std::string &filename)
Check if the file with given filename exists.
Definition: FileSystem.cc:32
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:31
Belle2::GeoMagneticField::readConstantBField
void readConstantBField(const GearDir &component)
Reads the parameters for a homogeneous magnetic field and adds the component to the global magnetic f...
Definition: GeoMagneticField.cc:173
Belle2::GeoMagneticField::add3dBField
void add3dBField(const GearDir &component, MagneticField &fielmap)
Add a 3D field component to a magnetic field configuration for the DB.
Definition: GeoMagneticField.cc:89
Belle2::GeoMagneticField::readBeamlineBField
void readBeamlineBField(const GearDir &component)
Reads the 3D Bfield map and parameters near beam pipes and adds the component to the global magnetic ...
Definition: GeoMagneticField.cc:238
Belle2::GeoMagneticField::readKlm1BField
void readKlm1BField(const GearDir &component)
Reads the 2D Bfield map and parameters outside of solenoid and adds the component to the global magne...
Definition: GeoMagneticField.cc:261
Belle2::GeoMagneticField::m_componentTypeMap
CompTypeMap m_componentTypeMap
Maps the name of the component to the function reading the parameters.
Definition: GeoMagneticField.h:70
Belle2::GeoMagneticField::~GeoMagneticField
virtual ~GeoMagneticField()
The destructor of the GeoMagneticField class.
Belle2::GeoMagneticField::readRadialBField
void readRadialBField(const GearDir &component)
Reads the parameters for a radial magnetic field and adds the component to the global magnetic field.
Definition: GeoMagneticField.cc:187
Belle2::GeoMagneticField::createPayloads
virtual void createPayloads(const GearDir &content, const IntervalOfValidity &iov) override
Function to create the geometry database.
Definition: GeoMagneticField.cc:138
Belle2::GeoMagneticField::GeoMagneticField
GeoMagneticField()
Constructor of the GeoMagneticField class.
Definition: GeoMagneticField.cc:46
Belle2::GeoMagneticField::addConstantBField
void addConstantBField(const GearDir &component, MagneticField &fieldmap)
Add a constant field component to a magnetic field configuration for the DB.
Definition: GeoMagneticField.cc:75
Belle2::GeoMagneticField::createConfiguration
MagneticField createConfiguration(const GearDir &content)
Create a Database configuration from Gearbox parameters.
Definition: GeoMagneticField.cc:60
Belle2::GeoMagneticField::read3dBField
void read3dBField(const GearDir &component)
Reads the parameters for 3d magnetic field (r,phi,z).
Definition: GeoMagneticField.cc:289
Belle2::GeoMagneticField::create
virtual void create(const GearDir &content, G4LogicalVolume &topVolume, geometry::GeometryTypes type) override
Creates the global ROOT objects and prepares everything for other creators.
Definition: GeoMagneticField.cc:145
Belle2::GeoMagneticField::readQuadBField
void readQuadBField(const GearDir &component)
Reads the parameters for a quadrupole magnetic field inside beam pipes and adds the component to the ...
Definition: GeoMagneticField.cc:217
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition: IntervalOfValidity.h:25
Belle2::MagneticFieldComponent3D
Describe one component of the Geometry.
Definition: MagneticFieldComponent3D.h:22
Belle2::MagneticFieldComponentConstant
Describe one component of the Geometry.
Definition: MagneticFieldComponentConstant.h:19
Belle2::MagneticField
Magnetic field map.
Definition: MagneticField.h:32
Belle2::MagneticField::addComponent
void addComponent(MagneticFieldComponent *component)
Add a new component to the magnetic field.
Definition: MagneticField.h:59
Belle2::Unit::T
static const double T
[tesla]
Definition: Unit.h:120
Belle2::geometry::CreatorBase
Pure virtual base class for all geometry creators.
Definition: CreatorBase.h:28
Belle2::BFieldMap::addBFieldComponent
BFIELDCOMP & addBFieldComponent()
Adds a new BField component to the Belle II magnetic field.
Definition: BFieldMap.h:98
Belle2::geometry::GeometryTypes
GeometryTypes
Flag indiciating the type of geometry to be used.
Definition: GeometryManager.h:36
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
std
STL namespace.
Belle2::geometry::CreatorFactory
Very simple class to provide an easy way to register creators with the CreatorManager.
Definition: CreatorFactory.h:31