3 #include <TDirectory.h>
7 #include <Geant4/G4SystemOfUnits.hh>
9 #include <boost/stacktrace.hpp>
10 #include <boost/tuple/tuple.hpp>
11 #include <boost/tuple/tuple_comparison.hpp>
22 typedef boost::tuple<double, double, double>
trio;
23 std::map<boost::tuple<double, double, double>, boost::tuple<double, double, double> >
fieldmap;
31 for (
int i = 0; i < 2; i++)
33 for (
int j = 0; j < 2; j++)
35 for (
int k = 0;
k < 2;
k++)
37 for (
int l = 0; l < 3; l++)
39 xyz[i][j][
k][l] = NAN;
45 std::cout <<
"\n================ Begin Construct Mag Field =====================" << std::endl;
46 std::cout <<
"\n-----------------------------------------------------------"
47 <<
"\n Magnetic field Module - Verbosity:"
48 <<
"\n-----------------------------------------------------------";
51 TFile *rootinput = TFile::Open(
filename.c_str());
54 std::cout <<
"\n could not open " <<
filename <<
" exiting now" << std::endl;
57 std::cout <<
"\n ---> "
58 "Reading the field grid from "
63 TNtuple *field_map = (TNtuple *) gDirectory->Get(
"fieldmap");
65 Float_t ROOT_BX, ROOT_BY, ROOT_BZ;
66 field_map->SetBranchAddress(
"x", &ROOT_X);
67 field_map->SetBranchAddress(
"y", &ROOT_Y);
68 field_map->SetBranchAddress(
"z", &ROOT_Z);
69 field_map->SetBranchAddress(
"bx", &ROOT_BX);
70 field_map->SetBranchAddress(
"by", &ROOT_BY);
71 field_map->SetBranchAddress(
"bz", &ROOT_BZ);
73 for (
int i = 0; i < field_map->GetEntries(); i++)
75 field_map->GetEntry(i);
76 trio coord_key(ROOT_X *
cm, ROOT_Y * cm, ROOT_Z * cm);
77 trio field_val(ROOT_BX *
tesla * magfield_rescale, ROOT_BY *
tesla * magfield_rescale, ROOT_BZ *
tesla * magfield_rescale);
78 xvals.insert(ROOT_X * cm);
79 yvals.insert(ROOT_Y * cm);
80 zvals.insert(ROOT_Z * cm);
81 fieldmap[coord_key] = field_val;
83 xmin = *(xvals.begin());
84 xmax = *(xvals.rbegin());
86 ymin = *(yvals.begin());
87 ymax = *(yvals.rbegin());
90 std::cout <<
"PHField3DCartesian: Compiler bug!!!!!!!! Do not use inlining!!!!!!" << std::endl;
91 std::cout <<
"exiting now - recompile with -fno-inline" << std::endl;
95 zmin = *(zvals.begin());
96 zmax = *(zvals.rbegin());
104 std::cout <<
"\n================= End Construct Mag Field ======================\n"
112 std::cout <<
"PHField3DCartesian: cache hits: " <<
cache_hits
120 static double xsav = -1000000.;
121 static double ysav = -1000000.;
122 static double zsav = -1000000.;
133 static int ifirst = 0;
136 std::cout <<
"PHField3DCartesian::GetFieldValue: "
137 <<
"Invalid coordinates: "
141 <<
" bailing out returning zero bfield"
143 std::cout <<
"previous point: "
144 <<
"x: " << xsav /
cm
145 <<
", y: " << ysav /
cm
146 <<
", z: " << zsav /
cm
148 std::cout <<
"Here is the stacktrace: " << std::endl;
149 std::cout << boost::stacktrace::stacktrace();
150 std::cout <<
"This is not a segfault. Check the stacktrace for the guilty party (typically #2)" << std::endl;
159 if (point[0] <
xmin || point[0] >
xmax ||
160 point[1] <
ymin || point[1] >
ymax ||
165 std::set<double>::const_iterator
it = xvals.lower_bound(x);
166 if (it == xvals.begin())
168 std::cout <<
"x too small - outside range: " << x /
cm << std::endl;
176 it = yvals.lower_bound(y);
177 if (it == yvals.begin())
179 std::cout <<
"y too small - outside range: " << y /
cm << std::endl;
187 it = zvals.lower_bound(z);
188 if (it == zvals.begin())
190 std::cout <<
"z too small - outside range: " << z /
cm << std::endl;
207 std::map<boost::tuple<double, double, double>, boost::tuple<double, double, double> >::const_iterator magval;
209 for (
int i = 0; i < 2; i++)
211 for (
int j = 0; j < 2; j++)
213 for (
int k = 0;
k < 2;
k++)
215 key = boost::make_tuple(xkey[i], ykey[j], zkey[
k]);
216 magval = fieldmap.find(key);
217 if (magval == fieldmap.end())
219 std::cout <<
"could not locate key, x: " << xkey[i] /
cm
220 <<
", y: " << ykey[j] /
cm
221 <<
", z: " << zkey[
k] /
cm << std::endl;
224 xyz[i][j][
k][0] = (magval->first).get<0>();
225 xyz[i][j][
k][1] = (magval->first).get<1>();
226 xyz[i][j][
k][2] = (magval->first).get<2>();
227 bf[i][j][
k][0] = (magval->second).get<0>();
228 bf[i][j][
k][1] = (magval->second).get<1>();
229 bf[i][j][
k][2] = (magval->second).get<2>();
232 std::cout <<
"read x/y/z: " <<
xyz[i][j][
k][0] /
cm <<
"/"
233 <<
xyz[i][j][
k][1] /
cm <<
"/"
234 <<
xyz[i][j][
k][2] /
cm <<
" bx/by/bz: "
237 <<
bf[i][j][
k][2] /
tesla << std::endl;
249 double xinblock = point[0] - xkey[1];
250 double yinblock = point[1] - ykey[1];
251 double zinblock = point[2] - zkey[1];
259 std::cout <<
"x/y/z inblock: " << xinblock /
cm <<
"/" << yinblock /
cm <<
"/" << zinblock /
cm << std::endl;
260 std::cout <<
"x/y/z fraction: " << fractionx <<
"/" << fractiony <<
"/" << fractionz << std::endl;
275 for (
int i = 0; i < 3; i++)
277 Bfield[i] =
bf[0][0][0][i] * fractionx * fractiony * fractionz +
278 bf[1][0][0][i] * (1. - fractionx) * fractiony * fractionz +
279 bf[0][1][0][i] * fractionx * (1. - fractiony) * fractionz +
280 bf[0][0][1][i] * fractionx * fractiony * (1. - fractionz) +
281 bf[1][0][1][i] * (1. - fractionx) * fractiony * (1. - fractionz) +
282 bf[0][1][1][i] * fractionx * (1. - fractiony) * (1. - fractionz) +
283 bf[1][1][0][i] * (1. - fractionx) * (1. - fractiony) * fractionz +
284 bf[1][1][1][i] * (1. - fractionx) * (1. - fractiony) * (1. - fractionz);