1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
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 | |
28 | |
29 | GlueXUniformMagField::GlueXUniformMagField(G4ThreeVector B, G4double unit, |
30 | const G4AffineTransform &xform) |
31 | : G4UniformMagField(G4ThreeVector()), |
32 | fXform(xform) |
33 | { |
34 | |
35 | |
36 | |
37 | |
38 | |
39 | xform.ApplyAxisTransform(B); |
40 | SetMagField(B, unit); |
41 | } |
42 | |
43 | GlueXUniformMagField::~GlueXUniformMagField() |
44 | { } |
45 | |
46 | GlueXUniformMagField::GlueXUniformMagField(const GlueXUniformMagField &src) |
47 | : G4UniformMagField(src.GetConstantFieldValue()), |
48 | fXform(src.fXform) |
49 | { } |
50 | |
51 | GlueXUniformMagField& |
52 | GlueXUniformMagField::operator=(const GlueXUniformMagField &src) |
53 | { |
54 | |
55 | |
56 | fXform = src.fXform; |
57 | SetMagField(src.GetConstantFieldValue(), 1); |
58 | return *this; |
59 | } |
60 | |
61 | void GlueXUniformMagField::SetMagFieldTransform(const G4AffineTransform &xform) |
62 | { |
63 | |
64 | |
65 | G4ThreeVector B = GetConstantFieldValue(); |
66 | fXform.Inverse().ApplyAxisTransform(B); |
67 | xform.ApplyAxisTransform(B); |
68 | SetFieldValue(B); |
69 | fXform = xform; |
70 | } |
71 | |
72 | const G4AffineTransform& |
73 | GlueXUniformMagField::GetMagFieldTransform() const |
74 | { |
75 | |
76 | |
77 | return fXform; |
78 | } |
79 | |
80 | void GlueXUniformMagField::SetMagField(G4ThreeVector B, G4double unit) |
81 | { |
82 | |
83 | |
84 | |
85 | |
86 | |
87 | fXform.ApplyAxisTransform(B); |
88 | SetFieldValue(B *= unit); |
89 | } |
90 | |
91 | G4ThreeVector GlueXUniformMagField::GetMagField(G4double unit) const |
92 | { |
93 | |
94 | |
95 | |
96 | |
97 | G4ThreeVector B = GetConstantFieldValue(); |
98 | if (B[2] == 0) |
99 | B[2] = 1e-99; |
100 | return B *= unit; |
101 | } |
102 | |
103 | |
104 | |
105 | |
106 | GlueXMappedMagField::GlueXMappedMagField(G4double Bmax, G4double unit, |
107 | const G4AffineTransform &xform) |
108 | : fUnit(unit), |
109 | fBmax(Bmax), |
110 | fXform(xform), |
111 | fGridtype(0) |
112 | { |
113 | |
114 | |
115 | |
116 | |
117 | |
118 | |
119 | |
120 | |
121 | |
122 | |
123 | |
124 | fXfinv = xform.Inverse(); |
125 | } |
126 | |
127 | GlueXMappedMagField::~GlueXMappedMagField() |
128 | { } |
129 | |
130 | GlueXMappedMagField::GlueXMappedMagField(const GlueXMappedMagField &src) |
131 | : G4MagneticField(src) |
132 | { |
133 | |
134 | |
135 | *this = src; |
136 | } |
137 | |
138 | GlueXMappedMagField& |
139 | GlueXMappedMagField::operator=(const GlueXMappedMagField &src) |
140 | { |
141 | |
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 | |
158 | void GlueXMappedMagField::SetMagFieldTransform(const G4AffineTransform &xform) |
159 | { |
160 | |
161 | |
162 | |
163 | fXform = xform; |
164 | fXfinv = xform.Inverse(); |
165 | } |
166 | |
167 | const G4AffineTransform& GlueXMappedMagField::GetMagFieldTransform() const |
168 | { |
169 | |
170 | |
171 | return fXform; |
172 | } |
173 | |
174 | int 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 | |
181 | |
182 | |
183 | |
184 | |
185 | |
186 | |
187 | |
188 | |
189 | |
190 | |
191 | |
192 | |
193 | |
194 | |
195 | |
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 | |
211 | int 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 | |
218 | |
219 | |
220 | |
221 | |
222 | |
223 | |
224 | |
225 | |
226 | |
227 | |
228 | |
229 | |
230 | |
231 | |
232 | |
233 | |
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 | |
249 | int GlueXMappedMagField::ReadMapFile(const char *mapS) |
250 | { |
251 | |
252 | |
253 | |
254 | |
255 | |
256 | |
257 | |
258 | |
259 | |
260 | |
261 | |
262 | |
263 | |
264 | |
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 | |
290 | G4ThreeVector GlueXMappedMagField::GetMagField(const G4double point[4], |
291 | G4double unit) const |
292 | { |
293 | |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | |
300 | |
301 | |
302 | G4ThreeVector p(point[0], point[1], point[2]); |
303 | fXfinv.ApplyPointTransform(p); |
304 | double u[3]; |
305 | if (fGridtype == 1) |
| |
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) |
| |
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 || |
| |
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]; |
| |
| |
| |
343 | ui1[dim] = (ui[dim] < fDim[dim] - 1)? ui[dim] + 1 : ui[dim]; |
| |
| |
| |
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); |
| 16 | | Calling 'GlueXMappedMagField::lookup_field' | |
|
| 20 | | Returning from 'GlueXMappedMagField::lookup_field' | |
|
349 | lookup_field(ui1[0], ui[1], ui[2], mapfield1); |
350 | grad[0][0] = mapfield1[0] - mapfield0[0]; |
| 21 | | The right 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; |
385 | return Bvec *= fUnit / unit; |
386 | } |
387 | |
388 | void GlueXMappedMagField::GetFieldValue(const G4double point[4], |
389 | G4double Bfield[3] ) const |
390 | { |
391 | |
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 | |
399 | int GlueXMappedMagField::lookup_field(int i1, int i2, int i3, |
400 | G4double mapvalue[3]) const |
401 | { |
402 | |
403 | |
404 | |
405 | |
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()) |
| |
| 18 | | Returning from 'vector::size' | |
|
| |
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 | |
423 | |
424 | |
425 | |
426 | |
427 | |
428 | |
429 | GlueXComputedMagField::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 | |
438 | |
439 | |
440 | |
441 | |
442 | |
443 | |
444 | |
445 | |
446 | } |
447 | |
448 | GlueXComputedMagField::~GlueXComputedMagField() |
449 | { |
450 | if (fJanaFieldMap) |
451 | delete fJanaFieldMap; |
452 | if (fJanaFieldMapPS) |
453 | delete fJanaFieldMapPS; |
454 | } |
455 | |
456 | GlueXComputedMagField::GlueXComputedMagField(const GlueXComputedMagField &src) |
457 | : G4MagneticField(src) |
458 | { |
459 | |
460 | |
461 | *this = src; |
462 | } |
463 | |
464 | GlueXComputedMagField& |
465 | GlueXComputedMagField::operator=(const GlueXComputedMagField &src) |
466 | { |
467 | |
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 | |
502 | void GlueXComputedMagField::SetFunction(std::string function) |
503 | { |
504 | |
505 | |
506 | |
507 | |
508 | |
509 | |
510 | |
511 | |
512 | |
513 | |
514 | |
515 | |
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 | |
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 | |
553 | |
554 | |
555 | if (map_opts.find(1) != map_opts.end()) { |
556 | fJanaFieldMap = new DMagneticFieldMapFineMesh(japp, run_number, |
557 | map_opts[1]); |
558 | } |
559 | |
560 | |
561 | |
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 | |
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 | |
620 | |
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 | |
629 | |
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 | |
671 | std::string GlueXComputedMagField::GetFunction() const |
672 | { |
673 | |
674 | |
675 | |
676 | return fFunction; |
677 | } |
678 | |
679 | void GlueXComputedMagField::SetMagFieldTransform(const G4AffineTransform &xform) |
680 | { |
681 | |
682 | |
683 | |
684 | fXform = xform; |
685 | fXfinv = xform.Inverse(); |
686 | } |
687 | |
688 | const G4AffineTransform& GlueXComputedMagField::GetMagFieldTransform() const |
689 | { |
690 | |
691 | |
692 | return fXform; |
693 | } |
694 | |
695 | G4ThreeVector GlueXComputedMagField::GetMagField(const G4double point[4], |
696 | G4double unit) const |
697 | { |
698 | |
699 | |
700 | |
701 | |
702 | |
703 | |
704 | |
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 | |
711 | DMagneticFieldMapPS *mapps = fJanaFieldMapPS; |
712 | |
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; |
721 | return Bvec *= fUnit / unit; |
722 | } |
723 | |
724 | void GlueXComputedMagField::GetFieldValue(const G4double point[4], |
725 | G4double *Bfield ) const |
726 | { |
727 | |
728 | |
729 | G4ThreeVector Bvec = GetMagField(point, 1); |
730 | Bfield[0] = Bvec[0]; |
731 | Bfield[1] = Bvec[1]; |
732 | Bfield[2] = Bvec[2]; |
733 | } |