LDMX Software
Public Member Functions | Private Attributes | List of all members
simcore::MagneticFieldMap3D Class Reference

A 3D B-field map defined as a grid of points with associated B-field values. More...

#include <MagneticFieldMap3D.h>

Public Member Functions

 MagneticFieldMap3D (const char *filename, double xOffset, double yOffset, double zOffset)
 Class constructor.
 
void GetFieldValue (const double point[4], double *bfield) const
 Implementation of primary virtual method from G4MagneticField interface.
 

Private Attributes

vector< vector< vector< double > > > xField_
 
vector< vector< vector< double > > > yField_
 
vector< vector< vector< double > > > zField_
 
int nx_
 
int ny_
 
int nz_
 
double minx_
 
double maxx_
 
double miny_
 
double maxy_
 
double minz_
 
double maxz_
 
double dx_
 
double dy_
 
double dz_
 
double xOffset_
 
double yOffset_
 
double zOffset_
 
bool invertX_
 
bool invertY_
 
bool invertZ_
 

Detailed Description

A 3D B-field map defined as a grid of points with associated B-field values.

Note
Values are interpolated to obtain B-field information at points in between.

The header format of the input data file is as follows:

1) Blank line

2) int int int

Number of X, Y, and Z grid points (N_x, N_y, N_z)

3-8) int String

A description of fields in the lines of the map. In our case it is literally:

1 X 2 Y 3 Z 4 BX 5 BY 6 BZ

9) int String

Header terminator

10+) float float float float float float

N_x*N_y*N_z+9

x y z B_x B_y B_z

Original PurgMagTabulatedField3D code developed by: S.Larsson and J. Generowicz.

Definition at line 62 of file MagneticFieldMap3D.h.

Constructor & Destructor Documentation

◆ MagneticFieldMap3D()

simcore::MagneticFieldMap3D::MagneticFieldMap3D ( const char *  filename,
double  xOffset,
double  yOffset,
double  zOffset 
)

Class constructor.

Parameters
[in]filenameThe name of the file defining the B-field grid.
[in]xOffset,yOffset,zOffsetThe offset of the grid's coordinate system.

Definition at line 18 of file MagneticFieldMap3D.cxx.

20 : nx_(0),
21 ny_(0),
22 nz_(0),
23 xOffset_(xOffset),
24 yOffset_(yOffset),
25 zOffset_(zOffset),
26 invertX_(false),
27 invertY_(false),
28 invertZ_(false) {
29 ifstream file(filename); // Open the file for reading.
30
31 // Throw an error if file does not exist.
32 if (!file.good()) {
33 EXCEPTION_RAISE("FileDNE", "The field map file '" + std::string(filename) +
34 "' does not exist!");
35 }
36
37 G4cout << "-----------------------------------------------------------"
38 << G4endl;
39 G4cout << " Magnetic Field Map 3D" << G4endl;
40 G4cout << "-----------------------------------------------------------"
41 << G4endl << G4endl;
42
43 G4cout << "Reading the field grid from " << filename << " ... " << endl;
44 G4cout << " Offsets: " << xOffset << " " << yOffset << " " << zOffset
45 << G4endl;
46
47 // Ignore first blank line
48 char buffer[256];
49 file.getline(buffer, 256);
50
51 // Read table dimensions
52 file >> nx_ >> ny_ >> nz_; // Note dodgy order
53
54 G4cout << " Number of values: " << nx_ << " " << ny_ << " " << nz_ << G4endl;
55
56 // Set up storage space for table
57 xField_.resize(nx_);
58 yField_.resize(nx_);
59 zField_.resize(nx_);
60 int ix, iy, iz;
61 for (ix = 0; ix < nx_; ix++) {
62 xField_[ix].resize(ny_);
63 yField_[ix].resize(ny_);
64 zField_[ix].resize(ny_);
65 for (iy = 0; iy < ny_; iy++) {
66 xField_[ix][iy].resize(nz_);
67 yField_[ix][iy].resize(nz_);
68 zField_[ix][iy].resize(nz_);
69 }
70 }
71
72 // Ignore other header information
73 // The first line whose second character is '0' is considered to
74 // be the last line of the header.
75 do {
76 file.getline(buffer, 256);
77 } while (buffer[1] != '0');
78
79 // Read in the data
80 double xval, yval, zval, bx, by, bz;
81 for (ix = 0; ix < nx_; ix++) {
82 for (iy = 0; iy < ny_; iy++) {
83 for (iz = 0; iz < nz_; iz++) {
84 file >> xval >> yval >> zval >> bx >> by >> bz;
85 if (ix == 0 && iy == 0 && iz == 0) {
86 minx_ = xval;
87 miny_ = yval;
88 minz_ = zval;
89 }
90 xField_[ix][iy][iz] = bx;
91 yField_[ix][iy][iz] = by;
92 zField_[ix][iy][iz] = bz;
93 }
94 }
95 }
96 file.close();
97
98 maxx_ = xval;
99 maxy_ = yval;
100 maxz_ = zval;
101
102 G4cout << " ... done reading " << G4endl << G4endl;
103 G4cout << "Read values of field from file " << filename << G4endl;
104 G4cout << " Assumed the order: x, y, z, Bx, By, Bz" << G4endl;
105 G4cout << " Min values: " << minx_ << " " << miny_ << " " << minz_ << " mm "
106 << G4endl;
107 G4cout << " Max values: " << maxx_ << " " << maxy_ << " " << maxz_ << " mm "
108 << G4endl;
109 G4cout << " Field offsets: " << xOffset_ << " " << yOffset_ << " "
110 << zOffset_ << " mm " << G4endl << G4endl;
111
112 // Should really check that the limits are not the wrong way around.
113 if (maxx_ < minx_) {
114 swap(maxx_, minx_);
115 invertX_ = true;
116 }
117 if (maxy_ < miny_) {
118 swap(maxy_, miny_);
119 invertY_ = true;
120 }
121 if (maxz_ < minz_) {
122 swap(maxz_, minz_);
123 invertZ_ = true;
124 }
125
126 G4cout << "After reordering if necessary" << G4endl;
127 G4cout << " Min values: " << minx_ << " " << miny_ << " " << minz_ << " mm "
128 << G4endl;
129 G4cout << " Max values: " << maxx_ << " " << maxy_ << " " << maxz_ << " mm "
130 << G4endl;
131 ;
132
133 dx_ = maxx_ - minx_;
134 dy_ = maxy_ - miny_;
135 dz_ = maxz_ - minz_;
136
137 G4cout << " Range of values: " << dx_ << " " << dy_ << " " << dz_ << " mm"
138 << G4endl << G4endl;
139 G4cout << "Done loading field map" << G4endl << G4endl;
140 G4cout << "-----------------------------------------------------------"
141 << G4endl << G4endl;
142}

Member Function Documentation

◆ GetFieldValue()

void simcore::MagneticFieldMap3D::GetFieldValue ( const double  point[4],
double *  bfield 
) const

Implementation of primary virtual method from G4MagneticField interface.

Parameters
[in]pointThe point in 3D space.
[out]bfieldThe output B-field data at the point.

Definition at line 144 of file MagneticFieldMap3D.cxx.

145 {
146 double x = point[0] - xOffset_;
147 double y = point[1] - yOffset_;
148 double z = point[2] - zOffset_;
149 double eps = 1E-6;
150
151 // Check that the point is within the defined region
152 if (x >= minx_ && x < maxx_ - eps && y >= miny_ && y < maxy_ - eps &&
153 z >= minz_ && z < maxz_ - eps) {
154 // Position of given point within region, normalized to the range
155 // [0,1]
156 double xfraction = (x - minx_) / dx_;
157 double yfraction = (y - miny_) / dy_;
158 double zfraction = (z - minz_) / dz_;
159
160 if (invertX_) {
161 xfraction = 1 - xfraction;
162 }
163 if (invertY_) {
164 yfraction = 1 - yfraction;
165 }
166 if (invertZ_) {
167 zfraction = 1 - zfraction;
168 }
169
170 // Need addresses of these to pass to modf below.
171 // modf uses its second argument as an OUTPUT argument.
172 double xdindex, ydindex, zdindex;
173
174 // Position of the point within the cuboid defined by the
175 // nearest surrounding tabulated points
176 double xlocal = (std::modf(xfraction * (nx_ - 1), &xdindex));
177 double ylocal = (std::modf(yfraction * (ny_ - 1), &ydindex));
178 double zlocal = (std::modf(zfraction * (nz_ - 1), &zdindex));
179
180 // The indices of the nearest tabulated point whose coordinates
181 // are all less than those of the given point
182 int xindex = static_cast<int>(xdindex);
183 int yindex = static_cast<int>(ydindex);
184 int zindex = static_cast<int>(zdindex);
185
186#ifdef DEBUG_INTERPOLATING_FIELD
187 G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal
188 << G4endl;
189 G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex
190 << G4endl;
191 double valx0z0, mulx0z0, valx1z0, mulx1z0;
192 double valx0z1, mulx0z1, valx1z1, mulx1z1;
193 valx0z0 = table[xindex][0][zindex];
194 mulx0z0 = (1 - xlocal) * (1 - zlocal);
195 valx1z0 = table[xindex + 1][0][zindex];
196 mulx1z0 = xlocal * (1 - zlocal);
197 valx0z1 = table[xindex][0][zindex + 1];
198 mulx0z1 = (1 - xlocal) * zlocal;
199 valx1z1 = table[xindex + 1][0][zindex + 1];
200 mulx1z1 = xlocal * zlocal;
201#endif
202
203 // Full 3-dimensional version
204 bfield[0] =
205 xField_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
206 (1 - zlocal) +
207 xField_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
208 zlocal +
209 xField_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
210 (1 - zlocal) +
211 xField_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
212 zlocal +
213 xField_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
214 (1 - zlocal) +
215 xField_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
216 zlocal +
217 xField_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
218 (1 - zlocal) +
219 xField_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
220 bfield[1] =
221 yField_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
222 (1 - zlocal) +
223 yField_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
224 zlocal +
225 yField_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
226 (1 - zlocal) +
227 yField_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
228 zlocal +
229 yField_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
230 (1 - zlocal) +
231 yField_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
232 zlocal +
233 yField_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
234 (1 - zlocal) +
235 yField_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
236 bfield[2] =
237 zField_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
238 (1 - zlocal) +
239 zField_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
240 zlocal +
241 zField_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
242 (1 - zlocal) +
243 zField_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
244 zlocal +
245 zField_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
246 (1 - zlocal) +
247 zField_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
248 zlocal +
249 zField_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
250 (1 - zlocal) +
251 zField_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
252
253 } else {
254 bfield[0] = 0.0;
255 bfield[1] = 0.0;
256 bfield[2] = 0.0;
257 }
258}

Member Data Documentation

◆ dx_

double simcore::MagneticFieldMap3D::dx_
private

Definition at line 101 of file MagneticFieldMap3D.h.

◆ dy_

double simcore::MagneticFieldMap3D::dy_
private

Definition at line 101 of file MagneticFieldMap3D.h.

◆ dz_

double simcore::MagneticFieldMap3D::dz_
private

Definition at line 101 of file MagneticFieldMap3D.h.

◆ invertX_

bool simcore::MagneticFieldMap3D::invertX_
private

Definition at line 113 of file MagneticFieldMap3D.h.

◆ invertY_

bool simcore::MagneticFieldMap3D::invertY_
private

Definition at line 113 of file MagneticFieldMap3D.h.

◆ invertZ_

bool simcore::MagneticFieldMap3D::invertZ_
private

Definition at line 113 of file MagneticFieldMap3D.h.

◆ maxx_

double simcore::MagneticFieldMap3D::maxx_
private

Definition at line 96 of file MagneticFieldMap3D.h.

◆ maxy_

double simcore::MagneticFieldMap3D::maxy_
private

Definition at line 96 of file MagneticFieldMap3D.h.

◆ maxz_

double simcore::MagneticFieldMap3D::maxz_
private

Definition at line 96 of file MagneticFieldMap3D.h.

◆ minx_

double simcore::MagneticFieldMap3D::minx_
private

Definition at line 96 of file MagneticFieldMap3D.h.

◆ miny_

double simcore::MagneticFieldMap3D::miny_
private

Definition at line 96 of file MagneticFieldMap3D.h.

◆ minz_

double simcore::MagneticFieldMap3D::minz_
private

Definition at line 96 of file MagneticFieldMap3D.h.

◆ nx_

int simcore::MagneticFieldMap3D::nx_
private

Definition at line 91 of file MagneticFieldMap3D.h.

◆ ny_

int simcore::MagneticFieldMap3D::ny_
private

Definition at line 91 of file MagneticFieldMap3D.h.

◆ nz_

int simcore::MagneticFieldMap3D::nz_
private

Definition at line 91 of file MagneticFieldMap3D.h.

◆ xField_

vector<vector<vector<double> > > simcore::MagneticFieldMap3D::xField_
private

Definition at line 84 of file MagneticFieldMap3D.h.

◆ xOffset_

double simcore::MagneticFieldMap3D::xOffset_
private

Definition at line 106 of file MagneticFieldMap3D.h.

◆ yField_

vector<vector<vector<double> > > simcore::MagneticFieldMap3D::yField_
private

Definition at line 85 of file MagneticFieldMap3D.h.

◆ yOffset_

double simcore::MagneticFieldMap3D::yOffset_
private

Definition at line 107 of file MagneticFieldMap3D.h.

◆ zField_

vector<vector<vector<double> > > simcore::MagneticFieldMap3D::zField_
private

Definition at line 86 of file MagneticFieldMap3D.h.

◆ zOffset_

double simcore::MagneticFieldMap3D::zOffset_
private

Definition at line 108 of file MagneticFieldMap3D.h.

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.8