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