| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bi
28
29 #include "HadrontherapyElectricTabulatedField3
30 #include "G4SystemOfUnits.hh"
31 #include "G4AutoLock.hh"
32
33 namespace{ G4Mutex MyHadrontherapyLockEField=
34
35 using namespace std;
36
37 HadrontherapyElectricTabulatedField3D::Hadront
38 :feXoffset(exOffset),feYoffset(eyOffset),feZ
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
49 G4AutoLock lock(&MyHadrontherapyLockEField);
50
51 ifstream file( filename ); // Open the file
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 ord
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
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 " << G4e
111
112 // G4cout << " Read values of field from fil
113 G4cout << " ---> assumed the order: x, y, z
114 << "\n ---> Min values x,y,z: "
115 << Eminx/cm << " " << Eminy/cm << " " << Em
116 << "\n ---> Max values x,y,z: "
117 << Emaxx/cm << " " << Emaxy/cm << " " << Em
118 << "\n ---> The field will be offset in x b
119 << "\n ---> The field will be offset
120 << "\n ---> The field will be offset
121
122 // Should really check that the limits are n
123 if (Emaxx < Eminx) {swap(Emaxx,Eminx); einve
124 if (Emaxy < Eminy) {swap(Emaxy,Eminy); einve
125 if (Emaxz < Eminz) {swap(Emaxz,Eminz); einve
126 G4cout << "\nAfter reordering if neccesary"
127 << "\n ---> Min values x,y,z: "
128 << Eminx/cm << " " << Eminy/cm << " " << Em
129 << " \n ---> Max values x,y,z: "
130 << Emaxx/cm << " " << Emaxy/cm << " " << Em
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
137 << "\n-------------------------------------
138 }
139
140 void HadrontherapyElectricTabulatedField3D::Ge
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,
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 - Exfractio
154 if (einvertY) { Eyfraction = 1 - Eyfractio
155 if (einvertZ) { Ezfraction = 1 - Ezfractio
156
157 // Need addresses of these to pass to modf
158 // modf uses its second argument as an OUT
159 G4double exdindex, eydindex, ezdindex;
160
161 // Position of the point within the cuboid
162 // nearest surrounding tabulated points
163 G4double exlocal = ( std::modf(Exfraction*
164 G4double eylocal = ( std::modf(Eyfraction*
165 G4double ezlocal = ( std::modf(Ezfraction*
166
167 // The indices of the nearest tabulated po
168 // are all less than those of the given po
169 G4int exindex = static_cast<G4int>(std::fl
170 G4int eyindex = static_cast<G4int>(std::fl
171 G4int ezindex = static_cast<G4int>(std::fl
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 << "
190 G4cout << "Index x,y,z: " << exindex << "
191 G4double valx0z0, mulx0z0, valx1z0, mulx1z
192 G4double valx0z1, mulx0z1, valx1z1, mulx1z
193 valx0z0= table[exindex ][0][ezindex]; mu
194 valx1z0= table[exindex+1][0][ezindex]; mu
195 valx0z1= table[exindex ][0][ezindex+1]; m
196 valx1z1= table[exindex+1][0][ezindex+1]; m
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 ][ezinde
207 xEField[exindex ][eyindex ][ezinde
208 xEField[exindex ][eyindex+1][ezinde
209 xEField[exindex ][eyindex+1][ezinde
210 xEField[exindex+1][eyindex ][ezinde
211 xEField[exindex+1][eyindex ][ezinde
212 xEField[exindex+1][eyindex+1][ezinde
213 xEField[exindex+1][eyindex+1][ezinde
214 Efield[4] =
215 yEField[exindex ][eyindex ][ezinde
216 yEField[exindex ][eyindex ][ezinde
217 yEField[exindex ][eyindex+1][ezinde
218 yEField[exindex ][eyindex+1][ezinde
219 yEField[exindex+1][eyindex ][ezinde
220 yEField[exindex+1][eyindex ][ezinde
221 yEField[exindex+1][eyindex+1][ezinde
222 yEField[exindex+1][eyindex+1][ezinde
223 Efield[5] =
224 zEField[exindex ][eyindex ][ezinde
225 zEField[exindex ][eyindex ][ezinde
226 zEField[exindex ][eyindex+1][ezinde
227 zEField[exindex ][eyindex+1][ezinde
228 zEField[exindex+1][eyindex ][ezinde
229 zEField[exindex+1][eyindex ][ezinde
230 zEField[exindex+1][eyindex+1][ezinde
231 zEField[exindex+1][eyindex+1][ezinde
232 }
233 //G4cout << "Getting electric field " << Efiel
234 //G4cout << "For coordinates: " << Epoint[0] <
235
236 /*std::ofstream WriteDataIn("ElectricFieldFC.o
237 WriteDataIn << Epoint[0]
238 << Epoint[1] << '\t' << "
239 << Epoint[2] << '\t' << "
240 << Efield[3]/(volt/m) << '\
241 << Efield[4]/(volt/m) << '\t
242 << Efield[5]/(volt/m) << '\t
243 << std::endl; */
244 }
245