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