ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PurgMagTabulatedField3D.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PurgMagTabulatedField3D.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Code developed by:
27 // S.Larsson and J. Generowicz.
28 //
29 // *************************************
30 // * *
31 // * PurgMagTabulatedField3D.cc *
32 // * *
33 // *************************************
34 //
35 //
36 
37 #include "PurgMagTabulatedField3D.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4AutoLock.hh"
40 
41 namespace{
42  G4Mutex myPurgMagTabulatedField3DLock = G4MUTEX_INITIALIZER;
43 }
44 
45 using namespace std;
46 
47 PurgMagTabulatedField3D::PurgMagTabulatedField3D(const char* filename,
48  double zOffset )
49  :fZoffset(zOffset),invertX(false),invertY(false),invertZ(false)
50 {
51 
52  double lenUnit= meter;
53  double fieldUnit= tesla;
54  G4cout << "\n-----------------------------------------------------------"
55  << "\n Magnetic field"
56  << "\n-----------------------------------------------------------";
57 
58  G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
59 
60  //
61  //This is a thread-local class and we have to avoid that all workers open the
62  //file at the same time
63  G4AutoLock lock(&myPurgMagTabulatedField3DLock);
64 
65  ifstream file( filename ); // Open the file for reading.
66 
67  if (!file.is_open())
68  {
69  G4ExceptionDescription ed;
70  ed << "Could not open input file " << filename << std::endl;
71  G4Exception("PurgMagTabulatedField3D::PurgMagTabulatedField3D",
72  "pugmag001",FatalException,ed);
73  }
74 
75  // Ignore first blank line
76  char buffer[256];
77  file.getline(buffer,256);
78 
79  // Read table dimensions
80  file >> nx >> ny >> nz; // Note dodgy order
81 
82  G4cout << " [ Number of values x,y,z: "
83  << nx << " " << ny << " " << nz << " ] "
84  << G4endl;
85 
86  // Set up storage space for table
87  xField.resize( nx );
88  yField.resize( nx );
89  zField.resize( nx );
90  int ix, iy, iz;
91  for (ix=0; ix<nx; ix++) {
92  xField[ix].resize(ny);
93  yField[ix].resize(ny);
94  zField[ix].resize(ny);
95  for (iy=0; iy<ny; iy++) {
96  xField[ix][iy].resize(nz);
97  yField[ix][iy].resize(nz);
98  zField[ix][iy].resize(nz);
99  }
100  }
101 
102  // Ignore other header information
103  // The first line whose second character is '0' is considered to
104  // be the last line of the header.
105  do {
106  file.getline(buffer,256);
107  } while ( buffer[1]!='0');
108 
109  // Read in the data
110  double xval,yval,zval,bx,by,bz;
111  double permeability; // Not used in this example.
112  for (ix=0; ix<nx; ix++) {
113  for (iy=0; iy<ny; iy++) {
114  for (iz=0; iz<nz; iz++) {
115  file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
116  if ( ix==0 && iy==0 && iz==0 ) {
117  minx = xval * lenUnit;
118  miny = yval * lenUnit;
119  minz = zval * lenUnit;
120  }
121  xField[ix][iy][iz] = bx * fieldUnit;
122  yField[ix][iy][iz] = by * fieldUnit;
123  zField[ix][iy][iz] = bz * fieldUnit;
124  }
125  }
126  }
127  file.close();
128 
129  lock.unlock();
130 
131  maxx = xval * lenUnit;
132  maxy = yval * lenUnit;
133  maxz = zval * lenUnit;
134 
135  G4cout << "\n ---> ... done reading " << G4endl;
136 
137  // G4cout << " Read values of field from file " << filename << G4endl;
138  G4cout << " ---> assumed the order: x, y, z, Bx, By, Bz "
139  << "\n ---> Min values x,y,z: "
140  << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
141  << "\n ---> Max values x,y,z: "
142  << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm "
143  << "\n ---> The field will be offset by " << zOffset/cm << " cm " << G4endl;
144 
145  // Should really check that the limits are not the wrong way around.
146  if (maxx < minx) {swap(maxx,minx); invertX = true;}
147  if (maxy < miny) {swap(maxy,miny); invertY = true;}
148  if (maxz < minz) {swap(maxz,minz); invertZ = true;}
149  G4cout << "\nAfter reordering if neccesary"
150  << "\n ---> Min values x,y,z: "
151  << minx/cm << " " << miny/cm << " " << minz/cm << " cm "
152  << " \n ---> Max values x,y,z: "
153  << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm ";
154 
155  dx = maxx - minx;
156  dy = maxy - miny;
157  dz = maxz - minz;
158  G4cout << "\n ---> Dif values x,y,z (range): "
159  << dx/cm << " " << dy/cm << " " << dz/cm << " cm in z "
160  << "\n-----------------------------------------------------------" << G4endl;
161 }
162 
163 void PurgMagTabulatedField3D::GetFieldValue(const double point[4],
164  double *Bfield ) const
165 {
166 
167  double x = point[0];
168  double y = point[1];
169  double z = point[2] + fZoffset;
170 
171  // Check that the point is within the defined region
172  if ( x>=minx && x<=maxx &&
173  y>=miny && y<=maxy &&
174  z>=minz && z<=maxz ) {
175 
176  // Position of given point within region, normalized to the range
177  // [0,1]
178  double xfraction = (x - minx) / dx;
179  double yfraction = (y - miny) / dy;
180  double zfraction = (z - minz) / dz;
181 
182  if (invertX) { xfraction = 1 - xfraction;}
183  if (invertY) { yfraction = 1 - yfraction;}
184  if (invertZ) { zfraction = 1 - zfraction;}
185 
186  // Need addresses of these to pass to modf below.
187  // modf uses its second argument as an OUTPUT argument.
188  double xdindex, ydindex, zdindex;
189 
190  // Position of the point within the cuboid defined by the
191  // nearest surrounding tabulated points
192  double xlocal = ( std::modf(xfraction*(nx-1), &xdindex));
193  double ylocal = ( std::modf(yfraction*(ny-1), &ydindex));
194  double zlocal = ( std::modf(zfraction*(nz-1), &zdindex));
195 
196  // The indices of the nearest tabulated point whose coordinates
197  // are all less than those of the given point
198  int xindex = static_cast<int>(xdindex);
199  int yindex = static_cast<int>(ydindex);
200  int zindex = static_cast<int>(zdindex);
201 
202 
203 #ifdef DEBUG_INTERPOLATING_FIELD
204  G4cout << "Local x,y,z: " << xlocal << " " << ylocal << " " << zlocal << G4endl;
205  G4cout << "Index x,y,z: " << xindex << " " << yindex << " " << zindex << G4endl;
206  double valx0z0, mulx0z0, valx1z0, mulx1z0;
207  double valx0z1, mulx0z1, valx1z1, mulx1z1;
208  valx0z0= table[xindex ][0][zindex]; mulx0z0= (1-xlocal) * (1-zlocal);
209  valx1z0= table[xindex+1][0][zindex]; mulx1z0= xlocal * (1-zlocal);
210  valx0z1= table[xindex ][0][zindex+1]; mulx0z1= (1-xlocal) * zlocal;
211  valx1z1= table[xindex+1][0][zindex+1]; mulx1z1= xlocal * zlocal;
212 #endif
213 
214  // Full 3-dimensional version
215  Bfield[0] =
216  xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
217  xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
218  xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
219  xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
220  xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
221  xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
222  xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
223  xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
224  Bfield[1] =
225  yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
226  yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
227  yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
228  yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
229  yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
230  yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
231  yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
232  yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
233  Bfield[2] =
234  zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
235  zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
236  zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
237  zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
238  zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
239  zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
240  zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
241  zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
242 
243  } else {
244  Bfield[0] = 0.0;
245  Bfield[1] = 0.0;
246  Bfield[2] = 0.0;
247  }
248 }
249