EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHField3DCartesian.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHField3DCartesian.cc
1 #include "PHField3DCartesian.h"
2 
3 #include <TDirectory.h> // for TDirectory, gDirectory
4 #include <TFile.h>
5 #include <TNtuple.h>
6 
7 #include <Geant4/G4SystemOfUnits.hh>
8 
9 #include <boost/stacktrace.hpp>
10 #include <boost/tuple/tuple.hpp>
11 #include <boost/tuple/tuple_comparison.hpp>
12 
13 #include <cassert>
14 #include <cmath>
15 #include <cstdlib>
16 #include <iostream>
17 #include <iterator>
18 #include <map>
19 #include <set>
20 #include <utility>
21 
22 typedef boost::tuple<double, double, double> trio;
23 std::map<boost::tuple<double, double, double>, boost::tuple<double, double, double> > fieldmap;
24 std::set<double> xvals;
25 std::set<double> yvals;
26 std::set<double> zvals;
27 
28 PHField3DCartesian::PHField3DCartesian(const std::string &fname, const float magfield_rescale)
29  : filename(fname)
30 {
31  for (int i = 0; i < 2; i++)
32  {
33  for (int j = 0; j < 2; j++)
34  {
35  for (int k = 0; k < 2; k++)
36  {
37  for (int l = 0; l < 3; l++)
38  {
39  xyz[i][j][k][l] = NAN;
40  bf[i][j][k][l] = NAN;
41  }
42  }
43  }
44  }
45  std::cout << "\n================ Begin Construct Mag Field =====================" << std::endl;
46  std::cout << "\n-----------------------------------------------------------"
47  << "\n Magnetic field Module - Verbosity:"
48  << "\n-----------------------------------------------------------";
49 
50  // open file
51  TFile *rootinput = TFile::Open(filename.c_str());
52  if (!rootinput)
53  {
54  std::cout << "\n could not open " << filename << " exiting now" << std::endl;
55  exit(1);
56  }
57  std::cout << "\n ---> "
58  "Reading the field grid from "
59  << filename << " ... " << std::endl;
60  rootinput->cd();
61 
62  // get root NTuple objects
63  TNtuple *field_map = (TNtuple *) gDirectory->Get("fieldmap");
64  Float_t ROOT_X, ROOT_Y, ROOT_Z;
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);
72 
73  for (int i = 0; i < field_map->GetEntries(); i++)
74  {
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;
82  }
83  xmin = *(xvals.begin());
84  xmax = *(xvals.rbegin());
85 
86  ymin = *(yvals.begin());
87  ymax = *(yvals.rbegin());
88  if (ymin != xmin || ymax != xmax)
89  {
90  std::cout << "PHField3DCartesian: Compiler bug!!!!!!!! Do not use inlining!!!!!!" << std::endl;
91  std::cout << "exiting now - recompile with -fno-inline" << std::endl;
92  exit(1);
93  }
94 
95  zmin = *(zvals.begin());
96  zmax = *(zvals.rbegin());
97 
98  xstepsize = (xmax - xmin) / (xvals.size() - 1);
99  ystepsize = (ymax - ymin) / (yvals.size() - 1);
100  zstepsize = (zmax - zmin) / (zvals.size() - 1);
101 
102  rootinput->Close();
103 
104  std::cout << "\n================= End Construct Mag Field ======================\n"
105  << std::endl;
106 }
107 
108 PHField3DCartesian::~PHField3DCartesian()
109 {
110  if (Verbosity() > 0)
111  {
112  std::cout << "PHField3DCartesian: cache hits: " << cache_hits
113  << " cache misses: " << cache_misses
114  << std::endl;
115  }
116 }
117 
118 void PHField3DCartesian::GetFieldValue(const double point[4], double *Bfield) const
119 {
120  static double xsav = -1000000.;
121  static double ysav = -1000000.;
122  static double zsav = -1000000.;
123 
124  double x = point[0];
125  double y = point[1];
126  double z = point[2];
127 
128  Bfield[0] = 0.0;
129  Bfield[1] = 0.0;
130  Bfield[2] = 0.0;
131  if (!std::isfinite(x) || !std::isfinite(y) || !std::isfinite(z))
132  {
133  static int ifirst = 0;
134  if (ifirst < 10)
135  {
136  std::cout << "PHField3DCartesian::GetFieldValue: "
137  << "Invalid coordinates: "
138  << "x: " << x / cm
139  << ", y: " << y / cm
140  << ", z: " << z / cm
141  << " bailing out returning zero bfield"
142  << std::endl;
143  std::cout << "previous point: "
144  << "x: " << xsav / cm
145  << ", y: " << ysav / cm
146  << ", z: " << zsav / cm
147  << std::endl;
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;
151  ifirst++;
152  }
153  return;
154  }
155  xsav = x;
156  ysav = y;
157  zsav = z;
158 
159  if (point[0] < xmin || point[0] > xmax ||
160  point[1] < ymin || point[1] > ymax ||
161  point[2] < zmin || point[2] > zmax)
162  {
163  return;
164  }
165  std::set<double>::const_iterator it = xvals.lower_bound(x);
166  if (it == xvals.begin())
167  {
168  std::cout << "x too small - outside range: " << x / cm << std::endl;
169  return;
170  }
171  double xkey[2];
172  xkey[0] = *it;
173  --it;
174  xkey[1] = *it;
175 
176  it = yvals.lower_bound(y);
177  if (it == yvals.begin())
178  {
179  std::cout << "y too small - outside range: " << y / cm << std::endl;
180  return;
181  }
182  double ykey[2];
183  ykey[0] = *it;
184  --it;
185  ykey[1] = *it;
186 
187  it = zvals.lower_bound(z);
188  if (it == zvals.begin())
189  {
190  std::cout << "z too small - outside range: " << z / cm << std::endl;
191  return;
192  }
193  double zkey[2];
194  zkey[0] = *it;
195  --it;
196  zkey[1] = *it;
197 
198  if (xkey_save != xkey[0] ||
199  ykey_save != ykey[0] ||
200  zkey_save != zkey[0])
201  {
202  cache_misses++;
203  xkey_save = xkey[0];
204  ykey_save = ykey[0];
205  zkey_save = zkey[0];
206 
207  std::map<boost::tuple<double, double, double>, boost::tuple<double, double, double> >::const_iterator magval;
208  trio key;
209  for (int i = 0; i < 2; i++)
210  {
211  for (int j = 0; j < 2; j++)
212  {
213  for (int k = 0; k < 2; k++)
214  {
215  key = boost::make_tuple(xkey[i], ykey[j], zkey[k]);
216  magval = fieldmap.find(key);
217  if (magval == fieldmap.end())
218  {
219  std::cout << "could not locate key, x: " << xkey[i] / cm
220  << ", y: " << ykey[j] / cm
221  << ", z: " << zkey[k] / cm << std::endl;
222  return;
223  }
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>();
230  if (Verbosity() > 0)
231  {
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: "
235  << bf[i][j][k][0] / tesla << "/"
236  << bf[i][j][k][1] / tesla << "/"
237  << bf[i][j][k][2] / tesla << std::endl;
238  }
239  }
240  }
241  }
242  }
243  else
244  {
245  cache_hits++;
246  }
247 
248  // how far are we away from the reference point
249  double xinblock = point[0] - xkey[1];
250  double yinblock = point[1] - ykey[1];
251  double zinblock = point[2] - zkey[1];
252  // normalize distance to step size
253  double fractionx = xinblock / xstepsize;
254  double fractiony = yinblock / ystepsize;
255  double fractionz = zinblock / zstepsize;
256  if (Verbosity() > 0)
257  {
258  std::cout << "x/y/z stepsize: " << xstepsize / cm << "/" << ystepsize / cm << "/" << zstepsize / cm << std::endl;
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;
261  }
262 
263  // linear extrapolation in cube:
264 
265  //Vxyz =
266  //V000 * x * y * z +
267  //V100 * (1 - x) * y * z +
268  //V010 * x * (1 - y) * z +
269  //V001 * x y * (1 - z) +
270  //V101 * (1 - x) * y * (1 - z) +
271  //V011 * x * (1 - y) * (1 - z) +
272  //V110 * (1 - x) * (1 - y) * z +
273  //V111 * (1 - x) * (1 - y) * (1 - z)
274 
275  for (int i = 0; i < 3; i++)
276  {
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);
285  }
286 
287  return;
288 }