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