Implementation of primary virtual method from G4MagneticField interface.
145 {
146 double x = point[0] - x_offset_;
147 double y = point[1] - y_offset_;
148 double z = point[2] - z_offset_;
149 double eps = 1E-6;
150
151
152 if ((x >= minx_ && (x < maxx_ - eps)) && (y >= miny_ && (y < maxy_ - eps)) &&
153 (z >= minz_ && (z < maxz_ - eps))) {
154
155
156 double xfraction = (x - minx_) / dx_;
157 double yfraction = (y - miny_) / dy_;
158 double zfraction = (z - minz_) / dz_;
159
160 if (invert_x_) {
161 xfraction = 1 - xfraction;
162 }
163 if (invert_y_) {
164 yfraction = 1 - yfraction;
165 }
166 if (invert_z_) {
167 zfraction = 1 - zfraction;
168 }
169
170
171
172 double xdindex, ydindex, zdindex;
173
174
175
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
181
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 ldmx_log(trace) << "Local x_,y_,z_: " << xlocal << " " << ylocal << " "
188 << zlocal;
189 ldmx_log(trace) << "Index x_,y_,z_: " << xindex << " " << yindex << " "
190 << zindex;
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
204 bfield[0] =
205 x_field_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
206 (1 - zlocal) +
207 x_field_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
208 zlocal +
209 x_field_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
210 (1 - zlocal) +
211 x_field_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
212 zlocal +
213 x_field_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
214 (1 - zlocal) +
215 x_field_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
216 zlocal +
217 x_field_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
218 (1 - zlocal) +
219 x_field_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
220 bfield[1] =
221 y_field_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
222 (1 - zlocal) +
223 y_field_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
224 zlocal +
225 y_field_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
226 (1 - zlocal) +
227 y_field_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
228 zlocal +
229 y_field_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
230 (1 - zlocal) +
231 y_field_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
232 zlocal +
233 y_field_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
234 (1 - zlocal) +
235 y_field_[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
236 bfield[2] =
237 z_field_[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) *
238 (1 - zlocal) +
239 z_field_[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) *
240 zlocal +
241 z_field_[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal *
242 (1 - zlocal) +
243 z_field_[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal *
244 zlocal +
245 z_field_[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) *
246 (1 - zlocal) +
247 z_field_[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) *
248 zlocal +
249 z_field_[xindex + 1][yindex + 1][zindex] * xlocal * ylocal *
250 (1 - zlocal) +
251 z_field_[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}