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; << 36 << 37 HadrontherapyElectricTabulatedField3D::Hadront 35 HadrontherapyElectricTabulatedField3D::HadrontherapyElectricTabulatedField3D( const char* filename, G4double exOffset, G4double eyOffset, G4double ezOffset) 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 147 // Position of given point within region, 145 // Position of given point within region, normalized to the range 148 // [0,1] 146 // [0,1] 149 G4double Exfraction = (x1 - Eminx) / dx1; 147 G4double Exfraction = (x1 - Eminx) / dx1; 150 G4double Eyfraction = (y1 - Eminy) / dy1; 148 G4double Eyfraction = (y1 - Eminy) / dy1; 151 G4double Ezfraction = (z1 - Eminz) / dz1; 149 G4double Ezfraction = (z1 - Eminz) / dz1; 152 150 153 if (einvertX) { Exfraction = 1 - Exfractio 151 if (einvertX) { Exfraction = 1 - Exfraction;} 154 if (einvertY) { Eyfraction = 1 - Eyfractio 152 if (einvertY) { Eyfraction = 1 - Eyfraction;} 155 if (einvertZ) { Ezfraction = 1 - Ezfractio 153 if (einvertZ) { Ezfraction = 1 - Ezfraction;} 156 154 157 // Need addresses of these to pass to modf 155 // Need addresses of these to pass to modf below. 158 // modf uses its second argument as an OUT 156 // modf uses its second argument as an OUTPUT argument. 159 G4double exdindex, eydindex, ezdindex; 157 G4double exdindex, eydindex, ezdindex; 160 158 161 // Position of the point within the cuboid 159 // Position of the point within the cuboid defined by the 162 // nearest surrounding tabulated points 160 // nearest surrounding tabulated points 163 G4double exlocal = ( std::modf(Exfraction* 161 G4double exlocal = ( std::modf(Exfraction*(Enx-1), &exdindex)); 164 G4double eylocal = ( std::modf(Eyfraction* 162 G4double eylocal = ( std::modf(Eyfraction*(Eny-1), &eydindex)); 165 G4double ezlocal = ( std::modf(Ezfraction* 163 G4double ezlocal = ( std::modf(Ezfraction*(Enz-1), &ezdindex)); 166 164 167 // The indices of the nearest tabulated po 165 // The indices of the nearest tabulated point whose coordinates 168 // are all less than those of the given po 166 // are all less than those of the given point 169 G4int exindex = static_cast<G4int>(std::fl 167 G4int exindex = static_cast<G4int>(std::floor(exdindex)); 170 G4int eyindex = static_cast<G4int>(std::fl 168 G4int eyindex = static_cast<G4int>(std::floor(eydindex)); 171 G4int ezindex = static_cast<G4int>(std::fl 169 G4int ezindex = static_cast<G4int>(std::floor(ezdindex)); 172 170 173 if ((exindex < 0) || (exindex >= Enx - 1) 171 if ((exindex < 0) || (exindex >= Enx - 1) || 174 (eyindex < 0) || (eyindex >= Eny - 1) 172 (eyindex < 0) || (eyindex >= Eny - 1) || 175 (ezindex < 0) || (ezindex >= Enz - 1)) 173 (ezindex < 0) || (ezindex >= Enz - 1)) 176 { 174 { 177 Efield[0] = 0.0; 175 Efield[0] = 0.0; 178 Efield[1] = 0.0; 176 Efield[1] = 0.0; 179 Efield[2] = 0.0; 177 Efield[2] = 0.0; 180 Efield[3] = 0.0; 178 Efield[3] = 0.0; 181 Efield[4] = 0.0; 179 Efield[4] = 0.0; 182 Efield[5] = 0.0; 180 Efield[5] = 0.0; 183 } 181 } 184 else 182 else 185 { 183 { 186 184 187 /* 185 /* 188 #ifdef DEBUG_G4intERPOLATING_FIELD 186 #ifdef DEBUG_G4intERPOLATING_FIELD 189 G4cout << "Local x,y,z: " << exlocal << " << 187 G4cout << "Local x,y,z: " << exlocal << " " << eylocal << " " << ezlocal << endl; 190 G4cout << "Index x,y,z: " << exindex << " << 188 G4cout << "Index x,y,z: " << exindex << " " << eyindex << " " << ezindex << endl; 191 G4double valx0z0, mulx0z0, valx1z0, mulx1z 189 G4double valx0z0, mulx0z0, valx1z0, mulx1z0; 192 G4double valx0z1, mulx0z1, valx1z1, mulx1z 190 G4double valx0z1, mulx0z1, valx1z1, mulx1z1; 193 valx0z0= table[exindex ][0][ezindex]; mu 191 valx0z0= table[exindex ][0][ezindex]; mulx0z0= (1-exlocal) * (1-ezlocal); 194 valx1z0= table[exindex+1][0][ezindex]; mu 192 valx1z0= table[exindex+1][0][ezindex]; mulx1z0= exlocal * (1-ezlocal); 195 valx0z1= table[exindex ][0][ezindex+1]; m 193 valx0z1= table[exindex ][0][ezindex+1]; mulx0z1= (1-exlocal) * ezlocal; 196 valx1z1= table[exindex+1][0][ezindex+1]; m 194 valx1z1= table[exindex+1][0][ezindex+1]; mulx1z1= exlocal * ezlocal; 197 #endif 195 #endif 198 */ 196 */ 199 // Full 3-dimensional version 197 // Full 3-dimensional version 200 198 201 Efield[0] = 0.0; 199 Efield[0] = 0.0; 202 Efield[1] = 0.0; 200 Efield[1] = 0.0; 203 Efield[2] = 0.0; 201 Efield[2] = 0.0; 204 202 205 Efield[3] = 203 Efield[3] = 206 xEField[exindex ][eyindex ][ezinde 204 xEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) + 207 xEField[exindex ][eyindex ][ezinde 205 xEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal + 208 xEField[exindex ][eyindex+1][ezinde 206 xEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) + 209 xEField[exindex ][eyindex+1][ezinde 207 xEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal + 210 xEField[exindex+1][eyindex ][ezinde 208 xEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) + 211 xEField[exindex+1][eyindex ][ezinde 209 xEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal + 212 xEField[exindex+1][eyindex+1][ezinde 210 xEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) + 213 xEField[exindex+1][eyindex+1][ezinde 211 xEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ; 214 Efield[4] = 212 Efield[4] = 215 yEField[exindex ][eyindex ][ezinde 213 yEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) + 216 yEField[exindex ][eyindex ][ezinde 214 yEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal + 217 yEField[exindex ][eyindex+1][ezinde 215 yEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) + 218 yEField[exindex ][eyindex+1][ezinde 216 yEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal + 219 yEField[exindex+1][eyindex ][ezinde 217 yEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) + 220 yEField[exindex+1][eyindex ][ezinde 218 yEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal + 221 yEField[exindex+1][eyindex+1][ezinde 219 yEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) + 222 yEField[exindex+1][eyindex+1][ezinde 220 yEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ; 223 Efield[5] = 221 Efield[5] = 224 zEField[exindex ][eyindex ][ezinde 222 zEField[exindex ][eyindex ][ezindex ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) + 225 zEField[exindex ][eyindex ][ezinde 223 zEField[exindex ][eyindex ][ezindex+1] * (1-exlocal) * (1-eylocal) * ezlocal + 226 zEField[exindex ][eyindex+1][ezinde 224 zEField[exindex ][eyindex+1][ezindex ] * (1-exlocal) * eylocal * (1-ezlocal) + 227 zEField[exindex ][eyindex+1][ezinde 225 zEField[exindex ][eyindex+1][ezindex+1] * (1-exlocal) * eylocal * ezlocal + 228 zEField[exindex+1][eyindex ][ezinde 226 zEField[exindex+1][eyindex ][ezindex ] * exlocal * (1-eylocal) * (1-ezlocal) + 229 zEField[exindex+1][eyindex ][ezinde 227 zEField[exindex+1][eyindex ][ezindex+1] * exlocal * (1-eylocal) * ezlocal + 230 zEField[exindex+1][eyindex+1][ezinde 228 zEField[exindex+1][eyindex+1][ezindex ] * exlocal * eylocal * (1-ezlocal) + 231 zEField[exindex+1][eyindex+1][ezinde 229 zEField[exindex+1][eyindex+1][ezindex+1] * exlocal * eylocal * ezlocal ; 232 } 230 } 233 //G4cout << "Getting electric field " << Efiel << 231 //G4cout << "Getting electric field " << Efield[3]/(volt/m) << " " << Efield[4]/(volt/m) << " " << Efield[5]/(volt/m) << endl; 234 //G4cout << "For coordinates: " << Epoint[0] < 232 //G4cout << "For coordinates: " << Epoint[0] << " " << Epoint[1] << " " << Epoint[2] << G4endl; 235 233 236 /*std::ofstream WriteDataIn("ElectricFieldFC.o 234 /*std::ofstream WriteDataIn("ElectricFieldFC.out", std::ios::app); 237 WriteDataIn << Epoint[0] 235 WriteDataIn << Epoint[0] << '\t' << " " 238 << Epoint[1] << '\t' << " 236 << Epoint[1] << '\t' << " " 239 << Epoint[2] << '\t' << " 237 << Epoint[2] << '\t' << " " 240 << Efield[3]/(volt/m) << '\ 238 << Efield[3]/(volt/m) << '\t' << " " 241 << Efield[4]/(volt/m) << '\t 239 << Efield[4]/(volt/m) << '\t' << " " 242 << Efield[5]/(volt/m) << '\t 240 << Efield[5]/(volt/m) << '\t' << " " 243 << std::endl; */ << 241 << G4endl; */ 244 } 242 } 245 243