Newer
Older
#include "GenericBFieldMapBrBz.h"
#include "FieldMapFileProvider.h"
#include <cmath>
#include "DD4hep/DD4hepUnits.h"
GenericBFieldMapBrBz::GenericBFieldMapBrBz()
: m_provider(nullptr) {
type = dd4hep::CartesianField::MAGNETIC;
}
void GenericBFieldMapBrBz::fieldComponents(const double* pos, double* field) {
if (!m_provider) {
std::string error_msg = "[ERROR] GenericBFieldMapBrBz: No provider! ";
throw std::runtime_error(error_msg);
}
if (m_length_unit<=0||m_bfield_unit<=0) {
std::string error_msg = "[ERROR] GenericBFieldMapBrBz: Not set units or error (<0)! ";
throw std::runtime_error(error_msg);
}
double x = pos[0] / m_length_unit; // convert to length unit from input
double y = pos[1] / m_length_unit;
double z = pos[2] / m_length_unit;
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
double r = sqrt(x*x+y*y);
double phi = atan2(y, x);
// std::cout << " r: " << r
// << " x: " << x
// << " y: " << y
// << " z: " << z
// << " field Bx: " << field[0]
// << " field By: " << field[1]
// << " field Bz: " << field[2]
// << std::endl;
// get the point r0/z0, r1/z1
// -- normalized r/z at the bin
double rn = 0.0;
double zn = 0.0;
int ir0 = 0;
int iz0 = 0;
ir0 = m_provider->rBinIdx(r, rn);
iz0 = m_provider->zBinIdx(z, zn);
// not in the valid return
if (ir0 < 0 || iz0 < 0) {
// std::cout << "SKIP due to "
// << " r: " << r
// << " x: " << x
// << " y: " << y
// << " z: " << z
// << " ir0: " << ir0
// << " iz0: " << iz0
// << " rn: " << rn
// << " zn: " << zn
// << std::endl;
return;
}
// std::cout << " r: " << r
// << " x: " << x
// << " y: " << y
// << " z: " << z
// << " ir0: " << ir0
// << " iz0: " << iz0
// << " rn: " << rn
// << " zn: " << zn
// << std::endl;
int ir1 = ir0 + 1;
int iz1 = iz0 + 1;
// get the Br/Bz at four points
double Br_r0z0 = 0.0; double Bz_r0z0 = 0.0;
double Br_r1z0 = 0.0; double Bz_r1z0 = 0.0;
double Br_r0z1 = 0.0; double Bz_r0z1 = 0.0;
double Br_r1z1 = 0.0; double Bz_r1z1 = 0.0;
// calculate the field at r/z
m_provider->access(ir0, iz0, Br_r0z0, Bz_r0z0);
m_provider->access(ir1, iz0, Br_r1z0, Bz_r1z0);
m_provider->access(ir0, iz1, Br_r0z1, Bz_r0z1);
m_provider->access(ir1, iz1, Br_r1z1, Bz_r1z1);
double Br = (1.0 - rn) * (1.0 - zn) * Br_r0z0
+ rn * (1.0 - zn) * Br_r1z0
+ (1.0 - rn) * zn * Br_r0z1
+ rn * zn * Br_r1z1;
double Bz = (1.0 - rn) * (1.0 - zn) * Bz_r0z0
+ rn * (1.0 - zn) * Bz_r1z0
+ (1.0 - rn) * zn * Bz_r0z1
+ rn * zn * Bz_r1z1;
field[0] += Br*cos(phi)*m_bfield_unit; // convert to input unit
field[1] += Br*sin(phi)*m_bfield_unit;
field[2] += Bz*m_bfield_unit;
return;
}
void GenericBFieldMapBrBz::init_provider(const std::string& provider, const std::string& url) {
if (provider == "file") {
std::cout << "Initialize provider with file. " << std::endl;
m_provider = new FieldMapFileProvider(url);
} else if (provider == "db") {
std::cout << "Initialize provider with file. " << std::endl;
} else {
std::string error_msg = "[ERROR] GenericBFieldMapBrBz: Unknown provider: " + provider;
throw std::runtime_error(error_msg);
}
}
void GenericBFieldMapBrBz::init_unit(double l, double b) {
m_length_unit = l;
m_bfield_unit = b;
std::cout << "Initialize units to l = " << m_length_unit/dd4hep::m << " m" << ", b = " << m_bfield_unit/dd4hep::tesla << " tesla" << std::endl;
}