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