Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/hadrontherapy/src/HadrontherapyElectricTabulatedField3D.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 
 29 #include "HadrontherapyElectricTabulatedField3D.hh"
 30 #include "G4SystemOfUnits.hh"
 31 #include "G4AutoLock.hh"
 32 
 33 namespace{  G4Mutex MyHadrontherapyLockEField=G4MUTEX_INITIALIZER;  }
 34 
 35 using namespace std;
 36 
 37 HadrontherapyElectricTabulatedField3D::HadrontherapyElectricTabulatedField3D( const char* filename, G4double exOffset, G4double eyOffset, G4double ezOffset)
 38   :feXoffset(exOffset),feYoffset(eyOffset),feZoffset(ezOffset),einvertX(false),einvertY(false),einvertZ(false)
 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 from " << filename << " ... " << G4endl;
 49   G4AutoLock lock(&MyHadrontherapyLockEField);
 50 
 51   ifstream file( filename ); // Open the file for reading.
 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 order
 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 >> Ey >> Ez;
 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 " << G4endl;
111 
112   // G4cout << " Read values of field from file " << filename << G4endl;
113   G4cout << " ---> assumed the order:  x, y, z, Ex, Ey, Ez "
114    << "\n ---> Min values x,y,z: "
115    << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
116    << "\n ---> Max values x,y,z: "
117    << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm "
118    << "\n ---> The field will be offset in x by " << exOffset/cm << " cm "
119          << "\n ---> The field will be offset in y by " << eyOffset/cm << " cm "
120          << "\n ---> The field will be offset in z by " << ezOffset/cm << " cm " << G4endl;
121 
122   // Should really check that the limits are not the wrong way around.
123   if (Emaxx < Eminx) {swap(Emaxx,Eminx); einvertX = true;}
124   if (Emaxy < Eminy) {swap(Emaxy,Eminy); einvertY = true;}
125   if (Emaxz < Eminz) {swap(Emaxz,Eminz); einvertZ = true;}
126   G4cout << "\nAfter reordering if neccesary"
127    << "\n ---> Min values x,y,z: "
128    << Eminx/cm << " " << Eminy/cm << " " << Eminz/cm << " cm "
129    << " \n ---> Max values x,y,z: "
130    << Emaxx/cm << " " << Emaxy/cm << " " << Emaxz/cm << " cm ";
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 << " cm  "
137    << "\n-----------------------------------------------------------" << G4endl;
138 }
139 
140 void HadrontherapyElectricTabulatedField3D::GetFieldValue(const G4double Epoint[4],
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, normalized to the range
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 - Exfraction;}
154     if (einvertY) { Eyfraction = 1 - Eyfraction;}
155     if (einvertZ) { Ezfraction = 1 - Ezfraction;}
156 
157     // Need addresses of these to pass to modf below.
158     // modf uses its second argument as an OUTPUT argument.
159     G4double exdindex, eydindex, ezdindex;
160 
161     // Position of the point within the cuboid defined by the
162     // nearest surrounding tabulated points
163     G4double exlocal = ( std::modf(Exfraction*(Enx-1), &exdindex));
164     G4double eylocal = ( std::modf(Eyfraction*(Eny-1), &eydindex));
165     G4double ezlocal = ( std::modf(Ezfraction*(Enz-1), &ezdindex));
166 
167     // The indices of the nearest tabulated point whose coordinates
168     // are all less than those of the given point
169     G4int exindex = static_cast<G4int>(std::floor(exdindex));
170     G4int eyindex = static_cast<G4int>(std::floor(eydindex));
171     G4int ezindex = static_cast<G4int>(std::floor(ezdindex));
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 << " " << eylocal << " " << ezlocal << G4endl;
190     G4cout << "Index x,y,z: " << exindex << " " << eyindex << " " << ezindex << G4endl;
191     G4double valx0z0, mulx0z0, valx1z0, mulx1z0;
192     G4double valx0z1, mulx0z1, valx1z1, mulx1z1;
193     valx0z0= table[exindex  ][0][ezindex];  mulx0z0=  (1-exlocal) * (1-ezlocal);
194     valx1z0= table[exindex+1][0][ezindex];  mulx1z0=   exlocal    * (1-ezlocal);
195     valx0z1= table[exindex  ][0][ezindex+1]; mulx0z1= (1-exlocal) * ezlocal;
196     valx1z1= table[exindex+1][0][ezindex+1]; mulx1z1=  exlocal    * ezlocal;
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  ][ezindex  ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
207           xEField[exindex  ][eyindex  ][ezindex+1] * (1-exlocal) * (1-eylocal) *    ezlocal  +
208           xEField[exindex  ][eyindex+1][ezindex  ] * (1-exlocal) *    eylocal  * (1-ezlocal) +
209           xEField[exindex  ][eyindex+1][ezindex+1] * (1-exlocal) *    eylocal  *    ezlocal  +
210           xEField[exindex+1][eyindex  ][ezindex  ] *    exlocal  * (1-eylocal) * (1-ezlocal) +
211           xEField[exindex+1][eyindex  ][ezindex+1] *    exlocal  * (1-eylocal) *    ezlocal  +
212           xEField[exindex+1][eyindex+1][ezindex  ] *    exlocal  *    eylocal  * (1-ezlocal) +
213           xEField[exindex+1][eyindex+1][ezindex+1] *    exlocal  *    eylocal  *    ezlocal ;
214         Efield[4] =
215           yEField[exindex  ][eyindex  ][ezindex  ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
216           yEField[exindex  ][eyindex  ][ezindex+1] * (1-exlocal) * (1-eylocal) *    ezlocal  +
217           yEField[exindex  ][eyindex+1][ezindex  ] * (1-exlocal) *    eylocal  * (1-ezlocal) +
218           yEField[exindex  ][eyindex+1][ezindex+1] * (1-exlocal) *    eylocal  *    ezlocal  +
219           yEField[exindex+1][eyindex  ][ezindex  ] *    exlocal  * (1-eylocal) * (1-ezlocal) +
220           yEField[exindex+1][eyindex  ][ezindex+1] *    exlocal  * (1-eylocal) *    ezlocal  +
221           yEField[exindex+1][eyindex+1][ezindex  ] *    exlocal  *    eylocal  * (1-ezlocal) +
222           yEField[exindex+1][eyindex+1][ezindex+1] *    exlocal  *    eylocal  *    ezlocal ;
223         Efield[5] =
224           zEField[exindex  ][eyindex  ][ezindex  ] * (1-exlocal) * (1-eylocal) * (1-ezlocal) +
225           zEField[exindex  ][eyindex  ][ezindex+1] * (1-exlocal) * (1-eylocal) *    ezlocal  +
226           zEField[exindex  ][eyindex+1][ezindex  ] * (1-exlocal) *    eylocal  * (1-ezlocal) +
227           zEField[exindex  ][eyindex+1][ezindex+1] * (1-exlocal) *    eylocal  *    ezlocal  +
228           zEField[exindex+1][eyindex  ][ezindex  ] *    exlocal  * (1-eylocal) * (1-ezlocal) +
229           zEField[exindex+1][eyindex  ][ezindex+1] *    exlocal  * (1-eylocal) *    ezlocal  +
230           zEField[exindex+1][eyindex+1][ezindex  ] *    exlocal  *    eylocal  * (1-ezlocal) +
231           zEField[exindex+1][eyindex+1][ezindex+1] *    exlocal  *    eylocal  *    ezlocal ;
232   }
233 //G4cout << "Getting electric field " << Efield[3]/(volt/m) << " " << Efield[4]/(volt/m) << " " << Efield[5]/(volt/m) << G4endl;
234 //G4cout << "For coordinates: " << Epoint[0] << " " << Epoint[1] << " " << Epoint[2] << G4endl;
235 
236 /*std::ofstream WriteDataIn("ElectricFieldFC.out", std::ios::app);
237        WriteDataIn  <<   Epoint[0]             << '\t' << "   "
238       <<   Epoint[1]           << '\t' << "   "
239       <<   Epoint[2]            << '\t' << "   "
240       <<   Efield[3]/(volt/m)            << '\t' << "   "
241       <<   Efield[4]/(volt/m)           << '\t' << "   "
242       <<   Efield[5]/(volt/m)           << '\t' << "   "
243       << std::endl;    */
244 }
245