| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Hadrontherapy advanced example for Geant4 26 // Hadrontherapy advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bi 27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
28 28
29 #include "HadrontherapyElectricTabulatedField3 29 #include "HadrontherapyElectricTabulatedField3D.hh"
30 #include "G4SystemOfUnits.hh" 30 #include "G4SystemOfUnits.hh"
31 #include "G4AutoLock.hh" 31 #include "G4AutoLock.hh"
32 32
33 namespace{ G4Mutex MyHadrontherapyLockEField= 33 namespace{ G4Mutex MyHadrontherapyLockEField=G4MUTEX_INITIALIZER; }
34 34
35 using namespace std; << 35 HadrontherapyElectricTabulatedField3D::HadrontherapyElectricTabulatedField3D( const char* filename, G4double exOffset, G4double eyOffset, G4double ezOffset)
36 <<
37 HadrontherapyElectricTabulatedField3D::Hadront <<
38 :feXoffset(exOffset),feYoffset(eyOffset),feZ 36 :feXoffset(exOffset),feYoffset(eyOffset),feZoffset(ezOffset),einvertX(false),einvertY(false),einvertZ(false)
39 { << 37 {
40 //The format file is: X Y Z Ex Ey Ez 38 //The format file is: X Y Z Ex Ey Ez
41 39
42 G4double ElenUnit= cm; 40 G4double ElenUnit= cm;
43 G4double EfieldUnit= volt/m; << 41 G4double EfieldUnit= volt/m;
44 G4cout << "\n------------------------------- 42 G4cout << "\n-----------------------------------------------------------"
45 << "\n Electric field" << 43 << "\n Electric field"
46 << "\n------------------------------- << 44 << "\n-----------------------------------------------------------";
47 << 45
48 G4cout << "\n ---> " "Reading the field grid << 46 G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << endl;
49 G4AutoLock lock(&MyHadrontherapyLockEField); 47 G4AutoLock lock(&MyHadrontherapyLockEField);
50 48
51 ifstream file( filename ); // Open the file 49 ifstream file( filename ); // Open the file for reading.
52 << 50
53 // Ignore first blank line 51 // Ignore first blank line
54 char ebuffer[256]; 52 char ebuffer[256];
55 file.getline(ebuffer,256); 53 file.getline(ebuffer,256);
56 << 54
57 // Read table dimensions << 55 // Read table dimensions
58 file >> Enx >> Eny >> Enz; // Note dodgy ord 56 file >> Enx >> Eny >> Enz; // Note dodgy order
59 57
60 G4cout << " [ Number of values x,y,z: " << 58 G4cout << " [ Number of values x,y,z: "
61 << Enx << " " << Eny << " " << Enz << " ] " 59 << Enx << " " << Eny << " " << Enz << " ] "
62 << G4endl; << 60 << endl;
63 61
64 // Set up storage space for table 62 // Set up storage space for table
65 xEField.resize( Enx ); 63 xEField.resize( Enx );
66 yEField.resize( Enx ); 64 yEField.resize( Enx );
67 zEField.resize( Enx ); 65 zEField.resize( Enx );
68 G4int ix, iy, iz; 66 G4int ix, iy, iz;
69 for (ix=0; ix<Enx; ix++) { 67 for (ix=0; ix<Enx; ix++) {
70 xEField[ix].resize(Eny); 68 xEField[ix].resize(Eny);
71 yEField[ix].resize(Eny); 69 yEField[ix].resize(Eny);
72 zEField[ix].resize(Eny); 70 zEField[ix].resize(Eny);
73 for (iy=0; iy<Eny; iy++) { 71 for (iy=0; iy<Eny; iy++) {
74 xEField[ix][iy].resize(Enz); 72 xEField[ix][iy].resize(Enz);
75 yEField[ix][iy].resize(Enz); 73 yEField[ix][iy].resize(Enz);
76 zEField[ix][iy].resize(Enz); 74 zEField[ix][iy].resize(Enz);
77 } 75 }
78 } 76 }
79 77
80 // Read in the data 78 // Read in the data
81 G4double Exval=0.; 79 G4double Exval=0.;
82 G4double Eyval=0.; 80 G4double Eyval=0.;
83 G4double Ezval=0.; 81 G4double Ezval=0.;
84 G4double Ex=0.; 82 G4double Ex=0.;
85 G4double Ey=0.; 83 G4double Ey=0.;
86 G4double Ez=0.; 84 G4double Ez=0.;
87 for (iz=0; iz<Enz; iz++) { 85 for (iz=0; iz<Enz; iz++) {
88 for (iy=0; iy<Eny; iy++) { 86 for (iy=0; iy<Eny; iy++) {
89 for (ix=0; ix<Enx; ix++) { 87 for (ix=0; ix<Enx; ix++) {
90 file >> Exval >> Eyval >> Ezval >> Ex 88 file >> Exval >> Eyval >> Ezval >> Ex >> Ey >> Ez;
91 << 89
92 if ( ix==0 && iy==0 && iz==0 ) { 90 if ( ix==0 && iy==0 && iz==0 ) {
93 Eminx = Exval * ElenUnit; 91 Eminx = Exval * ElenUnit;
94 Eminy = Eyval * ElenUnit; 92 Eminy = Eyval * ElenUnit;
95 Eminz = Ezval * ElenUnit; 93 Eminz = Ezval * ElenUnit;
96 } 94 }
97 xEField[ix][iy][iz] = Ex * EfieldUnit; 95 xEField[ix][iy][iz] = Ex * EfieldUnit;
98 yEField[ix][iy][iz] = Ey * EfieldUnit; 96 yEField[ix][iy][iz] = Ey * EfieldUnit;
99 zEField[ix][iy][iz] = Ez * EfieldUnit; 97 zEField[ix][iy][iz] = Ez * EfieldUnit;
100 } 98 }
101 } 99 }
102 } 100 }
103 file.close(); 101 file.close();
104 lock.unlock(); 102 lock.unlock();
105 103
106 Emaxx = Exval * ElenUnit; 104 Emaxx = Exval * ElenUnit;
107 Emaxy = Eyval * ElenUnit; 105 Emaxy = Eyval * ElenUnit;
108 Emaxz = Ezval * ElenUnit; 106 Emaxz = Ezval * ElenUnit;
109 107
110 G4cout << "\n ---> ... done reading " << G4e << 108 G4cout << "\n ---> ... done reading " << endl;
111 109
112 // G4cout << " Read values of field from fil << 110 // G4cout << " Read values of field from file " << filename << endl;
113 G4cout << " ---> assumed the order: x, y, z 111 G4cout << " ---> assumed the order: x, y, z, Ex, Ey, Ez "
114 << "\n ---> Min values x,y,z: " << 112 << "\n ---> Min values x,y,z: "
115 << Eminx/cm << " " << Eminy/cm << " " << Em 113 << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
116 << "\n ---> Max values x,y,z: " << 114 << "\n ---> Max values x,y,z: "
117 << Emaxx/cm << " " << Emaxy/cm << " " << Em 115 << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm "
118 << "\n ---> The field will be offset in x b << 116 << "\n ---> The field will be offset in x by " << exOffset/cm << " cm "
119 << "\n ---> The field will be offset << 117 << "\n ---> The field will be offset in y by " << eyOffset/cm << " cm "
120 << "\n ---> The field will be offset << 118 << "\n ---> The field will be offset in z by " << ezOffset/cm << " cm " << endl;
121 119
122 // Should really check that the limits are n 120 // Should really check that the limits are not the wrong way around.
123 if (Emaxx < Eminx) {swap(Emaxx,Eminx); einve << 121 if (Emaxx < Eminx) {swap(Emaxx,Eminx); einvertX = true;}
124 if (Emaxy < Eminy) {swap(Emaxy,Eminy); einve << 122 if (Emaxy < Eminy) {swap(Emaxy,Eminy); einvertY = true;}
125 if (Emaxz < Eminz) {swap(Emaxz,Eminz); einve << 123 if (Emaxz < Eminz) {swap(Emaxz,Eminz); einvertZ = true;}
126 G4cout << "\nAfter reordering if neccesary" << 124 G4cout << "\nAfter reordering if neccesary"
127 << "\n ---> Min values x,y,z: " << 125 << "\n ---> Min values x,y,z: "
128 << Eminx/cm << " " << Eminy/cm << " " << Em 126 << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
129 << " \n ---> Max values x,y,z: " << 127 << " \n ---> Max values x,y,z: "
130 << Emaxx/cm << " " << Emaxy/cm << " " << Em 128 << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm ";
131 129
132 dx1 = Emaxx - Eminx; 130 dx1 = Emaxx - Eminx;
133 dy1 = Emaxy - Eminy; 131 dy1 = Emaxy - Eminy;
134 dz1 = Emaxz - Eminz; 132 dz1 = Emaxz - Eminz;
135 G4cout << "\n ---> Dif values x,y,z (range): << 133 G4cout << "\n ---> Dif values x,y,z (range): "
136 << dx1/cm << " " << dy1/cm << " " << dz1/cm 134 << dx1/cm << " " << dy1/cm << " " << dz1/cm << " cm "
137 << "\n------------------------------------- << 135 << "\n-----------------------------------------------------------" << endl;
138 } 136 }
139 137
140 void HadrontherapyElectricTabulatedField3D::Ge 138 void HadrontherapyElectricTabulatedField3D::GetFieldValue(const G4double Epoint[4],
141 G4double *Efield ) const 139 G4double *Efield ) const
142 { 140 {
143 G4double x1 = Epoint[0] + feXoffset; << 141 G4double x1 = Epoint[0] + feXoffset;
144 G4double y1 = Epoint[1] + feYoffset; << 142 G4double y1 = Epoint[1] + feYoffset;
145 G4double z1 = Epoint[2] + feZoffset; << 143 G4double z1 = Epoint[2] + feZoffset;
146 << 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
147 // Position of given point within region, 150 // Position of given point within region, normalized to the range
148 // [0,1] 151 // [0,1]
149 G4double Exfraction = (x1 - Eminx) / dx1; 152 G4double Exfraction = (x1 - Eminx) / dx1;
150 G4double Eyfraction = (y1 - Eminy) / dy1; << 153 G4double Eyfraction = (y1 - Eminy) / dy1;
151 G4double Ezfraction = (z1 - Eminz) / dz1; 154 G4double Ezfraction = (z1 - Eminz) / dz1;
152 155
153 if (einvertX) { Exfraction = 1 - Exfractio 156 if (einvertX) { Exfraction = 1 - Exfraction;}
154 if (einvertY) { Eyfraction = 1 - Eyfractio 157 if (einvertY) { Eyfraction = 1 - Eyfraction;}
155 if (einvertZ) { Ezfraction = 1 - Ezfractio 158 if (einvertZ) { Ezfraction = 1 - Ezfraction;}
156 159
157 // Need addresses of these to pass to modf 160 // Need addresses of these to pass to modf below.
158 // modf uses its second argument as an OUT 161 // modf uses its second argument as an OUTPUT argument.
159 G4double exdindex, eydindex, ezdindex; 162 G4double exdindex, eydindex, ezdindex;
160 << 163
161 // Position of the point within the cuboid 164 // Position of the point within the cuboid defined by the
162 // nearest surrounding tabulated points 165 // nearest surrounding tabulated points
163 G4double exlocal = ( std::modf(Exfraction* 166 G4double exlocal = ( std::modf(Exfraction*(Enx-1), &exdindex));
164 G4double eylocal = ( std::modf(Eyfraction* 167 G4double eylocal = ( std::modf(Eyfraction*(Eny-1), &eydindex));
165 G4double ezlocal = ( std::modf(Ezfraction* 168 G4double ezlocal = ( std::modf(Ezfraction*(Enz-1), &ezdindex));
166 << 169
167 // The indices of the nearest tabulated po 170 // The indices of the nearest tabulated point whose coordinates
168 // are all less than those of the given po 171 // are all less than those of the given point
169 G4int exindex = static_cast<G4int>(std::fl << 172 G4int exindex = static_cast<G4int>(exdindex);
170 G4int eyindex = static_cast<G4int>(std::fl << 173 G4int eyindex = static_cast<G4int>(eydindex);
171 G4int ezindex = static_cast<G4int>(std::fl << 174 G4int ezindex = static_cast<G4int>(ezdindex);
172 << 175
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 /* 176 /*
188 #ifdef DEBUG_G4intERPOLATING_FIELD 177 #ifdef DEBUG_G4intERPOLATING_FIELD
189 G4cout << "Local x,y,z: " << exlocal << " << 178 G4cout << "Local x,y,z: " << exlocal << " " << eylocal << " " << ezlocal << endl;
190 G4cout << "Index x,y,z: " << exindex << " << 179 G4cout << "Index x,y,z: " << exindex << " " << eyindex << " " << ezindex << endl;
191 G4double valx0z0, mulx0z0, valx1z0, mulx1z 180 G4double valx0z0, mulx0z0, valx1z0, mulx1z0;
192 G4double valx0z1, mulx0z1, valx1z1, mulx1z 181 G4double valx0z1, mulx0z1, valx1z1, mulx1z1;
193 valx0z0= table[exindex ][0][ezindex]; mu 182 valx0z0= table[exindex ][0][ezindex]; mulx0z0= (1-exlocal) * (1-ezlocal);
194 valx1z0= table[exindex+1][0][ezindex]; mu 183 valx1z0= table[exindex+1][0][ezindex]; mulx1z0= exlocal * (1-ezlocal);
195 valx0z1= table[exindex ][0][ezindex+1]; m 184 valx0z1= table[exindex ][0][ezindex+1]; mulx0z1= (1-exlocal) * ezlocal;
196 valx1z1= table[exindex+1][0][ezindex+1]; m 185 valx1z1= table[exindex+1][0][ezindex+1]; mulx1z1= exlocal * ezlocal;
197 #endif 186 #endif
198 */ 187 */
199 // Full 3-dimensional version 188 // Full 3-dimensional version
200 189
201 Efield[0] = 0.0; << 190 Efield[0] = 0.0;
202 Efield[1] = 0.0; << 191 Efield[1] = 0.0;
203 Efield[2] = 0.0; << 192 Efield[2] = 0.0;
204 << 193
205 Efield[3] = << 194 Efield[3] =
206 xEField[exindex ][eyindex ][ezinde << 195 xEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
207 xEField[exindex ][eyindex ][ezinde << 196 xEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
208 xEField[exindex ][eyindex+1][ezinde << 197 xEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
209 xEField[exindex ][eyindex+1][ezinde << 198 xEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
210 xEField[exindex+1][eyindex ][ezinde << 199 xEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
211 xEField[exindex+1][eyindex ][ezinde << 200 xEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
212 xEField[exindex+1][eyindex+1][ezinde << 201 xEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
213 xEField[exindex+1][eyindex+1][ezinde << 202 xEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
214 Efield[4] = << 203 Efield[4] =
215 yEField[exindex ][eyindex ][ezinde << 204 yEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
216 yEField[exindex ][eyindex ][ezinde << 205 yEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
217 yEField[exindex ][eyindex+1][ezinde << 206 yEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
218 yEField[exindex ][eyindex+1][ezinde << 207 yEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
219 yEField[exindex+1][eyindex ][ezinde << 208 yEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
220 yEField[exindex+1][eyindex ][ezinde << 209 yEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
221 yEField[exindex+1][eyindex+1][ezinde << 210 yEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
222 yEField[exindex+1][eyindex+1][ezinde << 211 yEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ;
223 Efield[5] = << 212 Efield[5] =
224 zEField[exindex ][eyindex ][ezinde << 213 zEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
225 zEField[exindex ][eyindex ][ezinde << 214 zEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal +
226 zEField[exindex ][eyindex+1][ezinde << 215 zEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) +
227 zEField[exindex ][eyindex+1][ezinde << 216 zEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal +
228 zEField[exindex+1][eyindex ][ezinde << 217 zEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) +
229 zEField[exindex+1][eyindex ][ezinde << 218 zEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal +
230 zEField[exindex+1][eyindex+1][ezinde << 219 zEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) +
231 zEField[exindex+1][eyindex+1][ezinde << 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;
232 } 229 }
233 //G4cout << "Getting electric field " << Efiel << 230 //G4cout << "Getting electric field " << Efield[3]/(volt/m) << " " << Efield[4]/(volt/m) << " " << Efield[5]/(volt/m) << endl;
234 //G4cout << "For coordinates: " << Epoint[0] < 231 //G4cout << "For coordinates: " << Epoint[0] << " " << Epoint[1] << " " << Epoint[2] << G4endl;
235 232
236 /*std::ofstream WriteDataIn("ElectricFieldFC.o 233 /*std::ofstream WriteDataIn("ElectricFieldFC.out", std::ios::app);
237 WriteDataIn << Epoint[0] 234 WriteDataIn << Epoint[0] << '\t' << " "
238 << Epoint[1] << '\t' << " 235 << Epoint[1] << '\t' << " "
239 << Epoint[2] << '\t' << " 236 << Epoint[2] << '\t' << " "
240 << Efield[3]/(volt/m) << '\ 237 << Efield[3]/(volt/m) << '\t' << " "
241 << Efield[4]/(volt/m) << '\t 238 << Efield[4]/(volt/m) << '\t' << " "
242 << Efield[5]/(volt/m) << '\t 239 << Efield[5]/(volt/m) << '\t' << " "
243 << std::endl; */ << 240 << G4endl; */
244 } 241 }
>> 242
245 243