ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HadrontherapyElectricTabulatedField3D.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HadrontherapyElectricTabulatedField3D.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 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
28 
29 #include "HadrontherapyElectricTabulatedField3D.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4AutoLock.hh"
32 
33 namespace{ G4Mutex MyHadrontherapyLockEField=G4MUTEX_INITIALIZER; }
34 
35 using namespace std;
36 
37 HadrontherapyElectricTabulatedField3D::HadrontherapyElectricTabulatedField3D( const char* filename, G4double exOffset, G4double eyOffset, G4double ezOffset)
38  :feXoffset(exOffset),feYoffset(eyOffset),feZoffset(ezOffset),einvertX(false),einvertY(false),einvertZ(false)
39 {
40  //The format file is: X Y Z Ex Ey Ez
41 
42  G4double ElenUnit= cm;
43  G4double EfieldUnit= volt/m;
44  G4cout << "\n-----------------------------------------------------------"
45  << "\n Electric field"
46  << "\n-----------------------------------------------------------";
47 
48  G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
49  G4AutoLock lock(&MyHadrontherapyLockEField);
50 
51  ifstream file( filename ); // Open the file for reading.
52 
53  // Ignore first blank line
54  char ebuffer[256];
55  file.getline(ebuffer,256);
56 
57  // Read table dimensions
58  file >> Enx >> Eny >> Enz; // Note dodgy order
59 
60  G4cout << " [ Number of values x,y,z: "
61  << Enx << " " << Eny << " " << Enz << " ] "
62  << G4endl;
63 
64  // Set up storage space for table
65  xEField.resize( Enx );
66  yEField.resize( Enx );
67  zEField.resize( Enx );
68  G4int ix, iy, iz;
69  for (ix=0; ix<Enx; ix++) {
70  xEField[ix].resize(Eny);
71  yEField[ix].resize(Eny);
72  zEField[ix].resize(Eny);
73  for (iy=0; iy<Eny; iy++) {
74  xEField[ix][iy].resize(Enz);
75  yEField[ix][iy].resize(Enz);
76  zEField[ix][iy].resize(Enz);
77  }
78  }
79 
80  // Read in the data
81  G4double Exval=0.;
82  G4double Eyval=0.;
83  G4double Ezval=0.;
84  G4double Ex=0.;
85  G4double Ey=0.;
86  G4double Ez=0.;
87  for (iz=0; iz<Enz; iz++) {
88  for (iy=0; iy<Eny; iy++) {
89  for (ix=0; ix<Enx; ix++) {
90  file >> Exval >> Eyval >> Ezval >> Ex >> Ey >> Ez;
91 
92  if ( ix==0 && iy==0 && iz==0 ) {
93  Eminx = Exval * ElenUnit;
94  Eminy = Eyval * ElenUnit;
95  Eminz = Ezval * ElenUnit;
96  }
97  xEField[ix][iy][iz] = Ex * EfieldUnit;
98  yEField[ix][iy][iz] = Ey * EfieldUnit;
99  zEField[ix][iy][iz] = Ez * EfieldUnit;
100  }
101  }
102  }
103  file.close();
104  lock.unlock();
105 
106  Emaxx = Exval * ElenUnit;
107  Emaxy = Eyval * ElenUnit;
108  Emaxz = Ezval * ElenUnit;
109 
110  G4cout << "\n ---> ... done reading " << G4endl;
111 
112  // G4cout << " Read values of field from file " << filename << G4endl;
113  G4cout << " ---> assumed the order: x, y, z, Ex, Ey, Ez "
114  << "\n ---> Min values x,y,z: "
115  << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
116  << "\n ---> Max values x,y,z: "
117  << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm "
118  << "\n ---> The field will be offset in x by " << exOffset/cm << " cm "
119  << "\n ---> The field will be offset in y by " << eyOffset/cm << " cm "
120  << "\n ---> The field will be offset in z by " << ezOffset/cm << " cm " << G4endl;
121 
122  // Should really check that the limits are not the wrong way around.
123  if (Emaxx < Eminx) {swap(Emaxx,Eminx); einvertX = true;}
124  if (Emaxy < Eminy) {swap(Emaxy,Eminy); einvertY = true;}
125  if (Emaxz < Eminz) {swap(Emaxz,Eminz); einvertZ = true;}
126  G4cout << "\nAfter reordering if neccesary"
127  << "\n ---> Min values x,y,z: "
128  << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
129  << " \n ---> Max values x,y,z: "
130  << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm ";
131 
132  dx1 = Emaxx - Eminx;
133  dy1 = Emaxy - Eminy;
134  dz1 = Emaxz - Eminz;
135  G4cout << "\n ---> Dif values x,y,z (range): "
136  << dx1/cm << " " << dy1/cm << " " << dz1/cm << " cm "
137  << "\n-----------------------------------------------------------" << G4endl;
138 }
139 
140 void HadrontherapyElectricTabulatedField3D::GetFieldValue(const G4double Epoint[4],
141  G4double *Efield ) const
142 {
143  G4double x1 = Epoint[0] + feXoffset;
144  G4double y1 = Epoint[1] + feYoffset;
145  G4double z1 = Epoint[2] + feZoffset;
146 
147  // Position of given point within region, normalized to the range
148  // [0,1]
149  G4double Exfraction = (x1 - Eminx) / dx1;
150  G4double Eyfraction = (y1 - Eminy) / dy1;
151  G4double Ezfraction = (z1 - Eminz) / dz1;
152 
153  if (einvertX) { Exfraction = 1 - Exfraction;}
154  if (einvertY) { Eyfraction = 1 - Eyfraction;}
155  if (einvertZ) { Ezfraction = 1 - Ezfraction;}
156 
157  // Need addresses of these to pass to modf below.
158  // modf uses its second argument as an OUTPUT argument.
159  G4double exdindex, eydindex, ezdindex;
160 
161  // Position of the point within the cuboid defined by the
162  // nearest surrounding tabulated points
163  G4double exlocal = ( std::modf(Exfraction*(Enx-1), &exdindex));
164  G4double eylocal = ( std::modf(Eyfraction*(Eny-1), &eydindex));
165  G4double ezlocal = ( std::modf(Ezfraction*(Enz-1), &ezdindex));
166 
167  // The indices of the nearest tabulated point whose coordinates
168  // are all less than those of the given point
169  G4int exindex = static_cast<G4int>(std::floor(exdindex));
170  G4int eyindex = static_cast<G4int>(std::floor(eydindex));
171  G4int ezindex = static_cast<G4int>(std::floor(ezdindex));
172 
173  if ((exindex < 0) || (exindex >= Enx - 1) ||
174  (eyindex < 0) || (eyindex >= Eny - 1) ||
175  (ezindex < 0) || (ezindex >= Enz - 1))
176  {
177  Efield[0] = 0.0;
178  Efield[1] = 0.0;
179  Efield[2] = 0.0;
180  Efield[3] = 0.0;
181  Efield[4] = 0.0;
182  Efield[5] = 0.0;
183  }
184  else
185  {
186 
187 /*
188 #ifdef DEBUG_G4intERPOLATING_FIELD
189  G4cout << "Local x,y,z: " << exlocal << " " << eylocal << " " << ezlocal << G4endl;
190  G4cout << "Index x,y,z: " << exindex << " " << eyindex << " " << ezindex << G4endl;
191  G4double valx0z0, mulx0z0, valx1z0, mulx1z0;
192  G4double valx0z1, mulx0z1, valx1z1, mulx1z1;
193  valx0z0= table[exindex ][0][ezindex]; mulx0z0= (1-exlocal) * (1-ezlocal);
194  valx1z0= table[exindex+1][0][ezindex]; mulx1z0= exlocal * (1-ezlocal);
195  valx0z1= table[exindex ][0][ezindex+1]; mulx0z1= (1-exlocal) * ezlocal;
196  valx1z1= table[exindex+1][0][ezindex+1]; mulx1z1= exlocal * ezlocal;
197 #endif
198 */
199  // Full 3-dimensional version
200 
201  Efield[0] = 0.0;
202  Efield[1] = 0.0;
203  Efield[2] = 0.0;
204 
205  Efield[3] =
206  xEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
207  xEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
208  xEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
209  xEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
210  xEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
211  xEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
212  xEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
213  xEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
214  Efield[4] =
215  yEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
216  yEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
217  yEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
218  yEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
219  yEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
220  yEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
221  yEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
222  yEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
223  Efield[5] =
224  zEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
225  zEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
226  zEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
227  zEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
228  zEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
229  zEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
230  zEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
231  zEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
232  }
233 //G4cout << "Getting electric field " << Efield[3]/(volt/m) << " " << Efield[4]/(volt/m) << " " << Efield[5]/(volt/m) << G4endl;
234 //G4cout << "For coordinates: " << Epoint[0] << " " << Epoint[1] << " " << Epoint[2] << G4endl;
235 
236 /*std::ofstream WriteDataIn("ElectricFieldFC.out", std::ios::app);
237  WriteDataIn << Epoint[0] << '\t' << " "
238  << Epoint[1] << '\t' << " "
239  << Epoint[2] << '\t' << " "
240  << Efield[3]/(volt/m) << '\t' << " "
241  << Efield[4]/(volt/m) << '\t' << " "
242  << Efield[5]/(volt/m) << '\t' << " "
243  << std::endl; */
244 }