ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHField2D.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHField2D.cc
1 #include "PHField2D.h"
2 
3 //root framework
4 #include <TDirectory.h>
5 #include <TFile.h>
6 #include <TNtuple.h>
7 
8 #include <Geant4/G4SystemOfUnits.hh>
9 
10 #include <TSystem.h>
11 
12 #include <boost/tuple/tuple_comparison.hpp>
13 
14 #include <algorithm>
15 #include <cmath>
16 #include <cstdlib>
17 #include <iostream>
18 #include <iterator>
19 #include <set>
20 #include <utility>
21 
22 using namespace std;
23 
24 PHField2D::PHField2D(const string &filename, const int verb, const float magfield_rescale)
25  : PHField(verb)
26  , r_index0_cache(0)
27  , r_index1_cache(0)
28  , z_index0_cache(0)
29  , z_index1_cache(0)
30 {
31  if (Verbosity() > 0)
32  {
33  cout << " ------------- PHField2D::PHField2D() ------------------" << endl;
34  }
35  // open file
36  TFile *rootinput = TFile::Open(filename.c_str());
37  if (!rootinput)
38  {
39  cout << " could not open " << filename << " exiting now" << endl;
40  gSystem->Exit(1);
41  exit(1);
42  }
43  if (Verbosity() > 0) cout << " Field grid file: " << filename << endl;
44  rootinput->cd();
45 
46  Float_t ROOT_Z, ROOT_R;
47  Float_t ROOT_BZ, ROOT_BR;
48  // get root NTuple objects
49  TNtuple *field_map = (TNtuple *) gDirectory->Get("fieldmap");
50  if (field_map)
51  {
53  }
54  else
55  {
56  field_map = (TNtuple *) gDirectory->Get("map");
57  if (!field_map)
58  {
59  cout << "PHField2D: could not locate ntuple of name map or fieldmap, exiting now" << endl;
60  exit(1);
61  }
63  }
64  field_map->SetBranchAddress("z", &ROOT_Z);
65  field_map->SetBranchAddress("r", &ROOT_R);
66  field_map->SetBranchAddress("bz", &ROOT_BZ);
67  field_map->SetBranchAddress("br", &ROOT_BR);
68 
69  // get the number of entries in the tree
70  int nz, nr;
71  nz = field_map->GetEntries("z>-1e6");
72  nr = field_map->GetEntries("r>-1e6");
73  static const int NENTRIES = field_map->GetEntries();
74 
75  // run checks on entries
76  if (Verbosity() > 0) cout << " The field grid contained " << NENTRIES << " entries" << endl;
77  if (Verbosity() > 1)
78  {
79  cout << "\n NENTRIES should be the same as the following values:"
80  << "\n [ Number of values r,z: "
81  << nr << " " << nz << " ]! " << endl;
82  }
83 
84  if (nz != nr)
85  {
86  cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
87  << "\n The file you entered is not a \"table\" of values"
88  << "\n Something very likely went oh so wrong"
89  << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
90  }
91 
92  // Keep track of the unique z, r, phi values in the grid using sets
93  std::set<float> z_set, r_set;
94 
95  // Sort the entries to get rid of any stupid ordering problems...
96 
97  // We copy the TNtuple into a std::map (which is always sorted)
98  // using a 3-tuple of (z, r, phi) so it is sorted in z, then r, then
99  // phi.
100  if (Verbosity() > 1)
101  {
102  cout << " --> Sorting Entries..." << endl;
103  }
104  std::map<trio, trio> sorted_map;
105  for (int i = 0; i < field_map->GetEntries(); i++)
106  {
107  field_map->GetEntry(i);
108  trio coord_key(ROOT_Z, ROOT_R);
109  trio field_val(ROOT_BZ, ROOT_BR);
110  sorted_map[coord_key] = field_val;
111 
112  z_set.insert(ROOT_Z * cm);
113  r_set.insert(ROOT_R * cm);
114  }
115 
116  // couts for assurance
117  if (Verbosity() > 4)
118  {
119  map<trio, trio>::iterator it = sorted_map.begin();
120  print_map(it);
121  float last_z = it->first.get<0>();
122  for (it = sorted_map.begin(); it != sorted_map.end(); ++it)
123  {
124  if (it->first.get<0>() != last_z)
125  {
126  last_z = it->first.get<0>();
127  print_map(it);
128  }
129  }
130  }
131 
132  if (Verbosity() > 1)
133  {
134  cout << " --> Putting entries into containers... " << endl;
135  }
136 
137  // grab the minimum and maximum z values
138  minz_ = *(z_set.begin());
139  std::set<float>::iterator ziter = z_set.end();
140  --ziter;
141  maxz_ = *ziter;
142 
143  // initialize maps
144  nz = z_set.size();
145  nr = r_set.size();
146  r_map_.resize(nr, 0);
147  z_map_.resize(nz, 0);
148 
149  std::copy(z_set.begin(), z_set.end(), z_map_.begin());
150  std::copy(r_set.begin(), r_set.end(), r_map_.begin());
151 
152  // initialize the field map vectors to the correct sizes
153  BFieldR_.resize(nz, vector<float>(nr, 0));
154  BFieldZ_.resize(nz, vector<float>(nr, 0));
155 
156  // all of this assumes that z_prev < z , i.e. the table is ordered (as of right now)
157  unsigned int ir = 0, iz = 0; // useful indexes to keep track of
158  map<trio, trio>::iterator iter = sorted_map.begin();
159  for (; iter != sorted_map.end(); ++iter)
160  {
161  // equivalent to ->GetEntry(iter)
162  float z = iter->first.get<0>() * cm;
163  float r = iter->first.get<1>() * cm;
164  float Bz = iter->second.get<0>() * magfield_unit;
165  float Br = iter->second.get<1>() * magfield_unit;
166 
167  if (z > maxz_) maxz_ = z;
168  if (z < minz_) minz_ = z;
169 
170  // check for change in z value, when z changes we have a ton of updates to do
171  if (z != z_map_[iz])
172  {
173  ++iz;
174  ir = 0;
175  }
176  else if (r != r_map_[ir]) // check for change in r value
177  {
178  ++ir;
179  }
180 
181  // shouldn't happen
182  if (iz > 0 && z < z_map_[iz - 1])
183  {
184  cout << "!!!!!!!!! Your map isn't ordered.... z: " << z << " zprev: " << z_map_[iz - 1] << endl;
185  }
186 
187  BFieldR_[iz][ir] = Br * magfield_rescale;
188  BFieldZ_[iz][ir] = Bz * magfield_rescale;
189 
190  // you can change this to check table values for correctness
191  // print_map prints the values in the root table, and the
192  // couts print the values entered into the vectors
193  if (fabs(z) < 10 && ir < 10 /*&& iphi==2*/ && Verbosity() > 3)
194  {
195  print_map(iter);
196 
197  cout << " B("
198  << r_map_[ir] << ", "
199  << z_map_[iz] << "): ("
200  << BFieldR_[iz][ir] << ", "
201  << BFieldZ_[iz][ir] << ")" << endl;
202  }
203 
204  } // end loop over root field map file
205 
206  rootinput->Close();
207 
208  if (Verbosity() > 0) cout << " Mag field z boundaries (min,max): (" << minz_ / cm << ", " << maxz_ / cm << ") cm" << endl;
209  if (Verbosity() > 0) cout << " Mag field r max boundary: " << r_map_.back() / cm << " cm" << endl;
210 
211  if (Verbosity() > 0)
212  cout << " -----------------------------------------------------------" << endl;
213 }
214 
215 void PHField2D::GetFieldValue(const double point[4], double *Bfield) const
216 {
217  if (Verbosity() > 2)
218  cout << "\nPHField2D::GetFieldValue" << endl;
219  double x = point[0];
220  double y = point[1];
221  double z = point[2];
222  double r = sqrt(x * x + y * y);
223  double phi;
224  phi = atan2(y, x);
225  if (phi < 0) phi += 2 * M_PI; // normalize phi to be over the range [0,2*pi]
226 
227  // Check that the point is within the defined z region (check r in a second)
228  if ((z >= minz_) && (z <= maxz_))
229  {
230  double BFieldCyl[3];
231  double cylpoint[4] = {z, r, phi, 0};
232 
233  // take <z,r,phi> location and return a vector of <Bz, Br, Bphi>
234  GetFieldCyl(cylpoint, BFieldCyl);
235 
236  // X direction of B-field ( Bx = Br*cos(phi) - Bphi*sin(phi)
237  Bfield[0] = cos(phi) * BFieldCyl[1] - sin(phi) * BFieldCyl[2]; // unit vector transformations
238 
239  // Y direction of B-field ( By = Br*sin(phi) + Bphi*cos(phi)
240  Bfield[1] = sin(phi) * BFieldCyl[1] + cos(phi) * BFieldCyl[2];
241 
242  // Z direction of B-field
243  Bfield[2] = BFieldCyl[0];
244  }
245  else // x,y,z is outside of z range of the field map
246  {
247  Bfield[0] = 0.0;
248  Bfield[1] = 0.0;
249  Bfield[2] = 0.0;
250  if (Verbosity() > 2)
251  cout << "!!!!!!!!!! Field point not in defined region (outside of z bounds)" << endl;
252  }
253 
254  if (Verbosity() > 2)
255  {
256  cout << "END PHField2D::GetFieldValue\n"
257  << " ---> {Bx, By, Bz} : "
258  << "< " << Bfield[0] << ", " << Bfield[1] << ", " << Bfield[2] << " >" << endl;
259  }
260 
261  return;
262 }
263 
264 void PHField2D::GetFieldCyl(const double CylPoint[4], double *BfieldCyl) const
265 {
266  float z = CylPoint[0];
267  float r = CylPoint[1];
268 
269  BfieldCyl[0] = 0.0;
270  BfieldCyl[1] = 0.0;
271  BfieldCyl[2] = 0.0;
272 
273  if (Verbosity() > 2)
274  cout << "GetFieldCyl@ <z,r>: {" << z << "," << r << "}" << endl;
275 
276  if (z < z_map_[0] || z > z_map_[z_map_.size() - 1])
277  {
278  if (Verbosity() > 2)
279  cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
280  return;
281  }
282 
283  // since GEANT4 looks up the field ~95% of the time in the same voxel
284  // between subsequent calls, we can save on the expense of the upper_bound
285  // lookup (~10-15% of central event run time) with some caching between calls
286 
287  unsigned int r_index0 = r_index0_cache;
288  unsigned int r_index1 = r_index1_cache;
289 
290  if (!((r > r_map_[r_index0]) && (r < r_map_[r_index1])))
291  {
292  // if miss cached r values, search through the lookup table
293  vector<float>::const_iterator riter = upper_bound(r_map_.begin(), r_map_.end(), r);
294  r_index0 = distance(r_map_.begin(), riter) - 1;
295  if (r_index0 >= r_map_.size())
296  {
297  if (Verbosity() > 2)
298  cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
299  return;
300  }
301 
302  r_index1 = r_index0 + 1;
303  if (r_index1 >= r_map_.size())
304  {
305  if (Verbosity() > 2)
306  cout << "!!!! Point not in defined region (radius too large in specific z-plane)" << endl;
307  return;
308  }
309 
310  // update cache
311  r_index0_cache = r_index0;
312  r_index1_cache = r_index1;
313  }
314 
315  unsigned int z_index0 = z_index0_cache;
316  unsigned int z_index1 = z_index1_cache;
317 
318  if (!((z > z_map_[z_index0]) && (z < z_map_[z_index1])))
319  {
320  // if miss cached z values, search through the lookup table
321  vector<float>::const_iterator ziter = upper_bound(z_map_.begin(), z_map_.end(), z);
322  z_index0 = distance(z_map_.begin(), ziter) - 1;
323  z_index1 = z_index0 + 1;
324  if (z_index1 >= z_map_.size())
325  {
326  if (Verbosity() > 2)
327  cout << "!!!! Point not in defined region (z too large in specific r-plane)" << endl;
328  return;
329  }
330 
331  // update cache
332  z_index0_cache = z_index0;
333  z_index1_cache = z_index1;
334  }
335 
336  double Br000 = BFieldR_[z_index0][r_index0];
337  double Br010 = BFieldR_[z_index0][r_index1];
338  double Br100 = BFieldR_[z_index1][r_index0];
339  double Br110 = BFieldR_[z_index1][r_index1];
340 
341  double Bz000 = BFieldZ_[z_index0][r_index0];
342  double Bz100 = BFieldZ_[z_index1][r_index0];
343  double Bz010 = BFieldZ_[z_index0][r_index1];
344  double Bz110 = BFieldZ_[z_index1][r_index1];
345 
346  double zweight = z - z_map_[z_index0];
347  double zspacing = z_map_[z_index1] - z_map_[z_index0];
348  zweight /= zspacing;
349 
350  double rweight = r - r_map_[r_index0];
351  double rspacing = r_map_[r_index1] - r_map_[r_index0];
352  rweight /= rspacing;
353 
354  // Z direction of B-field
355  BfieldCyl[0] =
356  (1 - zweight) * ((1 - rweight) * Bz000 +
357  rweight * Bz010) +
358  zweight * ((1 - rweight) * Bz100 +
359  rweight * Bz110);
360 
361  // R direction of B-field
362  BfieldCyl[1] =
363  (1 - zweight) * ((1 - rweight) * Br000 +
364  rweight * Br010) +
365  zweight * ((1 - rweight) * Br100 +
366  rweight * Br110);
367 
368  // PHI Direction of B-field
369  BfieldCyl[2] = 0;
370 
371  if (Verbosity() > 2)
372  {
373  cout << "End GFCyl Call: <bz,br,bphi> : {"
374  << BfieldCyl[0] / gauss << "," << BfieldCyl[1] / gauss << "," << BfieldCyl[2] / gauss << "}"
375  << endl;
376  }
377 
378  return;
379 }
380 
381 // debug function to print key/value pairs in map
382 void PHField2D::print_map(map<trio, trio>::iterator &it) const
383 {
384  cout << " Key: <"
385  << it->first.get<0>() / cm << ","
386  << it->first.get<1>() / cm << ">"
387 
388  << " Value: <"
389  << it->second.get<0>() / magfield_unit << ","
390  << it->second.get<1>() / magfield_unit << ">\n";
391 }