Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/GlueXMagneticField.cc
Location:line 350, column 33
Description:The left operand of '-' is a garbage value

Annotated Source Code

1//
2// class implementation for
3// * GlueXUniformMagField
4// * GlueXMappedMagField
5// * GlueXComputedMagField
6//
7// author: richard.t.jones at uconn.edu
8// version: may 12, 2012
9
10#include <iostream>
11#include <fstream>
12
13#include "GlueXMagneticField.hh"
14#include "GlueXUserOptions.hh"
15#include "GlueXDetectorConstruction.hh"
16#include "HddmOutput.hh"
17
18#include <JANA/JApplication.h>
19#include <JANA/JCalibration.h>
20
21#include <HDGEOMETRY/DMagneticFieldMapFineMesh.h>
22#include <HDGEOMETRY/DMagneticFieldMapNoField.h>
23#include <HDGEOMETRY/DMagneticFieldMapConst.h>
24#include <HDGEOMETRY/DMagneticFieldMapPS2DMap.h>
25#include <HDGEOMETRY/DMagneticFieldMapPSConst.h>
26
27// Implementation code for class GlueXUniformMagField
28
29GlueXUniformMagField::GlueXUniformMagField(G4ThreeVector B, G4double unit,
30 const G4AffineTransform &xform)
31 : G4UniformMagField(G4ThreeVector()),
32 fXform(xform)
33{
34 // uniform magnetic field constructor, constant field value B required
35 // as input. Factor unit converts B to standard G4 magnetic field units.
36 // Transform xform performs immediate active rotation on B, and is not
37 // saved. Translation in xform is ignored.
38
39 xform.ApplyAxisTransform(B);
40 SetMagField(B, unit);
41}
42
43GlueXUniformMagField::~GlueXUniformMagField()
44{ }
45
46GlueXUniformMagField::GlueXUniformMagField(const GlueXUniformMagField &src)
47 : G4UniformMagField(src.GetConstantFieldValue()),
48 fXform(src.fXform)
49{ }
50
51GlueXUniformMagField&
52GlueXUniformMagField::operator=(const GlueXUniformMagField &src)
53{
54 // assignment operator
55
56 fXform = src.fXform;
57 SetMagField(src.GetConstantFieldValue(), 1);
58 return *this;
59}
60
61void GlueXUniformMagField::SetMagFieldTransform(const G4AffineTransform &xform)
62{
63 // change the rotation applied to the field vector
64
65 G4ThreeVector B = GetConstantFieldValue();
66 fXform.Inverse().ApplyAxisTransform(B);
67 xform.ApplyAxisTransform(B);
68 SetFieldValue(B);
69 fXform = xform;
70}
71
72const G4AffineTransform&
73GlueXUniformMagField::GetMagFieldTransform() const
74{
75 // return the rotation applied to the field vector
76
77 return fXform;
78}
79
80void GlueXUniformMagField::SetMagField(G4ThreeVector B, G4double unit)
81{
82 // supplementary magnetic field setter, allows user to specify factor
83 // unit that converts input field B into G4 standard field units.
84 // Standard magnetic field setter method SetFieldValue works as well
85 // for cases when B is already in G4 units.
86
87 fXform.ApplyAxisTransform(B);
88 SetFieldValue(B *= unit);
89}
90
91G4ThreeVector GlueXUniformMagField::GetMagField(G4double unit) const
92{
93 // supplementary magnetic field getter, allows user to specify factor
94 // unit that is used to convert internal G4 field values into any
95 // desired units, eg. tesla.
96
97 G4ThreeVector B = GetConstantFieldValue();
98 if (B[2] == 0)
99 B[2] = 1e-99; // avoid divide-by-zero by excluding |B|=0
100 return B *= unit;
101}
102
103
104// Implementation code for class GlueXMappedMagField
105
106GlueXMappedMagField::GlueXMappedMagField(G4double Bmax, G4double unit,
107 const G4AffineTransform &xform)
108 : fUnit(unit),
109 fBmax(Bmax),
110 fXform(xform),
111 fGridtype(0)
112{
113 // mapped magnetic field constructor, maximum field value Bmax required
114 // as input. Factor unit converts B (Bmax and field components that will
115 // be read from the map) to standard G4 magnetic field units. Transform
116 // xform provides 2 functions: its inverse is used to transform the user
117 // coordinates used to access the map into local field map coordinates,
118 // and then once the map components are extracted from the map, an active
119 // rotation from xform is performed on B before it is returned to the user.
120 // Note that one or more calls to either method AddCartesianGrid() or
121 // AddCylindricalGrid() followed by ReadMapFile() must be issued before
122 // the object is ready for calls to GetFieldValue() or GetMagField().
123
124 fXfinv = xform.Inverse();
125}
126
127GlueXMappedMagField::~GlueXMappedMagField()
128{ }
129
130GlueXMappedMagField::GlueXMappedMagField(const GlueXMappedMagField &src)
131 : G4MagneticField(src)
132{
133 // copy constructor
134
135 *this = src;
136}
137
138GlueXMappedMagField&
139GlueXMappedMagField::operator=(const GlueXMappedMagField &src)
140{
141 // assignment operator, map copies are deep so use sparingly
142
143 fUnit = src.fUnit;
144 fBmax = src.fBmax;
145 fXform = src.fXform;
146 fXfinv = src.fXfinv;
147 fGridtype = src.fGridtype;
148 for (int dim=0; dim < 3; ++dim)
149 {
150 fDim[dim] = src.fDim[dim];
151 fOrder[dim] = src.fOrder[dim];
152 }
153 fGrid = src.fGrid;
154 fFieldMap = src.fFieldMap;
155 return *this;
156}
157
158void GlueXMappedMagField::SetMagFieldTransform(const G4AffineTransform &xform)
159{
160 // provides a means to modify the placement of the mapped field in
161 // the geometry without having to reinitialize the map.
162
163 fXform = xform;
164 fXfinv = xform.Inverse();
165}
166
167const G4AffineTransform& GlueXMappedMagField::GetMagFieldTransform() const
168{
169 // returns a copy of the current map placement transform
170
171 return fXform;
172}
173
174int GlueXMappedMagField::AddCartesianGrid(const int axsamples[4],
175 const int axorder[4],
176 const int axsense[4],
177 const G4double axlower[4],
178 const G4double axupper[4])
179{
180 // specifies that the field map is specified on the nodes of a Cartesian
181 // grid, what the dimensions of the grid are (axsamples), how the 3D
182 // array of grid points is ordered into a linear list (axorder), whether
183 // or not to reverse the sign of the field coordinates when they are read
184 // from the map (axsense), and the lower (axlower) and upper (axupper)
185 // coordinates of the bounding box of the grid. Multiple calls are allowed
186 // to AddCartesianGrid, to enable the same map to be replicated several
187 // times (with different axsense inversions in each case, as needed) so
188 // that symmetries in the field geometry can be exploited to minimize the
189 // total size of the map and eliminate redundancies.
190 // Input arrays all have the same component structure: x=1, y=2, z=3
191 // and the 0'th component is unused. Units for axlower and axupper are cm.
192 // Values in axorder[1..3] should be in the range 1..3. Values in axsense
193 // should be either +1 or -1. Repeated calls to AddCartesianGrid must all
194 // agree in the contents of axsamples and axorder, but can differ in the
195 // remaining arguments. This is assumed but not checked.
196
197 fGridtype = 1;
198 struct field_map_grid_t grid;
199 for (int dim=0; dim < 3; ++dim)
200 {
201 fDim[dim] = axsamples[dim + 1];
202 fOrder[dim] = axorder[dim + 1];
203 grid.sense.push_back(axsense[dim + 1]);
204 grid.lower.push_back(axlower[dim + 1]);
205 grid.upper.push_back(axupper[dim + 1]);
206 }
207 fGrid.push_back(grid);
208 return 1;
209}
210
211int GlueXMappedMagField::AddCylindricalGrid(const int axsamples[4],
212 const int axorder[4],
213 const int axsense[4],
214 const G4double axlower[4],
215 const G4double axupper[4])
216{
217 // specifies that the field map is specified on the nodes of a cylindrical
218 // grid, what the dimensions of the grid are (axsamples), how the 3D
219 // array of grid points is ordered into a linear list (axorder), whether
220 // or not to reverse the sign of the field coordinates when they are read
221 // from the map (axsense), and the lower (axlower) and upper (axupper)
222 // coordinates of the bounding box of the grid. Multiple calls are allowed
223 // to AddCylindricalGrid, to enable the same map to be replicated several
224 // times (with different axsense inversions in each case, as needed) so
225 // that symmetries in the field geometry can be exploited to minimize the
226 // total size of the map and eliminate redundancies.
227 // Input arrays all have the same component structure: rho=1, phi=2, z=3
228 // and the 0'th component is unused. Units for axlower and axupper are cm
229 // for rho and z, and radians for phi. Values in axorder[1..3] should be in
230 // the range 1..3. Values in axsense should be either +1 or -1. Repeated
231 // calls to AddCylindricalGrid must all agree in the contents of axsamples
232 // and axorder, but can differ in the remaining arguments. This is assumed
233 // but not checked.
234
235 fGridtype = 2;
236 struct field_map_grid_t grid;
237 for (int dim=0; dim < 3; ++dim)
238 {
239 fDim[dim] = axsamples[dim + 1];
240 fOrder[dim] = axorder[dim + 1];
241 grid.sense.push_back(axsense[dim + 1]);
242 grid.lower.push_back(axlower[dim + 1]);
243 grid.upper.push_back(axupper[dim + 1]);
244 }
245 fGrid.push_back(grid);
246 return 1;
247}
248
249int GlueXMappedMagField::ReadMapFile(const char *mapS)
250{
251 // reads field values from an input text file and stores them in an
252 // internal table for interpolation according to a preconceived spatial
253 // grid. The grid coordinates are not contained in the file, and are
254 // specified separately through calls to methods AddCartesianGrid() or
255 // AddCylindricalGrid(). A map is assumed to be either Cartesian or
256 // cylindrical, and invoking both methods on the same map will produce
257 // nonsense. The input file should contain one triplet of floating point
258 // numbers on each line separated by spaces, each representing the
259 // magnetic field components at a given grid site ordered as Bx By Bz
260 // for a Cartesian grid or Brho Bphi Bz for a cylindrical grid. It is
261 // assumed (but not checked) that the number of lines in the file equals
262 // the product of the three dimensions of the grid. Note that at least
263 // one call to either AddCartesianGrid() or AddCylindricalGrid() must
264 // have occurred on this object prior to the invocation of ReadMapFile.
265
266 std::ifstream mapfile(mapS);
267 if (! mapfile.good())
268 {
269 G4cerr(*G4cerr_p) << "GlueXMappedMagField::ReadMapFile error - "
270 << "open failed on input map file " << mapS
271 << G4endlstd::endl;
272 return 0;
273 }
274
275 while (mapfile.good())
276 {
277 union field_map_entry_t entry;
278 mapfile >> entry.cart.x >> entry.cart.y >> entry.cart.z;
279 if (mapfile.good())
280 {
281 fFieldMap.push_back(entry);
282 }
283 else {
284 break;
285 }
286 }
287 return 1;
288}
289
290G4ThreeVector GlueXMappedMagField::GetMagField(const G4double point[4],
291 G4double unit) const
292{
293 // interpolates the magnetic field map at point, and returns the result
294 // scaled to the requested unit (eg. tesla). Position point[0..2] gives
295 // the Cartesian coordinates at which the map is to be evaluated. The
296 // time stored in point[3] is presently ignored. The position is transformed
297 // into local map coordinates, and a 3D linear interpolation based on
298 // central estimates for the local field gradient is used. Finally, the
299 // field is rotated back into user coordinates and returned in the form
300 // of a Cartesian 3-vector.
301
302 G4ThreeVector p(point[0], point[1], point[2]);
303 fXfinv.ApplyPointTransform(p);
304 double u[3];
305 if (fGridtype == 1)
2
←
Taking false branch
→
306 {
307 u[0] = p[0] / cm;
308 u[1] = p[1] / cm;
309 u[2] = p[2] / cm;
310 }
311 else if (fGridtype == 2)
3
←
Taking true branch
→
312 {
313 u[0] = sqrt(p[0] * p[0] + p[1] * p[1]) / cm;
314 u[1] = atan2(p[1], p[0]);
315 u[2] = p[2] / cm;
316 }
317 else
318 {
319 return G4ThreeVector(0, 0, 1e-99);
320 }
321
322 double B[3] = {0,0,0};
323 std::vector<struct field_map_grid_t>::const_iterator iter;
324 for (iter = fGrid.begin(); iter != fGrid.end(); ++iter)
4
←
Loop condition is true. Entering loop body
→
325 {
326 double unorm[3];
327 unorm[0] = (u[0] - iter->lower[0]) / (iter->upper[0] - iter->lower[0]);
328 unorm[1] = (u[1] - iter->lower[1]) / (iter->upper[1] - iter->lower[1]);
329 unorm[2] = (u[2] - iter->lower[2]) / (iter->upper[2] - iter->lower[2]);
330 if (unorm[0] < 0 || unorm[0] > 1 ||
5
←
Taking false branch
→
331 unorm[1] < 0 || unorm[1] > 1 ||
332 unorm[2] < 0 || unorm[2] > 1)
333 {
334 break;
335 }
336 double ur[3], dur[3];
337 int ui[3], ui0[3], ui1[3];
338 for (int dim = 0; dim < 3; ++dim)
6
←
Loop condition is true. Entering loop body
→
9
←
Loop condition is true. Entering loop body
→
12
←
Loop condition is true. Entering loop body
→
15
←
Loop condition is false. Execution continues on line 346
→
339 {
340 ur[dim] = unorm[dim] * (fDim[dim] - 1);
341 ui[dim] = ceil(ur[dim] - 0.5);
342 ui0[dim] = (ui[dim] > 0)? ui[dim] - 1 : ui[dim];
7
←
'?' condition is false
→
10
←
'?' condition is false
→
13
←
'?' condition is false
→
343 ui1[dim] = (ui[dim] < fDim[dim] - 1)? ui[dim] + 1 : ui[dim];
8
←
'?' condition is false
→
11
←
'?' condition is true
→
14
←
'?' condition is false
→
344 dur[dim] = (ur[dim] - ui[dim]) / (ui1[dim] - ui0[dim] + 1e-20);
345 }
346 double grad[3][3];
347 double mapfield0[3], mapfield1[3];
348 lookup_field(ui0[0], ui[1], ui[2], mapfield0);
349 lookup_field(ui1[0], ui[1], ui[2], mapfield1);
16
←
Calling 'GlueXMappedMagField::lookup_field'
→
20
←
Returning from 'GlueXMappedMagField::lookup_field'
→
350 grad[0][0] = mapfield1[0] - mapfield0[0];
21
←
The left operand of '-' is a garbage value
351 grad[1][0] = mapfield1[1] - mapfield0[1];
352 grad[2][0] = mapfield1[2] - mapfield0[2];
353 lookup_field(ui[0], ui0[1], ui[2], mapfield0);
354 lookup_field(ui[0], ui1[1], ui[2], mapfield1);
355 grad[0][1] = mapfield1[0] - mapfield0[0];
356 grad[1][1] = mapfield1[1] - mapfield0[1];
357 grad[2][1] = mapfield1[2] - mapfield0[2];
358 lookup_field(ui[0], ui[1], ui0[2], mapfield0);
359 lookup_field(ui[0], ui[1], ui1[2], mapfield1);
360 grad[0][2] = mapfield1[0] - mapfield0[0];
361 grad[1][2] = mapfield1[1] - mapfield0[1];
362 grad[2][2] = mapfield1[2] - mapfield0[2];
363
364 lookup_field(ui[0], ui[1], ui[2], B);
365 B[0] += grad[0][0] * dur[0] + grad[0][1] * dur[1] + grad[0][2] * dur[2];
366 B[1] += grad[1][0] * dur[0] + grad[1][1] * dur[1] + grad[1][2] * dur[2];
367 B[2] += grad[2][0] * dur[0] + grad[2][1] * dur[1] + grad[2][2] * dur[2];
368 B[0] *= iter->sense[0];
369 B[1] *= iter->sense[1];
370 B[2] *= iter->sense[2];
371 }
372
373 G4ThreeVector Bvec;
374 if (fGridtype == 1)
375 {
376 Bvec.set(B[0], B[1], B[2]);
377 }
378 else if (fGridtype == 2)
379 {
380 Bvec.setRhoPhiZ(B[0], B[1], B[2]);
381 }
382 fXform.ApplyAxisTransform(Bvec);
383 if (Bvec[2] == 0)
384 Bvec[2] = 1e-99; // avoid divide-by-zero by excluding |B|=0
385 return Bvec *= fUnit / unit;
386}
387
388void GlueXMappedMagField::GetFieldValue(const G4double point[4],
389 G4double Bfield[3] ) const
390{
391 // implements the standard G4MagneticField::GetFieldValue() interface
392
393 G4ThreeVector Bvec = GetMagField(point, 1);
1
Calling 'GlueXMappedMagField::GetMagField'
→
394 Bfield[0] = Bvec[0];
395 Bfield[1] = Bvec[1];
396 Bfield[2] = Bvec[2];
397}
398
399int GlueXMappedMagField::lookup_field(int i1, int i2, int i3,
400 G4double mapvalue[3]) const
401{
402 // private helper method to perform lookup of the field at a single
403 // site in the field map. Many such accesses are performed in a
404 // single interpolation of the field, so it should be reasonably
405 // efficient.
406
407 int ivec[3] = {i1, i2, i3};
408 int iord[3] = {ivec[fOrder[0]], ivec[fOrder[1]], ivec[fOrder[2]]};
409 int nord[3] = {fDim[fOrder[0]], fDim[fOrder[1]], fDim[fOrder[2]]};
410 int index = nord[1] * (nord[0] * iord[0] + iord[1]) + iord[2];
411 if (index < (int)fFieldMap.size())
17
←
Calling 'vector::size'
→
18
←
Returning from 'vector::size'
→
19
←
Taking false branch
→
412 {
413 mapvalue[0] = fFieldMap[index].cart.x;
414 mapvalue[1] = fFieldMap[index].cart.y;
415 mapvalue[2] = fFieldMap[index].cart.z;
416 return 1;
417 }
418 return 0;
419}
420
421
422// Implementation code for class GlueXComputedMagField
423// Here a "computed" field covers any case where the map from position to
424// field value is returned by an external function, whether that function
425// implements the field as a formula or looks it up in a database. This is
426// the choice to use for any case where the field is not uniform, and the
427// map is not stored in the specific format specified for HDDS field maps.
428
429GlueXComputedMagField::GlueXComputedMagField(G4double Bmax, G4double unit,
430 const G4AffineTransform &xform)
431 : fUnit(unit),
432 fBmax(Bmax),
433 fXform(xform),
434 fJanaFieldMap(0),
435 fJanaFieldMapPS(0)
436{
437 // computed magnetic field constructor, maximum field value Bmax required
438 // as input. Factor unit converts B (Bmax and field components that will
439 // be read from the map) to standard G4 magnetic field units. Transform
440 // xform provides 2 functions: its inverse is used to transform the user
441 // coordinates used to access the map into local field map coordinates,
442 // and then once the map components are extracted from the map, an active
443 // rotation from xform is performed on B before it is returned to the user.
444 // Note that method SetFunction() must be issued before the object is ready
445 // for calls to GetFieldValue() or GetMagField().
446}
447
448GlueXComputedMagField::~GlueXComputedMagField()
449{
450 if (fJanaFieldMap)
451 delete fJanaFieldMap;
452 if (fJanaFieldMapPS)
453 delete fJanaFieldMapPS;
454}
455
456GlueXComputedMagField::GlueXComputedMagField(const GlueXComputedMagField &src)
457 : G4MagneticField(src)
458{
459 // copy constructor
460
461 *this = src;
462}
463
464GlueXComputedMagField&
465GlueXComputedMagField::operator=(const GlueXComputedMagField &src)
466{
467 // assignment operator
468
469 fBmax = src.fBmax;
470 fUnit = src.fUnit;
471 fXform = src.fXform;
472 fFunction = src.fFunction;
473 fJanaFieldMap = 0;
474 DMagneticFieldMapFineMesh *mapFineMesh = 0;
475 DMagneticFieldMapNoField *mapNoField = 0;
476 DMagneticFieldMapConst *mapConstField = 0;
477 if (src.fJanaFieldMap) {
478 mapFineMesh = dynamic_cast<DMagneticFieldMapFineMesh*>(src.fJanaFieldMap);
479 mapNoField = dynamic_cast<DMagneticFieldMapNoField*>(src.fJanaFieldMap);
480 mapConstField = dynamic_cast<DMagneticFieldMapConst*>(src.fJanaFieldMap);
481 if (mapFineMesh)
482 fJanaFieldMap = new DMagneticFieldMapFineMesh(*mapFineMesh);
483 else if (mapNoField)
484 fJanaFieldMap = new DMagneticFieldMapNoField(*mapNoField);
485 else if (mapConstField)
486 fJanaFieldMap = new DMagneticFieldMapConst(*mapConstField);
487 }
488 fJanaFieldMapPS = 0;
489 DMagneticFieldMapPS2DMap *mapPS2D = 0;
490 DMagneticFieldMapPSConst *mapPSConst = 0;
491 if (src.fJanaFieldMapPS) {
492 mapPS2D = dynamic_cast<DMagneticFieldMapPS2DMap*>(src.fJanaFieldMapPS);
493 mapPSConst = dynamic_cast<DMagneticFieldMapPSConst*>(src.fJanaFieldMapPS);
494 if (mapPS2D)
495 fJanaFieldMapPS = new DMagneticFieldMapPS2DMap(*mapPS2D);
496 else if (mapPSConst)
497 fJanaFieldMapPS = new DMagneticFieldMapPSConst(*mapPSConst);
498 }
499 return *this;
500}
501
502void GlueXComputedMagField::SetFunction(std::string function)
503{
504 // Pull a DMagneticField object from the JANA framework
505 // and store a pointer to it in this object. Note that
506 // GlueXComputedMagField never actually owns the object.
507 // The function argument supports a few special values:
508 // 1) gufld_db(r,B) - refers to the GlueX solenoid field
509 // 2) gufld_ps(r,B) - refers to the GlueX pair spectrometer field
510 // Lookup of the appropriate magnetic field map in the ccdb
511 // is done using a database key whose value is fetched from the
512 // global GlueXUserOptions object. Hard-coded defaults are also
513 // provided below, but normally the user should specify the
514 // correct field map key as an input option, eg. through the
515 // control.in file.
516
517 extern jana::JApplication *japp;
518 if (japp == 0) {
519 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
520 << "jana global DApplication object not set, "
521 << "cannot continue." << G4endlstd::endl;
522 exit(-1);
523 }
524 jana::JParameterManager *jpars = japp->GetJParameterManager();
525 if (jpars == 0) {
526 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
527 << "jana global DApplication object returns null "
528 << "JParameterManager object, cannot continue." << G4endlstd::endl;
529 exit(-1);
530 }
531 GlueXUserOptions *user_opts = GlueXUserOptions::GetInstance();
532 if (user_opts == 0) {
533 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
534 << "GlueXUserOptions::GetInstance() returns null, "
535 << "cannot continue." << G4endlstd::endl;
536 exit(-1);
537 }
538
539 int run_number = HddmOutput::getRunNo();
540 jana::JCalibration *jcalib = japp->GetJCalibration(run_number);
541
542 if (function == "gufld_db(r,B)") {
543
544 // load the solenoid field map
545
546 std::map<int, std::string> type_opts;
547 std::map<int, std::string> map_opts;
548 GlueXUserOptions::GetInstance()->Find("BFIELDTYPE", type_opts);
549 GlueXUserOptions::GetInstance()->Find("BFIELDMAP", map_opts);
550 if (type_opts.find(1) == type_opts.end() || type_opts[1] == "CalibDB") {
551
552 // if the magnetic field is specified in control.in
553 // then use that value instead of the CCDB values
554
555 if (map_opts.find(1) != map_opts.end()) {
556 fJanaFieldMap = new DMagneticFieldMapFineMesh(japp, run_number,
557 map_opts[1]);
558 }
559
560 // otherwise see if we can load the name of the
561 // magnetic field map to use from the ccdb itself
562
563 else {
564 std::map<std::string, std::string> map_name;
565 const char *map_key = "/Magnets/Solenoid/solenoid_map";
566 if (jcalib->GetCalib(map_key, map_name)) {
567 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
568 << "failed to figure out which magnetic field map "
569 << "to use for the solenoid in this simulation, "
570 << "cannot continue." << G4endlstd::endl;
571 exit(-1);
572 }
573 else if (map_name.find("map_name") != map_name.end()) {
574 if (map_name["map_name"] == "NoField")
575 fJanaFieldMap = new DMagneticFieldMapNoField(japp);
576 else
577 fJanaFieldMap = new DMagneticFieldMapFineMesh(japp,
578 run_number, map_name["map_name"]);
579 }
580 else {
581 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
582 << "no solenoid magnetic field map specified for "
583 << "this run in either the simulation options "
584 << "or in ccdb, cannot continue." << G4endlstd::endl;
585 exit(-1);
586 }
587 }
588 }
589 else if (type_opts[1] =="NoField") {
590 fJanaFieldMap = new DMagneticFieldMapNoField(japp);
591 }
592 else if (type_opts[1] == "Const") {
593 const GlueXDetectorConstruction *geom =
594 GlueXDetectorConstruction::GetInstance();
595 double field_T = (geom)? geom->GetUniformField(tesla) : 1.9;
596 if (type_opts.find(2) != type_opts.end()) {
597 field_T = std::atof(type_opts[2].c_str());
598 }
599 fJanaFieldMap = new DMagneticFieldMapConst(0.0, 0.0, field_T);
600 }
601 else {
602 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
603 << "unknown DMagneticFieldMap type " << type_opts[1]
604 << ", cannot continue." << G4endlstd::endl;
605 exit(-1);
606 }
607 }
608
609 else if (function == "gufld_ps(r,B)") {
610
611 // load the pair spectrometer field map
612
613 std::map<int, std::string> type_opts;
614 std::map<int, std::string> map_opts;
615 GlueXUserOptions::GetInstance()->Find("PSBFIELDTYPE", type_opts);
616 GlueXUserOptions::GetInstance()->Find("PSBFIELDMAP", map_opts);
617 if (type_opts.find(1) == type_opts.end() || type_opts[1] == "CalibDB") {
618
619 // if the magnetic field is specified in control.in
620 // then use that value instead of the CCDB values
621
622 if (map_opts.find(1) != map_opts.end()) {
623 fJanaFieldMapPS = new DMagneticFieldMapPS2DMap(japp, run_number,
624 map_opts[1]);
625 }
626 else {
627
628 // see if we can load the name of the magnetic
629 // field map to use from the calib DB
630
631 map<string,string> map_name;
632 const char *map_key = "/Magnets/PairSpectrometer/ps_magnet_map";
633 if (jcalib->GetCalib(map_key, map_name)) {
634 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
635 << "failed to figure out which magnetic field map "
636 << "to use for the pair spectrometer, "
637 << "cannot continue." << G4endlstd::endl;
638 exit(-1);
639 }
640 else if (map_name.find("map_name") != map_name.end()) {
641 fJanaFieldMapPS = new DMagneticFieldMapPS2DMap(japp,
642 run_number, map_name["map_name"]);
643 }
644 else {
645 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
646 << "no pair spectrometer magnetic field map "
647 << "specified for this run either in the options "
648 << "or in ccdb, cannot continue." << G4endlstd::endl;
649 exit(-1);
650 }
651 }
652 }
653 else if (type_opts[1] == "Const") {
654 double field_T = 1.64;
655 if (type_opts.find(2) != type_opts.end()) {
656 field_T = std::atof(type_opts[2].c_str());
657 }
658 fJanaFieldMapPS = new DMagneticFieldMapPSConst(0.0, field_T, 0.0);
659 }
660 else {
661 G4cerr(*G4cerr_p) << "Error in GlueXComputedMagField::SetFunction - "
662 << "unknown DMagneticFieldMapPS type " << type_opts[1]
663 << ", cannot continue." << G4endlstd::endl;
664 exit(-1);
665 }
666 }
667
668 fFunction = function;
669}
670
671std::string GlueXComputedMagField::GetFunction() const
672{
673 // returns the string specified in a prior call to SetFunction,
674 // or the null string
675
676 return fFunction;
677}
678
679void GlueXComputedMagField::SetMagFieldTransform(const G4AffineTransform &xform)
680{
681 // provides a means to modify the placement of the mapped field in
682 // the geometry without having to reinitialize the map.
683
684 fXform = xform;
685 fXfinv = xform.Inverse();
686}
687
688const G4AffineTransform& GlueXComputedMagField::GetMagFieldTransform() const
689{
690 // returns a copy of the current map placement transform
691
692 return fXform;
693}
694
695G4ThreeVector GlueXComputedMagField::GetMagField(const G4double point[4],
696 G4double unit) const
697{
698 // looks up the magnetic field at the specified point, and returns
699 // the result scaled to the requested unit (eg. gauss). The first 3
700 // elements of point are the Cartesian coordinates in the G4 world
701 // coordinate system at which the map is to be evaluated. The time
702 // stored in point[3] is presently ignored. The position is transformed
703 // into local map coordinates for lookup, and then the field is rotated
704 // back into user coordinates and returned in the requested units.
705
706 G4ThreeVector p(point[0] / cm, point[1] / cm, point[2] / cm);
707 fXfinv.ApplyPointTransform(p);
708 double B[3] = {0,0,0};
709 DMagneticFieldMap *mapso = fJanaFieldMap;
710 //dynamic_cast <DMagneticFieldMap* const> (fJanaFieldMap);
711 DMagneticFieldMapPS *mapps = fJanaFieldMapPS;
712 //dynamic_cast <DMagneticFieldMapPS* const> (fJanaFieldMapPS);
713 if (mapso)
714 mapso->GetField(p[0], p[1], p[2], B[0], B[1], B[2]);
715 else if (mapps)
716 mapps->GetField(p[0], p[1], p[2], B[0], B[1], B[2]);
717 G4ThreeVector Bvec(B[0], B[1], B[2]);
718 fXform.ApplyAxisTransform(Bvec);
719 if (Bvec[2] == 0)
720 Bvec[2] = 1e-99; // avoid divide-by-zero by excluding |B|=0
721 return Bvec *= fUnit / unit;
722}
723
724void GlueXComputedMagField::GetFieldValue(const G4double point[4],
725 G4double *Bfield ) const
726{
727 // implements the standard G4MagneticField::GetFieldValue() interface
728
729 G4ThreeVector Bvec = GetMagField(point, 1);
730 Bfield[0] = Bvec[0];
731 Bfield[1] = Bvec[1];
732 Bfield[2] = Bvec[2];
733}