ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHFieldCleo.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHFieldCleo.cc
1 #include "PHFieldCleo.h"
2 
3 #include <Geant4/G4SystemOfUnits.hh>
4 
5 #include <TSystem.h>
6 
7 #include <cmath> // for modf, fabs, isfinite, NAN
8 #include <cstdlib> // for exit
9 #include <fstream>
10 #include <iostream>
11 #include <memory> // for allocator_traits<>::value_type
12 
13 using namespace std;
14 
15 PHFieldCleo::PHFieldCleo(const string &filename, const int verb, const float magfield_rescale)
16  : PHField(verb)
17  , m_MagFieldScale(magfield_rescale)
18 {
19  ifstream file(filename);
20  if (!file.is_open())
21  {
22  cout << "Could not open " << filename << " exiting now" << endl;
23  gSystem->Exit(1);
24  exit(1);
25  }
26  char buffer[256];
27  // Read table dimensions
28  int tmp;
29  file >> nz >> ny >> nx >> tmp; // first z, then y and x (y,x are 50, it does not matter)
30  xField.resize(nx);
31  yField.resize(nx);
32  zField.resize(nx);
33  int ix, iy, iz;
34  for (ix = 0; ix < nx; ix++)
35  {
36  xField[ix].resize(ny);
37  yField[ix].resize(ny);
38  zField[ix].resize(ny);
39  for (iy = 0; iy < ny; iy++)
40  {
41  xField[ix][iy].resize(nz);
42  yField[ix][iy].resize(nz);
43  zField[ix][iy].resize(nz);
44  }
45  }
46 
47  // Ignore other header information
48  // The first line whose second character is '0' is considered to
49  // be the last line of the header.
50  do
51  {
52  file.getline(buffer, 256);
53  } while (buffer[1] != '0');
54 
55  // Read in the data
56  double xval, yval, zval, bx, by, bz;
57  double save_x = NAN;
58  double save_y = NAN;
59  for (ix = 0; ix < nx; ix++)
60  {
61  for (iy = 0; iy < ny; iy++)
62  {
63  for (iz = 0; iz < nz; iz++)
64  {
65  file >> xval >> yval >> zval >> bx >> by >> bz;
66  if (!isfinite(save_x))
67  {
68  save_x = xval;
69  }
70  else
71  {
72  if (save_x != xval)
73  {
74  cout << "mismatch expected x coordinate: " << save_x
75  << " read " << xval << endl;
76  gSystem->Exit(1);
77  }
78  }
79  if (!isfinite(save_y))
80  {
81  save_y = yval;
82  }
83  else
84  {
85  if (save_y != yval)
86  {
87  cout << "mismatch expected y coordinate: " << save_y
88  << " read " << yval << endl;
89  gSystem->Exit(1);
90  }
91  }
92  if (ix == 0 && iy == 0 && iz == 0)
93  {
94  minx = xval * cm;
95  miny = yval * cm;
96  minz = zval * cm;
97  }
98  xField[ix][iy][iz] = bx * gauss * m_MagFieldScale;
99  yField[ix][iy][iz] = by * gauss * m_MagFieldScale;
100  zField[ix][iy][iz] = bz * gauss * m_MagFieldScale;
101  }
102  save_y = NAN;
103  }
104  save_x = NAN;
105  }
106  file.close();
107  maxx = xval * cm;
108  maxy = yval * cm;
109  maxz = zval * cm;
110  dx = maxx - minx;
111  dy = maxy - miny;
112  dz = maxz - minz;
113 }
114 
115 void PHFieldCleo::GetFieldValue(const double point[4], double *Bfield) const
116 {
117  double x = fabs(point[0]);
118  double y = fabs(point[1]);
119  double z = point[2];
120 
121  // Check that the point is within the defined region
122  if (x >= minx && x <= maxx && y >= miny && y <= maxy && z >= minz && z <= maxz)
123  {
124  // Position of given point within region, normalized to the range
125  // [0,1]
126  double xfraction = (x - minx) / dx;
127  double yfraction = (y - miny) / dy;
128  double zfraction = (z - minz) / dz;
129 
130  // Need addresses of these to pass to modf below.
131  // modf uses its second argument as an OUTPUT argument.
132  double xdindex, ydindex, zdindex;
133 
134  // Position of the point within the cuboid defined by the
135  // nearest surrounding tabulated points
136  double xlocal = (std::modf(xfraction * (nx - 1), &xdindex));
137  double ylocal = (std::modf(yfraction * (ny - 1), &ydindex));
138  double zlocal = (std::modf(zfraction * (nz - 1), &zdindex));
139 
140  // The indices of the nearest tabulated point whose coordinates
141  // are all less than those of the given point
142  int xindex = static_cast<int>(xdindex);
143  int yindex = static_cast<int>(ydindex);
144  int zindex = static_cast<int>(zdindex);
145 
146  // Full 3-dimensional version
147  Bfield[0] = xField[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) * (1 - zlocal) + xField[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) * zlocal + xField[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal * (1 - zlocal) + xField[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal * zlocal + xField[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) * (1 - zlocal) + xField[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) * zlocal + xField[xindex + 1][yindex + 1][zindex] * xlocal * ylocal * (1 - zlocal) + xField[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
148  Bfield[1] = yField[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) * (1 - zlocal) + yField[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) * zlocal + yField[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal * (1 - zlocal) + yField[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal * zlocal + yField[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) * (1 - zlocal) + yField[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) * zlocal + yField[xindex + 1][yindex + 1][zindex] * xlocal * ylocal * (1 - zlocal) + yField[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
149  Bfield[2] = zField[xindex][yindex][zindex] * (1 - xlocal) * (1 - ylocal) * (1 - zlocal) + zField[xindex][yindex][zindex + 1] * (1 - xlocal) * (1 - ylocal) * zlocal + zField[xindex][yindex + 1][zindex] * (1 - xlocal) * ylocal * (1 - zlocal) + zField[xindex][yindex + 1][zindex + 1] * (1 - xlocal) * ylocal * zlocal + zField[xindex + 1][yindex][zindex] * xlocal * (1 - ylocal) * (1 - zlocal) + zField[xindex + 1][yindex][zindex + 1] * xlocal * (1 - ylocal) * zlocal + zField[xindex + 1][yindex + 1][zindex] * xlocal * ylocal * (1 - zlocal) + zField[xindex + 1][yindex + 1][zindex + 1] * xlocal * ylocal * zlocal;
150  }
151  else
152  {
153  Bfield[0] = 0.0;
154  Bfield[1] = 0.0;
155  Bfield[2] = 0.0;
156  }
157 }