Belle II Software light-2406-ragdoll
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
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
60
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
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
140{
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
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
167}
168
169
170//============================================================
171// Protected methods
172//============================================================
173
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
184 bComp.setMagneticFieldValues(xValue, yValue, zValue, RmaxValue, ZminValue, ZmaxValue);
185}
186
187
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
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
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
233 bComp.setMapFilename(mapFilenameHER, mapFilenameLER, mapFilenameHERleak);
234 bComp.setApertFilename(apertFilenameHER, apertFilenameLER);
235 bComp.setMapSize(mapSizeHER, mapSizeLER, mapSizeHERleak);
236 bComp.setApertSize(apertSizeHER, apertSizeLER);
237}
238
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
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
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
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
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
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}
The BFieldComponent3d class.
void setMapFilename(const std::string &filename)
Sets the filename of the magnetic field map.
The BFieldComponentBeamline class.
void setMapFilename(const std::string &filename_her, const std::string &filename_ler)
Sets the filename of the magnetic field map.
The BFieldComponentConstant class.
void setMagneticFieldValues(double x, double y, double z, double rmax, double zmin, double zmax)
Sets the values for the homogeneous magnetic field vector.
The Bfieldcomponentklm1 class.
void setMapFilename(const std::string &filename)
Sets the filename of the magnetic field map.
The BFieldComponentQuad class.
void setMapFilename(const std::string &filenameHER, const std::string &filenameLER, const std::string &filenameHERleak)
Sets the filename of the magnetic field map.
The BFieldComponentRadial class.
void setMapFilename(const std::string &filename)
Sets the filename of the magnetic field map.
void initialize()
Initialize the magnetic field after adding all components.
Definition: BFieldMap.cc:21
void clear()
Clear the existing components.
Definition: BFieldMap.cc:33
static BFieldMap & Instance()
Static method to get a reference to the BFieldMap instance.
Definition: BFieldMap.cc:15
bool import(const IntervalOfValidity &iov)
Import the object to database.
Definition: DBImportBase.cc:36
Class for importing a single object to the database.
void construct(Args &&... params)
Construct an object of type T in this DBImportObjPtr using the provided constructor arguments.
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
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:31
virtual std::string getString(const std::string &path="") const noexcept(false) override
Get the parameter path as a string.
Definition: GearDir.h:69
void readConstantBField(const GearDir &component)
Reads the parameters for a homogeneous magnetic field and adds the component to the global magnetic f...
void add3dBField(const GearDir &component, MagneticField &fielmap)
Add a 3D field component to a magnetic field configuration for the DB.
void readBeamlineBField(const GearDir &component)
Reads the 3D Bfield map and parameters near beam pipes and adds the component to the global magnetic ...
void readKlm1BField(const GearDir &component)
Reads the 2D Bfield map and parameters outside of solenoid and adds the component to the global magne...
CompTypeMap m_componentTypeMap
Maps the name of the component to the function reading the parameters.
virtual ~GeoMagneticField()
The destructor of the GeoMagneticField class.
void readRadialBField(const GearDir &component)
Reads the parameters for a radial magnetic field and adds the component to the global magnetic field.
virtual void createPayloads(const GearDir &content, const IntervalOfValidity &iov) override
Function to create the geometry database.
GeoMagneticField()
Constructor of the GeoMagneticField class.
void addConstantBField(const GearDir &component, MagneticField &fieldmap)
Add a constant field component to a magnetic field configuration for the DB.
MagneticField createConfiguration(const GearDir &content)
Create a Database configuration from Gearbox parameters.
void read3dBField(const GearDir &component)
Reads the parameters for 3d magnetic field (r,phi,z).
virtual void create(const GearDir &content, G4LogicalVolume &topVolume, geometry::GeometryTypes type) override
Creates the global ROOT objects and prepares everything for other creators.
void readQuadBField(const GearDir &component)
Reads the parameters for a quadrupole magnetic field inside beam pipes and adds the component to the ...
A class that describes the interval of experiments/runs for which an object in the database is valid.
Describe one component of the Geometry.
Describe one component of the Geometry.
Magnetic field map.
Definition: MagneticField.h:32
void addComponent(MagneticFieldComponent *component)
Add a new component to the magnetic field.
Definition: MagneticField.h:59
static const double T
[tesla]
Definition: Unit.h:120
double getDouble(const std::string &path="") const noexcept(false)
Get the parameter path as a double.
Definition: Interface.cc:41
double getLength(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard length unit.
Definition: Interface.h:259
double getWithUnit(const std::string &path) const noexcept(false)
Get the parameter path as a double converted to the standard unit.
Definition: Interface.cc:97
int getInt(const std::string &path="") const noexcept(false)
Get the parameter path as a int.
Definition: Interface.cc:60
BFIELDCOMP & addBFieldComponent()
Adds a new BField component to the Belle II magnetic field.
Definition: BFieldMap.h:98
GeometryTypes
Flag indiciating the type of geometry to be used.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
Very simple class to provide an easy way to register creators with the CreatorManager.