ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHField3DCylindrical.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHField3DCylindrical.cc
1 #include "PHField3DCylindrical.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/tuple/tuple_comparison.hpp>
10 
11 #include <algorithm>
12 #include <cassert>
13 #include <cmath>
14 #include <cstdlib>
15 #include <iostream>
16 #include <iterator>
17 #include <set>
18 #include <utility>
19 
20 using namespace std;
21 
22 PHField3DCylindrical::PHField3DCylindrical(const string &filename, const int verb, const float magfield_rescale)
23  : PHField(verb)
24 {
25  cout << "\n================ Begin Construct Mag Field =====================" << endl;
26  cout << "\n-----------------------------------------------------------"
27  << "\n Magnetic field Module - Verbosity:" << Verbosity()
28  << "\n-----------------------------------------------------------";
29 
30  // open file
31  TFile *rootinput = TFile::Open(filename.c_str());
32  if (!rootinput)
33  {
34  cout << "\n could not open " << filename << " exiting now" << endl;
35  exit(1);
36  }
37  cout << "\n ---> "
38  "Reading the field grid from "
39  << filename << " ... " << endl;
40  rootinput->cd();
41 
42  // get root NTuple objects
43  TNtuple *field_map = (TNtuple *) gDirectory->Get("map");
44  Float_t ROOT_Z, ROOT_R, ROOT_PHI;
45  Float_t ROOT_BZ, ROOT_BR, ROOT_BPHI;
46  field_map->SetBranchAddress("z", &ROOT_Z);
47  field_map->SetBranchAddress("r", &ROOT_R);
48  field_map->SetBranchAddress("phi", &ROOT_PHI);
49  field_map->SetBranchAddress("bz", &ROOT_BZ);
50  field_map->SetBranchAddress("br", &ROOT_BR);
51  field_map->SetBranchAddress("bphi", &ROOT_BPHI);
52 
53  // get the number of entries in the tree
54  int nz, nr, nphi;
55  nz = field_map->GetEntries("z>-1e6");
56  nr = field_map->GetEntries("r>-1e6");
57  nphi = field_map->GetEntries("phi>-1e6");
58  static const int NENTRIES = field_map->GetEntries();
59 
60  // run checks on entries
61  cout << " ---> The field grid contained " << NENTRIES << " entries" << endl;
62  if (Verbosity() > 0)
63  {
64  cout << "\n NENTRIES should be the same as the following values:"
65  << "\n [ Number of values r,phi,z: "
66  << nr << " " << nphi << " " << nz << " ]! " << endl;
67  }
68 
69  if (nz != nr || nz != nphi || nr != nphi)
70  {
71  cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
72  << "\n The file you entered is not a \"table\" of values"
73  << "\n Something very likely went oh so wrong"
74  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
75  }
76 
77  // Keep track of the unique z, r, phi values in the grid using sets
78  std::set<float> z_set, r_set, phi_set;
79 
80  // Sort the entries to get rid of any stupid ordering problems...
81 
82  // We copy the TNtuple into a std::map (which is always sorted)
83  // using a 3-tuple of (z, r, phi) so it is sorted in z, then r, then
84  // phi.
85  if (Verbosity() > 0)
86  {
87  cout << " --> Sorting Entries..." << endl;
88  }
89  std::map<trio, trio> sorted_map;
90  for (int i = 0; i < field_map->GetEntries(); i++)
91  {
92  field_map->GetEntry(i);
93  trio coord_key(ROOT_Z, ROOT_R, ROOT_PHI);
94  trio field_val(ROOT_BZ, ROOT_BR, ROOT_BPHI);
95  sorted_map[coord_key] = field_val;
96 
97  z_set.insert(ROOT_Z * cm);
98  r_set.insert(ROOT_R * cm);
99  phi_set.insert(ROOT_PHI * deg);
100  }
101 
102  // couts for assurance
103  if (Verbosity() > 4)
104  {
105  map<trio, trio>::iterator it = sorted_map.begin();
106  print_map(it);
107  float last_z = it->first.get<0>();
108  for (it = sorted_map.begin(); it != sorted_map.end(); ++it)
109  {
110  if (it->first.get<0>() != last_z)
111  {
112  last_z = it->first.get<0>();
113  print_map(it);
114  }
115  }
116  }
117 
118  if (Verbosity() > 0)
119  {
120  cout << " --> Putting entries into containers... " << endl;
121  }
122 
123  // grab the minimum and maximum z values
124  minz_ = *(z_set.begin());
125  std::set<float>::iterator ziter = z_set.end();
126  --ziter;
127  maxz_ = *ziter;
128 
129  // initialize maps
130  nz = z_set.size();
131  nr = r_set.size();
132  nphi = phi_set.size();
133  r_map_.resize(nr, 0);
134  z_map_.resize(nz, 0);
135  phi_map_.resize(nphi, 0);
136 
137  std::copy(z_set.begin(), z_set.end(), z_map_.begin());
138  std::copy(phi_set.begin(), phi_set.end(), phi_map_.begin());
139  std::copy(r_set.begin(), r_set.end(), r_map_.begin());
140 
141  // initialize the field map vectors to the correct sizes
142  BFieldR_.resize(nz, vector<vector<float> >(nr, vector<float>(nphi, 0)));
143  BFieldPHI_.resize(nz, vector<vector<float> >(nr, vector<float>(nphi, 0)));
144  BFieldZ_.resize(nz, vector<vector<float> >(nr, vector<float>(nphi, 0)));
145 
146  // all of this assumes that z_prev < z , i.e. the table is ordered (as of right now)
147  unsigned int ir = 0, iphi = 0, iz = 0; // useful indexes to keep track of
148  map<trio, trio>::iterator iter = sorted_map.begin();
149  for (; iter != sorted_map.end(); ++iter)
150  {
151  // equivalent to ->GetEntry(iter)
152  float z = iter->first.get<0>() * cm;
153  float r = iter->first.get<1>() * cm;
154  float phi = iter->first.get<2>() * deg;
155  float Bz = iter->second.get<0>() * gauss;
156  float Br = iter->second.get<1>() * gauss;
157  float Bphi = iter->second.get<2>() * gauss;
158 
159  if (z > maxz_) maxz_ = z;
160  if (z < minz_) minz_ = z;
161 
162  // check for change in z value, when z changes we have a ton of updates to do
163  if (z != z_map_[iz])
164  {
165  ++iz;
166  ir = 0;
167  iphi = 0; // reset indices
168  }
169  else if (r != r_map_[ir])
170  { // check for change in r value
171  ++ir;
172  iphi = 0;
173  }
174  else if (phi != phi_map_[iphi])
175  { // change in phi value? (should be every time)
176  ++iphi;
177  }
178 
179  // shouldn't happen
180  if (iz > 0 && z < z_map_[iz - 1])
181  {
182  cout << "!!!!!!!!! Your map isn't ordered.... z: " << z << " zprev: " << z_map_[iz - 1] << endl;
183  }
184 
185  BFieldR_[iz][ir][iphi] = Br * magfield_rescale;
186  BFieldPHI_[iz][ir][iphi] = Bphi * magfield_rescale;
187  BFieldZ_[iz][ir][iphi] = Bz * magfield_rescale;
188 
189  // you can change this to check table values for correctness
190  // print_map prints the values in the root table, and the
191  // couts print the values entered into the vectors
192  if (fabs(z) < 10 && ir < 10 /*&& iphi==2*/ && Verbosity() > 3)
193  {
194  print_map(iter);
195 
196  cout << " B("
197  << r_map_[ir] << ", "
198  << phi_map_[iphi] << ", "
199  << z_map_[iz] << "): ("
200  << BFieldR_[iz][ir][iphi] << ", "
201  << BFieldPHI_[iz][ir][iphi] << ", "
202  << BFieldZ_[iz][ir][iphi] << ")" << endl;
203  }
204 
205  } // end loop over root field map file
206 
207  rootinput->Close();
208 
209  cout << "\n ---> ... read file successfully "
210  << "\n ---> Z Boundaries ~ zlow, zhigh: "
211  << minz_ / cm << "," << maxz_ / cm << " cm " << endl;
212 
213  cout << "\n================= End Construct Mag Field ======================\n"
214  << endl;
215 }
216 
217 void PHField3DCylindrical::GetFieldValue(const double point[4], double *Bfield) const
218 {
219  if (Verbosity() > 2)
220  cout << "\nPHField3DCylindrical::GetFieldValue" << endl;
221  double x = point[0];
222  double y = point[1];
223  double z = point[2];
224  double r = sqrt(x * x + y * y);
225  double phi;
226  if (x == 0)
227  {
228  phi = atan2(y, 0.00000000001);
229  }
230  else
231  {
232  phi = atan2(y, x);
233  }
234  if (phi < 0) phi += 2 * M_PI; // normalize phi to be over the range [0,2*pi]
235 
236  // Check that the point is within the defined z region (check r in a second)
237  if ((z >= minz_) && (z <= maxz_))
238  {
239  double BFieldCyl[3];
240  double cylpoint[4] = {z, r, phi, 0};
241 
242  // take <z,r,phi> location and return a vector of <Bz, Br, Bphi>
243  GetFieldCyl(cylpoint, BFieldCyl);
244 
245  // X direction of B-field ( Bx = Br*cos(phi) - Bphi*sin(phi)
246  Bfield[0] = cos(phi) * BFieldCyl[1] - sin(phi) * BFieldCyl[2]; // unit vector transformations
247 
248  // Y direction of B-field ( By = Br*sin(phi) + Bphi*cos(phi)
249  Bfield[1] = sin(phi) * BFieldCyl[1] + cos(phi) * BFieldCyl[2];
250 
251  // Z direction of B-field
252  Bfield[2] = BFieldCyl[0];
253  }
254  else
255  { // x,y,z is outside of z range of the field map
256 
257  Bfield[0] = 0.0;
258  Bfield[1] = 0.0;
259  Bfield[2] = 0.0;
260  if (Verbosity() > 2)
261  cout << "!!!!!!!!!! Field point not in defined region (outside of z bounds)" << endl;
262  }
263 
264  if (Verbosity() > 2)
265  {
266  cout << "END PHField3DCylindrical::GetFieldValue\n"
267  << " ---> {Bx, By, Bz} : "
268  << "< " << Bfield[0] << ", " << Bfield[1] << ", " << Bfield[2] << " >" << endl;
269  }
270 
271  return;
272 }
273 
274 void PHField3DCylindrical::GetFieldCyl(const double CylPoint[4], double *BfieldCyl) const
275 {
276  float z = CylPoint[0];
277  float r = CylPoint[1];
278  float phi = CylPoint[2];
279 
280  BfieldCyl[0] = 0.0;
281  BfieldCyl[1] = 0.0;
282  BfieldCyl[2] = 0.0;
283 
284  if (Verbosity() > 2)
285  cout << "GetFieldCyl@ <z,r,phi>: {" << z << "," << r << "," << phi << "}" << endl;
286 
287  if (z <= z_map_[0] || z >= z_map_[z_map_.size() - 1])
288  {
289  if (Verbosity() > 2)
290  cout << "!!!! Point not in defined region (|z| too large)" << endl;
291  return;
292  }
293  if (r < r_map_[0])
294  {
295  r = r_map_[0];
296  if (Verbosity() > 2)
297  cout << "!!!! Point not in defined region (radius too small in specific z-plane). Use min radius" << endl;
298  // return;
299  }
300  if (r > r_map_[r_map_.size() - 1])
301  {
302  if (Verbosity() > 2)
303  cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
304  return;
305  }
306 
307  vector<float>::const_iterator ziter = upper_bound(z_map_.begin(), z_map_.end(), z);
308  int z_index0 = distance(z_map_.begin(), ziter) - 1;
309  int z_index1 = z_index0 + 1;
310 
311  assert(z_index0 >= 0);
312  assert(z_index1 >= 0);
313  assert(z_index0 < (int) z_map_.size());
314  assert(z_index1 < (int) z_map_.size());
315 
316  vector<float>::const_iterator riter = upper_bound(r_map_.begin(), r_map_.end(), r);
317  int r_index0 = distance(r_map_.begin(), riter) - 1;
318  if (r_index0 >= (int) r_map_.size())
319  {
320  if (Verbosity() > 2)
321  cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
322  return;
323  }
324 
325  int r_index1 = r_index0 + 1;
326  if (r_index1 >= (int) r_map_.size())
327  {
328  if (Verbosity() > 2)
329  cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
330  return;
331  }
332 
333  assert(r_index0 >= 0);
334  assert(r_index1 >= 0);
335 
336  vector<float>::const_iterator phiiter = upper_bound(phi_map_.begin(), phi_map_.end(), phi);
337  int phi_index0 = distance(phi_map_.begin(), phiiter) - 1;
338  int phi_index1 = phi_index0 + 1;
339  if (phi_index1 >= (int) phi_map_.size())
340  phi_index1 = 0;
341 
342  assert(phi_index0 >= 0);
343  assert(phi_index0 < (int) phi_map_.size());
344  assert(phi_index1 >= 0);
345 
346  double Br000 = BFieldR_[z_index0][r_index0][phi_index0];
347  double Br001 = BFieldR_[z_index0][r_index0][phi_index1];
348  double Br010 = BFieldR_[z_index0][r_index1][phi_index0];
349  double Br011 = BFieldR_[z_index0][r_index1][phi_index1];
350  double Br100 = BFieldR_[z_index1][r_index0][phi_index0];
351  double Br101 = BFieldR_[z_index1][r_index0][phi_index1];
352  double Br110 = BFieldR_[z_index1][r_index1][phi_index0];
353  double Br111 = BFieldR_[z_index1][r_index1][phi_index1];
354 
355  double Bphi000 = BFieldPHI_[z_index0][r_index0][phi_index0];
356  double Bphi001 = BFieldPHI_[z_index0][r_index0][phi_index1];
357  double Bphi010 = BFieldPHI_[z_index0][r_index1][phi_index0];
358  double Bphi011 = BFieldPHI_[z_index0][r_index1][phi_index1];
359  double Bphi100 = BFieldPHI_[z_index1][r_index0][phi_index0];
360  double Bphi101 = BFieldPHI_[z_index1][r_index0][phi_index1];
361  double Bphi110 = BFieldPHI_[z_index1][r_index1][phi_index0];
362  double Bphi111 = BFieldPHI_[z_index1][r_index1][phi_index1];
363 
364  double Bz000 = BFieldZ_[z_index0][r_index0][phi_index0];
365  double Bz001 = BFieldZ_[z_index0][r_index0][phi_index1];
366  double Bz100 = BFieldZ_[z_index1][r_index0][phi_index0];
367  double Bz101 = BFieldZ_[z_index1][r_index0][phi_index1];
368  double Bz010 = BFieldZ_[z_index0][r_index1][phi_index0];
369  double Bz110 = BFieldZ_[z_index1][r_index1][phi_index0];
370  double Bz011 = BFieldZ_[z_index0][r_index1][phi_index1];
371  double Bz111 = BFieldZ_[z_index1][r_index1][phi_index1];
372 
373  double zweight = z - z_map_[z_index0];
374  double zspacing = z_map_[z_index1] - z_map_[z_index0];
375  zweight /= zspacing;
376 
377  double rweight = r - r_map_[r_index0];
378  double rspacing = r_map_[r_index1] - r_map_[r_index0];
379  rweight /= rspacing;
380 
381  double phiweight = phi - phi_map_[phi_index0];
382  double phispacing = phi_map_[phi_index1] - phi_map_[phi_index0];
383  if (phi_index1 == 0)
384  phispacing += 2 * M_PI;
385  phiweight /= phispacing;
386 
387  // Z direction of B-field
388  BfieldCyl[0] =
389  (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Bz000 + phiweight * Bz001) +
390  rweight * ((1 - phiweight) * Bz010 + phiweight * Bz011)) +
391  zweight * ((1 - rweight) * ((1 - phiweight) * Bz100 + phiweight * Bz101) +
392  rweight * ((1 - phiweight) * Bz110 + phiweight * Bz111));
393 
394  // R direction of B-field
395  BfieldCyl[1] =
396  (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Br000 + phiweight * Br001) +
397  rweight * ((1 - phiweight) * Br010 + phiweight * Br011)) +
398  zweight * ((1 - rweight) * ((1 - phiweight) * Br100 + phiweight * Br101) +
399  rweight * ((1 - phiweight) * Br110 + phiweight * Br111));
400 
401  // PHI Direction of B-field
402  BfieldCyl[2] =
403  (1 - zweight) * ((1 - rweight) * ((1 - phiweight) * Bphi000 + phiweight * Bphi001) +
404  rweight * ((1 - phiweight) * Bphi010 + phiweight * Bphi011)) +
405  zweight * ((1 - rweight) * ((1 - phiweight) * Bphi100 + phiweight * Bphi101) +
406  rweight * ((1 - phiweight) * Bphi110 + phiweight * Bphi111));
407 
408  // cout << "wr: " << rweight << " wz: " << zweight << " wphi: " << phiweight << endl;
409  // cout << "Bz000: " << Bz000 << endl
410  // << "Bz001: " << Bz001 << endl
411  // << "Bz010: " << Bz010 << endl
412  // << "Bz011: " << Bz011 << endl
413  // << "Bz100: " << Bz100 << endl
414  // << "Bz101: " << Bz101 << endl
415  // << "Bz110: " << Bz110 << endl
416  // << "Bz111: " << Bz111 << endl
417  // << "Bz: " << BfieldCyl[0] << endl << endl;
418 
419  if (Verbosity() > 2)
420  {
421  cout << "End GFCyl Call: <bz,br,bphi> : {"
422  << BfieldCyl[0] / gauss << "," << BfieldCyl[1] / gauss << "," << BfieldCyl[2] / gauss << "}"
423  << endl;
424  }
425 
426  return;
427 }
428 
429 // a binary search algorithm that puts the location that "key" would be, into index...
430 // it returns true if key was found, and false if not.
431 bool PHField3DCylindrical::bin_search(const vector<float> &vec, unsigned start, unsigned end, const float &key, unsigned &index) const
432 {
433  // Termination condition: start index greater than end index
434  if (start > end)
435  {
436  index = start;
437  return false;
438  }
439 
440  // Find the middle element of the vector and use that for splitting
441  // the array into two pieces.
442  unsigned middle = start + ((end - start) / 2);
443 
444  if (vec[middle] == key)
445  {
446  index = middle;
447  return true;
448  }
449  else if (vec[middle] > key)
450  {
451  return bin_search(vec, start, middle - 1, key, index);
452  }
453  return bin_search(vec, middle + 1, end, key, index);
454 }
455 
456 // debug function to print key/value pairs in map
457 void PHField3DCylindrical::print_map(map<trio, trio>::iterator &it) const
458 {
459  cout << " Key: <"
460  << it->first.get<0>() * cm << ","
461  << it->first.get<1>() * cm << ","
462  << it->first.get<2>() * deg << ">"
463 
464  << " Value: <"
465  << it->second.get<0>() * gauss << ","
466  << it->second.get<1>() * gauss << ","
467  << it->second.get<2>() * gauss << ">\n";
468 }