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